LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
AstrometryTransform.cc
Go to the documentation of this file.
1// -*- LSST-C++ -*-
2/*
3 * This file is part of jointcal.
4 *
5 * Developed for the LSST Data Management System.
6 * This product includes software developed by the LSST Project
7 * (https://www.lsst.org).
8 * See the COPYRIGHT file at the top-level directory of this distribution
9 * for details of code ownership.
10 *
11 * This program is free software: you can redistribute it and/or modify
12 * it under the terms of the GNU General Public License as published by
13 * the Free Software Foundation, either version 3 of the License, or
14 * (at your option) any later version.
15 *
16 * This program is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 * GNU General Public License for more details.
20 *
21 * You should have received a copy of the GNU General Public License
22 * along with this program. If not, see <https://www.gnu.org/licenses/>.
23 */
24
25#include <cassert>
26#include <cmath>
27#include <fstream>
28#include <iomanip>
29#include <iostream>
30#include <iterator> /* for ostream_iterator */
31#include <limits>
32#include <memory>
33#include <sstream>
34#include <utility>
35
36#include "Eigen/Core"
37
38#include "lsst/log/Log.h"
39#include "lsst/geom/Point.h"
41#include "lsst/jointcal/Frame.h"
43#include "lsst/pex/exceptions.h"
44#include "Eigen/Cholesky"
45
47
48using namespace std;
49
50namespace {
51LOG_LOGGER _log = LOG_GET("lsst.jointcal.AstrometryTransform");
52}
53
54namespace lsst {
55namespace jointcal {
56
57bool isIntegerShift(const AstrometryTransform *transform) {
58 const auto *shift =
59 dynamic_cast<const AstrometryTransformPolynomial *>(transform);
60 if (shift == nullptr) return false;
61
62 static const double eps = 1e-5;
63
64 double dx = shift->getCoefficient(0, 0, 0);
65 double dy = shift->getCoefficient(0, 0, 1);
66
67 static Point dumb(4000, 4000);
68 if (fabs(dx - int(floor(dx + 0.5))) < eps && fabs(dy - int(floor(dy + 0.5))) < eps &&
69 fabs(dumb.x + dx - shift->apply(dumb).x) < eps && fabs(dumb.y + dy - shift->apply(dumb).y) < eps)
70 return true;
71
72 return false;
73}
74
75/********* AstrometryTransform ***********************/
76
77Frame AstrometryTransform::apply(Frame const &inputframe, bool inscribed) const {
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}
94
99
100double AstrometryTransform::getJacobian(const double x, const double y) const {
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}
115
120 const double step) const {
121 double x = where.x;
122 double y = where.y;
123 double xp0, yp0;
124 apply(x, y, xp0, yp0);
125
126 double xp, yp;
127 apply(x + step, y, xp, yp);
128 derivative.a11() = (xp - xp0) / step;
129 derivative.a21() = (yp - yp0) / step;
130 apply(x, y + step, xp, yp);
131 derivative.a12() = (xp - xp0) / step;
132 derivative.a22() = (yp - yp0) / step;
133 derivative.dx() = 0;
134 derivative.dy() = 0;
135}
136
138 const double step) const {
139 Point outwhere = apply(where);
141 computeDerivative(where, der, step);
142 return AstrometryTransformLinearShift(outwhere.x, outwhere.y) * der *
143 AstrometryTransformLinearShift(-where.x, -where.y);
144}
145
147 FatPoint res; // in case in and out are the same address...
148 res = apply(in);
150 // could save a call here, since Derivative needs the transform of where that we already have
151 // 0.01 may not be a very good idea in all cases. May be we should provide a way of altering that.
152 computeDerivative(in, der, 0.01);
153 double a11 = der.A11();
154 double a22 = der.A22();
155 double a21 = der.A21();
156 double a12 = der.A12();
157 res.vx = a11 * (a11 * in.vx + 2 * a12 * in.vxy) + a12 * a12 * in.vy;
158 res.vy = a21 * a21 * in.vx + a22 * a22 * in.vy + 2. * a21 * a22 * in.vxy;
159 res.vxy = a21 * a11 * in.vx + a22 * a12 * in.vy + (a21 * a12 + a11 * a22) * in.vxy;
160 out = res;
161}
162
163void AstrometryTransform::transformErrors(Point const &where, const double *vIn, double *vOut) const {
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}
195
197 // "in" and "out" refer to the inverse direction.
198 Point centerOut = region.getCenter();
199 Point centerIn = apply(centerOut);
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);
206}
207
208/* implement one in AstrometryTransform, so that all derived
209 classes do not need to provide one... */
210
211/* the routines that follow are used for ea generic parameter
212 transformation serialization, used e.g. for fits. Enables
213 to manipulate transformation parameters as vectors.
214*/
215
216// not dummy : what it does is virtual because paramRef is virtual.
217void AstrometryTransform::getParams(double *params) const {
218 std::size_t npar = getNpar();
219 for (std::size_t i = 0; i < npar; ++i) params[i] = paramRef(i);
220}
221
222void AstrometryTransform::offsetParams(Eigen::VectorXd const &delta) {
223 std::size_t npar = getNpar();
224 for (std::size_t i = 0; i < npar; ++i) paramRef(i) += delta[i];
225}
226
227double AstrometryTransform::paramRef(Eigen::Index const) const {
229 std::string("AstrometryTransform::paramRef should never be called "));
230}
231
232double &AstrometryTransform::paramRef(Eigen::Index const) {
234 "AstrometryTransform::paramRef should never be called ");
235}
236
237void AstrometryTransform::paramDerivatives(Point const &, double *, double *) const {
239 "AstrometryTransform::paramDerivatives() should never be called ");
240}
241
242ostream &operator<<(ostream &stream, AstrometryTransform const &transform) {
243 transform.print(stream);
244 return stream;
245}
246
247void AstrometryTransform::write(const std::string &fileName) const {
248 ofstream s(fileName.c_str());
249 write(s);
250 bool ok = !s.fail();
251 s.close();
252 if (!ok)
254 "AstrometryTransform::write, something went wrong for file " + fileName);
255}
256
258 throw LSST_EXCEPT(
260 "AstrometryTransform::write(ostream), should never be called. MEans that it is missing in some "
261 "derived class ");
262}
263
264/******************* GTransformInverse ****************/
265/* inverse transformation, solved by iterations. Before using
266 it (probably via AstrometryTransform::inverseTransform), consider
267 seriously StarMatchList::inverseTransform */
269private:
272 double precision2;
273
274public:
275 AstrometryTransformInverse(const AstrometryTransform *direct, double precision,
276 const Frame &region);
277
280 void apply(double xIn, double yIn, double &xOut, double &yOut) const;
281
282 void print(ostream &out) const;
283
284 double fit(StarMatchList const &starMatchList);
285
287
289
291 std::unique_ptr<AstrometryTransform> roughInverse(const Frame &) const { return _direct->clone(); }
292
295 return _direct->clone();
296 }
297
299
300private:
301 void operator=(AstrometryTransformInverse const &);
302};
303
305 const Frame &region) const {
306 return std::unique_ptr<AstrometryTransform>(new AstrometryTransformInverse(this, precision, region));
307}
308
310 const double precision, const Frame &region) {
311 _direct = direct->clone();
312 _roughInverse = _direct->roughInverse(region);
313 precision2 = precision * precision;
314}
315
318 _direct = model._direct->clone();
319 _roughInverse = model._roughInverse->clone();
320 precision2 = model.precision2;
321}
322
324
325void AstrometryTransformInverse::operator=(AstrometryTransformInverse const &model) {
326 _direct = model._direct->clone();
327 _roughInverse = model._roughInverse->clone();
328 precision2 = model.precision2;
329}
330
331void AstrometryTransformInverse::apply(const double xIn, const double yIn, double &xOut, double &yOut) const {
332 Point in(xIn, yIn);
333 Point outGuess = _roughInverse->apply(in);
334 AstrometryTransformLinear directDer, reverseDer;
335 int loop = 0;
336 int maxloop = 20;
337 double move2;
338 do {
339 loop++;
340 Point inGuess = _direct->apply(outGuess);
341 _direct->computeDerivative(outGuess, directDer);
342 reverseDer = directDer.inverted();
343 double xShift, yShift;
344 reverseDer.apply(xIn - inGuess.x, yIn - inGuess.y, xShift, yShift);
345 outGuess.x += xShift;
346 outGuess.y += yShift;
347 move2 = xShift * xShift + yShift * yShift;
348 } while ((move2 > precision2) && (loop < maxloop));
349 if (loop == maxloop) LOGLS_WARN(_log, "Problems applying AstrometryTransformInverse at " << in);
350 xOut = outGuess.x;
351 yOut = outGuess.y;
352}
353
355 stream << " AstrometryTransformInverse of :" << endl << *_direct << endl;
356}
357
360 "Cannot fit a AstrometryTransformInverse. Use StarMatchList::inverseTransform instead.");
361}
362
366
367/************* AstrometryTransformComposition **************/
368
369// This class was done to allow composition of AstrometryTransform's, without specifications of their types.
370// does not need to be public. Invoked by compose(left,right)
371
375private:
377
378public:
380
382 void apply(double xIn, double yIn, double &xOut, double &yOut) const;
383 void print(ostream &stream) const;
384
386 double fit(StarMatchList const &starMatchList);
387
390};
391
393 AstrometryTransform const &first) {
394 _first = first.clone();
395 _second = second.clone();
396}
397
398void AstrometryTransformComposition::apply(const double xIn, const double yIn, double &xOut,
399 double &yOut) const {
400 double xout, yout;
401 _first->apply(xIn, yIn, xout, yout);
402 _second->apply(xout, yout, xOut, yOut);
403}
404
406 stream << "Composed AstrometryTransform consisting of:" << std::endl;
407 _first->print(stream);
408 _second->print(stream);
409}
410
412 /* fits only one of them. could check that first can actually be fitted... */
413 return _first->fit(starMatchList);
414}
415
417 return std::make_unique<AstrometryTransformComposition>(*_second, *_first);
418}
419
421
426
428 AstrometryTransform const &right) {
429 // Try to use the composeAndReduce method from left. If absent, AstrometryTransform::composeAndReduce
430 // returns NULL. composeAndReduce is non trivial for polynomials.
431 std::unique_ptr<AstrometryTransform> composition(left.composeAndReduce(right));
432 // composition == NULL means no reduction: just build a Composition that pipelines "left" and "right".
433 if (composition == nullptr)
434 return std::make_unique<AstrometryTransformComposition>(left, right);
435 else
436 return composition;
437}
438
439// just a speed up, to avoid useless numerical derivation.
441 const double) const {
442 derivative = AstrometryTransformLinear();
443}
444
446 const double) const {
448 return result; // rely on default AstrometryTransformlin constructor;
449}
450
452 return std::make_shared<ast::UnitMap>(2); // a AstrometryTransformIdentity is identically ast::UnitMap(2)
453}
454
456 stream << "AstrometryTransformIdentity 1" << endl;
457}
458
460 int format;
461 stream >> format;
462 if (format != 1)
464 " AstrometryTransformIdentity::read : format is not 1 ");
465}
466
467/*************** AstrometryTransformPolynomial **************************************/
468
470
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}
482
483//#ifdef TO_BE_FIXED
485 const Frame &frame, std::size_t order,
486 std::size_t nPoint) {
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}
504//#endif
505
508 std::size_t order, std::size_t nSteps) {
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}
528
529void AstrometryTransformPolynomial::computeMonomials(double xIn, double yIn, double *monomial) const {
530 /* The ordering of monomials is implemented here.
531 You may not change it without updating the "mapping" routines
532 getCoefficient(std::size_t, std::size_t, std::size_t).
533 I (P.A.) did not find a clever way to loop over monomials.
534 Improvements welcome.
535 This routine is used also by the fit to fill monomials.
536 We could certainly be more elegant.
537 */
538
539 double xx = 1;
540 for (std::size_t ix = 0; ix <= _order; ++ix) {
541 double yy = 1;
542 std::size_t k = ix * (ix + 1) / 2;
543 for (std::size_t iy = 0; iy <= _order - ix; ++iy) {
544 monomial[k] = xx * yy;
545 yy *= yIn;
546 k += ix + iy + 2;
547 }
548 xx *= xIn;
549 }
550}
551
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}
571
572/* this is reasonably fast, when optimized */
573void AstrometryTransformPolynomial::apply(const double xIn, const double yIn, double &xOut,
574 double &yOut) const {
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}
596
598 AstrometryTransformLinear &derivative,
599 const double step)
600 const { /* 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}
656
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}
736
737/* The coefficient ordering is defined both here *AND* in the
738 AstrometryTransformPolynomial::apply, AstrometryTransformPolynomial::Derivative, ... routines
739 Change all or none ! */
740
742 std::size_t whichCoord) const {
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}
749
751 std::size_t whichCoord) {
752 assert((degX + degY <= _order) && whichCoord < 2);
753 return _coeffs[(degX + degY) * (degX + degY + 1) / 2 + degY + whichCoord * _nterms];
754}
755
757 std::size_t whichCoord) const {
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}
764
765/* parameter serialization for "virtual" fits */
766double AstrometryTransformPolynomial::paramRef(Eigen::Index const i) const {
767 assert(i < 2 * Eigen::Index(_nterms));
768 return _coeffs[i];
769}
770
771double &AstrometryTransformPolynomial::paramRef(Eigen::Index const i) {
772 assert(i < 2 * Eigen::Index(_nterms));
773 return _coeffs[i];
774}
775
776void AstrometryTransformPolynomial::paramDerivatives(Point const &where, double *dx, double *dy)
777 const { /* 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}
784
785/* utility for the print(ostream&) routine */
786static string monomialString(std::size_t powX, std::size_t powY) {
787 stringstream ss;
788 if (powX + powY) ss << "*";
789 if (powX > 0) ss << "x";
790 if (powX > 1) ss << "^" << powX;
791 if (powY > 0) ss << "y";
792 if (powY > 1) ss << "^" << powY;
793 return ss.str();
794}
795
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}
825
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}
832
835 Point center = frame.getCenter();
836 return AstrometryTransformLinearScale(2. / frame.getWidth(), 2. / frame.getHeight()) *
837 AstrometryTransformLinearShift(-center.x, -center.y);
838}
839
840/*utility for the AstrometryTransformPolynomial::fit() routine */
841static AstrometryTransformLinear shiftAndNormalize(StarMatchList const &starMatchList) {
842 double xav = 0;
843 double x2 = 0;
844 double yav = 0;
845 double y2 = 0;
846 double count = 0;
847 for (const auto & a_match : starMatchList) {
848 Point const &point1 = a_match.point1;
849 xav += point1.x;
850 yav += point1.y;
851 x2 += std::pow(point1.x, 2);
852 y2 += std::pow(point1.y, 2);
853 count++;
854 }
855 if (count == 0) return AstrometryTransformLinear();
856 xav /= count;
857 yav /= count;
858 // 3.5 stands for sqrt(12).
859 double xspan = 3.5 * std::sqrt(x2 / count - std::pow(xav, 2));
860 double yspan = 3.5 * std::sqrt(y2 / count - std::pow(yav, 2));
861 return AstrometryTransformLinearScale(2. / xspan, 2. / yspan) *
862 AstrometryTransformLinearShift(-xav, -yav);
863}
864
865static double sq(double x) { return x * x; }
866
867double AstrometryTransformPolynomial::computeFit(StarMatchList const &starMatchList,
868 AstrometryTransform const &shiftToCenter,
869 const bool useErrors) {
870 Eigen::MatrixXd A(2 * _nterms, 2 * _nterms);
871 A.setZero();
872 Eigen::VectorXd B(2 * _nterms);
873 B.setZero();
874 double sumr2 = 0;
875 double monomials[_nterms];
876 for (const auto & a_match : starMatchList) {
877 Point tmp = shiftToCenter.apply(a_match.point1);
878 FatPoint point1(tmp, a_match.point1.vx, a_match.point1.vy, a_match.point1.vxy);
879 FatPoint const &point2 = a_match.point2;
880 double wxx, wyy, wxy;
881 FatPoint tr1;
882 computeMonomials(point1.x, point1.y, monomials);
883 if (useErrors) {
884 transformPosAndErrors(point1, tr1); // we might consider recycling the monomials
885 double vxx = (tr1.vx + point2.vx);
886 double vyy = (tr1.vy + point2.vy);
887 double vxy = (tr1.vxy + point2.vxy);
888 double det = vxx * vyy - vxy * vxy;
889 wxx = vyy / det;
890 wyy = vxx / det;
891 wxy = -vxy / det;
892 } else {
893 wxx = wyy = 1;
894 wxy = 0;
895 apply(point1.x, point1.y, tr1.x, tr1.y);
896 }
897 double resx = point2.x - tr1.x;
898 double resy = point2.y - tr1.y;
899 sumr2 += wxx * sq(resx) + wyy * sq(resy) + 2 * wxy * resx * resy;
900
901 double bxcoeff = wxx * resx + wxy * resy;
902 double bycoeff = wyy * resy + wxy * resx;
903 for (std::size_t j = 0; j < _nterms; ++j) {
904 for (std::size_t i = j; i < _nterms; ++i) {
905 A(i, j) += wxx * monomials[i] * monomials[j];
906 A(i + _nterms, j + _nterms) += wyy * monomials[i] * monomials[j];
907 A(j, i + _nterms) = A(i, j + _nterms) += wxy * monomials[i] * monomials[j];
908 }
909 B(j) += bxcoeff * monomials[j];
910 B(j + _nterms) += bycoeff * monomials[j];
911 }
912 } // end loop on points
913 Eigen::LDLT<Eigen::MatrixXd, Eigen::Lower> factor(A);
914 // should probably throw
915 if (factor.info() != Eigen::Success) {
916 LOGL_ERROR(_log, "AstrometryTransformPolynomial::fit could not factorize");
917 return -1;
918 }
919
920 Eigen::VectorXd sol = factor.solve(B);
921 for (std::size_t k = 0; k < 2 * _nterms; ++k) _coeffs[k] += sol(k);
922 if (starMatchList.size() == _nterms) return 0;
923 return (sumr2 - B.dot(sol));
924}
925
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}
943
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}
951
952/* PolyXY the class used to perform polynomial algebra (and in
953 particular composition) at the coefficient level. This class
954 handles a single polynomial, while a AstrometryTransformPolynomial is a couple of
955 polynomials. This class does not have any routine to evaluate
956 polynomials. Efficiency is not a concern since these routines are
957 seldom used. There is no need to expose this tool class to
958 AstrometryTransform users.
959*/
960
961class PolyXY {
962 std::size_t order;
963 std::size_t nterms;
964 vector<long double> coeffs;
965
966public:
967 PolyXY(const int order) : order(order), nterms((order + 1) * (order + 2) / 2) {
968 coeffs.reserve(nterms);
969 coeffs.insert(coeffs.begin(), nterms, 0L); // fill & initialize to 0.
970 }
971
972 std::size_t getOrder() const { return order; }
973
975 : order(transform.getOrder()), nterms((order + 1) * (order + 2) / 2), coeffs(nterms, 0L) {
976 for (std::size_t px = 0; px <= order; ++px)
977 for (std::size_t py = 0; py <= order - px; ++py) {
978 getCoefficient(px, py) = transform.getCoefficient(px, py, whichCoord);
979 }
980 }
981
982 long double getCoefficient(std::size_t powX, std::size_t powY) const {
983 assert(powX + powY <= order);
984 return coeffs.at((powX + powY) * (powX + powY + 1) / 2 + powY);
985 }
986
987 long double &getCoefficient(std::size_t powX, std::size_t powY) {
988 assert(powX + powY <= order);
989 return coeffs.at((powX + powY) * (powX + powY + 1) / 2 + powY);
990 }
991};
992
993/* ===================== PolyXY Algebra routines ================== */
994
995static void operator+=(PolyXY &left, const PolyXY &right) {
996 std::size_t rdeg = right.getOrder();
997 assert(left.getOrder() >= rdeg);
998 for (std::size_t i = 0; i <= rdeg; ++i)
999 for (std::size_t j = 0; j <= rdeg - i; ++j) left.getCoefficient(i, j) += right.getCoefficient(i, j);
1000}
1001
1002/* multiplication by a scalar */
1003static PolyXY operator*(const long double &a, const PolyXY &polyXY) {
1004 PolyXY result(polyXY);
1005 // no direct access to coefficients: do it the soft way
1006 std::size_t order = polyXY.getOrder();
1007 for (std::size_t i = 0; i <= order; ++i)
1008 for (std::size_t j = 0; j <= order - i; ++j) result.getCoefficient(i, j) *= a;
1009 return result;
1010}
1011
1013static PolyXY product(const PolyXY &p1, const PolyXY &p2) {
1014 std::size_t deg1 = p1.getOrder();
1015 std::size_t deg2 = p2.getOrder();
1016 PolyXY result(deg1 + deg2);
1017 for (std::size_t i1 = 0; i1 <= deg1; ++i1)
1018 for (std::size_t j1 = 0; j1 <= deg1 - i1; ++j1)
1019 for (std::size_t i2 = 0; i2 <= deg2; ++i2)
1020 for (std::size_t j2 = 0; j2 <= deg2 - i2; ++j2)
1021 result.getCoefficient(i1 + i2, j1 + j2) +=
1022 p1.getCoefficient(i1, j1) * p2.getCoefficient(i2, j2);
1023 return result;
1024}
1025
1026/* powers[k](x,y) = polyXY(x,y)**k, 0 <= k <= maxP */
1027static void computePowers(const PolyXY &polyXY, std::size_t maxP, vector<PolyXY> &powers) {
1028 powers.reserve(maxP + 1);
1029 powers.emplace_back(0);
1030 powers[0].getCoefficient(0, 0) = 1L;
1031 for (std::size_t k = 1; k <= maxP; ++k) powers.push_back(product(powers[k - 1], polyXY));
1032}
1033
1035static PolyXY composition(const PolyXY &polyXY, const PolyXY &polyX, const PolyXY &polyY) {
1036 std::size_t pdeg = polyXY.getOrder();
1037 PolyXY result(pdeg * max(polyX.getOrder(), polyY.getOrder()));
1038 vector<PolyXY> pXPowers;
1039 vector<PolyXY> pYPowers;
1040 computePowers(polyX, pdeg, pXPowers);
1041 computePowers(polyY, pdeg, pYPowers);
1042 for (std::size_t px = 0; px <= pdeg; ++px)
1043 for (std::size_t py = 0; py <= pdeg - px; ++py)
1044 result += polyXY.getCoefficient(px, py) * product(pXPowers.at(px), pYPowers.at(py));
1045 return result;
1046}
1047
1048/* ===================== end of PolyXY Algebra routines ============= */
1049
1050/* reducing polynomial composition is the reason for PolyXY stuff : */
1051
1053 AstrometryTransformPolynomial const &right) const {
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}
1073
1075 AstrometryTransformPolynomial const &right) const {
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}
1087
1089 AstrometryTransformPolynomial const &right) const {
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}
1098
1100 auto inverse = inversePolyTransform(*this, domain, 1e-7, _order + 2, 100);
1101 return std::make_shared<ast::PolyMap>(toAstPolyMapCoefficients(), inverse->toAstPolyMapCoefficients());
1102}
1103
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}
1113
1115 int format;
1116 s >> format;
1117 if (format != 1)
1119 " AstrometryTransformPolynomial::read : format is not 1 ");
1120
1121 string order;
1122 s >> order >> _order;
1123 if (order != "order")
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}
1129
1130ndarray::Array<double, 2, 2> AstrometryTransformPolynomial::toAstPolyMapCoefficients() const {
1131 std::size_t nCoeffs = _coeffs.size();
1132 ndarray::Array<double, 2, 2> result = ndarray::allocate(ndarray::makeVector(nCoeffs, std::size_t(4)));
1133
1134 ndarray::Size k = 0;
1135 for (std::size_t iCoord = 0; iCoord < 2; ++iCoord) {
1136 for (std::size_t p = 0; p <= _order; ++p) {
1137 for (std::size_t py = 0; py <= p; ++py, ++k) {
1138 result[k][0] = getCoefficient(p - py, py, iCoord);
1139 result[k][1] = iCoord + 1;
1140 result[k][2] = p - py;
1141 result[k][3] = py;
1142 }
1143 }
1144 }
1145
1146 return result;
1147}
1148
1150 Frame const &domain,
1151 double const precision,
1152 std::size_t maxOrder,
1153 std::size_t nSteps) {
1154 StarMatchList sm;
1155 double xStart = domain.xMin;
1156 double yStart = domain.yMin;
1157 double xStep = domain.getWidth() / (nSteps - 1);
1158 double yStep = domain.getHeight() / (nSteps - 1);
1159 for (std::size_t i = 0; i < nSteps; ++i) {
1160 for (std::size_t j = 0; j < nSteps; ++j) {
1161 Point in(xStart + i * xStep, yStart + j * yStep);
1162 Point out(forward.apply(in));
1163 sm.push_back(StarMatch(out, in, nullptr, nullptr));
1164 }
1165 }
1166 std::size_t npairs = sm.size();
1170 double chi2 = 0;
1171 double oldChi2 = std::numeric_limits<double>::infinity();
1172 for (order = 1; order <= maxOrder; ++order) {
1174 auto success = poly->fit(sm);
1175 if (success == -1) {
1176 std::stringstream errMsg;
1177 errMsg << "Cannot fit a polynomial of order " << order << " with " << nSteps << "^2 points";
1178 throw pexExcept::RuntimeError(errMsg.str());
1179 }
1180 // compute the chi2 ignoring errors:
1181 chi2 = 0;
1182 for (auto const &i : sm) chi2 += i.point2.computeDist2(poly->apply((i.point1)));
1183 LOGLS_TRACE(_log, "inversePoly order " << order << ": " << chi2 << " / " << npairs << " = "
1184 << chi2 / npairs << " < " << precision * precision);
1185
1186 if (chi2 / npairs < precision * precision) break;
1187
1188 // If this triggers, we know we did not reach the required precision.
1189 if (chi2 > oldChi2) {
1190 LOGLS_WARN(_log, "inversePolyTransform: chi2 increases ("
1191 << chi2 << " > " << oldChi2 << "); ending fit with order: " << order);
1192 LOGLS_WARN(_log, "inversePolyTransform: requested precision not reached: "
1193 << chi2 << " / " << npairs << " = " << chi2 / npairs << " < "
1194 << precision * precision);
1195 poly = std::move(oldPoly);
1196 order--;
1197 break;
1198 } else {
1199 oldChi2 = chi2;
1200 // Clone it so we don't lose it in the next iteration.
1201 oldPoly = dynamic_pointer_cast<AstrometryTransformPolynomial>(
1203 }
1204 }
1205 if (order > maxOrder)
1206 LOGLS_WARN(_log, "inversePolyTransform: Reached max order without reaching requested precision: "
1207 << chi2 << " / " << npairs << " = " << chi2 / npairs << " < "
1208 << precision * precision);
1209 return poly;
1210}
1211
1212/**************** AstrometryTransformLinear ***************************************/
1213/* AstrometryTransformLinear is a specialized constructor of AstrometryTransformPolynomial
1214 May be it could just disappear ??
1215*/
1216
1217AstrometryTransformLinear::AstrometryTransformLinear(const double Dx, const double Dy, const double A11,
1218 const double A12, const double A21, const double A22)
1220 dx() = Dx;
1221 a11() = A11;
1222 a12() = A12;
1223 dy() = Dy;
1224 a21() = A21;
1225 a22() = A22;
1226}
1227
1230 if (transform.getOrder() != 1)
1232 "Trying to build a AstrometryTransformLinear from a higher order transform. Aborting. ");
1234}
1235
1237 // There is a general routine in AstrometryTransformPolynomial that would do the job:
1238 // return AstrometryTransformLinear(AstrometryTransformPolynomial::operator*(right));
1239 // however, we are using this composition of linear stuff heavily in
1240 // AstrometryTransform::linearApproximation, itself used in inverseTransform::apply.
1241 // So, for sake of efficiency, and since it is easy, we take a shortcut:
1243 apply(right.Dx(), right.Dy(), result.dx(), result.dy());
1244 result.a11() = A11() * right.A11() + A12() * right.A21();
1245 result.a12() = A11() * right.A12() + A12() * right.A22();
1246 result.a21() = A21() * right.A11() + A22() * right.A21();
1247 result.a22() = A21() * right.A12() + A22() * right.A22();
1248 return result;
1249}
1250
1252 const double) const {
1253 derivative = *this;
1254 derivative.getCoefficient(0, 0, 0) = 0;
1255 derivative.getCoefficient(0, 0, 1) = 0;
1256}
1257
1259 return *this;
1260}
1261
1263 //
1264 // (T1,M1) * (T2,M2) = 1 i.e (0,1) implies
1265 // T1 = -M1 * T2
1266 // M1 = M2^(-1)
1267 //
1268
1269 double a11 = A11();
1270 double a12 = A12();
1271 double a21 = A21();
1272 double a22 = A22();
1273 double d = (a11 * a22 - a12 * a21);
1274 if (d == 0) {
1275 LOGL_FATAL(_log,
1276 "AstrometryTransformLinear::invert singular transformation: transform contents will be "
1277 "dumped to stderr.");
1278 print(cerr);
1279 }
1280
1281 AstrometryTransformLinear result(0, 0, a22 / d, -a12 / d, -a21 / d, a11 / d);
1282 double rdx, rdy;
1283 result.apply(Dx(), Dy(), rdx, rdy);
1284 result.dx() = -rdx;
1285 result.dy() = -rdy;
1286 return result;
1287}
1288
1293
1295 auto oldPrecision = stream.precision();
1296 stream.precision(12);
1297 stream << "Linear AstrometryTransform:" << std::endl;
1298 stream << A11() << " " << A12() << " + " << Dx() << std::endl;
1299 stream << A21() << " " << A22() << " + " << Dy() << std::endl;
1300 stream.precision(oldPrecision);
1301 stream << "determinant = " << determinant();
1302 stream.precision(oldPrecision);
1303}
1304
1306 throw pexExcept::NotFoundError("AstrometryTransformLinearRot::fit not implemented! aborting");
1307}
1308
1310 std::size_t npairs = starMatchList.size();
1311 if (npairs < 3) {
1312 LOGLS_FATAL(_log, "AstrometryTransformLinearShift::fit trying to fit a linear transform with only "
1313 << npairs << " matches.");
1314 return -1;
1315 }
1316
1317 double sumr2 = 0; /* used to compute chi2 without relooping */
1318 /* loop on pairs and fill */
1319 Eigen::VectorXd B(2);
1320 B.setZero();
1321 Eigen::MatrixXd A(2, 2);
1322 A.setZero();
1323
1324 for (auto const &it : starMatchList) {
1325 FatPoint const &point1 = it.point1;
1326 FatPoint const &point2 = it.point2;
1327 double deltax = point2.x - point1.x;
1328 double deltay = point2.y - point1.y;
1329 double vxx = point1.vx + point2.vx;
1330 double vyy = point1.vy + point2.vy;
1331 double vxy = point1.vxy + point2.vxy;
1332 double det = vxx * vyy - vxy * vxy;
1333 double wxx = vyy / det;
1334 double wyy = vxx / det;
1335 double wxy = -vxy / det;
1336 B(0) += deltax * wxx + wxy * deltay;
1337 B(1) += deltay * wyy + wxy * deltax;
1338 A(0, 0) += wxx;
1339 A(1, 1) += wyy;
1340 A(0, 1) += wxy;
1341 sumr2 += deltax * deltax * wxx + deltay * deltay * wyy + 2. * wxy * deltax * deltay;
1342 }
1343 double det = A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
1344 if (det <= 0) return -1;
1345 double tmp = A(0, 0);
1346 A(0, 0) = A(1, 1) / det;
1347 A(1, 1) = tmp / det;
1348 A(0, 1) = A(1, 0) = -A(0, 1) / det;
1349 Eigen::VectorXd sol = A * B;
1350 (*this) = AstrometryTransformLinearShift(sol(0), sol(1));
1351 return (sumr2 - sol.dot(B)); // chi2 compact form
1352}
1353
1355 const double scaleFactor) {
1356 double c = scaleFactor * std::cos(angleRad);
1357 double s = scaleFactor * std::sin(angleRad);
1358 a11() = a22() = c;
1359 a21() = s;
1360 a12() = -s;
1361
1362 // we want that the center does not move : transform+M*C = C ==> transform = C - M*C
1363 Point a_point(0., 0.);
1364 if (center) a_point = *center;
1365
1366 dx() = dy() = 0;
1367 AstrometryTransformPolynomial::apply(a_point.x, a_point.y, dx(), dy()); // compute M*C
1368 dx() = a_point.x - Dx();
1369 dy() = a_point.y - dy();
1370}
1371
1372static double deg2rad(double degree) { return degree * M_PI / 180.; }
1373
1374static double rad2deg(double rad) { return rad * 180. / M_PI; }
1375
1376/************* WCS transform ******************/
1377/************** LinPixelToTan *******************/
1378
1379/* Implementation note : it seemed wise to incorporate
1380 the radians to degreess convertion into the linPixelToTan
1381 part (right in the constructor), and to do the
1382 opposite operation in the LinPart routine.
1383 When I was coding the fit, I realized that it was a
1384 bad idea. then I realized that the fitting routine
1385 itself was probably useless. Finaly, for a persistor,
1386 it seems bizarre that the stored data is different
1387 from what convention (angles in degrees for WCS's)
1388 would expect.
1389 So, no more "automatic" degrees to radians and
1390 radians to degrees conversion. They are explicitely
1391 done in apply (for TanPixelToRaDec and TanRaDecToPixel).
1392 This is a minor concern though....
1393*/
1394BaseTanWcs::BaseTanWcs(AstrometryTransformLinear const &pixToTan, Point const &tangentPoint,
1395 const AstrometryTransformPolynomial *corrections) {
1396 // the angles returned by linPixelToTan should be in degrees.
1397 linPixelToTan = pixToTan;
1398 ra0 = deg2rad(tangentPoint.x);
1399 dec0 = deg2rad(tangentPoint.y);
1400 cos0 = std::cos(dec0);
1401 sin0 = std::sin(dec0);
1402 corr = nullptr;
1403 if (corrections) corr = std::make_unique<AstrometryTransformPolynomial>(*corrections);
1404}
1405
1406/* with some sort of smart pointer ro handle "corr", we could remove the
1407 copy constructor, the operator = and the destructor */
1408
1409// ": AstrometryTransform" suppresses a warning
1411 corr = nullptr;
1412 *this = original;
1413}
1414
1416 linPixelToTan = original.linPixelToTan;
1417 ra0 = original.ra0;
1418 dec0 = original.dec0;
1419 cos0 = std::cos(dec0);
1420 sin0 = std::sin(dec0);
1421 corr = nullptr;
1422 if (original.corr) corr = std::make_unique<AstrometryTransformPolynomial>(*original.corr);
1423 return *this;
1424}
1425
1426void BaseTanWcs::apply(const double xIn, const double yIn, double &xOut, double &yOut) const {
1427 double l, m; // radians in the tangent plane
1428 pixToTangentPlane(xIn, yIn, l, m); // l, m in degrees.
1429 l = deg2rad(l);
1430 m = deg2rad(m); // now in radians
1431 // Code inspired from worldpos.c in wcssubs (ancestor of the wcslib)
1432 /* At variance with wcslib, it collapses the projection to a plane
1433 and expression of sidereal cooordinates into a single set of
1434 operations. */
1435 double dect = cos0 - m * sin0;
1436 if (dect == 0) {
1437 LOGL_WARN(_log, "No sidereal coordinates at pole!");
1438 xOut = 0;
1439 yOut = 0;
1440 return;
1441 }
1442 double rat = ra0 + atan2(l, dect);
1443 dect = atan(std::cos(rat - ra0) * (m * cos0 + sin0) / dect);
1444 if (rat - ra0 > M_PI) rat -= (2. * M_PI);
1445 if (rat - ra0 < -M_PI) rat += (2. * M_PI);
1446 if (rat < 0.0) rat += (2. * M_PI);
1447 // convert to degree
1448 xOut = rad2deg(rat);
1449 yOut = rad2deg(dect);
1450}
1451
1452Point BaseTanWcs::getTangentPoint() const { return Point(rad2deg(ra0), rad2deg(dec0)); }
1453
1455
1459
1461 /* CRPIX's are defined by:
1462 ( CD1_1 CD1_2 ) (x - crpix1)
1463 transformed = ( ) * ( )
1464 ( CD2_1 CD2_2 ) (y - crpix2)
1465
1466 so that CrPix is the point which transforms to (0,0)
1467 */
1469 return Point(inverse.Dx(), inverse.Dy());
1470}
1471
1472BaseTanWcs::~BaseTanWcs() = default;
1473
1476
1477void AstrometryTransformSkyWcs::apply(const double xIn, const double yIn, double &xOut, double &yOut) const {
1478 auto const outCoord = _skyWcs->pixelToSky(geom::Point2D(xIn, yIn));
1479 xOut = outCoord[0].asDegrees();
1480 yOut = outCoord[1].asDegrees();
1481}
1482
1484 out << "AstrometryTransformSkyWcs(" << *_skyWcs << ")";
1485}
1486
1487double AstrometryTransformSkyWcs::fit(const StarMatchList &starMatchList) {
1488 throw LSST_EXCEPT(pex::exceptions::LogicError, "Not implemented");
1489}
1490
1492 return std::make_unique<AstrometryTransformSkyWcs>(getSkyWcs());
1493}
1494
1495/*************************** TanPixelToRaDec ***************/
1496
1498 const AstrometryTransformPolynomial *corrections)
1499 : BaseTanWcs(pixToTan, tangentPoint, corrections) {}
1500
1501// ": AstrometryTransform" suppresses a warning
1503
1505 AstrometryTransformLinear const &right) const {
1506 if (right.getOrder() == 1) {
1507 return std::make_unique<TanPixelToRaDec>((*this) * (right));
1508 } else {
1510 }
1511}
1512
1514 TanPixelToRaDec result(*this);
1515 result.linPixelToTan = result.linPixelToTan * right;
1516 return result;
1517}
1518
1520 if (corr != nullptr) {
1521 LOGL_WARN(_log, "You are inverting a TanPixelToRaDec with corrections.");
1522 LOGL_WARN(_log, "The inverse you get ignores the corrections!");
1523 }
1525}
1526
1531
1533 const Frame &region) const {
1534 if (!corr)
1537 else
1538 return std::unique_ptr<AstrometryTransform>(new AstrometryTransformInverse(this, precision, region));
1539}
1540
1542 if (corr)
1543 return (*corr) * linPixelToTan;
1544 else
1545 return linPixelToTan;
1546}
1547
1548void TanPixelToRaDec::pixToTangentPlane(double xPixel, double yPixel, double &xTangentPlane,
1549 double &yTangentPlane) const {
1550 // xTangentPlane, yTangentPlane in degrees.
1551 linPixelToTan.apply(xPixel, yPixel, xTangentPlane, yTangentPlane);
1552 if (corr) {
1553 double xtmp = xTangentPlane;
1554 double ytmp = yTangentPlane;
1555 corr->apply(xtmp, ytmp, xTangentPlane, yTangentPlane); // still in degrees.
1556 }
1557}
1558
1563
1564void TanPixelToRaDec::print(ostream &stream) const {
1565 stream << " TanPixelToRaDec, lin part :" << endl << linPixelToTan;
1566 Point tp = getTangentPoint();
1567 stream << " tangent point " << tp.x << ' ' << tp.y << endl;
1568 Point crpix = getCrPix();
1569 stream << " crpix " << crpix.x << ' ' << crpix.y << endl;
1570 if (corr) stream << "PV correction: " << endl << *corr;
1571}
1572
1574 /* OK we could implement this routine, but it is
1575 probably useless since to do the match, we have to
1576 project from sky to tangent plane. When a match is
1577 found, it is easier to carry out the fit in the
1578 tangent plane, rather than going back to the celestial
1579 sphere (and reproject to fit...). Anyway if this
1580 message shows up, we'll think about it.
1581 */
1583 "TanPixelToRaDec::fit is NOT implemented (although it is doable)) ");
1584 return -1;
1585}
1586
1587/*************************** TanSipPixelToRaDec ***************/
1588
1590 const AstrometryTransformPolynomial *corrections)
1591 : BaseTanWcs(pixToTan, tangentPoint, corrections) {}
1592
1593// ": AstrometryTransform" suppresses a warning
1595
1596/* Would require some checks before cooking up something more efficient
1597 than just a linear approximation */
1598#if 0
1600{
1601 if (&region) {}
1603}
1604#endif
1605
1607 const Frame &region)
1608 const { /* We have not implemented (yet) the reverse corrections available in SIP */
1609 return std::unique_ptr<AstrometryTransform>(new AstrometryTransformInverse(this, precision, region));
1610}
1611
1618
1619void TanSipPixelToRaDec::pixToTangentPlane(double xPixel, double yPixel, double &xTangentPlane,
1620 double &yTangentPlane) const {
1621 // xTangentPlane, yTangentPlane returned in degrees
1622 if (corr) {
1623 double xtmp, ytmp;
1624 corr->apply(xPixel, yPixel, xtmp, ytmp);
1625 linPixelToTan.apply(xtmp, ytmp, xTangentPlane, yTangentPlane);
1626 } else
1627 linPixelToTan.apply(xPixel, yPixel, xTangentPlane, yTangentPlane);
1628}
1629
1634
1636 stream << " TanSipPixelToRaDec, lin part :" << endl << linPixelToTan;
1637 Point tp = getTangentPoint();
1638 stream << " tangent point " << tp.x << ' ' << tp.y << endl;
1639 Point crpix = getCrPix();
1640 stream << " crpix " << crpix.x << ' ' << crpix.y << endl;
1641 if (corr) stream << "PV correction: " << endl << *corr;
1642}
1643
1645 /* OK we could implement this routine, but it is
1646 probably useless since to do the match, we have to
1647 project from sky to tangent plane. When a match is
1648 found, it is easier to carry out the fit in the
1649 tangent plane, rather than going back to the celestial
1650 sphere (and reproject to fit...). Anyway if this
1651 message shows up, we'll think about it.
1652 */
1654 "TanSipPixelToRaDec::fit is NOT implemented (although it is doable)) ");
1655 return -1;
1656}
1657
1658/*************** reverse transform of TanPixelToRaDec: TanRaDecToPixel ********/
1659
1661 : linTan2Pix(std::move(tan2Pix)) {
1662 setTangentPoint(tangentPoint);
1663}
1664
1665void TanRaDecToPixel::setTangentPoint(Point const &tangentPoint) {
1666 /* the radian to degrees conversion after projection
1667 is handled in apply */
1668 ra0 = deg2rad(tangentPoint.x);
1669 dec0 = deg2rad(tangentPoint.y);
1670 cos0 = std::cos(dec0);
1671 sin0 = std::sin(dec0);
1672}
1673
1675 ra0 = dec0 = 0;
1676 cos0 = 1;
1677 sin0 = 0;
1678}
1679
1680Point TanRaDecToPixel::getTangentPoint() const { return Point(rad2deg(ra0), rad2deg(dec0)); }
1681
1683
1684// Use analytic derivatives, computed at the same time as the transform itself
1686 /* this routine is very similar to apply, but also propagates errors.
1687 The deg2rad and rad2deg are ignored for errors because they act as
1688 2 global scalings that cancel each other.
1689 Derivatives were computed using maple:
1690
1691 l1 := sin(a - a0)*cos(d);
1692 m1 := sin(d)*sin(d0)+cos(d)*cos(d0)*cos(a-a0);
1693 l2 := sin(d)*cos(d0)-cos(d)*sin(d0)*cos(a-a0);
1694 simplify(diff(l1/m1,a));
1695 simplify(diff(l1/m1,d));
1696 simplify(diff(l2/m1,a));
1697 simplify(diff(l2/m1,d));
1698
1699 Checked against AstrometryTransform::transformPosAndErrors (dec 09)
1700 */
1701 double ra = deg2rad(in.x);
1702 double dec = deg2rad(in.y);
1703 if (ra - ra0 > M_PI) ra -= (2. * M_PI);
1704 if (ra - ra0 < -M_PI) ra += (2. * M_PI);
1705 // Code inspired from worldpos.c in wcssubs (ancestor of the wcslib)
1706 // The same code is copied in ::apply()
1707
1708 double coss = std::cos(dec);
1709 double sins = std::sin(dec);
1710 double sinda = std::sin(ra - ra0);
1711 double cosda = std::cos(ra - ra0);
1712 double l = sinda * coss;
1713 double m = sins * sin0 + coss * cos0 * cosda;
1714 l = l / m;
1715 m = (sins * cos0 - coss * sin0 * cosda) / m;
1716
1717 // derivatives
1718 double deno =
1719 sq(sin0) - sq(coss) + sq(coss * cos0) * (1 + sq(cosda)) + 2 * sins * sin0 * coss * cos0 * cosda;
1720 double a11 = coss * (cosda * sins * sin0 + coss * cos0) / deno;
1721 double a12 = -sinda * sin0 / deno;
1722 double a21 = coss * sinda * sins / deno;
1723 double a22 = cosda / deno;
1724
1725 FatPoint tmp;
1726 tmp.vx = a11 * (a11 * in.vx + 2 * a12 * in.vxy) + a12 * a12 * in.vy;
1727 tmp.vy = a21 * a21 * in.vx + a22 * a22 * in.vy + 2. * a21 * a22 * in.vxy;
1728 tmp.vxy = a21 * a11 * in.vx + a22 * a12 * in.vy + (a21 * a12 + a11 * a22) * in.vxy;
1729
1730 // l and m are now coordinates in the tangent plane, in radians.
1731 tmp.x = rad2deg(l);
1732 tmp.y = rad2deg(m);
1733
1734 linTan2Pix.transformPosAndErrors(tmp, out);
1735}
1736
1737void TanRaDecToPixel::apply(const double xIn, const double yIn, double &xOut, double &yOut) const {
1738 double ra = deg2rad(xIn);
1739 double dec = deg2rad(yIn);
1740 if (ra - ra0 > M_PI) ra -= (2. * M_PI);
1741 if (ra - ra0 < -M_PI) ra += (2. * M_PI);
1742 // Code inspired from worldpos.c in wcssubs (ancestor of the wcslib)
1743 // The same code is copied in ::transformPosAndErrors()
1744 double coss = std::cos(dec);
1745 double sins = std::sin(dec);
1746 double l = std::sin(ra - ra0) * coss;
1747 double m = sins * sin0 + coss * cos0 * std::cos(ra - ra0);
1748 l = l / m;
1749 m = (sins * cos0 - coss * sin0 * std::cos(ra - ra0)) / m;
1750 // l and m are now coordinates in the tangent plane, in radians.
1751 l = rad2deg(l);
1752 m = rad2deg(m);
1753 linTan2Pix.apply(l, m, xOut, yOut);
1754}
1755
1759
1760void TanRaDecToPixel::print(ostream &stream) const {
1761 Point tp = getTangentPoint();
1762 stream << "tan2pix " << linTan2Pix << std::endl;
1763 stream << "tangent point " << tp.x << ' ' << tp.y << endl;
1764}
1765
1770
1775
1779
1782 "TanRaDecToPixel::fit is NOT implemented (although it is doable)) ");
1783 return -1;
1784}
1785
1786/************* a "run-time" transform, that does not require to modify this file */
1787
1789 : _userFun(userFun), _userData(userData) {}
1790
1791void UserTransform::apply(const double xIn, const double yIn, double &xOut, double &yOut) const {
1792 _userFun(xIn, yIn, xOut, yOut, _userData);
1793}
1794
1795void UserTransform::print(ostream &stream) const {
1796 stream << "UserTransform with user function @ " << _userFun << "and userData@ " << _userData << endl;
1797}
1798
1801 "UserTransform::fit is NOT implemented (and will never be)) ");
1802 return -1;
1803}
1804
1808
1809/*************************************************************/
1810
1812 ifstream s(fileName.c_str());
1813 if (!s)
1815 " astrometryTransformRead : cannot open " + fileName);
1816 try {
1818 s.close();
1819 return res;
1822 std::string(e.what()) + " in file " + fileName);
1823 }
1824}
1825
1828 s >> type;
1829 if (s.fail())
1831 "astrometryTransformRead : could not find a AstrometryTransformtype");
1832 if (type == "AstrometryTransformIdentity") {
1834 res->read(s);
1835 return std::move(res);
1836 } else if (type == "AstrometryTransformPolynomial") {
1838 res->read(s);
1839 return std::move(res);
1840 } else
1842 " astrometryTransformRead : No reader for AstrometryTransform type " + type);
1843}
1844} // namespace jointcal
1845} // namespace lsst
py::object result
Definition _schema.cc:429
int const step
int min
int max
table::Key< int > type
Definition Detector.cc:163
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
#define M_PI
Definition ListMatch.cc:31
LSST DM logging module built on log4cxx.
#define LOGLS_WARN(logger, message)
Log a warn-level message using an iostream-based interface.
Definition Log.h:659
#define LOGL_WARN(logger, message...)
Log a warn-level message using a varargs/printf style interface.
Definition Log.h:547
#define LOGLS_FATAL(logger, message)
Log a fatal-level message using an iostream-based interface.
Definition Log.h:699
#define LOG_GET(logger)
Returns a Log object associated with logger.
Definition Log.h:75
#define LOG_LOGGER
Definition Log.h:714
#define LOGL_ERROR(logger, message...)
Log a error-level message using a varargs/printf style interface.
Definition Log.h:563
#define LOGLS_TRACE(logger, message)
Log a trace-level message using an iostream-based interface.
Definition Log.h:599
#define LOGL_FATAL(logger, message...)
Log a fatal-level message using a varargs/printf style interface.
Definition Log.h:579
double dec
Definition Match.cc:41
table::PointKey< double > crpix
Definition OldWcs.cc:129
int y
Definition SpanSet.cc:48
int m
Definition SpanSet.cc:48
pairs of points
table::Key< int > transform
table::Key< int > a
T at(T... args)
T atan2(T... args)
T atan(T... args)
T begin(T... args)
T c_str(T... args)
Private class to handle AstrometryTransform compositions (i.e.
void apply(double xIn, double yIn, double &xOut, double &yOut) const
return second(first(xIn,yIn))
double fit(StarMatchList const &starMatchList)
fits a transform to a std::list of Point pairs (p1,p2, the Point fields in StarMatch).
void print(ostream &stream) const
prints the transform coefficients to stream.
std::unique_ptr< AstrometryTransform > clone() const
returns a copy (allocated by new) of the transformation.
AstrometryTransformComposition(AstrometryTransform const &second, AstrometryTransform const &first)
a virtual (interface) class for geometric transformations.
virtual void paramDerivatives(Point const &where, double *dx, double *dy) const
Derivative w.r.t parameters.
virtual std::unique_ptr< AstrometryTransform > roughInverse(const Frame &region) const
Rough inverse.
virtual std::size_t getNpar() const
returns the number of parameters (to compute chi2's)
virtual double getJacobian(Point const &point) const
returns the local jacobian.
void write(const std::string &fileName) const
virtual AstrometryTransformLinear linearApproximation(Point const &where, double step=0.01) const
linear (local) approximation.
virtual double paramRef(Eigen::Index i) const
virtual void apply(double xIn, double yIn, double &xOut, double &yOut) const =0
virtual void computeDerivative(Point const &where, AstrometryTransformLinear &derivative, double step=0.01) const
Computes the local Derivative of a transform, w.r.t.
virtual void transformPosAndErrors(const FatPoint &in, FatPoint &out) const
virtual std::unique_ptr< AstrometryTransform > composeAndReduce(AstrometryTransform const &right) const
Return a reduced composition of newTransform = this(right()), or nullptr if it cannot be reduced.
void offsetParams(Eigen::VectorXd const &delta)
virtual std::unique_ptr< AstrometryTransform > inverseTransform(double precision, const Frame &region) const
returns an inverse transform. Numerical if not overloaded.
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))
void getParams(double *params) const
params should be at least Npar() long
virtual std::unique_ptr< AstrometryTransform > clone() const =0
returns a copy (allocated by new) of the transformation.
A do-nothing transformation. It anyway has dummy routines to mimick a AstrometryTransform.
void computeDerivative(Point const &where, AstrometryTransformLinear &derivative, double step=0.01) const override
Computes the local Derivative of a transform, w.r.t.
virtual AstrometryTransformLinear linearApproximation(Point const &where, double step=0.01) const override
linear approximation.
void write(std::ostream &s) const override
std::shared_ptr< ast::Mapping > toAstMap(jointcal::Frame const &domain) const override
Create an equivalent AST mapping for this transformation, including an analytic inverse if possible.
double fit(StarMatchList const &starMatchList)
fits a transform to a std::list of Point pairs (p1,p2, the Point fields in StarMatch).
void apply(double xIn, double yIn, double &xOut, double &yOut) const
implements an iterative (Gauss-Newton) solver.
void print(ostream &out) const
prints the transform coefficients to stream.
AstrometryTransformInverse(const AstrometryTransform *direct, double precision, const Frame &region)
virtual std::unique_ptr< AstrometryTransform > clone() const
returns a copy (allocated by new) of the transformation.
std::unique_ptr< AstrometryTransform > inverseTransform(double, const Frame &) const
Inverse transform: returns the direct one!
std::unique_ptr< AstrometryTransform > roughInverse(const Frame &) const
Overload the "generic routine".
implements the linear transformations (6 real coefficients).
void computeDerivative(Point const &where, AstrometryTransformLinear &derivative, double step=0.01) const override
Computes the local Derivative of a transform, w.r.t.
AstrometryTransformLinear linearApproximation(Point const &where, double step=0.01) const override
linear (local) approximation.
void print(std::ostream &out) const override
prints the transform coefficients to stream.
AstrometryTransformLinear operator*(AstrometryTransformLinear const &right) const
enables to combine linear tranformations: T1=T2*T3 is legal.
AstrometryTransformLinear inverted() const
returns the inverse: T1 = T2.inverted();
AstrometryTransformLinear()
the default constructor constructs the do-nothing transformation.
void apply(double xIn, double yIn, double &xOut, double &yOut) const override
std::unique_ptr< AstrometryTransform > inverseTransform(double precision, const Frame &region) const override
returns an inverse transform. Numerical if not overloaded.
double fit(StarMatchList const &starMatchList)
fits a transform to a std::list of Point pairs (p1,p2, the Point fields in StarMatch).
just here to provide specialized constructors. AstrometryTransformLinear fit routine.
just here to provide a specialized constructor, and fit.
double fit(StarMatchList const &starMatchList)
fits a transform to a std::list of Point pairs (p1,p2, the Point fields in StarMatch).
AstrometryTransformLinearShift(double ox=0., double oy=0.)
Add ox and oy.
double coeffOrZero(std::size_t powX, std::size_t powY, std::size_t whichCoord) const
read access, zero if beyond order
std::shared_ptr< ast::Mapping > toAstMap(jointcal::Frame const &domain) const override
Create an equivalent AST mapping for this transformation, including an analytic inverse if possible.
AstrometryTransformPolynomial operator-(AstrometryTransformPolynomial const &right) const
Subtraction.
AstrometryTransformPolynomial operator*(AstrometryTransformPolynomial const &right) const
Composition (internal stuff in quadruple precision)
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 paramRef(Eigen::Index i) const override
void setOrder(std::size_t order)
Sets the polynomial order (the highest sum of exponents of the largest monomial).
double fit(StarMatchList const &starMatchList) override
guess what
void print(std::ostream &out) const override
print out of coefficients in a readable form.
void paramDerivatives(Point const &where, double *dx, double *dy) const override
Derivative w.r.t parameters.
void write(std::ostream &s) const override
virtual void transformPosAndErrors(const FatPoint &in, FatPoint &out) const override
a mix of apply and Derivative
std::size_t getOrder() const
Returns the polynomial order.
AstrometryTransformPolynomial(std::size_t order=1)
Default transform : identity for all orders (>=1 ).
AstrometryTransformPolynomial operator+(AstrometryTransformPolynomial const &right) const
Addition.
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
std::unique_ptr< AstrometryTransform > composeAndReduce(AstrometryTransformPolynomial const &right) const
Return a reduced composition of newTransform = this(right()), or nullptr if it cannot be reduced.
void print(std::ostream &out) const override
prints the transform coefficients to stream.
std::unique_ptr< AstrometryTransform > clone() const override
returns a copy (allocated by new) of the transformation.
std::shared_ptr< afw::geom::SkyWcs > getSkyWcs() const
void apply(double xIn, double yIn, double &xOut, double &yOut) const override
double fit(const StarMatchList &starMatchList) override
Not implemented; throws pex::exceptions::LogicError.
AstrometryTransformSkyWcs(std::shared_ptr< afw::geom::SkyWcs > skyWcs)
std::unique_ptr< AstrometryTransformPolynomial > corr
virtual void pixToTangentPlane(double xPixel, double yPixel, double &xTangentPlane, double &yTangentPlane) const =0
Transform from pixels to tangent plane (degrees)
Point getCrPix() const
Get the pixel origin of the WCS (CRPIX in FITS WCS terminology, but zero-based)
BaseTanWcs & operator=(const BaseTanWcs &original)
BaseTanWcs(AstrometryTransformLinear const &pixToTan, Point const &tangentPoint, const AstrometryTransformPolynomial *corrections=nullptr)
AstrometryTransformLinear getLinPart() const
The Linear part (corresponding to CD's and CRPIX's)
Point getTangentPoint() const
Get the sky origin (CRVAL in FITS WCS terminology) in degrees.
void apply(double xIn, double yIn, double &xOut, double &yOut) const
Transform pixels to ICRS RA, Dec in degrees.
AstrometryTransformLinear linPixelToTan
void setCorrections(std::unique_ptr< AstrometryTransformPolynomial > corrections)
Assign the correction polynomial (what it means is left to derived classes)
A Point with uncertainties.
Definition FatPoint.h:34
rectangle with sides parallel to axes.
Definition Frame.h:38
Point getCenter() const
Center of the frame.
Definition Frame.h:60
double getHeight() const
size along y axis
Definition Frame.h:57
double xMin
coordinate of boundary.
Definition Frame.h:41
double getArea() const
Definition Frame.cc:118
double getWidth() const
size along x axis
Definition Frame.h:54
A point in a plane.
Definition Point.h:37
double x
coordinate
Definition Point.h:42
PolyXY(AstrometryTransformPolynomial const &transform, std::size_t whichCoord)
long double getCoefficient(std::size_t powX, std::size_t powY) const
std::size_t getOrder() const
long double & getCoefficient(std::size_t powX, std::size_t powY)
A hanger for star associations.
Definition StarMatch.h:55
The transformation that handles pixels to sideral transformations (Gnomonic, possibly with polynomial...
TanPixelToRaDec operator*(AstrometryTransformLinear const &right) const
composition with AstrometryTransformLinear
std::unique_ptr< AstrometryTransform > clone() const
returns a copy (allocated by new) of the transformation.
TanRaDecToPixel inverted() const
approximate inverse : it ignores corrections;
double fit(StarMatchList const &starMatchList)
Not implemented yet, because we do it otherwise.
virtual void pixToTangentPlane(double xPixel, double yPixel, double &xTangentPlane, double &yTangentPlane) const
transforms from pixel space to tangent plane (degrees)
void print(std::ostream &out) const
prints the transform coefficients to stream.
std::unique_ptr< AstrometryTransform > roughInverse(const Frame &region) const
Overload the "generic routine" (available for all AstrometryTransform types.
AstrometryTransformPolynomial getPixelToTangentPlane() const
the transformation from pixels to tangent plane (degrees)
std::unique_ptr< AstrometryTransform > inverseTransform(double precision, const Frame &region) const
Inverse transform: returns a TanRaDecToPixel if there are no corrections, or the iterative solver if ...
std::unique_ptr< AstrometryTransform > composeAndReduce(AstrometryTransformLinear const &right) const
Return a reduced composition of newTransform = this(right()), or nullptr if it cannot be reduced.
This one is the Tangent Plane (called gnomonic) projection (from celestial sphere to tangent plane)
void transformPosAndErrors(const FatPoint &in, FatPoint &out) const
transform with analytical derivatives
Point getTangentPoint() const
tangent point coordinates (degrees)
double fit(StarMatchList const &starMatchList)
fits a transform to a std::list of Point pairs (p1,p2, the Point fields in StarMatch).
void setTangentPoint(Point const &tangentPoint)
Resets the projection (or tangent) point.
std::unique_ptr< AstrometryTransform > inverseTransform(double precision, const Frame &region) const
Inverse transform: returns a TanPixelToRaDec.
TanPixelToRaDec inverted() const
exact typed inverse:
void apply(double xIn, double yIn, double &xOut, double &yOut) const
std::unique_ptr< AstrometryTransform > clone() const
returns a copy (allocated by new) of the transformation.
AstrometryTransformLinear getLinPart() const
The Linear part (corresponding to CD's and CRPIX's)
std::unique_ptr< AstrometryTransform > roughInverse(const Frame &region) const
Overload the "generic routine" (available for all AstrometryTransform types.
void print(std::ostream &out) const
prints the transform coefficients to stream.
std::unique_ptr< AstrometryTransform > inverseTransform(double precision, const Frame &region) const
Inverse transform: returns a TanRaDecToPixel if there are no corrections, or the iterative solver if ...
void print(std::ostream &out) const
prints the transform coefficients to stream.
double fit(StarMatchList const &starMatchList)
Not implemented yet, because we do it otherwise.
virtual void pixToTangentPlane(double xPixel, double yPixel, double &xTangentPlane, double &yTangentPlane) const
transforms from pixel space to tangent plane (degrees)
std::unique_ptr< AstrometryTransform > clone() const
returns a copy (allocated by new) of the transformation.
AstrometryTransformPolynomial getPixelToTangentPlane() const
the transformation from pixels to tangent plane (degrees)
void print(std::ostream &out) const
prints the transform coefficients to stream.
std::unique_ptr< AstrometryTransform > clone() const
returns a copy (allocated by new) of the transformation.
UserTransform(AstrometryTransformFun &fun, const void *userData)
the transform routine and extra data that it may need.
double fit(StarMatchList const &starMatchList)
fits a transform to a std::list of Point pairs (p1,p2, the Point fields in StarMatch).
void apply(double xIn, double yIn, double &xOut, double &yOut) const
virtual char const * what(void) const noexcept
Return a character string summarizing this exception.
Definition Exception.cc:99
Reports invalid arguments.
Definition Runtime.h:66
Reports errors in the logical structure of the program.
Definition Runtime.h:46
Reports attempts to access elements using an invalid key.
Definition Runtime.h:151
Reports errors that are due to events beyond the control of the program.
Definition Runtime.h:104
T cos(T... args)
T count(T... args)
T emplace_back(T... args)
T endl(T... args)
T fabs(T... args)
T floor(T... args)
T forward(T... args)
T infinity(T... args)
T insert(T... args)
T isnan(T... args)
T left(T... args)
T max(T... args)
T min(T... args)
T move(T... args)
Low-level polynomials (including special polynomials) in C++.
AstrometryTransformLinear normalizeCoordinatesTransform(const Frame &frame)
Returns the transformation that maps the input frame along both axes to [-1,1].
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.
bool isIntegerShift(const AstrometryTransform *transform)
Shorthand test to tell if a transform is a simple integer shift.
std::unique_ptr< AstrometryTransform > compose(AstrometryTransform const &left, AstrometryTransform const &right)
Returns a pointer to a composition of transforms, representing left(right()).
std::unique_ptr< AstrometryTransform > astrometryTransformRead(const std::string &fileName)
The virtual constructor from a file.
void(const double, const double, double &, double &, const void *) AstrometryTransformFun
signature of the user-provided routine that actually does the coordinate transform for UserTransform.
std::ostream & operator<<(std::ostream &stream, AstrometryMapping const &mapping)
STL namespace.
T pow(T... args)
T precision(T... args)
T push_back(T... args)
T reserve(T... args)
T resize(T... args)
T setprecision(T... args)
T sin(T... args)
T size(T... args)
T sqrt(T... args)
table::Key< table::Array< int > > degree
Definition PsfexPsf.cc:360
T str(T... args)
table::Key< int > order