LSST Applications g04c3c9f7ca+2075667efa,g1e125bf412+5f448d5fcf,g2079a07aa2+3e9fd84d81,g2305ad1205+b635cf1488,g2bbee38e9b+6c6beb4891,g337abbeb29+6c6beb4891,g33d1c0ed96+6c6beb4891,g3a166c0a6a+6c6beb4891,g50ff169b8f+96c6868917,g52b1c1532d+585e252eca,g591dd9f2cf+42f171e1e6,g5c3423f6d4+d536b04327,g607f77f49a+d536b04327,g6f43f06aed+ca1339dc19,g858d7b2824+d536b04327,g8ee334c5b4+d7f9608c2f,g9963eaa53e+b3dc1655d3,g998f4353bf+d536b04327,g99cad8db69+8ef2408349,g9ddcbc5298+9a081db1e4,ga1e77700b3+2cbb763275,gadfd92a7e4+aec2f3b930,gae0086650b+585e252eca,gb0e22166c9+0e73c8378f,gb3b7280ab2+cb5fdb229e,gbb8dafda3b+a327199e22,gc120e1dc64+88074880ea,gc28159a63d+6c6beb4891,gcdd4ae20e8+bd241b2308,gcde1bda545+903e937d91,gcf0d15dbbd+bd241b2308,gdaeeff99f8+f9a426f77a,gddc38dedce+585e252eca,ge79ae78c31+6c6beb4891,gfbcc870c63+b310236976,w.2024.23
LSST Data Management Base Package
Loading...
Searching...
No Matches
Public Member Functions | List of all members
lsst::jointcal::AstrometryTransformPolynomial Class Reference

Polynomial transformation class. More...

#include <AstrometryTransform.h>

Inheritance diagram for lsst::jointcal::AstrometryTransformPolynomial:
lsst::jointcal::AstrometryTransform lsst::jointcal::AstrometryTransformLinear lsst::jointcal::AstrometryTransformLinearRot lsst::jointcal::AstrometryTransformLinearScale lsst::jointcal::AstrometryTransformLinearShift

Public Member Functions

 AstrometryTransformPolynomial (std::size_t order=1)
 Default transform : identity for all orders (>=1 ).
 
 AstrometryTransformPolynomial (const AstrometryTransform *transform, const Frame &frame, std::size_t order, std::size_t nPoint=1000)
 Constructs a "polynomial image" from an existing transform, over a specified domain.
 
 AstrometryTransformPolynomial (const std::shared_ptr< afw::geom::TransformPoint2ToPoint2 > &transform, jointcal::Frame const &domain, std::size_t order, std::size_t nSteps=50)
 Constructs a polynomial approximation to an afw::geom::TransformPoint2ToPoint2.
 
void setOrder (std::size_t order)
 Sets the polynomial order (the highest sum of exponents of the largest monomial).
 
std::size_t getOrder () const
 Returns the polynomial order.
 
void apply (double xIn, double yIn, double &xOut, double &yOut) const override
 
void computeDerivative (Point const &where, AstrometryTransformLinear &derivative, double step=0.01) const override
 specialised analytic routine
 
virtual void transformPosAndErrors (const FatPoint &in, FatPoint &out) const override
 a mix of apply and Derivative
 
std::size_t getNpar () const override
 total number of parameters
 
void print (std::ostream &out) const override
 print out of coefficients in a readable form.
 
double fit (StarMatchList const &starMatchList) override
 guess what
 
AstrometryTransformPolynomial operator* (AstrometryTransformPolynomial const &right) const
 Composition (internal stuff in quadruple precision)
 
AstrometryTransformPolynomial operator+ (AstrometryTransformPolynomial const &right) const
 Addition.
 
AstrometryTransformPolynomial operator- (AstrometryTransformPolynomial const &right) const
 Subtraction.
 
std::unique_ptr< AstrometryTransformcomposeAndReduce (AstrometryTransformPolynomial const &right) const
 Return a reduced composition of newTransform = this(right()), or nullptr if it cannot be reduced.
 
std::unique_ptr< AstrometryTransformclone () const override
 returns a copy (allocated by new) of the transformation.
 
double coeffOrZero (std::size_t powX, std::size_t powY, std::size_t whichCoord) const
 read access, zero if beyond order
 
double determinant () const
 
double paramRef (Eigen::Index i) const override
 
double & paramRef (Eigen::Index i) override
 
void paramDerivatives (Point const &where, double *dx, double *dy) const override
 Derivative w.r.t parameters.
 
std::shared_ptr< ast::MappingtoAstMap (jointcal::Frame const &domain) const override
 Create an equivalent AST mapping for this transformation, including an analytic inverse if possible.
 
void write (std::ostream &s) const override
 
void read (std::istream &s)
 
virtual void apply (double xIn, double yIn, double &xOut, double &yOut) const=0
 
void apply (Point const &in, Point &out) const
 applies the tranfo to in and writes into out. Is indeed virtual.
 
Point apply (Point const &in) const
 All these apply(..) shadow the virtual one in derived classes, unless one writes "using AstrometryTransform::apply".
 
Frame apply (Frame const &inputframe, bool inscribed) const
 Transform a bounding box, taking either the inscribed or circumscribed box.
 
virtual std::unique_ptr< AstrometryTransformcomposeAndReduce (AstrometryTransform const &right) const
 Return a reduced composition of newTransform = this(right()), or nullptr if it cannot be reduced.
 
std::string __str__ () const
 
void transformStar (FatPoint &in) const
 
virtual double getJacobian (Point const &point) const
 returns the local jacobian.
 
virtual double getJacobian (double x, double y) const
 returns the local jacobian.
 
virtual AstrometryTransformLinear linearApproximation (Point const &where, double step=0.01) const
 linear (local) approximation.
 
virtual void transformErrors (Point const &where, const double *vIn, double *vOut) const
 transform errors (represented as double[3] in order V(xx),V(yy),Cov(xy))
 
virtual std::unique_ptr< AstrometryTransforminverseTransform (double precision, const Frame &region) const
 returns an inverse transform. Numerical if not overloaded.
 
void getParams (double *params) const
 params should be at least Npar() long
 
void offsetParams (Eigen::VectorXd const &delta)
 
virtual std::unique_ptr< AstrometryTransformroughInverse (const Frame &region) const
 Rough inverse.
 
void write (const std::string &fileName) const
 
double getCoefficient (std::size_t powX, std::size_t powY, std::size_t whichCoord) const
 Get the coefficient of a given power in x and y, for either the x or y coordinate.
 
double & getCoefficient (std::size_t powX, std::size_t powY, std::size_t whichCoord)
 Get the coefficient of a given power in x and y, for either the x or y coordinate.
 

Detailed Description

Polynomial transformation class.

Definition at line 280 of file AstrometryTransform.h.

Constructor & Destructor Documentation

◆ AstrometryTransformPolynomial() [1/3]

lsst::jointcal::AstrometryTransformPolynomial::AstrometryTransformPolynomial ( std::size_t order = 1)

Default transform : identity for all orders (>=1 ).

Default transform : identity for all orders (>=1 )

Parameters
orderThe highest total power (x+y) of monomials of this polynomial.

Definition at line 471 of file AstrometryTransform.cc.

471 : _order(order) {
472 _nterms = (order + 1) * (order + 2) / 2;
473
474 // allocate and fill coefficients
475 _coeffs.resize(2 * _nterms, 0.);
476 // the default is supposed to be the identity, (for order>=1).
477 if (_order >= 1) {
478 getCoefficient(1, 0, 0) = 1;
479 getCoefficient(0, 1, 1) = 1;
480 }
481}
double getCoefficient(std::size_t powX, std::size_t powY, std::size_t whichCoord) const
Get the coefficient of a given power in x and y, for either the x or y coordinate.
T resize(T... args)
table::Key< int > order

◆ AstrometryTransformPolynomial() [2/3]

lsst::jointcal::AstrometryTransformPolynomial::AstrometryTransformPolynomial ( const AstrometryTransform * transform,
const Frame & frame,
std::size_t order,
std::size_t nPoint = 1000 )

Constructs a "polynomial image" from an existing transform, over a specified domain.

Definition at line 484 of file AstrometryTransform.cc.

486 {
487 StarMatchList sm;
488
489 double step = std::sqrt(fabs(frame.getArea()) / double(nPoint));
490 for (double x = frame.xMin + step / 2; x <= frame.xMax; x += step)
491 for (double y = frame.yMin + step / 2; y <= frame.yMax; y += step) {
492 auto pix = std::make_shared<BaseStar>(x, y, 0, 0);
493 double xtr, ytr;
494 transform->apply(x, y, xtr, ytr);
495 auto tp = std::make_shared<BaseStar>(xtr, ytr, 0, 0);
496 /* These are fake stars so no need to transform fake errors.
497 all errors (and weights) will be equal : */
498 sm.push_back(StarMatch(*pix, *tp, pix, tp));
499 }
501 ret.fit(sm);
502 *this = ret;
503}
int const step
int y
Definition SpanSet.cc:48
AstrometryTransformPolynomial(std::size_t order=1)
Default transform : identity for all orders (>=1 ).
T fabs(T... args)
T sqrt(T... args)

◆ AstrometryTransformPolynomial() [3/3]

lsst::jointcal::AstrometryTransformPolynomial::AstrometryTransformPolynomial ( const std::shared_ptr< afw::geom::TransformPoint2ToPoint2 > & transform,
jointcal::Frame const & domain,
std::size_t order,
std::size_t nSteps = 50 )

Constructs a polynomial approximation to an afw::geom::TransformPoint2ToPoint2.

Parameters
[in]transformThe transform to be approximated.
[in]domainThe valid domain of the transform.
[in]orderThe polynomial order to use when approximating.
[in]nStepsThe number of sample points per axis (nSteps^2 total points).

Definition at line 506 of file AstrometryTransform.cc.

508 {
509 jointcal::StarMatchList starMatchList;
510 double xStart = domain.xMin;
511 double yStart = domain.yMin;
512 double xStep = domain.getWidth() / (nSteps + 1);
513 double yStep = domain.getHeight() / (nSteps + 1);
514 for (std::size_t i = 0; i < nSteps; ++i) {
515 for (std::size_t j = 0; j < nSteps; ++j) {
516 // TODO: once DM-4044 is done, we can remove the redundancy in `Point`/`Point2D` here
517 jointcal::Point in(xStart + i * xStep, yStart + j * yStep);
518 geom::Point2D inAfw(in.x, in.y);
519 geom::Point2D outAfw = transform->applyForward(inAfw);
520 jointcal::Point out(outAfw.getX(), outAfw.getY());
521 starMatchList.emplace_back(in, out, nullptr, nullptr);
522 }
523 }
525 poly.fit(starMatchList);
526 *this = poly;
527}
A point in a plane.
Definition Point.h:37
T emplace_back(T... args)
Low-level polynomials (including special polynomials) in C++.

Member Function Documentation

◆ __str__()

std::string lsst::jointcal::AstrometryTransform::__str__ ( ) const
inlineinherited

Definition at line 94 of file AstrometryTransform.h.

94 {
96 print(s);
97 return s.str();
98 }
virtual void print(std::ostream &out) const =0
prints the transform coefficients to stream.

◆ apply() [1/5]

void lsst::jointcal::AstrometryTransformPolynomial::apply ( double xIn,
double yIn,
double & xOut,
double & yOut ) const
overridevirtual

Implements lsst::jointcal::AstrometryTransform.

Definition at line 573 of file AstrometryTransform.cc.

574 {
575 /*
576 This routine computes the monomials only once for both
577 polynomials. This is why AstrometryTransformPolynomial does not use an auxilary
578 class (such as PolyXY) to handle each polynomial.
579
580 The code works even if &xIn == &xOut (or &yIn == &yOut)
581 It uses Variable Length Allocation (VLA) rather than a vector<double>
582 because allocating the later costs about 50 ns. All VLA uses are tagged.
583 */
584 double monomials[_nterms]; // this is VLA, which is (perhaps) not casher C++
585 computeMonomials(xIn, yIn, monomials);
586
587 xOut = 0;
588 yOut = 0;
589 const double *c = &_coeffs[0];
590 const double *pm = &monomials[0];
591 // the ordering of the coefficients and the monomials are identical.
592 for (int k = _nterms; k--;) xOut += (*(pm++)) * (*(c++));
593 pm = &monomials[0];
594 for (int k = _nterms; k--;) yOut += (*(pm++)) * (*(c++));
595}

◆ apply() [2/5]

virtual void lsst::jointcal::AstrometryTransform::apply ( double xIn,
double yIn,
double & xOut,
double & yOut ) const
virtual

◆ apply() [3/5]

Frame lsst::jointcal::AstrometryTransform::apply ( Frame const & inputframe,
bool inscribed ) const

Transform a bounding box, taking either the inscribed or circumscribed box.

Parameters
[in]inputframeThe frame to be transformed.
[in]inscribedReturn the inscribed (true) or circumscribed (false) box.
Returns
The transformed frame.

Definition at line 89 of file AstrometryTransform.cc.

77 {
78 // 2 opposite corners
79 double xtmin1, xtmax1, ytmin1, ytmax1;
80 apply(inputframe.xMin, inputframe.yMin, xtmin1, ytmin1);
81 apply(inputframe.xMax, inputframe.yMax, xtmax1, ytmax1);
82 Frame fr1(std::min(xtmin1, xtmax1), std::min(ytmin1, ytmax1), std::max(xtmin1, xtmax1),
83 std::max(ytmin1, ytmax1));
84 // 2 other corners
85 double xtmin2, xtmax2, ytmin2, ytmax2;
86 apply(inputframe.xMin, inputframe.yMax, xtmin2, ytmax2);
87 apply(inputframe.xMax, inputframe.yMin, xtmax2, ytmin2);
88 Frame fr2(std::min(xtmin2, xtmax2), std::min(ytmin2, ytmax2), std::max(xtmin2, xtmax2),
89 std::max(ytmin2, ytmax2));
90
91 if (inscribed) return fr1 * fr2;
92 return fr1 + fr2;
93}
void apply(double xIn, double yIn, double &xOut, double &yOut) const override
T max(T... args)
T min(T... args)

◆ apply() [4/5]

Point lsst::jointcal::AstrometryTransform::apply ( Point const & in) const
inline

All these apply(..) shadow the virtual one in derived classes, unless one writes "using AstrometryTransform::apply".

Definition at line 75 of file AstrometryTransform.h.

75 {
76 double xout, yout;
77 apply(in.x, in.y, xout, yout);
78 return Point(xout, yout);
79 }

◆ apply() [5/5]

void lsst::jointcal::AstrometryTransform::apply ( Point const & in,
Point & out ) const
inline

applies the tranfo to in and writes into out. Is indeed virtual.

Definition at line 71 of file AstrometryTransform.h.

71{ apply(in.x, in.y, out.x, out.y); }

◆ clone()

std::unique_ptr< AstrometryTransform > lsst::jointcal::AstrometryTransformPolynomial::clone ( ) const
inlineoverridevirtual

returns a copy (allocated by new) of the transformation.

Implements lsst::jointcal::AstrometryTransform.

Definition at line 344 of file AstrometryTransform.h.

◆ coeffOrZero()

double lsst::jointcal::AstrometryTransformPolynomial::coeffOrZero ( std::size_t powX,
std::size_t powY,
std::size_t whichCoord ) const

read access, zero if beyond order

Definition at line 756 of file AstrometryTransform.cc.

757 {
758 // assert((degX+degY<=order) && whichCoord<2);
759 assert(whichCoord < 2);
760 if (degX + degY <= _order)
761 return _coeffs[(degX + degY) * (degX + degY + 1) / 2 + degY + whichCoord * _nterms];
762 return 0;
763}

◆ composeAndReduce() [1/2]

std::unique_ptr< AstrometryTransform > lsst::jointcal::AstrometryTransform::composeAndReduce ( AstrometryTransform const & right) const
virtual

Return a reduced composition of newTransform = this(right()), or nullptr if it cannot be reduced.

"Reduced" in this context means that they are capable of being merged into a single transform, for example, for two polynomials:

\[ f(x) = 1 + x^2, g(x) = -1 + 3x \]

we would have h = f.composeAndReduce(g) == 2 - 6x + 9x^2.

To be overloaded by derived classes if they can properly reduce the composition.

Parameters
rightThe transform to apply first.
Returns
The new reduced and composed AstrometryTransform, or nullptr if no such reduction is possible.

Reimplemented from lsst::jointcal::AstrometryTransform.

Definition at line 131 of file AstrometryTransform.cc.

96 { // by default no way to compose
98}

◆ composeAndReduce() [2/2]

std::unique_ptr< AstrometryTransform > lsst::jointcal::AstrometryTransformPolynomial::composeAndReduce ( AstrometryTransformPolynomial const & right) const

Return a reduced composition of newTransform = this(right()), or nullptr if it cannot be reduced.

"Reduced" in this context means that they are capable of being merged into a single transform, for example, for two polynomials:

\[ f(x) = 1 + x^2, g(x) = -1 + 3x \]

we would have h = f.composeAndReduce(g) == 2 - 6x + 9x^2.

To be overloaded by derived classes if they can properly reduce the composition.

Parameters
rightThe transform to apply first.
Returns
The new reduced and composed AstrometryTransform, or nullptr if no such reduction is possible.

Definition at line 944 of file AstrometryTransform.cc.

945 {
946 if (getOrder() == 1 && right.getOrder() == 1)
947 return std::make_unique<AstrometryTransformLinear>((*this) * (right)); // does the composition
948 else
949 return std::make_unique<AstrometryTransformPolynomial>((*this) * (right)); // does the composition
950}
std::size_t getOrder() const
Returns the polynomial order.
T right(T... args)

◆ computeDerivative()

void lsst::jointcal::AstrometryTransformPolynomial::computeDerivative ( Point const & where,
AstrometryTransformLinear & derivative,
double step = 0.01 ) const
overridevirtual

specialised analytic routine

Reimplemented from lsst::jointcal::AstrometryTransform.

Definition at line 597 of file AstrometryTransform.cc.

600 { /* routine checked against numerical derivatives from AstrometryTransform::Derivative */
601 if (_order == 1) {
602 derivative = AstrometryTransformLinear(*this);
603 derivative.dx() = derivative.dy() = 0;
604 return;
605 }
606
607 double dermx[2 * _nterms]; // VLA
608 double *dermy = dermx + _nterms;
609 double xin = where.x;
610 double yin = where.y;
611
612 double xx = 1;
613 double xxm1 = 1; // xx^(ix-1)
614 for (std::size_t ix = 0; ix <= _order; ++ix) {
615 std::size_t k = (ix) * (ix + 1) / 2;
616 // iy = 0
617 dermx[k] = ix * xxm1;
618 dermy[k] = 0;
619 k += ix + 2;
620 double yym1 = 1; // yy^(iy-1)
621 for (std::size_t iy = 1; iy <= _order - ix; ++iy) {
622 dermx[k] = ix * xxm1 * yym1 * yin;
623 dermy[k] = iy * xx * yym1;
624 yym1 *= yin;
625 k += ix + iy + 2;
626 }
627 xx *= xin;
628 if (ix >= 1) xxm1 *= xin;
629 }
630
631 derivative.dx() = 0;
632 derivative.dy() = 0;
633
634 const double *mx = &dermx[0];
635 const double *my = &dermy[0];
636 const double *c = &_coeffs[0];
637 // dx'
638 double a11 = 0, a12 = 0;
639 for (int k = _nterms; k--;) {
640 a11 += (*(mx++)) * (*c);
641 a12 += (*(my++)) * (*(c++));
642 }
643 derivative.a11() = a11;
644 derivative.a12() = a12;
645 // dy'
646 double a21 = 0, a22 = 0;
647 mx = &dermx[0];
648 my = &dermy[0];
649 for (int k = _nterms; k--;) {
650 a21 += (*(mx++)) * (*c);
651 a22 += (*(my++)) * (*(c++));
652 }
653 derivative.a21() = a21;
654 derivative.a22() = a22;
655}

◆ determinant()

double lsst::jointcal::AstrometryTransformPolynomial::determinant ( ) const

Definition at line 826 of file AstrometryTransform.cc.

826 {
827 if (_order >= 1)
828 return getCoefficient(1, 0, 0) * getCoefficient(0, 1, 1) -
829 getCoefficient(0, 1, 0) * getCoefficient(1, 0, 1);
830 return 0;
831}

◆ fit()

double lsst::jointcal::AstrometryTransformPolynomial::fit ( StarMatchList const & starMatchList)
overridevirtual

guess what

Implements lsst::jointcal::AstrometryTransform.

Definition at line 926 of file AstrometryTransform.cc.

926 {
927 if (starMatchList.size() < _nterms) {
928 LOGLS_FATAL(_log, "AstrometryTransformPolynomial::fit trying to fit a polynomial transform of order "
929 << _order << " with only " << starMatchList.size() << " matches.");
930 return -1;
931 }
932
933 AstrometryTransformPolynomial conditionner = shiftAndNormalize(starMatchList);
934
935 computeFit(starMatchList, conditionner, false); // get a rough solution
936 computeFit(starMatchList, conditionner, true); // weight with it
937 double chi2 = computeFit(starMatchList, conditionner, true); // once more
938
939 (*this) = (*this) * conditionner;
940 if (starMatchList.size() == _nterms) return 0;
941 return chi2;
942}
#define LOGLS_FATAL(logger, message)
Log a fatal-level message using an iostream-based interface.
Definition Log.h:699

◆ getCoefficient() [1/2]

double & lsst::jointcal::AstrometryTransformPolynomial::getCoefficient ( std::size_t powX,
std::size_t powY,
std::size_t whichCoord )

Get the coefficient of a given power in x and y, for either the x or y coordinate.

Definition at line 750 of file AstrometryTransform.cc.

751 {
752 assert((degX + degY <= _order) && whichCoord < 2);
753 return _coeffs[(degX + degY) * (degX + degY + 1) / 2 + degY + whichCoord * _nterms];
754}

◆ getCoefficient() [2/2]

double lsst::jointcal::AstrometryTransformPolynomial::getCoefficient ( std::size_t powX,
std::size_t powY,
std::size_t whichCoord ) const

Get the coefficient of a given power in x and y, for either the x or y coordinate.

Definition at line 741 of file AstrometryTransform.cc.

742 {
743 assert((degX + degY <= _order) && whichCoord < 2);
744 /* this assertion above is enough to ensure that the index used just
745 below is within bounds since the reserved length is
746 2*_nterms=(order+1)*(order+2) */
747 return _coeffs[(degX + degY) * (degX + degY + 1) / 2 + degY + whichCoord * _nterms];
748}

◆ getJacobian() [1/2]

double lsst::jointcal::AstrometryTransform::getJacobian ( double x,
double y ) const
virtualinherited

returns the local jacobian.

Definition at line 100 of file AstrometryTransform.cc.

100 {
101 double x2, y2;
102 double eps = x * 0.01;
103 if (eps == 0) eps = 0.01;
104 apply(x, y, x2, y2);
105 double dxdx, dydx;
106 apply(x + eps, y, dxdx, dydx);
107 dxdx -= x2;
108 dydx -= y2;
109 double dxdy, dydy;
110 apply(x, y + eps, dxdy, dydy);
111 dxdy -= x2;
112 dydy -= y2;
113 return ((dxdx * dydy - dxdy * dydx) / (eps * eps));
114}
virtual void apply(double xIn, double yIn, double &xOut, double &yOut) const =0

◆ getJacobian() [2/2]

virtual double lsst::jointcal::AstrometryTransform::getJacobian ( Point const & point) const
inlinevirtualinherited

returns the local jacobian.

Definition at line 110 of file AstrometryTransform.h.

110{ return getJacobian(point.x, point.y); }
virtual double getJacobian(Point const &point) const
returns the local jacobian.

◆ getNpar()

std::size_t lsst::jointcal::AstrometryTransformPolynomial::getNpar ( ) const
inlineoverridevirtual

total number of parameters

Reimplemented from lsst::jointcal::AstrometryTransform.

Definition at line 321 of file AstrometryTransform.h.

321{ return 2 * _nterms; }

◆ getOrder()

std::size_t lsst::jointcal::AstrometryTransformPolynomial::getOrder ( ) const
inline

Returns the polynomial order.

Definition at line 307 of file AstrometryTransform.h.

307{ return _order; }

◆ getParams()

void lsst::jointcal::AstrometryTransform::getParams ( double * params) const
inherited

params should be at least Npar() long

Definition at line 217 of file AstrometryTransform.cc.

217 {
218 std::size_t npar = getNpar();
219 for (std::size_t i = 0; i < npar; ++i) params[i] = paramRef(i);
220}
virtual std::size_t getNpar() const
returns the number of parameters (to compute chi2's)
virtual double paramRef(Eigen::Index i) const

◆ inverseTransform()

std::unique_ptr< AstrometryTransform > lsst::jointcal::AstrometryTransform::inverseTransform ( double precision,
const Frame & region ) const
virtualinherited

returns an inverse transform. Numerical if not overloaded.

precision and region refer to the "input" side of this, and hence to the output side of the returned AstrometryTransform.

Reimplemented in lsst::jointcal::TanPixelToRaDec, lsst::jointcal::TanSipPixelToRaDec, lsst::jointcal::TanRaDecToPixel, lsst::jointcal::AstrometryTransformLinear, and lsst::jointcal::AstrometryTransformInverse.

Definition at line 304 of file AstrometryTransform.cc.

305 {
306 return std::unique_ptr<AstrometryTransform>(new AstrometryTransformInverse(this, precision, region));
307}

◆ linearApproximation()

AstrometryTransformLinear lsst::jointcal::AstrometryTransform::linearApproximation ( Point const & where,
double step = 0.01 ) const
virtualinherited

linear (local) approximation.

Reimplemented in lsst::jointcal::AstrometryTransformIdentity, and lsst::jointcal::AstrometryTransformLinear.

Definition at line 137 of file AstrometryTransform.cc.

138 {
139 Point outwhere = apply(where);
140 AstrometryTransformLinear der;
141 computeDerivative(where, der, step);
142 return AstrometryTransformLinearShift(outwhere.x, outwhere.y) * der *
143 AstrometryTransformLinearShift(-where.x, -where.y);
144}
virtual void computeDerivative(Point const &where, AstrometryTransformLinear &derivative, double step=0.01) const
Computes the local Derivative of a transform, w.r.t.

◆ offsetParams()

void lsst::jointcal::AstrometryTransform::offsetParams ( Eigen::VectorXd const & delta)
inherited

Definition at line 222 of file AstrometryTransform.cc.

222 {
223 std::size_t npar = getNpar();
224 for (std::size_t i = 0; i < npar; ++i) paramRef(i) += delta[i];
225}

◆ operator*()

AstrometryTransformPolynomial lsst::jointcal::AstrometryTransformPolynomial::operator* ( AstrometryTransformPolynomial const & right) const

Composition (internal stuff in quadruple precision)

Definition at line 1052 of file AstrometryTransform.cc.

1053 {
1054 // split each transform into 2d polynomials
1055 PolyXY plx(*this, 0);
1056 PolyXY ply(*this, 1);
1057 PolyXY prx(right, 0);
1058 PolyXY pry(right, 1);
1059
1060 // compute the compositions
1061 PolyXY rx(composition(plx, prx, pry));
1062 PolyXY ry(composition(ply, prx, pry));
1063
1064 // copy the results the hard way.
1066 for (std::size_t px = 0; px <= result._order; ++px)
1067 for (std::size_t py = 0; py <= result._order - px; ++py) {
1068 result.getCoefficient(px, py, 0) = rx.getCoefficient(px, py);
1069 result.getCoefficient(px, py, 1) = ry.getCoefficient(px, py);
1070 }
1071 return result;
1072}
py::object result
Definition _schema.cc:429

◆ operator+()

AstrometryTransformPolynomial lsst::jointcal::AstrometryTransformPolynomial::operator+ ( AstrometryTransformPolynomial const & right) const

Addition.

Definition at line 1074 of file AstrometryTransform.cc.

1075 {
1076 if (_order >= right._order) {
1078 for (std::size_t i = 0; i <= right._order; ++i)
1079 for (std::size_t j = 0; j <= right._order - i; ++j) {
1080 res.getCoefficient(i, j, 0) += right.getCoefficient(i, j, 0);
1081 res.getCoefficient(i, j, 1) += right.getCoefficient(i, j, 1);
1082 }
1083 return res;
1084 } else
1085 return (right + (*this));
1086}

◆ operator-()

AstrometryTransformPolynomial lsst::jointcal::AstrometryTransformPolynomial::operator- ( AstrometryTransformPolynomial const & right) const

Subtraction.

Definition at line 1088 of file AstrometryTransform.cc.

1089 {
1090 AstrometryTransformPolynomial res(std::max(_order, right._order));
1091 for (std::size_t i = 0; i <= res._order; ++i)
1092 for (std::size_t j = 0; j <= res._order - i; ++j) {
1093 res.getCoefficient(i, j, 0) = coeffOrZero(i, j, 0) - right.coeffOrZero(i, j, 0);
1094 res.getCoefficient(i, j, 1) = coeffOrZero(i, j, 1) - right.coeffOrZero(i, j, 1);
1095 }
1096 return res;
1097}
double coeffOrZero(std::size_t powX, std::size_t powY, std::size_t whichCoord) const
read access, zero if beyond order

◆ paramDerivatives()

void lsst::jointcal::AstrometryTransformPolynomial::paramDerivatives ( Point const & where,
double * dx,
double * dy ) const
overridevirtual

Derivative w.r.t parameters.

Derivatives should be al least 2*NPar long. first Npar, for x, last Npar for y.

Reimplemented from lsst::jointcal::AstrometryTransform.

Definition at line 776 of file AstrometryTransform.cc.

777 { /* first half : dxout/dpar, second half : dyout/dpar */
778 computeMonomials(where.x, where.y, dx);
779 for (std::size_t k = 0; k < _nterms; ++k) {
780 dy[_nterms + k] = dx[k];
781 dx[_nterms + k] = dy[k] = 0;
782 }
783}

◆ paramRef() [1/2]

double lsst::jointcal::AstrometryTransformPolynomial::paramRef ( Eigen::Index i) const
overridevirtual

Reimplemented from lsst::jointcal::AstrometryTransform.

Definition at line 766 of file AstrometryTransform.cc.

766 {
767 assert(i < 2 * Eigen::Index(_nterms));
768 return _coeffs[i];
769}

◆ paramRef() [2/2]

double & lsst::jointcal::AstrometryTransformPolynomial::paramRef ( Eigen::Index i)
overridevirtual

Reimplemented from lsst::jointcal::AstrometryTransform.

Definition at line 771 of file AstrometryTransform.cc.

771 {
772 assert(i < 2 * Eigen::Index(_nterms));
773 return _coeffs[i];
774}

◆ print()

void lsst::jointcal::AstrometryTransformPolynomial::print ( std::ostream & out) const
overridevirtual

print out of coefficients in a readable form.

Implements lsst::jointcal::AstrometryTransform.

Definition at line 796 of file AstrometryTransform.cc.

796 {
797 auto oldPrecision = stream.precision();
798 stream.precision(12);
799 stream << "AstrometryTransformPolynomial: order=" << getOrder() << std::endl;
800 for (std::size_t ic = 0; ic < 2; ++ic) {
801 if (ic == 0)
802 stream << "newx = ";
803 else
804 stream << "newy = ";
805 bool printed = false; // whether we've printed any monomials
806 for (std::size_t p = 0; p <= _order; ++p)
807 for (std::size_t py = 0; py <= p; ++py) {
808 double coefficient = getCoefficient(p - py, py, ic);
809 if (coefficient != 0 || isnan(coefficient)) {
810 // Only print "interesting" coefficients.
811 if (printed && p + py != 0) stream << " + ";
812 printed = true;
813 stream << coefficient << monomialString(p - py, py);
814 }
815 }
816 // if all components are zero, nothing is output by the loops above
817 if (!printed) {
818 stream << 0;
819 }
820 stream << endl;
821 }
822 if (_order > 0) stream << "Linear determinant = " << determinant();
823 stream.precision(oldPrecision);
824}
T endl(T... args)
T isnan(T... args)

◆ read()

void lsst::jointcal::AstrometryTransformPolynomial::read ( std::istream & s)

Definition at line 1114 of file AstrometryTransform.cc.

1114 {
1115 int format;
1116 s >> format;
1117 if (format != 1)
1118 throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
1119 " AstrometryTransformPolynomial::read : format is not 1 ");
1120
1121 string order;
1122 s >> order >> _order;
1123 if (order != "order")
1124 throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
1125 " AstrometryTransformPolynomial::read : expecting \"order\" and found " + order);
1126 setOrder(_order);
1127 for (std::size_t k = 0; k < 2 * _nterms; ++k) s >> _coeffs[k];
1128}
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
void setOrder(std::size_t order)
Sets the polynomial order (the highest sum of exponents of the largest monomial).
format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition history.py:173

◆ roughInverse()

std::unique_ptr< AstrometryTransform > lsst::jointcal::AstrometryTransform::roughInverse ( const Frame & region) const
virtualinherited

Rough inverse.

Stored by the numerical inverter to guess starting point for the trials. Just here to enable overloading.

Reimplemented in lsst::jointcal::AstrometryTransformInverse, lsst::jointcal::TanPixelToRaDec, and lsst::jointcal::TanRaDecToPixel.

Definition at line 196 of file AstrometryTransform.cc.

196 {
197 // "in" and "out" refer to the inverse direction.
198 Point centerOut = region.getCenter();
199 Point centerIn = apply(centerOut);
200 AstrometryTransformLinear der;
201 computeDerivative(centerOut, der, std::sqrt(region.getArea()) / 5.);
202 der = der.inverted();
203 der = AstrometryTransformLinearShift(centerOut.x, centerOut.y) * der *
204 AstrometryTransformLinearShift(-centerIn.x, -centerIn.y);
205 return std::unique_ptr<AstrometryTransform>(new AstrometryTransformLinear(der));
206}

◆ setOrder()

void lsst::jointcal::AstrometryTransformPolynomial::setOrder ( std::size_t order)

Sets the polynomial order (the highest sum of exponents of the largest monomial).

Definition at line 552 of file AstrometryTransform.cc.

552 {
553 _order = order;
554 std::size_t old_nterms = _nterms;
555 _nterms = (_order + 1) * (_order + 2) / 2;
556
557 // temporarily save coefficients
558 vector<double> old_coeffs = _coeffs;
559 // reallocate enough size
560 _coeffs.resize(2 * _nterms);
561 // reassign to zero (this is necessary because ycoeffs
562 // are after xcoeffs and so their meaning changes
563 for (std::size_t k = 0; k < _nterms; ++k) _coeffs[k] = 0;
564 // put back what we had before
565 std::size_t kmax = min(old_nterms, _nterms);
566 for (std::size_t k = 0; k < kmax; ++k) {
567 _coeffs[k] = old_coeffs[k]; // x terms
568 _coeffs[k + _nterms] = old_coeffs[k + old_nterms]; // y terms
569 }
570}
int min

◆ toAstMap()

std::shared_ptr< ast::Mapping > lsst::jointcal::AstrometryTransformPolynomial::toAstMap ( jointcal::Frame const & domain) const
overridevirtual

Create an equivalent AST mapping for this transformation, including an analytic inverse if possible.

Parameters
domainThe domain of the transform, to help find an inverse.
Returns
An AST Mapping that represents this transformation.

Reimplemented from lsst::jointcal::AstrometryTransform.

Definition at line 1099 of file AstrometryTransform.cc.

1099 {
1100 auto inverse = inversePolyTransform(*this, domain, 1e-7, _order + 2, 100);
1101 return std::make_shared<ast::PolyMap>(toAstPolyMapCoefficients(), inverse->toAstPolyMapCoefficients());
1102}
std::shared_ptr< AstrometryTransformPolynomial > inversePolyTransform(AstrometryTransform const &forward, Frame const &domain, double precision, std::size_t maxOrder=9, std::size_t nSteps=50)
Approximate the inverse by a polynomial, to some precision.

◆ transformErrors()

void lsst::jointcal::AstrometryTransform::transformErrors ( Point const & where,
const double * vIn,
double * vOut ) const
virtualinherited

transform errors (represented as double[3] in order V(xx),V(yy),Cov(xy))

Definition at line 163 of file AstrometryTransform.cc.

163 {
164 AstrometryTransformLinear der;
165 computeDerivative(where, der, 0.01);
166 double a11 = der.A11();
167 double a22 = der.A22();
168 double a21 = der.A21();
169 double a12 = der.A12();
170
171 /* (a11 a12) (vxx vxy)
172 M = ( ) and V = ( )
173 (a21 a22) (xvy vyy)
174
175 Vxx = Vin[0], vyy = Vin[1], Vxy = Vin[2];
176 we want to compute M*V*tp(M)
177 A lin alg light package would be perfect...
178 */
179 int xx = 0;
180 int yy = 1;
181 int xy = 2;
182 // M*V :
183
184 double b11 = a11 * vIn[xx] + a12 * vIn[xy];
185 double b22 = a21 * vIn[xy] + a22 * vIn[yy];
186 double b12 = a11 * vIn[xy] + a12 * vIn[yy];
187 double b21 = a21 * vIn[xx] + a22 * vIn[xy];
188
189 // (M*V) * tp(M)
190
191 vOut[xx] = b11 * a11 + b12 * a12;
192 vOut[xy] = b11 * a21 + b12 * a22;
193 vOut[yy] = b21 * a21 + b22 * a22;
194}

◆ transformPosAndErrors()

void lsst::jointcal::AstrometryTransformPolynomial::transformPosAndErrors ( const FatPoint & in,
FatPoint & out ) const
overridevirtual

a mix of apply and Derivative

Reimplemented from lsst::jointcal::AstrometryTransform.

Definition at line 657 of file AstrometryTransform.cc.

657 {
658 /*
659 The results from this routine were compared to what comes out
660 from apply and transformErrors. The Derivative routine was
661 checked against numerical derivatives from
662 AstrometryTransform::Derivative. (P.A dec 2009).
663
664 This routine could be made much simpler by calling apply and
665 Derivative (i.e. you just suppress it, and the fallback is the
666 generic version in AstrometryTransform). BTW, I checked that both routines
667 provide the same result. This version is however faster
668 (monomials get recycled).
669 */
670 double monomials[_nterms]; // VLA
671
672 FatPoint res; // to store the result, because nothing forbids &in == &out.
673
674 double dermx[2 * _nterms]; // monomials for derivative w.r.t. x (VLA)
675 double *dermy = dermx + _nterms; // same for y
676 double xin = in.x;
677 double yin = in.y;
678
679 double xx = 1;
680 double xxm1 = 1; // xx^(ix-1)
681 for (std::size_t ix = 0; ix <= _order; ++ix) {
682 std::size_t k = (ix) * (ix + 1) / 2;
683 // iy = 0
684 dermx[k] = ix * xxm1;
685 dermy[k] = 0;
686 monomials[k] = xx;
687 k += ix + 2;
688 double yy = yin;
689 double yym1 = 1; // yy^(iy-1)
690 for (std::size_t iy = 1; iy <= _order - ix; ++iy) {
691 monomials[k] = xx * yy;
692 dermx[k] = ix * xxm1 * yy;
693 dermy[k] = iy * xx * yym1;
694 yym1 *= yin;
695 yy *= yin;
696 k += ix + iy + 2;
697 }
698 xx *= xin;
699 if (ix >= 1) xxm1 *= xin;
700 }
701
702 // output position
703 double xout = 0, yout = 0;
704 const double *c = &_coeffs[0];
705 const double *pm = &monomials[0];
706 for (int k = _nterms; k--;) xout += (*(pm++)) * (*(c++));
707 pm = &monomials[0];
708 for (int k = _nterms; k--;) yout += (*(pm++)) * (*(c++));
709 res.x = xout;
710 res.y = yout;
711
712 // derivatives
713 c = &_coeffs[0];
714 const double *mx = &dermx[0];
715 const double *my = &dermy[0];
716 double a11 = 0, a12 = 0;
717 for (int k = _nterms; k--;) {
718 a11 += (*(mx++)) * (*c);
719 a12 += (*(my++)) * (*(c++));
720 }
721
722 double a21 = 0, a22 = 0;
723 mx = &dermx[0];
724 my = &dermy[0];
725 for (int k = _nterms; k--;) {
726 a21 += (*(mx++)) * (*c);
727 a22 += (*(my++)) * (*(c++));
728 }
729
730 // output co-variance
731 res.vx = a11 * (a11 * in.vx + 2 * a12 * in.vxy) + a12 * a12 * in.vy;
732 res.vy = a21 * a21 * in.vx + a22 * a22 * in.vy + 2. * a21 * a22 * in.vxy;
733 res.vxy = a21 * a11 * in.vx + a22 * a12 * in.vy + (a21 * a12 + a11 * a22) * in.vxy;
734 out = res;
735}

◆ transformStar()

void lsst::jointcal::AstrometryTransform::transformStar ( FatPoint & in) const
inlineinherited

Definition at line 107 of file AstrometryTransform.h.

107{ transformPosAndErrors(in, in); }
virtual void transformPosAndErrors(const FatPoint &in, FatPoint &out) const

◆ write() [1/2]

void lsst::jointcal::AstrometryTransform::write ( const std::string & fileName) const
inherited

Definition at line 247 of file AstrometryTransform.cc.

247 {
248 ofstream s(fileName.c_str());
249 write(s);
250 bool ok = !s.fail();
251 s.close();
252 if (!ok)
253 throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
254 "AstrometryTransform::write, something went wrong for file " + fileName);
255}
T c_str(T... args)
void write(const std::string &fileName) const

◆ write() [2/2]

void lsst::jointcal::AstrometryTransformPolynomial::write ( std::ostream & s) const
overridevirtual

Reimplemented from lsst::jointcal::AstrometryTransform.

Definition at line 1104 of file AstrometryTransform.cc.

1104 {
1105 s << " AstrometryTransformPolynomial 1" << endl;
1106 s << "order " << _order << endl;
1107 int oldprec = s.precision();
1108 s << setprecision(12);
1109 for (std::size_t k = 0; k < 2 * _nterms; ++k) s << _coeffs[k] << ' ';
1110 s << endl;
1111 s << setprecision(oldprec);
1112}
T setprecision(T... args)

The documentation for this class was generated from the following files: