LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
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 <iostream>
26 #include <iomanip>
27 #include <iterator> /* for ostream_iterator */
28 #include <limits>
29 #include <cmath>
30 #include <math.h>
31 #include <fstream>
32 #include "assert.h"
33 #include <sstream>
34 
35 #include "Eigen/Core"
36 
37 #include "lsst/log/Log.h"
38 #include "lsst/geom/Point.h"
40 #include "lsst/jointcal/Frame.h"
42 #include "lsst/pex/exceptions.h"
43 #include "Eigen/Cholesky"
44 
46 
47 using namespace std;
48 
49 namespace {
50 LOG_LOGGER _log = LOG_GET("lsst.jointcal.AstrometryTransform");
51 }
52 
53 namespace lsst {
54 namespace jointcal {
55 
57  const AstrometryTransformPolynomial *shift =
58  dynamic_cast<const AstrometryTransformPolynomial *>(transform);
59  if (shift == nullptr) return false;
60 
61  static const double eps = 1e-5;
62 
63  double dx = shift->getCoefficient(0, 0, 0);
64  double dy = shift->getCoefficient(0, 0, 1);
65 
66  static Point dumb(4000, 4000);
67  if (fabs(dx - int(floor(dx + 0.5))) < eps && fabs(dy - int(floor(dy + 0.5))) < eps &&
68  fabs(dumb.x + dx - shift->apply(dumb).x) < eps && fabs(dumb.y + dy - shift->apply(dumb).y) < eps)
69  return true;
70 
71  return false;
72 }
73 
74 /********* AstrometryTransform ***********************/
75 
76 Frame AstrometryTransform::apply(Frame const &inputframe, bool inscribed) const {
77  // 2 opposite corners
78  double xtmin1, xtmax1, ytmin1, ytmax1;
79  apply(inputframe.xMin, inputframe.yMin, xtmin1, ytmin1);
80  apply(inputframe.xMax, inputframe.yMax, xtmax1, ytmax1);
81  Frame fr1(std::min(xtmin1, xtmax1), std::min(ytmin1, ytmax1), std::max(xtmin1, xtmax1),
82  std::max(ytmin1, ytmax1));
83  // 2 other corners
84  double xtmin2, xtmax2, ytmin2, ytmax2;
85  apply(inputframe.xMin, inputframe.yMax, xtmin2, ytmax2);
86  apply(inputframe.xMax, inputframe.yMin, xtmax2, ytmin2);
87  Frame fr2(std::min(xtmin2, xtmax2), std::min(ytmin2, ytmax2), std::max(xtmin2, xtmax2),
88  std::max(ytmin2, ytmax2));
89 
90  if (inscribed) return fr1 * fr2;
91  return fr1 + fr2;
92 }
93 
94 std::unique_ptr<AstrometryTransform> AstrometryTransform::composeAndReduce(
95  AstrometryTransform const &) const { // by default no way to compose
97 }
98 
99 double AstrometryTransform::getJacobian(const double x, const double y) const {
100  double x2, y2;
101  double eps = x * 0.01;
102  if (eps == 0) eps = 0.01;
103  apply(x, y, x2, y2);
104  double dxdx, dydx;
105  apply(x + eps, y, dxdx, dydx);
106  dxdx -= x2;
107  dydx -= y2;
108  double dxdy, dydy;
109  apply(x, y + eps, dxdy, dydy);
110  dxdy -= x2;
111  dydy -= y2;
112  return ((dxdx * dydy - dxdy * dydx) / (eps * eps));
113 }
114 
118 void AstrometryTransform::computeDerivative(Point const &where, AstrometryTransformLinear &derivative,
119  const double step) const {
120  double x = where.x;
121  double y = where.y;
122  double xp0, yp0;
123  apply(x, y, xp0, yp0);
124 
125  double xp, yp;
126  apply(x + step, y, xp, yp);
127  derivative.a11() = (xp - xp0) / step;
128  derivative.a21() = (yp - yp0) / step;
129  apply(x, y + step, xp, yp);
130  derivative.a12() = (xp - xp0) / step;
131  derivative.a22() = (yp - yp0) / step;
132  derivative.dx() = 0;
133  derivative.dy() = 0;
134 }
135 
136 AstrometryTransformLinear AstrometryTransform::linearApproximation(Point const &where,
137  const double step) const {
138  Point outwhere = apply(where);
140  computeDerivative(where, der, step);
141  return AstrometryTransformLinearShift(outwhere.x, outwhere.y) * der *
142  AstrometryTransformLinearShift(-where.x, -where.y);
143 }
144 
145 void AstrometryTransform::transformPosAndErrors(FatPoint const &in, FatPoint &out) const {
146  FatPoint res; // in case in and out are the same address...
147  res = apply(in);
149  // could save a call here, since Derivative needs the transform of where that we already have
150  // 0.01 may not be a very good idea in all cases. May be we should provide a way of altering that.
151  computeDerivative(in, der, 0.01);
152  double a11 = der.A11();
153  double a22 = der.A22();
154  double a21 = der.A21();
155  double a12 = der.A12();
156  res.vx = a11 * (a11 * in.vx + 2 * a12 * in.vxy) + a12 * a12 * in.vy;
157  res.vy = a21 * a21 * in.vx + a22 * a22 * in.vy + 2. * a21 * a22 * in.vxy;
158  res.vxy = a21 * a11 * in.vx + a22 * a12 * in.vy + (a21 * a12 + a11 * a22) * in.vxy;
159  out = res;
160 }
161 
162 void AstrometryTransform::transformErrors(Point const &where, const double *vIn, double *vOut) const {
164  computeDerivative(where, der, 0.01);
165  double a11 = der.A11();
166  double a22 = der.A22();
167  double a21 = der.A21();
168  double a12 = der.A12();
169 
170  /* (a11 a12) (vxx vxy)
171  M = ( ) and V = ( )
172  (a21 a22) (xvy vyy)
173 
174  Vxx = Vin[0], vyy = Vin[1], Vxy = Vin[2];
175  we want to compute M*V*tp(M)
176  A lin alg light package would be perfect...
177  */
178  int xx = 0;
179  int yy = 1;
180  int xy = 2;
181  // M*V :
182 
183  double b11 = a11 * vIn[xx] + a12 * vIn[xy];
184  double b22 = a21 * vIn[xy] + a22 * vIn[yy];
185  double b12 = a11 * vIn[xy] + a12 * vIn[yy];
186  double b21 = a21 * vIn[xx] + a22 * vIn[xy];
187 
188  // (M*V) * tp(M)
189 
190  vOut[xx] = b11 * a11 + b12 * a12;
191  vOut[xy] = b11 * a21 + b12 * a22;
192  vOut[yy] = b21 * a21 + b22 * a22;
193 }
194 
195 std::unique_ptr<AstrometryTransform> AstrometryTransform::roughInverse(const Frame &region) const {
196  // "in" and "out" refer to the inverse direction.
197  Point centerOut = region.getCenter();
198  Point centerIn = apply(centerOut);
200  computeDerivative(centerOut, der, std::sqrt(region.getArea()) / 5.);
201  der = der.inverted();
202  der = AstrometryTransformLinearShift(centerOut.x, centerOut.y) * der *
203  AstrometryTransformLinearShift(-centerIn.x, -centerIn.y);
205 }
206 
207 /* implement one in AstrometryTransform, so that all derived
208  classes do not need to provide one... */
209 
210 /* the routines that follow are used for ea generic parameter
211  transformation serialization, used e.g. for fits. Enables
212  to manipulate transformation parameters as vectors.
213 */
214 
215 // not dummy : what it does is virtual because paramRef is virtual.
216 void AstrometryTransform::getParams(double *params) const {
217  std::size_t npar = getNpar();
218  for (std::size_t i = 0; i < npar; ++i) params[i] = paramRef(i);
219 }
220 
221 void AstrometryTransform::offsetParams(Eigen::VectorXd const &delta) {
222  std::size_t npar = getNpar();
223  for (std::size_t i = 0; i < npar; ++i) paramRef(i) += delta[i];
224 }
225 
226 double AstrometryTransform::paramRef(Eigen::Index const) const {
228  std::string("AstrometryTransform::paramRef should never be called "));
229 }
230 
231 double &AstrometryTransform::paramRef(Eigen::Index const) {
233  "AstrometryTransform::paramRef should never be called ");
234 }
235 
236 void AstrometryTransform::paramDerivatives(Point const &, double *, double *) const {
238  "AstrometryTransform::paramDerivatives() should never be called ");
239 }
240 
242  transform.print(stream);
243  return stream;
244 }
245 
246 void AstrometryTransform::write(const std::string &fileName) const {
247  ofstream s(fileName.c_str());
248  write(s);
249  bool ok = !s.fail();
250  s.close();
251  if (!ok)
253  "AstrometryTransform::write, something went wrong for file " + fileName);
254 }
255 
256 void AstrometryTransform::write(ostream &stream) const {
257  throw LSST_EXCEPT(
259  "AstrometryTransform::write(ostream), should never be called. MEans that it is missing in some "
260  "derived class ");
261 }
262 
263 /******************* GTransformInverse ****************/
264 /* inverse transformation, solved by iterations. Before using
265  it (probably via AstrometryTransform::inverseTransform), consider
266  seriously StarMatchList::inverseTransform */
268 private:
271  double precision2;
272 
273 public:
274  AstrometryTransformInverse(const AstrometryTransform *direct, const double precision,
275  const Frame &region);
276 
279  void apply(const double xIn, const double yIn, double &xOut, double &yOut) const;
280 
281  void print(ostream &out) const;
282 
283  double fit(StarMatchList const &starMatchList);
284 
286 
288 
290  std::unique_ptr<AstrometryTransform> roughInverse(const Frame &) const { return _direct->clone(); }
291 
294  return _direct->clone();
295  }
296 
298 
299 private:
301 };
302 
303 std::unique_ptr<AstrometryTransform> AstrometryTransform::inverseTransform(const double precision,
304  const Frame &region) const {
306 }
307 
308 AstrometryTransformInverse::AstrometryTransformInverse(const AstrometryTransform *direct,
309  const double precision, const Frame &region) {
310  _direct = direct->clone();
311  _roughInverse = _direct->roughInverse(region);
312  precision2 = precision * precision;
313 }
314 
315 AstrometryTransformInverse::AstrometryTransformInverse(AstrometryTransformInverse const &model)
316  : AstrometryTransform() {
317  _direct = model._direct->clone();
318  _roughInverse = model._roughInverse->clone();
319  precision2 = model.precision2;
320 }
321 
323 
324 void AstrometryTransformInverse::operator=(AstrometryTransformInverse const &model) {
325  _direct = model._direct->clone();
326  _roughInverse = model._roughInverse->clone();
327  precision2 = model.precision2;
328 }
329 
330 void AstrometryTransformInverse::apply(const double xIn, const double yIn, double &xOut, double &yOut) const {
331  Point in(xIn, yIn);
332  Point outGuess = _roughInverse->apply(in);
333  AstrometryTransformLinear directDer, reverseDer;
334  int loop = 0;
335  int maxloop = 20;
336  double move2;
337  do {
338  loop++;
339  Point inGuess = _direct->apply(outGuess);
340  _direct->computeDerivative(outGuess, directDer);
341  reverseDer = directDer.inverted();
342  double xShift, yShift;
343  reverseDer.apply(xIn - inGuess.x, yIn - inGuess.y, xShift, yShift);
344  outGuess.x += xShift;
345  outGuess.y += yShift;
346  move2 = xShift * xShift + yShift * yShift;
347  } while ((move2 > precision2) && (loop < maxloop));
348  if (loop == maxloop) LOGLS_WARN(_log, "Problems applying AstrometryTransformInverse at " << in);
349  xOut = outGuess.x;
350  yOut = outGuess.y;
351 }
352 
354  stream << " AstrometryTransformInverse of :" << endl << *_direct << endl;
355 }
356 
359  "Cannot fit a AstrometryTransformInverse. Use StarMatchList::inverseTransform instead.");
360 }
361 
364 }
365 
366 /************* AstrometryTransformComposition **************/
367 
368 // This class was done to allow composition of AstrometryTransform's, without specifications of their types.
369 // does not need to be public. Invoked by compose(left,right)
370 
374 private:
375  std::unique_ptr<AstrometryTransform> _first, _second;
376 
377 public:
379 
381  void apply(const double xIn, const double yIn, double &xOut, double &yOut) const;
382  void print(ostream &stream) const;
383 
385  double fit(StarMatchList const &starMatchList);
386 
389 };
390 
392  AstrometryTransform const &first) {
393  _first = first.clone();
394  _second = second.clone();
395 }
396 
397 void AstrometryTransformComposition::apply(const double xIn, const double yIn, double &xOut,
398  double &yOut) const {
399  double xout, yout;
400  _first->apply(xIn, yIn, xout, yout);
401  _second->apply(xout, yout, xOut, yOut);
402 }
403 
405  stream << "Composed AstrometryTransform consisting of:" << std::endl;
406  _first->print(stream);
407  _second->print(stream);
408 }
409 
411  /* fits only one of them. could check that first can actually be fitted... */
412  return _first->fit(starMatchList);
413 }
414 
416  return std::make_unique<AstrometryTransformComposition>(*_second, *_first);
417 }
418 
420 
422  AstrometryTransformIdentity const &right) {
423  return left.clone();
424 }
425 
427  AstrometryTransform const &right) {
428  // Try to use the composeAndReduce method from left. If absent, AstrometryTransform::composeAndReduce
429  // returns NULL. composeAndReduce is non trivial for polynomials.
430  std::unique_ptr<AstrometryTransform> composition(left.composeAndReduce(right));
431  // composition == NULL means no reduction: just build a Composition that pipelines "left" and "right".
432  if (composition == nullptr)
433  return std::make_unique<AstrometryTransformComposition>(left, right);
434  else
435  return composition;
436 }
437 
438 // just a speed up, to avoid useless numerical derivation.
440  const double) const {
441  derivative = AstrometryTransformLinear();
442 }
443 
445  const double) const {
447  return result; // rely on default AstrometryTransformlin constructor;
448 }
449 
451  return std::make_shared<ast::UnitMap>(2); // a AstrometryTransformIdentity is identically ast::UnitMap(2)
452 }
453 
455  stream << "AstrometryTransformIdentity 1" << endl;
456 }
457 
459  int format;
460  stream >> format;
461  if (format != 1)
463  " AstrometryTransformIdentity::read : format is not 1 ");
464 }
465 
466 /*************** AstrometryTransformPolynomial **************************************/
467 
469 
471  _nterms = (order + 1) * (order + 2) / 2;
472 
473  // allocate and fill coefficients
474  _coeffs.resize(2 * _nterms, 0.);
475  // the default is supposed to be the identity, (for order>=1).
476  if (_order >= 1) {
477  getCoefficient(1, 0, 0) = 1;
478  getCoefficient(0, 1, 1) = 1;
479  }
480 }
481 
482 //#ifdef TO_BE_FIXED
484  const Frame &frame, std::size_t order,
485  std::size_t nPoint) {
486  StarMatchList sm;
487 
488  double step = std::sqrt(fabs(frame.getArea()) / double(nPoint));
489  for (double x = frame.xMin + step / 2; x <= frame.xMax; x += step)
490  for (double y = frame.yMin + step / 2; y <= frame.yMax; y += step) {
491  auto pix = std::make_shared<BaseStar>(x, y, 0, 0);
492  double xtr, ytr;
493  transform->apply(x, y, xtr, ytr);
494  auto tp = std::make_shared<BaseStar>(xtr, ytr, 0, 0);
495  /* These are fake stars so no need to transform fake errors.
496  all errors (and weights) will be equal : */
497  sm.push_back(StarMatch(*pix, *tp, pix, tp));
498  }
500  ret.fit(sm);
501  *this = ret;
502 }
503 //#endif
504 
507  std::size_t order, std::size_t nSteps) {
508  jointcal::StarMatchList starMatchList;
509  double xStart = domain.xMin;
510  double yStart = domain.yMin;
511  double xStep = domain.getWidth() / (nSteps + 1);
512  double yStep = domain.getHeight() / (nSteps + 1);
513  for (std::size_t i = 0; i < nSteps; ++i) {
514  for (std::size_t j = 0; j < nSteps; ++j) {
515  // TODO: once DM-4044 is done, we can remove the redundancy in `Point`/`Point2D` here
516  jointcal::Point in(xStart + i * xStep, yStart + j * yStep);
517  geom::Point2D inAfw(in.x, in.y);
518  geom::Point2D outAfw = transform->applyForward(inAfw);
519  jointcal::Point out(outAfw.getX(), outAfw.getY());
520  starMatchList.emplace_back(in, out, nullptr, nullptr);
521  }
522  }
524  poly.fit(starMatchList);
525  *this = poly;
526 }
527 
528 void AstrometryTransformPolynomial::computeMonomials(double xIn, double yIn, double *monomial) const {
529  /* The ordering of monomials is implemented here.
530  You may not change it without updating the "mapping" routines
531  getCoefficient(std::size_t, std::size_t, std::size_t).
532  I (P.A.) did not find a clever way to loop over monomials.
533  Improvements welcome.
534  This routine is used also by the fit to fill monomials.
535  We could certainly be more elegant.
536  */
537 
538  double xx = 1;
539  for (std::size_t ix = 0; ix <= _order; ++ix) {
540  double yy = 1;
541  std::size_t k = ix * (ix + 1) / 2;
542  for (std::size_t iy = 0; iy <= _order - ix; ++iy) {
543  monomial[k] = xx * yy;
544  yy *= yIn;
545  k += ix + iy + 2;
546  }
547  xx *= xIn;
548  }
549 }
550 
552  _order = order;
553  std::size_t old_nterms = _nterms;
554  _nterms = (_order + 1) * (_order + 2) / 2;
555 
556  // temporarily save coefficients
557  vector<double> old_coeffs = _coeffs;
558  // reallocate enough size
559  _coeffs.resize(2 * _nterms);
560  // reassign to zero (this is necessary because ycoeffs
561  // are after xcoeffs and so their meaning changes
562  for (std::size_t k = 0; k < _nterms; ++k) _coeffs[k] = 0;
563  // put back what we had before
564  std::size_t kmax = min(old_nterms, _nterms);
565  for (std::size_t k = 0; k < kmax; ++k) {
566  _coeffs[k] = old_coeffs[k]; // x terms
567  _coeffs[k + _nterms] = old_coeffs[k + old_nterms]; // y terms
568  }
569 }
570 
571 /* this is reasonably fast, when optimized */
572 void AstrometryTransformPolynomial::apply(const double xIn, const double yIn, double &xOut,
573  double &yOut) const {
574  /*
575  This routine computes the monomials only once for both
576  polynomials. This is why AstrometryTransformPolynomial does not use an auxilary
577  class (such as PolyXY) to handle each polynomial.
578 
579  The code works even if &xIn == &xOut (or &yIn == &yOut)
580  It uses Variable Length Allocation (VLA) rather than a vector<double>
581  because allocating the later costs about 50 ns. All VLA uses are tagged.
582  */
583  double monomials[_nterms]; // this is VLA, which is (perhaps) not casher C++
584  computeMonomials(xIn, yIn, monomials);
585 
586  xOut = 0;
587  yOut = 0;
588  const double *c = &_coeffs[0];
589  const double *pm = &monomials[0];
590  // the ordering of the coefficients and the monomials are identical.
591  for (int k = _nterms; k--;) xOut += (*(pm++)) * (*(c++));
592  pm = &monomials[0];
593  for (int k = _nterms; k--;) yOut += (*(pm++)) * (*(c++));
594 }
595 
597  AstrometryTransformLinear &derivative,
598  const double step)
599  const { /* routine checked against numerical derivatives from AstrometryTransform::Derivative */
600  if (_order == 1) {
601  derivative = AstrometryTransformLinear(*this);
602  derivative.dx() = derivative.dy() = 0;
603  return;
604  }
605 
606  double dermx[2 * _nterms]; // VLA
607  double *dermy = dermx + _nterms;
608  double xin = where.x;
609  double yin = where.y;
610 
611  double xx = 1;
612  double xxm1 = 1; // xx^(ix-1)
613  for (std::size_t ix = 0; ix <= _order; ++ix) {
614  std::size_t k = (ix) * (ix + 1) / 2;
615  // iy = 0
616  dermx[k] = ix * xxm1;
617  dermy[k] = 0;
618  k += ix + 2;
619  double yym1 = 1; // yy^(iy-1)
620  for (std::size_t iy = 1; iy <= _order - ix; ++iy) {
621  dermx[k] = ix * xxm1 * yym1 * yin;
622  dermy[k] = iy * xx * yym1;
623  yym1 *= yin;
624  k += ix + iy + 2;
625  }
626  xx *= xin;
627  if (ix >= 1) xxm1 *= xin;
628  }
629 
630  derivative.dx() = 0;
631  derivative.dy() = 0;
632 
633  const double *mx = &dermx[0];
634  const double *my = &dermy[0];
635  const double *c = &_coeffs[0];
636  // dx'
637  double a11 = 0, a12 = 0;
638  for (int k = _nterms; k--;) {
639  a11 += (*(mx++)) * (*c);
640  a12 += (*(my++)) * (*(c++));
641  }
642  derivative.a11() = a11;
643  derivative.a12() = a12;
644  // dy'
645  double a21 = 0, a22 = 0;
646  mx = &dermx[0];
647  my = &dermy[0];
648  for (int k = _nterms; k--;) {
649  a21 += (*(mx++)) * (*c);
650  a22 += (*(my++)) * (*(c++));
651  }
652  derivative.a21() = a21;
653  derivative.a22() = a22;
654 }
655 
657  /*
658  The results from this routine were compared to what comes out
659  from apply and transformErrors. The Derivative routine was
660  checked against numerical derivatives from
661  AstrometryTransform::Derivative. (P.A dec 2009).
662 
663  This routine could be made much simpler by calling apply and
664  Derivative (i.e. you just suppress it, and the fallback is the
665  generic version in AstrometryTransform). BTW, I checked that both routines
666  provide the same result. This version is however faster
667  (monomials get recycled).
668  */
669  double monomials[_nterms]; // VLA
670 
671  FatPoint res; // to store the result, because nothing forbids &in == &out.
672 
673  double dermx[2 * _nterms]; // monomials for derivative w.r.t. x (VLA)
674  double *dermy = dermx + _nterms; // same for y
675  double xin = in.x;
676  double yin = in.y;
677 
678  double xx = 1;
679  double xxm1 = 1; // xx^(ix-1)
680  for (std::size_t ix = 0; ix <= _order; ++ix) {
681  std::size_t k = (ix) * (ix + 1) / 2;
682  // iy = 0
683  dermx[k] = ix * xxm1;
684  dermy[k] = 0;
685  monomials[k] = xx;
686  k += ix + 2;
687  double yy = yin;
688  double yym1 = 1; // yy^(iy-1)
689  for (std::size_t iy = 1; iy <= _order - ix; ++iy) {
690  monomials[k] = xx * yy;
691  dermx[k] = ix * xxm1 * yy;
692  dermy[k] = iy * xx * yym1;
693  yym1 *= yin;
694  yy *= yin;
695  k += ix + iy + 2;
696  }
697  xx *= xin;
698  if (ix >= 1) xxm1 *= xin;
699  }
700 
701  // output position
702  double xout = 0, yout = 0;
703  const double *c = &_coeffs[0];
704  const double *pm = &monomials[0];
705  for (int k = _nterms; k--;) xout += (*(pm++)) * (*(c++));
706  pm = &monomials[0];
707  for (int k = _nterms; k--;) yout += (*(pm++)) * (*(c++));
708  res.x = xout;
709  res.y = yout;
710 
711  // derivatives
712  c = &_coeffs[0];
713  const double *mx = &dermx[0];
714  const double *my = &dermy[0];
715  double a11 = 0, a12 = 0;
716  for (int k = _nterms; k--;) {
717  a11 += (*(mx++)) * (*c);
718  a12 += (*(my++)) * (*(c++));
719  }
720 
721  double a21 = 0, a22 = 0;
722  mx = &dermx[0];
723  my = &dermy[0];
724  for (int k = _nterms; k--;) {
725  a21 += (*(mx++)) * (*c);
726  a22 += (*(my++)) * (*(c++));
727  }
728 
729  // output co-variance
730  res.vx = a11 * (a11 * in.vx + 2 * a12 * in.vxy) + a12 * a12 * in.vy;
731  res.vy = a21 * a21 * in.vx + a22 * a22 * in.vy + 2. * a21 * a22 * in.vxy;
732  res.vxy = a21 * a11 * in.vx + a22 * a12 * in.vy + (a21 * a12 + a11 * a22) * in.vxy;
733  out = res;
734 }
735 
736 /* The coefficient ordering is defined both here *AND* in the
737  AstrometryTransformPolynomial::apply, AstrometryTransformPolynomial::Derivative, ... routines
738  Change all or none ! */
739 
741  std::size_t whichCoord) const {
742  assert((degX + degY <= _order) && whichCoord < 2);
743  /* this assertion above is enough to ensure that the index used just
744  below is within bounds since the reserved length is
745  2*_nterms=(order+1)*(order+2) */
746  return _coeffs[(degX + degY) * (degX + degY + 1) / 2 + degY + whichCoord * _nterms];
747 }
748 
750  std::size_t whichCoord) {
751  assert((degX + degY <= _order) && whichCoord < 2);
752  return _coeffs[(degX + degY) * (degX + degY + 1) / 2 + degY + whichCoord * _nterms];
753 }
754 
756  std::size_t whichCoord) const {
757  // assert((degX+degY<=order) && whichCoord<2);
758  assert(whichCoord < 2);
759  if (degX + degY <= _order)
760  return _coeffs[(degX + degY) * (degX + degY + 1) / 2 + degY + whichCoord * _nterms];
761  return 0;
762 }
763 
764 /* parameter serialization for "virtual" fits */
765 double AstrometryTransformPolynomial::paramRef(Eigen::Index const i) const {
766  assert(i < 2 * Eigen::Index(_nterms));
767  return _coeffs[i];
768 }
769 
770 double &AstrometryTransformPolynomial::paramRef(Eigen::Index const i) {
771  assert(i < 2 * Eigen::Index(_nterms));
772  return _coeffs[i];
773 }
774 
775 void AstrometryTransformPolynomial::paramDerivatives(Point const &where, double *dx, double *dy)
776  const { /* first half : dxout/dpar, second half : dyout/dpar */
777  computeMonomials(where.x, where.y, dx);
778  for (std::size_t k = 0; k < _nterms; ++k) {
779  dy[_nterms + k] = dx[k];
780  dx[_nterms + k] = dy[k] = 0;
781  }
782 }
783 
784 /* utility for the print(ostream&) routine */
785 static string monomialString(std::size_t powX, std::size_t powY) {
786  stringstream ss;
787  if (powX + powY) ss << "*";
788  if (powX > 0) ss << "x";
789  if (powX > 1) ss << "^" << powX;
790  if (powY > 0) ss << "y";
791  if (powY > 1) ss << "^" << powY;
792  return ss.str();
793 }
794 
796  auto oldPrecision = stream.precision();
797  stream.precision(12);
798  stream << "AstrometryTransformPolynomial: order=" << getOrder() << std::endl;
799  for (std::size_t ic = 0; ic < 2; ++ic) {
800  if (ic == 0)
801  stream << "newx = ";
802  else
803  stream << "newy = ";
804  bool printed = false; // whether we've printed any monomials
805  for (std::size_t p = 0; p <= _order; ++p)
806  for (std::size_t py = 0; py <= p; ++py) {
807  double coefficient = getCoefficient(p - py, py, ic);
808  if (coefficient != 0 || isnan(coefficient)) {
809  // Only print "interesting" coefficients.
810  if (printed && p + py != 0) stream << " + ";
811  printed = true;
812  stream << coefficient << monomialString(p - py, py);
813  }
814  }
815  // if all components are zero, nothing is output by the loops above
816  if (!printed) {
817  stream << 0;
818  }
819  stream << endl;
820  }
821  if (_order > 0) stream << "Linear determinant = " << determinant();
822  stream.precision(oldPrecision);
823 }
824 
826  if (_order >= 1)
827  return getCoefficient(1, 0, 0) * getCoefficient(0, 1, 1) -
828  getCoefficient(0, 1, 0) * getCoefficient(1, 0, 1);
829  return 0;
830 }
831 
834  Point center = frame.getCenter();
835  return AstrometryTransformLinearScale(2. / frame.getWidth(), 2. / frame.getHeight()) *
836  AstrometryTransformLinearShift(-center.x, -center.y);
837 }
838 
839 /*utility for the AstrometryTransformPolynomial::fit() routine */
840 static AstrometryTransformLinear shiftAndNormalize(StarMatchList const &starMatchList) {
841  double xav = 0;
842  double x2 = 0;
843  double yav = 0;
844  double y2 = 0;
845  double count = 0;
846  for (auto it = starMatchList.begin(); it != starMatchList.end(); ++it) {
847  const StarMatch &a_match = *it;
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 
865 static double sq(double x) { return x * x; }
866 
867 double 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 (auto it = starMatchList.begin(); it != starMatchList.end(); ++it) {
877  const StarMatch &a_match = *it;
878  Point tmp = shiftToCenter.apply(a_match.point1);
879  FatPoint point1(tmp, a_match.point1.vx, a_match.point1.vy, a_match.point1.vxy);
880  FatPoint const &point2 = a_match.point2;
881  double wxx, wyy, wxy;
882  FatPoint tr1;
883  computeMonomials(point1.x, point1.y, monomials);
884  if (useErrors) {
885  transformPosAndErrors(point1, tr1); // we might consider recycling the monomials
886  double vxx = (tr1.vx + point2.vx);
887  double vyy = (tr1.vy + point2.vy);
888  double vxy = (tr1.vxy + point2.vxy);
889  double det = vxx * vyy - vxy * vxy;
890  wxx = vyy / det;
891  wyy = vxx / det;
892  wxy = -vxy / det;
893  } else {
894  wxx = wyy = 1;
895  wxy = 0;
896  apply(point1.x, point1.y, tr1.x, tr1.y);
897  }
898  double resx = point2.x - tr1.x;
899  double resy = point2.y - tr1.y;
900  sumr2 += wxx * sq(resx) + wyy * sq(resy) + 2 * wxy * resx * resy;
901 
902  double bxcoeff = wxx * resx + wxy * resy;
903  double bycoeff = wyy * resy + wxy * resx;
904  for (std::size_t j = 0; j < _nterms; ++j) {
905  for (std::size_t i = j; i < _nterms; ++i) {
906  A(i, j) += wxx * monomials[i] * monomials[j];
907  A(i + _nterms, j + _nterms) += wyy * monomials[i] * monomials[j];
908  A(j, i + _nterms) = A(i, j + _nterms) += wxy * monomials[i] * monomials[j];
909  }
910  B(j) += bxcoeff * monomials[j];
911  B(j + _nterms) += bycoeff * monomials[j];
912  }
913  } // end loop on points
914  Eigen::LDLT<Eigen::MatrixXd, Eigen::Lower> factor(A);
915  // should probably throw
916  if (factor.info() != Eigen::Success) {
917  LOGL_ERROR(_log, "AstrometryTransformPolynomial::fit could not factorize");
918  return -1;
919  }
920 
921  Eigen::VectorXd sol = factor.solve(B);
922  for (std::size_t k = 0; k < 2 * _nterms; ++k) _coeffs[k] += sol(k);
923  if (starMatchList.size() == _nterms) return 0;
924  return (sumr2 - B.dot(sol));
925 }
926 
927 double AstrometryTransformPolynomial::fit(StarMatchList const &starMatchList) {
928  if (starMatchList.size() < _nterms) {
929  LOGLS_FATAL(_log, "AstrometryTransformPolynomial::fit trying to fit a polynomial transform of order "
930  << _order << " with only " << starMatchList.size() << " matches.");
931  return -1;
932  }
933 
934  AstrometryTransformPolynomial conditionner = shiftAndNormalize(starMatchList);
935 
936  computeFit(starMatchList, conditionner, false); // get a rough solution
937  computeFit(starMatchList, conditionner, true); // weight with it
938  double chi2 = computeFit(starMatchList, conditionner, true); // once more
939 
940  (*this) = (*this) * conditionner;
941  if (starMatchList.size() == _nterms) return 0;
942  return chi2;
943 }
944 
946  AstrometryTransformPolynomial const &right) const {
947  if (getOrder() == 1 && right.getOrder() == 1)
948  return std::make_unique<AstrometryTransformLinear>((*this) * (right)); // does the composition
949  else
950  return std::make_unique<AstrometryTransformPolynomial>((*this) * (right)); // does the composition
951 }
952 
953 /* PolyXY the class used to perform polynomial algebra (and in
954  particular composition) at the coefficient level. This class
955  handles a single polynomial, while a AstrometryTransformPolynomial is a couple of
956  polynomials. This class does not have any routine to evaluate
957  polynomials. Efficiency is not a concern since these routines are
958  seldom used. There is no need to expose this tool class to
959  AstrometryTransform users.
960 */
961 
962 class PolyXY {
963  std::size_t order;
964  std::size_t nterms;
965  vector<long double> coeffs;
966 
967 public:
968  PolyXY(const int order) : order(order), nterms((order + 1) * (order + 2) / 2) {
969  coeffs.reserve(nterms);
970  coeffs.insert(coeffs.begin(), nterms, 0L); // fill & initialize to 0.
971  }
972 
973  std::size_t getOrder() const { return order; }
974 
976  : order(transform.getOrder()), nterms((order + 1) * (order + 2) / 2), coeffs(nterms, 0L) {
977  for (std::size_t px = 0; px <= order; ++px)
978  for (std::size_t py = 0; py <= order - px; ++py) {
979  getCoefficient(px, py) = transform.getCoefficient(px, py, whichCoord);
980  }
981  }
982 
983  long double getCoefficient(std::size_t powX, std::size_t powY) const {
984  assert(powX + powY <= order);
985  return coeffs.at((powX + powY) * (powX + powY + 1) / 2 + powY);
986  }
987 
988  long double &getCoefficient(std::size_t powX, std::size_t powY) {
989  assert(powX + powY <= order);
990  return coeffs.at((powX + powY) * (powX + powY + 1) / 2 + powY);
991  }
992 };
993 
994 /* ===================== PolyXY Algebra routines ================== */
995 
996 static void operator+=(PolyXY &left, const PolyXY &right) {
997  std::size_t rdeg = right.getOrder();
998  assert(left.getOrder() >= rdeg);
999  for (std::size_t i = 0; i <= rdeg; ++i)
1000  for (std::size_t j = 0; j <= rdeg - i; ++j) left.getCoefficient(i, j) += right.getCoefficient(i, j);
1001 }
1002 
1003 /* multiplication by a scalar */
1004 static PolyXY operator*(const long double &a, const PolyXY &polyXY) {
1005  PolyXY result(polyXY);
1006  // no direct access to coefficients: do it the soft way
1007  std::size_t order = polyXY.getOrder();
1008  for (std::size_t i = 0; i <= order; ++i)
1009  for (std::size_t j = 0; j <= order - i; ++j) result.getCoefficient(i, j) *= a;
1010  return result;
1011 }
1012 
1014 static PolyXY product(const PolyXY &p1, const PolyXY &p2) {
1015  std::size_t deg1 = p1.getOrder();
1016  std::size_t deg2 = p2.getOrder();
1017  PolyXY result(deg1 + deg2);
1018  for (std::size_t i1 = 0; i1 <= deg1; ++i1)
1019  for (std::size_t j1 = 0; j1 <= deg1 - i1; ++j1)
1020  for (std::size_t i2 = 0; i2 <= deg2; ++i2)
1021  for (std::size_t j2 = 0; j2 <= deg2 - i2; ++j2)
1022  result.getCoefficient(i1 + i2, j1 + j2) +=
1023  p1.getCoefficient(i1, j1) * p2.getCoefficient(i2, j2);
1024  return result;
1025 }
1026 
1027 /* powers[k](x,y) = polyXY(x,y)**k, 0 <= k <= maxP */
1028 static void computePowers(const PolyXY &polyXY, std::size_t maxP, vector<PolyXY> &powers) {
1029  powers.reserve(maxP + 1);
1030  powers.push_back(PolyXY(0));
1031  powers[0].getCoefficient(0, 0) = 1L;
1032  for (std::size_t k = 1; k <= maxP; ++k) powers.push_back(product(powers[k - 1], polyXY));
1033 }
1034 
1036 static PolyXY composition(const PolyXY &polyXY, const PolyXY &polyX, const PolyXY &polyY) {
1037  std::size_t pdeg = polyXY.getOrder();
1038  PolyXY result(pdeg * max(polyX.getOrder(), polyY.getOrder()));
1039  vector<PolyXY> pXPowers;
1040  vector<PolyXY> pYPowers;
1041  computePowers(polyX, pdeg, pXPowers);
1042  computePowers(polyY, pdeg, pYPowers);
1043  for (std::size_t px = 0; px <= pdeg; ++px)
1044  for (std::size_t py = 0; py <= pdeg - px; ++py)
1045  result += polyXY.getCoefficient(px, py) * product(pXPowers.at(px), pYPowers.at(py));
1046  return result;
1047 }
1048 
1049 /* ===================== end of PolyXY Algebra routines ============= */
1050 
1051 /* reducing polynomial composition is the reason for PolyXY stuff : */
1052 
1054  AstrometryTransformPolynomial const &right) const {
1055  // split each transform into 2d polynomials
1056  PolyXY plx(*this, 0);
1057  PolyXY ply(*this, 1);
1058  PolyXY prx(right, 0);
1059  PolyXY pry(right, 1);
1060 
1061  // compute the compositions
1062  PolyXY rx(composition(plx, prx, pry));
1063  PolyXY ry(composition(ply, prx, pry));
1064 
1065  // copy the results the hard way.
1066  AstrometryTransformPolynomial result(_order * right._order);
1067  for (std::size_t px = 0; px <= result._order; ++px)
1068  for (std::size_t py = 0; py <= result._order - px; ++py) {
1069  result.getCoefficient(px, py, 0) = rx.getCoefficient(px, py);
1070  result.getCoefficient(px, py, 1) = ry.getCoefficient(px, py);
1071  }
1072  return result;
1073 }
1074 
1076  AstrometryTransformPolynomial const &right) const {
1077  if (_order >= right._order) {
1078  AstrometryTransformPolynomial res(*this);
1079  for (std::size_t i = 0; i <= right._order; ++i)
1080  for (std::size_t j = 0; j <= right._order - i; ++j) {
1081  res.getCoefficient(i, j, 0) += right.getCoefficient(i, j, 0);
1082  res.getCoefficient(i, j, 1) += right.getCoefficient(i, j, 1);
1083  }
1084  return res;
1085  } else
1086  return (right + (*this));
1087 }
1088 
1090  AstrometryTransformPolynomial const &right) const {
1091  AstrometryTransformPolynomial res(std::max(_order, right._order));
1092  for (std::size_t i = 0; i <= res._order; ++i)
1093  for (std::size_t j = 0; j <= res._order - i; ++j) {
1094  res.getCoefficient(i, j, 0) = coeffOrZero(i, j, 0) - right.coeffOrZero(i, j, 0);
1095  res.getCoefficient(i, j, 1) = coeffOrZero(i, j, 1) - right.coeffOrZero(i, j, 1);
1096  }
1097  return res;
1098 }
1099 
1101  auto inverse = inversePolyTransform(*this, domain, 1e-7, _order + 2, 100);
1102  return std::make_shared<ast::PolyMap>(toAstPolyMapCoefficients(), inverse->toAstPolyMapCoefficients());
1103 }
1104 
1106  s << " AstrometryTransformPolynomial 1" << endl;
1107  s << "order " << _order << endl;
1108  int oldprec = s.precision();
1109  s << setprecision(12);
1110  for (std::size_t k = 0; k < 2 * _nterms; ++k) s << _coeffs[k] << ' ';
1111  s << endl;
1112  s << setprecision(oldprec);
1113 }
1114 
1116  int format;
1117  s >> format;
1118  if (format != 1)
1120  " AstrometryTransformPolynomial::read : format is not 1 ");
1121 
1122  string order;
1123  s >> order >> _order;
1124  if (order != "order")
1126  " AstrometryTransformPolynomial::read : expecting \"order\" and found " + order);
1127  setOrder(_order);
1128  for (std::size_t k = 0; k < 2 * _nterms; ++k) s >> _coeffs[k];
1129 }
1130 
1131 ndarray::Array<double, 2, 2> AstrometryTransformPolynomial::toAstPolyMapCoefficients() const {
1132  std::size_t nCoeffs = _coeffs.size();
1133  ndarray::Array<double, 2, 2> result = ndarray::allocate(ndarray::makeVector(nCoeffs, std::size_t(4)));
1134 
1135  ndarray::Size k = 0;
1136  for (std::size_t iCoord = 0; iCoord < 2; ++iCoord) {
1137  for (std::size_t p = 0; p <= _order; ++p) {
1138  for (std::size_t py = 0; py <= p; ++py, ++k) {
1139  result[k][0] = getCoefficient(p - py, py, iCoord);
1140  result[k][1] = iCoord + 1;
1141  result[k][2] = p - py;
1142  result[k][3] = py;
1143  }
1144  }
1145  }
1146 
1147  return result;
1148 }
1149 
1151  Frame const &domain,
1152  double const precision,
1153  std::size_t maxOrder,
1154  std::size_t nSteps) {
1155  StarMatchList sm;
1156  double xStart = domain.xMin;
1157  double yStart = domain.yMin;
1158  double xStep = domain.getWidth() / (nSteps - 1);
1159  double yStep = domain.getHeight() / (nSteps - 1);
1160  for (std::size_t i = 0; i < nSteps; ++i) {
1161  for (std::size_t j = 0; j < nSteps; ++j) {
1162  Point in(xStart + i * xStep, yStart + j * yStep);
1163  Point out(forward.apply(in));
1164  sm.push_back(StarMatch(out, in, nullptr, nullptr));
1165  }
1166  }
1167  std::size_t npairs = sm.size();
1171  double chi2 = 0;
1172  double oldChi2 = std::numeric_limits<double>::infinity();
1173  for (order = 1; order <= maxOrder; ++order) {
1175  auto success = poly->fit(sm);
1176  if (success == -1) {
1177  std::stringstream errMsg;
1178  errMsg << "Cannot fit a polynomial of order " << order << " with " << nSteps << "^2 points";
1179  throw pexExcept::RuntimeError(errMsg.str());
1180  }
1181  // compute the chi2 ignoring errors:
1182  chi2 = 0;
1183  for (auto const &i : sm) chi2 += i.point2.computeDist2(poly->apply((i.point1)));
1184  LOGLS_TRACE(_log, "inversePoly order " << order << ": " << chi2 << " / " << npairs << " = "
1185  << chi2 / npairs << " < " << precision * precision);
1186 
1187  if (chi2 / npairs < precision * precision) break;
1188 
1189  // If this triggers, we know we did not reach the required precision.
1190  if (chi2 > oldChi2) {
1191  LOGLS_WARN(_log, "inversePolyTransform: chi2 increases ("
1192  << chi2 << " > " << oldChi2 << "); ending fit with order: " << order);
1193  LOGLS_WARN(_log, "inversePolyTransform: requested precision not reached: "
1194  << chi2 << " / " << npairs << " = " << chi2 / npairs << " < "
1195  << precision * precision);
1196  poly = std::move(oldPoly);
1197  order--;
1198  break;
1199  } else {
1200  oldChi2 = chi2;
1201  // Clone it so we don't lose it in the next iteration.
1202  oldPoly = dynamic_pointer_cast<AstrometryTransformPolynomial>(
1204  }
1205  }
1206  if (order > maxOrder)
1207  LOGLS_WARN(_log, "inversePolyTransform: Reached max order without reaching requested precision: "
1208  << chi2 << " / " << npairs << " = " << chi2 / npairs << " < "
1209  << precision * precision);
1210  return poly;
1211 }
1212 
1213 /**************** AstrometryTransformLinear ***************************************/
1214 /* AstrometryTransformLinear is a specialized constructor of AstrometryTransformPolynomial
1215  May be it could just disappear ??
1216 */
1217 
1218 AstrometryTransformLinear::AstrometryTransformLinear(const double Dx, const double Dy, const double A11,
1219  const double A12, const double A21, const double A22)
1221  dx() = Dx;
1222  a11() = A11;
1223  a12() = A12;
1224  dy() = Dy;
1225  a21() = A21;
1226  a22() = A22;
1227 }
1228 
1231  if (transform.getOrder() != 1)
1233  "Trying to build a AstrometryTransformLinear from a higher order transform. Aborting. ");
1235 }
1236 
1238  // There is a general routine in AstrometryTransformPolynomial that would do the job:
1239  // return AstrometryTransformLinear(AstrometryTransformPolynomial::operator*(right));
1240  // however, we are using this composition of linear stuff heavily in
1241  // AstrometryTransform::linearApproximation, itself used in inverseTransform::apply.
1242  // So, for sake of efficiency, and since it is easy, we take a shortcut:
1244  apply(right.Dx(), right.Dy(), result.dx(), result.dy());
1245  result.a11() = A11() * right.A11() + A12() * right.A21();
1246  result.a12() = A11() * right.A12() + A12() * right.A22();
1247  result.a21() = A21() * right.A11() + A22() * right.A21();
1248  result.a22() = A21() * right.A12() + A22() * right.A22();
1249  return result;
1250 }
1251 
1253  const double) const {
1254  derivative = *this;
1255  derivative.getCoefficient(0, 0, 0) = 0;
1256  derivative.getCoefficient(0, 0, 1) = 0;
1257 }
1258 
1260  return *this;
1261 }
1262 
1264  //
1265  // (T1,M1) * (T2,M2) = 1 i.e (0,1) implies
1266  // T1 = -M1 * T2
1267  // M1 = M2^(-1)
1268  //
1269 
1270  double a11 = A11();
1271  double a12 = A12();
1272  double a21 = A21();
1273  double a22 = A22();
1274  double d = (a11 * a22 - a12 * a21);
1275  if (d == 0) {
1276  LOGL_FATAL(_log,
1277  "AstrometryTransformLinear::invert singular transformation: transform contents will be "
1278  "dumped to stderr.");
1279  print(cerr);
1280  }
1281 
1282  AstrometryTransformLinear result(0, 0, a22 / d, -a12 / d, -a21 / d, a11 / d);
1283  double rdx, rdy;
1284  result.apply(Dx(), Dy(), rdx, rdy);
1285  result.dx() = -rdx;
1286  result.dy() = -rdy;
1287  return result;
1288 }
1289 
1291  const Frame &) const {
1293 }
1294 
1296  auto oldPrecision = stream.precision();
1297  stream.precision(12);
1298  stream << "Linear AstrometryTransform:" << std::endl;
1299  stream << A11() << " " << A12() << " + " << Dx() << std::endl;
1300  stream << A21() << " " << A22() << " + " << Dy() << std::endl;
1301  stream.precision(oldPrecision);
1302  stream << "determinant = " << determinant();
1303  stream.precision(oldPrecision);
1304 }
1305 
1307  throw pexExcept::NotFoundError("AstrometryTransformLinearRot::fit not implemented! aborting");
1308 }
1309 
1311  std::size_t npairs = starMatchList.size();
1312  if (npairs < 3) {
1313  LOGLS_FATAL(_log, "AstrometryTransformLinearShift::fit trying to fit a linear transform with only "
1314  << npairs << " matches.");
1315  return -1;
1316  }
1317 
1318  double sumr2 = 0; /* used to compute chi2 without relooping */
1319  /* loop on pairs and fill */
1320  Eigen::VectorXd B(2);
1321  B.setZero();
1322  Eigen::MatrixXd A(2, 2);
1323  A.setZero();
1324 
1325  for (auto const &it : starMatchList) {
1326  FatPoint const &point1 = it.point1;
1327  FatPoint const &point2 = it.point2;
1328  double deltax = point2.x - point1.x;
1329  double deltay = point2.y - point1.y;
1330  double vxx = point1.vx + point2.vx;
1331  double vyy = point1.vy + point2.vy;
1332  double vxy = point1.vxy + point2.vxy;
1333  double det = vxx * vyy - vxy * vxy;
1334  double wxx = vyy / det;
1335  double wyy = vxx / det;
1336  double wxy = -vxy / det;
1337  B(0) += deltax * wxx + wxy * deltay;
1338  B(1) += deltay * wyy + wxy * deltax;
1339  A(0, 0) += wxx;
1340  A(1, 1) += wyy;
1341  A(0, 1) += wxy;
1342  sumr2 += deltax * deltax * wxx + deltay * deltay * wyy + 2. * wxy * deltax * deltay;
1343  }
1344  double det = A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
1345  if (det <= 0) return -1;
1346  double tmp = A(0, 0);
1347  A(0, 0) = A(1, 1) / det;
1348  A(1, 1) = tmp / det;
1349  A(0, 1) = A(1, 0) = -A(0, 1) / det;
1350  Eigen::VectorXd sol = A * B;
1351  (*this) = AstrometryTransformLinearShift(sol(0), sol(1));
1352  return (sumr2 - sol.dot(B)); // chi2 compact form
1353 }
1354 
1356  const double scaleFactor) {
1357  double c = scaleFactor * std::cos(angleRad);
1358  double s = scaleFactor * std::sin(angleRad);
1359  a11() = a22() = c;
1360  a21() = s;
1361  a12() = -s;
1362 
1363  // we want that the center does not move : transform+M*C = C ==> transform = C - M*C
1364  Point a_point(0., 0.);
1365  if (center) a_point = *center;
1366 
1367  dx() = dy() = 0;
1368  AstrometryTransformPolynomial::apply(a_point.x, a_point.y, dx(), dy()); // compute M*C
1369  dx() = a_point.x - Dx();
1370  dy() = a_point.y - dy();
1371 }
1372 
1373 static double deg2rad(double degree) { return degree * M_PI / 180.; }
1374 
1375 static double rad2deg(double rad) { return rad * 180. / M_PI; }
1376 
1377 /************* WCS transform ******************/
1378 /************** LinPixelToTan *******************/
1379 
1380 /* Implementation note : it seemed wise to incorporate
1381  the radians to degreess convertion into the linPixelToTan
1382  part (right in the constructor), and to do the
1383  opposite operation in the LinPart routine.
1384  When I was coding the fit, I realized that it was a
1385  bad idea. then I realized that the fitting routine
1386  itself was probably useless. Finaly, for a persistor,
1387  it seems bizarre that the stored data is different
1388  from what convention (angles in degrees for WCS's)
1389  would expect.
1390  So, no more "automatic" degrees to radians and
1391  radians to degrees conversion. They are explicitely
1392  done in apply (for TanPixelToRaDec and TanRaDecToPixel).
1393  This is a minor concern though....
1394 */
1395 BaseTanWcs::BaseTanWcs(AstrometryTransformLinear const &pixToTan, Point const &tangentPoint,
1396  const AstrometryTransformPolynomial *corrections) {
1397  // the angles returned by linPixelToTan should be in degrees.
1398  linPixelToTan = pixToTan;
1399  ra0 = deg2rad(tangentPoint.x);
1400  dec0 = deg2rad(tangentPoint.y);
1401  cos0 = std::cos(dec0);
1402  sin0 = std::sin(dec0);
1403  corr = nullptr;
1404  if (corrections) corr.reset(new AstrometryTransformPolynomial(*corrections));
1405 }
1406 
1407 /* with some sort of smart pointer ro handle "corr", we could remove the
1408  copy constructor, the operator = and the destructor */
1409 
1410 // ": AstrometryTransform" suppresses a warning
1412  corr = nullptr;
1413  *this = original;
1414 }
1415 
1416 void BaseTanWcs::operator=(const BaseTanWcs &original) {
1417  linPixelToTan = original.linPixelToTan;
1418  ra0 = original.ra0;
1419  dec0 = original.dec0;
1420  cos0 = std::cos(dec0);
1421  sin0 = std::sin(dec0);
1422  corr = nullptr;
1423  if (original.corr) corr.reset(new AstrometryTransformPolynomial(*original.corr));
1424 }
1425 
1426 void 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 
1452 Point BaseTanWcs::getTangentPoint() const { return Point(rad2deg(ra0), rad2deg(dec0)); }
1453 
1455 
1457  corr = std::move(corrections);
1458 }
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 
1473 
1475  : _skyWcs(skyWcs) {}
1476 
1477 void 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 
1487 double AstrometryTransformSkyWcs::fit(const StarMatchList &starMatchList) {
1488  throw LSST_EXCEPT(pex::exceptions::LogicError, "Not implemented");
1489 }
1490 
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 {
1509  return std::unique_ptr<AstrometryTransform>(nullptr);
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 
1530 }
1531 
1533  const Frame &region) const {
1534  if (!corr)
1537  else
1539 }
1540 
1542  if (corr)
1543  return (*corr) * linPixelToTan;
1544  else
1545  return linPixelToTan;
1546 }
1547 
1548 void 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 
1561  new TanPixelToRaDec(getLinPart(), getTangentPoint(), corr.get()));
1562 }
1563 
1564 void 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) {}
1602  return std::unique_ptr<AstrometryTransform>(new TanRaDecToPixel(getLinPart().inverted(),getTangentPoint()));
1603 }
1604 #endif
1605 
1607  const Frame &region)
1608  const { /* We have not implemented (yet) the reverse corrections available in SIP */
1610 }
1611 
1613  if (corr)
1614  return AstrometryTransformPolynomial(linPixelToTan) * (*corr);
1615  else
1616  return linPixelToTan;
1617 }
1618 
1619 void 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 
1633 }
1634 
1635 void TanSipPixelToRaDec::print(ostream &stream) const {
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(tan2Pix) {
1662  setTangentPoint(tangentPoint);
1663 }
1664 
1665 void 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 
1680 Point 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 
1737 void 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 
1758 }
1759 
1760 void 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 
1769 }
1770 
1774 }
1775 
1778 }
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 
1791 void UserTransform::apply(const double xIn, const double yIn, double &xOut, double &yOut) const {
1792  _userFun(xIn, yIn, xOut, yOut, _userData);
1793 }
1794 
1795 void 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 
1807 }
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 
1827  std::string type;
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
double x
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.
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.
void apply(const double xIn, const double yIn, double &xOut, double &yOut) const
return second(first(xIn,yIn))
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 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.
virtual AstrometryTransformLinear linearApproximation(Point const &where, const double step=0.01) const override
linear approximation.
void computeDerivative(Point const &where, AstrometryTransformLinear &derivative, const double step=0.01) const override
Computes the local Derivative of a transform, w.r.t.
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.
AstrometryTransformInverse(const AstrometryTransform *direct, const double precision, const Frame &region)
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 &out) const
prints the transform coefficients to stream.
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!
void apply(const double xIn, const double yIn, double &xOut, double &yOut) const
implements an iterative (Gauss-Newton) solver.
std::unique_ptr< AstrometryTransform > roughInverse(const Frame &) const
Overload the "generic routine".
implements the linear transformations (6 real coefficients).
AstrometryTransformLinear linearApproximation(Point const &where, const double step=0.01) const override
linear (local) approximation.
void computeDerivative(Point const &where, AstrometryTransformLinear &derivative, const double step=0.01) const override
specialised analytic routine
void print(std::ostream &out) const override
print out of coefficients in a readable form.
std::unique_ptr< AstrometryTransform > inverseTransform(const double precision, const Frame &region) const override
returns an inverse transform. Numerical if not overloaded.
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(const double xIn, const double yIn, double &xOut, double &yOut) const override
double fit(StarMatchList const &starMatchList)
guess what
just here to provide specialized constructors. AstrometryTransformLinear fit routine.
just here to provide a specialized constructor, and fit.
double fit(StarMatchList const &starMatchList)
guess what
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
void computeDerivative(Point const &where, AstrometryTransformLinear &derivative, const double step=0.01) const override
specialised analytic routine
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 const 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 std::unique_ptr< AstrometryTransform > composeAndReduce(AstrometryTransform const &right) const
Return a reduced composition of newTransform = this(right()), or nullptr if it cannot be reduced.
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 ).
virtual void apply(const double xIn, const double yIn, double &xOut, double &yOut) const=0
AstrometryTransformPolynomial operator+(AstrometryTransformPolynomial const &right) const
Addition.
void apply(const double xIn, const double yIn, double &xOut, double &yOut) const override
void print(std::ostream &out) const override
prints the transform coefficients to stream.
std::shared_ptr< afw::geom::SkyWcs > getSkyWcs() const
std::unique_ptr< AstrometryTransform > clone() const override
returns a copy (allocated by new) of the transformation.
virtual void apply(const double xIn, const double yIn, double &xOut, double &yOut) const=0
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(AstrometryTransformLinear const &pixToTan, Point const &tangentPoint, const AstrometryTransformPolynomial *corrections=nullptr)
AstrometryTransformLinear getLinPart() const
The Linear part (corresponding to CD's and CRPIX's)
void operator=(const BaseTanWcs &original)
Point getTangentPoint() const
Get the sky origin (CRVAL in FITS WCS terminology) in degrees.
virtual void apply(const double xIn, const double yIn, double &xOut, double &yOut) const=0
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:54
The transformation that handles pixels to sideral transformations (Gnomonic, possibly with polynomial...
std::unique_ptr< AstrometryTransform > inverseTransform(const double precision, const Frame &region) const
Inverse transform: returns a TanRaDecToPixel if there are no corrections, or the iterative solver if ...
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;
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.
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)
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.
TanPixelToRaDec inverted() const
exact typed inverse:
std::unique_ptr< AstrometryTransform > clone() const
returns a copy (allocated by new) of the transformation.
virtual void apply(const double xIn, const double yIn, double &xOut, double &yOut) const=0
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.
std::unique_ptr< AstrometryTransform > inverseTransform(const double precision, const Frame &region) const
Inverse transform: returns a TanPixelToRaDec.
void print(std::ostream &out) const
prints the transform coefficients to stream.
std::unique_ptr< AstrometryTransform > inverseTransform(const 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.
virtual void apply(const double xIn, const double yIn, double &xOut, double &yOut) const=0
double fit(StarMatchList const &starMatchList)
fits a transform to a std::list of Point pairs (p1,p2, the Point fields in StarMatch).
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 close(T... args)
T cos(T... args)
T count(T... args)
T emplace_back(T... args)
T endl(T... args)
T fabs(T... args)
T fail(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)
std::ostream & operator<<(std::ostream &os, CameraSysPrefix const &detSysPrefix)
Definition: CameraSys.cc:47
void write(OutputArchiveHandle &handle) const override
FilterProperty & operator=(FilterProperty const &)=default
Low-level polynomials (including special polynomials) in C++.
Definition: Basis1d.h:26
AstrometryTransformLinear normalizeCoordinatesTransform(const Frame &frame)
Returns the transformation that maps the input frame along both axes to [-1,1].
void() AstrometryTransformFun(const double, const double, double &, double &, const void *)
signature of the user-provided routine that actually does the coordinate transform for UserTransform.
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.
std::shared_ptr< AstrometryTransformPolynomial > inversePolyTransform(AstrometryTransform const &forward, Frame const &domain, double const precision, std::size_t maxOrder=9, std::size_t nSteps=50)
Approximate the inverse by a polynomial, to some precision.
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
A base class for image defects.
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)
T transform(T... args)
table::Key< int > order