LSSTApplications  17.0+124,17.0+14,17.0+73,18.0.0+37,18.0.0+80,18.0.0-4-g68ffd23+4,18.1.0-1-g0001055+12,18.1.0-1-g03d53ef+5,18.1.0-1-g1349e88+55,18.1.0-1-g2505f39+44,18.1.0-1-g5315e5e+4,18.1.0-1-g5e4b7ea+14,18.1.0-1-g7e8fceb+4,18.1.0-1-g85f8cd4+48,18.1.0-1-g8ff0b9f+4,18.1.0-1-ga2c679d+1,18.1.0-1-gd55f500+35,18.1.0-10-gb58edde+2,18.1.0-11-g0997b02+4,18.1.0-13-gfe4edf0b+12,18.1.0-14-g259bd21+21,18.1.0-19-gdb69f3f+2,18.1.0-2-g5f9922c+24,18.1.0-2-gd3b74e5+11,18.1.0-2-gfbf3545+32,18.1.0-26-g728bddb4+5,18.1.0-27-g6ff7ca9+2,18.1.0-3-g52aa583+25,18.1.0-3-g8ea57af+9,18.1.0-3-gb69f684+42,18.1.0-3-gfcaddf3+6,18.1.0-32-gd8786685a,18.1.0-4-gf3f9b77+6,18.1.0-5-g1dd662b+2,18.1.0-5-g6dbcb01+41,18.1.0-6-gae77429+3,18.1.0-7-g9d75d83+9,18.1.0-7-gae09a6d+30,18.1.0-9-gc381ef5+4,w.2019.45
LSSTDataManagementBasePackage
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 <fstream>
31 #include "assert.h"
32 #include <sstream>
33 
34 #include "Eigen/Core"
35 
36 #include "lsst/log/Log.h"
37 #include "lsst/geom/Point.h"
39 #include "lsst/jointcal/Frame.h"
41 #include "lsst/pex/exceptions.h"
42 #include "Eigen/Cholesky"
43 
45 
46 using namespace std;
47 
48 namespace {
49 LOG_LOGGER _log = LOG_GET("jointcal.AstrometryTransform");
50 }
51 
52 namespace lsst {
53 namespace jointcal {
54 
56  const AstrometryTransformPolynomial *shift =
57  dynamic_cast<const AstrometryTransformPolynomial *>(transform);
58  if (shift == nullptr) return false;
59 
60  static const double eps = 1e-5;
61 
62  double dx = shift->coeff(0, 0, 0);
63  double dy = shift->coeff(0, 0, 1);
64 
65  static Point dumb(4000, 4000);
66  if (fabs(dx - int(floor(dx + 0.5))) < eps && fabs(dy - int(floor(dy + 0.5))) < eps &&
67  fabs(dumb.x + dx - shift->apply(dumb).x) < eps && fabs(dumb.y + dy - shift->apply(dumb).y) < eps)
68  return true;
69 
70  return false;
71 }
72 
73 /********* AstrometryTransform ***********************/
74 
75 Frame AstrometryTransform::apply(Frame const &inputframe, bool inscribed) const {
76  // 2 opposite corners
77  double xtmin1, xtmax1, ytmin1, ytmax1;
78  apply(inputframe.xMin, inputframe.yMin, xtmin1, ytmin1);
79  apply(inputframe.xMax, inputframe.yMax, xtmax1, ytmax1);
80  Frame fr1(std::min(xtmin1, xtmax1), std::min(ytmin1, ytmax1), std::max(xtmin1, xtmax1),
81  std::max(ytmin1, ytmax1));
82  // 2 other corners
83  double xtmin2, xtmax2, ytmin2, ytmax2;
84  apply(inputframe.xMin, inputframe.yMax, xtmin2, ytmax2);
85  apply(inputframe.xMax, inputframe.yMin, xtmax2, ytmin2);
86  Frame fr2(std::min(xtmin2, xtmax2), std::min(ytmin2, ytmax2), std::max(xtmin2, xtmax2),
87  std::max(ytmin2, ytmax2));
88 
89  if (inscribed) return fr1 * fr2;
90  return fr1 + fr2;
91 }
92 
93 std::unique_ptr<AstrometryTransform> AstrometryTransform::composeAndReduce(
94  AstrometryTransform const &) const { // by default no way to compose
96 }
97 
98 double AstrometryTransform::getJacobian(const double x, const double y) const {
99  double x2, y2;
100  double eps = x * 0.01;
101  if (eps == 0) eps = 0.01;
102  apply(x, y, x2, y2);
103  double dxdx, dydx;
104  apply(x + eps, y, dxdx, dydx);
105  dxdx -= x2;
106  dydx -= y2;
107  double dxdy, dydy;
108  apply(x, y + eps, dxdy, dydy);
109  dxdy -= x2;
110  dydy -= y2;
111  return ((dxdx * dydy - dxdy * dydx) / (eps * eps));
112 }
113 
117 void AstrometryTransform::computeDerivative(Point const &where, AstrometryTransformLinear &derivative,
118  const double step) const {
119  double x = where.x;
120  double y = where.y;
121  double xp0, yp0;
122  apply(x, y, xp0, yp0);
123 
124  double xp, yp;
125  apply(x + step, y, xp, yp);
126  derivative.a11() = (xp - xp0) / step;
127  derivative.a21() = (yp - yp0) / step;
128  apply(x, y + step, xp, yp);
129  derivative.a12() = (xp - xp0) / step;
130  derivative.a22() = (yp - yp0) / step;
131  derivative.dx() = 0;
132  derivative.dy() = 0;
133 }
134 
135 AstrometryTransformLinear AstrometryTransform::linearApproximation(Point const &where,
136  const double step) const {
137  Point outwhere = apply(where);
139  computeDerivative(where, der, step);
140  return AstrometryTransformLinearShift(outwhere.x, outwhere.y) * der *
141  AstrometryTransformLinearShift(-where.x, -where.y);
142 }
143 
144 void AstrometryTransform::transformPosAndErrors(FatPoint const &in, FatPoint &out) const {
145  FatPoint res; // in case in and out are the same address...
146  res = apply(in);
148  // could save a call here, since Derivative needs the transform of where that we already have
149  // 0.01 may not be a very good idea in all cases. May be we should provide a way of altering that.
150  computeDerivative(in, der, 0.01);
151  double a11 = der.A11();
152  double a22 = der.A22();
153  double a21 = der.A21();
154  double a12 = der.A12();
155  res.vx = a11 * (a11 * in.vx + 2 * a12 * in.vxy) + a12 * a12 * in.vy;
156  res.vy = a21 * a21 * in.vx + a22 * a22 * in.vy + 2. * a21 * a22 * in.vxy;
157  res.vxy = a21 * a11 * in.vx + a22 * a12 * in.vy + (a21 * a12 + a11 * a22) * in.vxy;
158  out = res;
159 }
160 
161 void AstrometryTransform::transformErrors(Point const &where, const double *vIn, double *vOut) const {
163  computeDerivative(where, der, 0.01);
164  double a11 = der.A11();
165  double a22 = der.A22();
166  double a21 = der.A21();
167  double a12 = der.A12();
168 
169  /* (a11 a12) (vxx vxy)
170  M = ( ) and V = ( )
171  (a21 a22) (xvy vyy)
172 
173  Vxx = Vin[0], vyy = Vin[1], Vxy = Vin[2];
174  we want to compute M*V*tp(M)
175  A lin alg light package would be perfect...
176  */
177  int xx = 0;
178  int yy = 1;
179  int xy = 2;
180  // M*V :
181 
182  double b11 = a11 * vIn[xx] + a12 * vIn[xy];
183  double b22 = a21 * vIn[xy] + a22 * vIn[yy];
184  double b12 = a11 * vIn[xy] + a12 * vIn[yy];
185  double b21 = a21 * vIn[xx] + a22 * vIn[xy];
186 
187  // (M*V) * tp(M)
188 
189  vOut[xx] = b11 * a11 + b12 * a12;
190  vOut[xy] = b11 * a21 + b12 * a22;
191  vOut[yy] = b21 * a21 + b22 * a22;
192 }
193 
194 std::unique_ptr<AstrometryTransform> AstrometryTransform::roughInverse(const Frame &region) const {
195  // "in" and "out" refer to the inverse direction.
196  Point centerOut = region.getCenter();
197  Point centerIn = apply(centerOut);
199  computeDerivative(centerOut, der, std::sqrt(region.getArea()) / 5.);
200  der = der.inverted();
201  der = AstrometryTransformLinearShift(centerOut.x, centerOut.y) * der *
202  AstrometryTransformLinearShift(-centerIn.x, -centerIn.y);
204 }
205 
206 /* implement one in AstrometryTransform, so that all derived
207  classes do not need to provide one... */
208 
209 /* the routines that follow are used for ea generic parameter
210  transformation serialization, used e.g. for fits. Enables
211  to manipulate transformation parameters as vectors.
212 */
213 
214 // not dummy : what it does is virtual because paramRef is virtual.
215 void AstrometryTransform::getParams(double *params) const {
216  std::size_t npar = getNpar();
217  for (std::size_t i = 0; i < npar; ++i) params[i] = paramRef(i);
218 }
219 
220 void AstrometryTransform::offsetParams(Eigen::VectorXd const &delta) {
221  std::size_t npar = getNpar();
222  for (std::size_t i = 0; i < npar; ++i) paramRef(i) += delta[i];
223 }
224 
225 double AstrometryTransform::paramRef(Eigen::Index const) const {
227  std::string("AstrometryTransform::paramRef should never be called "));
228 }
229 
230 double &AstrometryTransform::paramRef(Eigen::Index const) {
232  "AstrometryTransform::paramRef should never be called ");
233 }
234 
235 void AstrometryTransform::paramDerivatives(Point const &, double *, double *) const {
237  "AstrometryTransform::paramDerivatives() should never be called ");
238 }
239 
241  transform.dump(stream);
242  return stream;
243 }
244 
245 void AstrometryTransform::write(const std::string &fileName) const {
246  ofstream s(fileName.c_str());
247  write(s);
248  bool ok = !s.fail();
249  s.close();
250  if (!ok)
252  "AstrometryTransform::write, something went wrong for file " + fileName);
253 }
254 
256  throw LSST_EXCEPT(
258  "AstrometryTransform::write(ostream), should never be called. MEans that it is missing in some "
259  "derived class ");
260 }
261 
262 /******************* GTransformInverse ****************/
263 /* inverse transformation, solved by iterations. Before using
264  it (probably via AstrometryTransform::inverseTransform), consider
265  seriously StarMatchList::inverseTransform */
267 private:
270  double precision2;
271 
272 public:
273  AstrometryTransformInverse(const AstrometryTransform *direct, const double precision,
274  const Frame &region);
275 
278  void apply(const double xIn, const double yIn, double &xOut, double &yOut) const;
279 
280  void dump(ostream &stream) const;
281 
282  double fit(StarMatchList const &starMatchList);
283 
285 
287 
289  std::unique_ptr<AstrometryTransform> roughInverse(const Frame &) const { return _direct->clone(); }
290 
293  return _direct->clone();
294  }
295 
297 
298 private:
299  void operator=(AstrometryTransformInverse const &);
300 };
301 
302 std::unique_ptr<AstrometryTransform> AstrometryTransform::inverseTransform(const double precision,
303  const Frame &region) const {
304  return std::unique_ptr<AstrometryTransform>(new AstrometryTransformInverse(this, precision, region));
305 }
306 
307 AstrometryTransformInverse::AstrometryTransformInverse(const AstrometryTransform *direct,
308  const double precision, const Frame &region) {
309  _direct = direct->clone();
310  _roughInverse = _direct->roughInverse(region);
311  precision2 = precision * precision;
312 }
313 
314 AstrometryTransformInverse::AstrometryTransformInverse(AstrometryTransformInverse const &model)
315  : AstrometryTransform() {
316  _direct = model._direct->clone();
317  _roughInverse = model._roughInverse->clone();
318  precision2 = model.precision2;
319 }
320 
322 
323 void AstrometryTransformInverse::operator=(AstrometryTransformInverse const &model) {
324  _direct = model._direct->clone();
325  _roughInverse = model._roughInverse->clone();
326  precision2 = model.precision2;
327 }
328 
329 void AstrometryTransformInverse::apply(const double xIn, const double yIn, double &xOut, double &yOut) const {
330  Point in(xIn, yIn);
331  Point outGuess = _roughInverse->apply(in);
332  AstrometryTransformLinear directDer, reverseDer;
333  int loop = 0;
334  int maxloop = 20;
335  double move2;
336  do {
337  loop++;
338  Point inGuess = _direct->apply(outGuess);
339  _direct->computeDerivative(outGuess, directDer);
340  reverseDer = directDer.inverted();
341  double xShift, yShift;
342  reverseDer.apply(xIn - inGuess.x, yIn - inGuess.y, xShift, yShift);
343  outGuess.x += xShift;
344  outGuess.y += yShift;
345  move2 = xShift * xShift + yShift * yShift;
346  } while ((move2 > precision2) && (loop < maxloop));
347  if (loop == maxloop) LOGLS_WARN(_log, "Problems applying AstrometryTransformInverse at " << in);
348  xOut = outGuess.x;
349  yOut = outGuess.y;
350 }
351 
353  stream << " AstrometryTransformInverse of :" << endl << *_direct << endl;
354 }
355 
358  "Cannot fit a AstrometryTransformInverse. Use StarMatchList::inverseTransform instead.");
359 }
360 
363 }
364 
365 /************* AstrometryTransformComposition **************/
366 
367 // This class was done to allow composition of AstrometryTransform's, without specifications of their types.
368 // does not need to be public. Invoked by compose(left,right)
369 
373 private:
374  std::unique_ptr<AstrometryTransform> _first, _second;
375 
376 public:
378 
380  void apply(const double xIn, const double yIn, double &xOut, double &yOut) const;
381  void dump(ostream &stream = cout) const;
382 
384  double fit(StarMatchList const &starMatchList);
385 
388 };
389 
391  AstrometryTransform const &first) {
392  _first = first.clone();
393  _second = second.clone();
394 }
395 
396 void AstrometryTransformComposition::apply(const double xIn, const double yIn, double &xOut,
397  double &yOut) const {
398  double xout, yout;
399  _first->apply(xIn, yIn, xout, yout);
400  _second->apply(xout, yout, xOut, yOut);
401 }
402 
404  _first->dump(stream);
405  _second->dump(stream);
406 }
407 
409  /* fits only one of them. could check that first can actually be fitted... */
410  return _first->fit(starMatchList);
411 }
412 
414  return std::make_unique<AstrometryTransformComposition>(*_second, *_first);
415 }
416 
418 
420  AstrometryTransformIdentity const &right) {
421  return left.clone();
422 }
423 
425  AstrometryTransform const &right) {
426  // Try to use the composeAndReduce method from left. If absent, AstrometryTransform::composeAndReduce
427  // returns NULL. composeAndReduce is non trivial for polynomials.
429  // composition == NULL means no reduction: just build a Composition that pipelines "left" and "right".
430  if (composition == nullptr)
431  return std::make_unique<AstrometryTransformComposition>(left, right);
432  else
433  return composition;
434 }
435 
436 // just a speed up, to avoid useless numerical derivation.
438  const double) const {
439  derivative = AstrometryTransformLinear();
440 }
441 
443  const double) const {
445  return result; // rely on default AstrometryTransformlin constructor;
446 }
447 
449  return std::make_shared<ast::UnitMap>(2); // a AstrometryTransformIdentity is identically ast::UnitMap(2)
450 }
451 
453  stream << "AstrometryTransformIdentity 1" << endl;
454 }
455 
457  int format;
458  stream >> format;
459  if (format != 1)
461  " AstrometryTransformIdentity::read : format is not 1 ");
462 }
463 
464 /*************** AstrometryTransformPolynomial **************************************/
465 
467 
469  _nterms = (order + 1) * (order + 2) / 2;
470 
471  // allocate and fill coefficients
472  _coeffs.resize(2 * _nterms, 0.);
473  // the default is supposed to be the identity, (for order>=1).
474  if (_order >= 1) {
475  coeff(1, 0, 0) = 1;
476  coeff(0, 1, 1) = 1;
477  }
478 }
479 
480 //#ifdef TO_BE_FIXED
482  const Frame &frame, std::size_t order,
483  std::size_t nPoint) {
484  StarMatchList sm;
485 
486  double step = std::sqrt(fabs(frame.getArea()) / double(nPoint));
487  for (double x = frame.xMin + step / 2; x <= frame.xMax; x += step)
488  for (double y = frame.yMin + step / 2; y <= frame.yMax; y += step) {
489  auto pix = std::make_shared<BaseStar>(x, y, 0, 0);
490  double xtr, ytr;
491  transform->apply(x, y, xtr, ytr);
492  auto tp = std::make_shared<BaseStar>(xtr, ytr, 0, 0);
493  /* These are fake stars so no need to transform fake errors.
494  all errors (and weights) will be equal : */
495  sm.push_back(StarMatch(*pix, *tp, pix, tp));
496  }
498  ret.fit(sm);
499  *this = ret;
500 }
501 //#endif
502 
505  std::size_t order, std::size_t nSteps) {
506  jointcal::StarMatchList starMatchList;
507  double xStart = domain.xMin;
508  double yStart = domain.yMin;
509  double xStep = domain.getWidth() / (nSteps + 1);
510  double yStep = domain.getHeight() / (nSteps + 1);
511  for (std::size_t i = 0; i < nSteps; ++i) {
512  for (std::size_t j = 0; j < nSteps; ++j) {
513  // TODO: once DM-4044 is done, we can remove the redundancy in `Point`/`Point2D` here
514  jointcal::Point in(xStart + i * xStep, yStart + j * yStep);
515  geom::Point2D inAfw(in.x, in.y);
516  geom::Point2D outAfw = transform->applyForward(inAfw);
517  jointcal::Point out(outAfw.getX(), outAfw.getY());
518  starMatchList.emplace_back(in, out, nullptr, nullptr);
519  }
520  }
522  poly.fit(starMatchList);
523  *this = poly;
524 }
525 
526 void AstrometryTransformPolynomial::computeMonomials(double xIn, double yIn, double *monomial) const {
527  /* The ordering of monomials is implemented here.
528  You may not change it without updating the "mapping" routines
529  coeff(std::size_t, std::size_t, std::size_t).
530  I (P.A.) did not find a clever way to loop over monomials.
531  Improvements welcome.
532  This routine is used also by the fit to fill monomials.
533  We could certainly be more elegant.
534  */
535 
536  double xx = 1;
537  for (std::size_t ix = 0; ix <= _order; ++ix) {
538  double yy = 1;
539  std::size_t k = ix * (ix + 1) / 2;
540  for (std::size_t iy = 0; iy <= _order - ix; ++iy) {
541  monomial[k] = xx * yy;
542  yy *= yIn;
543  k += ix + iy + 2;
544  }
545  xx *= xIn;
546  }
547 }
548 
550  _order = order;
551  std::size_t old_nterms = _nterms;
552  _nterms = (_order + 1) * (_order + 2) / 2;
553 
554  // temporarily save coefficients
555  vector<double> old_coeffs = _coeffs;
556  // reallocate enough size
557  _coeffs.resize(2 * _nterms);
558  // reassign to zero (this is necessary because ycoeffs
559  // are after xcoeffs and so their meaning changes
560  for (std::size_t k = 0; k < _nterms; ++k) _coeffs[k] = 0;
561  // put back what we had before
562  std::size_t kmax = min(old_nterms, _nterms);
563  for (std::size_t k = 0; k < kmax; ++k) {
564  _coeffs[k] = old_coeffs[k]; // x terms
565  _coeffs[k + _nterms] = old_coeffs[k + old_nterms]; // y terms
566  }
567 }
568 
569 /* this is reasonably fast, when optimized */
570 void AstrometryTransformPolynomial::apply(const double xIn, const double yIn, double &xOut,
571  double &yOut) const {
572  /*
573  This routine computes the monomials only once for both
574  polynomials. This is why AstrometryTransformPolynomial does not use an auxilary
575  class (such as PolyXY) to handle each polynomial.
576 
577  The code works even if &xIn == &xOut (or &yIn == &yOut)
578  It uses Variable Length Allocation (VLA) rather than a vector<double>
579  because allocating the later costs about 50 ns. All VLA uses are tagged.
580  */
581  double monomials[_nterms]; // this is VLA, which is (perhaps) not casher C++
582  computeMonomials(xIn, yIn, monomials);
583 
584  xOut = 0;
585  yOut = 0;
586  const double *c = &_coeffs[0];
587  const double *pm = &monomials[0];
588  // the ordering of the coefficients and the monomials are identical.
589  for (int k = _nterms; k--;) xOut += (*(pm++)) * (*(c++));
590  pm = &monomials[0];
591  for (int k = _nterms; k--;) yOut += (*(pm++)) * (*(c++));
592 }
593 
595  AstrometryTransformLinear &derivative,
596  const double step)
597  const { /* routine checked against numerical derivatives from AstrometryTransform::Derivative */
598  if (_order == 1) {
599  derivative = AstrometryTransformLinear(*this);
600  derivative.dx() = derivative.dy() = 0;
601  return;
602  }
603 
604  double dermx[2 * _nterms]; // VLA
605  double *dermy = dermx + _nterms;
606  double xin = where.x;
607  double yin = where.y;
608 
609  double xx = 1;
610  double xxm1 = 1; // xx^(ix-1)
611  for (std::size_t ix = 0; ix <= _order; ++ix) {
612  std::size_t k = (ix) * (ix + 1) / 2;
613  // iy = 0
614  dermx[k] = ix * xxm1;
615  dermy[k] = 0;
616  k += ix + 2;
617  double yym1 = 1; // yy^(iy-1)
618  for (std::size_t iy = 1; iy <= _order - ix; ++iy) {
619  dermx[k] = ix * xxm1 * yym1 * yin;
620  dermy[k] = iy * xx * yym1;
621  yym1 *= yin;
622  k += ix + iy + 2;
623  }
624  xx *= xin;
625  if (ix >= 1) xxm1 *= xin;
626  }
627 
628  derivative.dx() = 0;
629  derivative.dy() = 0;
630 
631  const double *mx = &dermx[0];
632  const double *my = &dermy[0];
633  const double *c = &_coeffs[0];
634  // dx'
635  double a11 = 0, a12 = 0;
636  for (int k = _nterms; k--;) {
637  a11 += (*(mx++)) * (*c);
638  a12 += (*(my++)) * (*(c++));
639  }
640  derivative.a11() = a11;
641  derivative.a12() = a12;
642  // dy'
643  double a21 = 0, a22 = 0;
644  mx = &dermx[0];
645  my = &dermy[0];
646  for (int k = _nterms; k--;) {
647  a21 += (*(mx++)) * (*c);
648  a22 += (*(my++)) * (*(c++));
649  }
650  derivative.a21() = a21;
651  derivative.a22() = a22;
652 }
653 
655  /*
656  The results from this routine were compared to what comes out
657  from apply and transformErrors. The Derivative routine was
658  checked against numerical derivatives from
659  AstrometryTransform::Derivative. (P.A dec 2009).
660 
661  This routine could be made much simpler by calling apply and
662  Derivative (i.e. you just suppress it, and the fallback is the
663  generic version in AstrometryTransform). BTW, I checked that both routines
664  provide the same result. This version is however faster
665  (monomials get recycled).
666  */
667  double monomials[_nterms]; // VLA
668 
669  FatPoint res; // to store the result, because nothing forbids &in == &out.
670 
671  double dermx[2 * _nterms]; // monomials for derivative w.r.t. x (VLA)
672  double *dermy = dermx + _nterms; // same for y
673  double xin = in.x;
674  double yin = in.y;
675 
676  double xx = 1;
677  double xxm1 = 1; // xx^(ix-1)
678  for (std::size_t ix = 0; ix <= _order; ++ix) {
679  std::size_t k = (ix) * (ix + 1) / 2;
680  // iy = 0
681  dermx[k] = ix * xxm1;
682  dermy[k] = 0;
683  monomials[k] = xx;
684  k += ix + 2;
685  double yy = yin;
686  double yym1 = 1; // yy^(iy-1)
687  for (std::size_t iy = 1; iy <= _order - ix; ++iy) {
688  monomials[k] = xx * yy;
689  dermx[k] = ix * xxm1 * yy;
690  dermy[k] = iy * xx * yym1;
691  yym1 *= yin;
692  yy *= yin;
693  k += ix + iy + 2;
694  }
695  xx *= xin;
696  if (ix >= 1) xxm1 *= xin;
697  }
698 
699  // output position
700  double xout = 0, yout = 0;
701  const double *c = &_coeffs[0];
702  const double *pm = &monomials[0];
703  for (int k = _nterms; k--;) xout += (*(pm++)) * (*(c++));
704  pm = &monomials[0];
705  for (int k = _nterms; k--;) yout += (*(pm++)) * (*(c++));
706  res.x = xout;
707  res.y = yout;
708 
709  // derivatives
710  c = &_coeffs[0];
711  const double *mx = &dermx[0];
712  const double *my = &dermy[0];
713  double a11 = 0, a12 = 0;
714  for (int k = _nterms; k--;) {
715  a11 += (*(mx++)) * (*c);
716  a12 += (*(my++)) * (*(c++));
717  }
718 
719  double a21 = 0, a22 = 0;
720  mx = &dermx[0];
721  my = &dermy[0];
722  for (int k = _nterms; k--;) {
723  a21 += (*(mx++)) * (*c);
724  a22 += (*(my++)) * (*(c++));
725  }
726 
727  // output co-variance
728  res.vx = a11 * (a11 * in.vx + 2 * a12 * in.vxy) + a12 * a12 * in.vy;
729  res.vy = a21 * a21 * in.vx + a22 * a22 * in.vy + 2. * a21 * a22 * in.vxy;
730  res.vxy = a21 * a11 * in.vx + a22 * a12 * in.vy + (a21 * a12 + a11 * a22) * in.vxy;
731  out = res;
732 }
733 
734 /* The coefficient ordering is defined both here *AND* in the
735  AstrometryTransformPolynomial::apply, AstrometryTransformPolynomial::Derivative, ... routines
736  Change all or none ! */
737 
739  std::size_t whichCoord) const {
740  assert((degX + degY <= _order) && whichCoord < 2);
741  /* this assertion above is enough to ensure that the index used just
742  below is within bounds since the reserved length is
743  2*_nterms=(order+1)*(order+2) */
744  return _coeffs[(degX + degY) * (degX + degY + 1) / 2 + degY + whichCoord * _nterms];
745 }
746 
748  std::size_t whichCoord) {
749  assert((degX + degY <= _order) && whichCoord < 2);
750  return _coeffs[(degX + degY) * (degX + degY + 1) / 2 + degY + whichCoord * _nterms];
751 }
752 
754  std::size_t whichCoord) const {
755  // assert((degX+degY<=order) && whichCoord<2);
756  assert(whichCoord < 2);
757  if (degX + degY <= _order)
758  return _coeffs[(degX + degY) * (degX + degY + 1) / 2 + degY + whichCoord * _nterms];
759  return 0;
760 }
761 
762 /* parameter serialization for "virtual" fits */
763 double AstrometryTransformPolynomial::paramRef(Eigen::Index const i) const {
764  assert(i < 2 * Eigen::Index(_nterms));
765  return _coeffs[i];
766 }
767 
768 double &AstrometryTransformPolynomial::paramRef(Eigen::Index const i) {
769  assert(i < 2 * Eigen::Index(_nterms));
770  return _coeffs[i];
771 }
772 
773 void AstrometryTransformPolynomial::paramDerivatives(Point const &where, double *dx, double *dy)
774  const { /* first half : dxout/dpar, second half : dyout/dpar */
775  computeMonomials(where.x, where.y, dx);
776  for (std::size_t k = 0; k < _nterms; ++k) {
777  dy[_nterms + k] = dx[k];
778  dx[_nterms + k] = dy[k] = 0;
779  }
780 }
781 
782 /* utility for the dump(ostream&) routine */
783 static string monomialString(std::size_t powX, std::size_t powY) {
784  stringstream ss;
785  if (powX + powY) ss << "*";
786  if (powX > 0) ss << "x";
787  if (powX > 1) ss << "^" << powX;
788  if (powY > 0) ss << "y";
789  if (powY > 1) ss << "^" << powY;
790  return ss.str();
791 }
792 
794  auto oldPrecision = stream.precision();
795  stream.precision(12);
796  for (std::size_t ic = 0; ic < 2; ++ic) {
797  if (ic == 0)
798  stream << "newx = ";
799  else
800  stream << "newy = ";
801  for (std::size_t p = 0; p <= _order; ++p)
802  for (std::size_t py = 0; py <= p; ++py) {
803  if (p + py != 0) stream << " + ";
804  stream << coeff(p - py, py, ic) << monomialString(p - py, py);
805  }
806  stream << endl;
807  }
808  if (_order > 0) stream << " Linear determinant = " << determinant() << endl;
809  stream.precision(oldPrecision);
810 }
811 
813  if (_order >= 1) return coeff(1, 0, 0) * coeff(0, 1, 1) - coeff(0, 1, 0) * coeff(1, 0, 1);
814  return 0;
815 }
816 
819  Point center = frame.getCenter();
820  return AstrometryTransformLinearScale(2. / frame.getWidth(), 2. / frame.getHeight()) *
821  AstrometryTransformLinearShift(-center.x, -center.y);
822 }
823 
824 /*utility for the AstrometryTransformPolynomial::fit() routine */
825 static AstrometryTransformLinear shiftAndNormalize(StarMatchList const &starMatchList) {
826  double xav = 0;
827  double x2 = 0;
828  double yav = 0;
829  double y2 = 0;
830  double count = 0;
831  for (auto it = starMatchList.begin(); it != starMatchList.end(); ++it) {
832  const StarMatch &a_match = *it;
833  Point const &point1 = a_match.point1;
834  xav += point1.x;
835  yav += point1.y;
836  x2 += std::pow(point1.x, 2);
837  y2 += std::pow(point1.y, 2);
838  count++;
839  }
840  if (count == 0) return AstrometryTransformLinear();
841  xav /= count;
842  yav /= count;
843  // 3.5 stands for sqrt(12).
844  double xspan = 3.5 * std::sqrt(x2 / count - std::pow(xav, 2));
845  double yspan = 3.5 * std::sqrt(y2 / count - std::pow(yav, 2));
846  return AstrometryTransformLinearScale(2. / xspan, 2. / yspan) *
847  AstrometryTransformLinearShift(-xav, -yav);
848 }
849 
850 static double sq(double x) { return x * x; }
851 
852 double AstrometryTransformPolynomial::computeFit(StarMatchList const &starMatchList,
853  AstrometryTransform const &shiftToCenter,
854  const bool useErrors) {
855  Eigen::MatrixXd A(2 * _nterms, 2 * _nterms);
856  A.setZero();
857  Eigen::VectorXd B(2 * _nterms);
858  B.setZero();
859  double sumr2 = 0;
860  double monomials[_nterms];
861  for (auto it = starMatchList.begin(); it != starMatchList.end(); ++it) {
862  const StarMatch &a_match = *it;
863  Point tmp = shiftToCenter.apply(a_match.point1);
864  FatPoint point1(tmp, a_match.point1.vx, a_match.point1.vy, a_match.point1.vxy);
865  FatPoint const &point2 = a_match.point2;
866  double wxx, wyy, wxy;
867  FatPoint tr1;
868  computeMonomials(point1.x, point1.y, monomials);
869  if (useErrors) {
870  transformPosAndErrors(point1, tr1); // we might consider recycling the monomials
871  double vxx = (tr1.vx + point2.vx);
872  double vyy = (tr1.vy + point2.vy);
873  double vxy = (tr1.vxy + point2.vxy);
874  double det = vxx * vyy - vxy * vxy;
875  wxx = vyy / det;
876  wyy = vxx / det;
877  wxy = -vxy / det;
878  } else {
879  wxx = wyy = 1;
880  wxy = 0;
881  apply(point1.x, point1.y, tr1.x, tr1.y);
882  }
883  double resx = point2.x - tr1.x;
884  double resy = point2.y - tr1.y;
885  sumr2 += wxx * sq(resx) + wyy * sq(resy) + 2 * wxy * resx * resy;
886 
887  double bxcoeff = wxx * resx + wxy * resy;
888  double bycoeff = wyy * resy + wxy * resx;
889  for (std::size_t j = 0; j < _nterms; ++j) {
890  for (std::size_t i = j; i < _nterms; ++i) {
891  A(i, j) += wxx * monomials[i] * monomials[j];
892  A(i + _nterms, j + _nterms) += wyy * monomials[i] * monomials[j];
893  A(j, i + _nterms) = A(i, j + _nterms) += wxy * monomials[i] * monomials[j];
894  }
895  B(j) += bxcoeff * monomials[j];
896  B(j + _nterms) += bycoeff * monomials[j];
897  }
898  } // end loop on points
899  Eigen::LDLT<Eigen::MatrixXd, Eigen::Lower> factor(A);
900  // should probably throw
901  if (factor.info() != Eigen::Success) {
902  LOGL_ERROR(_log, "AstrometryTransformPolynomial::fit could not factorize");
903  return -1;
904  }
905 
906  Eigen::VectorXd sol = factor.solve(B);
907  for (std::size_t k = 0; k < 2 * _nterms; ++k) _coeffs[k] += sol(k);
908  if (starMatchList.size() == _nterms) return 0;
909  return (sumr2 - B.dot(sol));
910 }
911 
912 double AstrometryTransformPolynomial::fit(StarMatchList const &starMatchList) {
913  if (starMatchList.size() < _nterms) {
914  LOGLS_FATAL(_log, "AstrometryTransformPolynomial::fit trying to fit a polynomial transform of order "
915  << _order << " with only " << starMatchList.size() << " matches.");
916  return -1;
917  }
918 
919  AstrometryTransformPolynomial conditionner = shiftAndNormalize(starMatchList);
920 
921  computeFit(starMatchList, conditionner, false); // get a rough solution
922  computeFit(starMatchList, conditionner, true); // weight with it
923  double chi2 = computeFit(starMatchList, conditionner, true); // once more
924 
925  (*this) = (*this) * conditionner;
926  if (starMatchList.size() == _nterms) return 0;
927  return chi2;
928 }
929 
931  AstrometryTransformPolynomial const &right) const {
932  if (getOrder() == 1 && right.getOrder() == 1)
933  return std::make_unique<AstrometryTransformLinear>((*this) * (right)); // does the composition
934  else
935  return std::make_unique<AstrometryTransformPolynomial>((*this) * (right)); // does the composition
936 }
937 
938 /* PolyXY the class used to perform polynomial algebra (and in
939  particular composition) at the coefficient level. This class
940  handles a single polynomial, while a AstrometryTransformPolynomial is a couple of
941  polynomials. This class does not have any routine to evaluate
942  polynomials. Efficiency is not a concern since these routines are
943  seldom used. There is no need to expose this tool class to
944  AstrometryTransform users.
945 */
946 
947 class PolyXY {
948  std::size_t order;
949  std::size_t nterms;
950  vector<long double> coeffs;
951 
952 public:
953  PolyXY(const int order) : order(order), nterms((order + 1) * (order + 2) / 2) {
954  coeffs.reserve(nterms);
955  coeffs.insert(coeffs.begin(), nterms, 0L); // fill & initialize to 0.
956  }
957 
958  std::size_t getOrder() const { return order; }
959 
961  : order(transform.getOrder()), nterms((order + 1) * (order + 2) / 2), coeffs(nterms, 0L) {
962  for (std::size_t px = 0; px <= order; ++px)
963  for (std::size_t py = 0; py <= order - px; ++py) coeff(px, py) = transform.coeff(px, py, whichCoord);
964  }
965 
966  long double coeff(std::size_t powX, std::size_t powY) const {
967  assert(powX + powY <= order);
968  return coeffs.at((powX + powY) * (powX + powY + 1) / 2 + powY);
969  }
970 
971  long double &coeff(std::size_t powX, std::size_t powY) {
972  assert(powX + powY <= order);
973  return coeffs.at((powX + powY) * (powX + powY + 1) / 2 + powY);
974  }
975 };
976 
977 /* ===================== PolyXY Algebra routines ================== */
978 
979 static void operator+=(PolyXY &left, const PolyXY &right) {
980  std::size_t rdeg = right.getOrder();
981  assert(left.getOrder() >= rdeg);
982  for (std::size_t i = 0; i <= rdeg; ++i)
983  for (std::size_t j = 0; j <= rdeg - i; ++j) left.coeff(i, j) += right.coeff(i, j);
984 }
985 
986 /* multiplication by a scalar */
987 static PolyXY operator*(const long double &a, const PolyXY &polyXY) {
988  PolyXY result(polyXY);
989  // no direct access to coefficients: do it the soft way
990  std::size_t order = polyXY.getOrder();
991  for (std::size_t i = 0; i <= order; ++i)
992  for (std::size_t j = 0; j <= order - i; ++j) result.coeff(i, j) *= a;
993  return result;
994 }
995 
997 static PolyXY product(const PolyXY &p1, const PolyXY &p2) {
998  std::size_t deg1 = p1.getOrder();
999  std::size_t deg2 = p2.getOrder();
1000  PolyXY result(deg1 + deg2);
1001  for (std::size_t i1 = 0; i1 <= deg1; ++i1)
1002  for (std::size_t j1 = 0; j1 <= deg1 - i1; ++j1)
1003  for (std::size_t i2 = 0; i2 <= deg2; ++i2)
1004  for (std::size_t j2 = 0; j2 <= deg2 - i2; ++j2)
1005  result.coeff(i1 + i2, j1 + j2) += p1.coeff(i1, j1) * p2.coeff(i2, j2);
1006  return result;
1007 }
1008 
1009 /* powers[k](x,y) = polyXY(x,y)**k, 0 <= k <= maxP */
1010 static void computePowers(const PolyXY &polyXY, std::size_t maxP, vector<PolyXY> &powers) {
1011  powers.reserve(maxP + 1);
1012  powers.push_back(PolyXY(0));
1013  powers[0].coeff(0, 0) = 1L;
1014  for (std::size_t k = 1; k <= maxP; ++k) powers.push_back(product(powers[k - 1], polyXY));
1015 }
1016 
1018 static PolyXY composition(const PolyXY &polyXY, const PolyXY &polyX, const PolyXY &polyY) {
1019  std::size_t pdeg = polyXY.getOrder();
1020  PolyXY result(pdeg * max(polyX.getOrder(), polyY.getOrder()));
1021  vector<PolyXY> pXPowers;
1022  vector<PolyXY> pYPowers;
1023  computePowers(polyX, pdeg, pXPowers);
1024  computePowers(polyY, pdeg, pYPowers);
1025  for (std::size_t px = 0; px <= pdeg; ++px)
1026  for (std::size_t py = 0; py <= pdeg - px; ++py)
1027  result += polyXY.coeff(px, py) * product(pXPowers.at(px), pYPowers.at(py));
1028  return result;
1029 }
1030 
1031 /* ===================== end of PolyXY Algebra routines ============= */
1032 
1033 /* reducing polynomial composition is the reason for PolyXY stuff : */
1034 
1036  AstrometryTransformPolynomial const &right) const {
1037  // split each transform into 2d polynomials
1038  PolyXY plx(*this, 0);
1039  PolyXY ply(*this, 1);
1040  PolyXY prx(right, 0);
1041  PolyXY pry(right, 1);
1042 
1043  // compute the compositions
1044  PolyXY rx(composition(plx, prx, pry));
1045  PolyXY ry(composition(ply, prx, pry));
1046 
1047  // copy the results the hard way.
1048  AstrometryTransformPolynomial result(_order * right._order);
1049  for (std::size_t px = 0; px <= result._order; ++px)
1050  for (std::size_t py = 0; py <= result._order - px; ++py) {
1051  result.coeff(px, py, 0) = rx.coeff(px, py);
1052  result.coeff(px, py, 1) = ry.coeff(px, py);
1053  }
1054  return result;
1055 }
1056 
1058  AstrometryTransformPolynomial const &right) const {
1059  if (_order >= right._order) {
1060  AstrometryTransformPolynomial res(*this);
1061  for (std::size_t i = 0; i <= right._order; ++i)
1062  for (std::size_t j = 0; j <= right._order - i; ++j) {
1063  res.coeff(i, j, 0) += right.coeff(i, j, 0);
1064  res.coeff(i, j, 1) += right.coeff(i, j, 1);
1065  }
1066  return res;
1067  } else
1068  return (right + (*this));
1069 }
1070 
1072  AstrometryTransformPolynomial const &right) const {
1073  AstrometryTransformPolynomial res(std::max(_order, right._order));
1074  for (std::size_t i = 0; i <= res._order; ++i)
1075  for (std::size_t j = 0; j <= res._order - i; ++j) {
1076  res.coeff(i, j, 0) = coeffOrZero(i, j, 0) - right.coeffOrZero(i, j, 0);
1077  res.coeff(i, j, 1) = coeffOrZero(i, j, 1) - right.coeffOrZero(i, j, 1);
1078  }
1079  return res;
1080 }
1081 
1083  auto inverse = inversePolyTransform(*this, domain, 1e-7, _order + 2, 100);
1084  return std::make_shared<ast::PolyMap>(toAstPolyMapCoefficients(), inverse->toAstPolyMapCoefficients());
1085 }
1086 
1088  s << " AstrometryTransformPolynomial 1" << endl;
1089  s << "order " << _order << endl;
1090  int oldprec = s.precision();
1091  s << setprecision(12);
1092  for (std::size_t k = 0; k < 2 * _nterms; ++k) s << _coeffs[k] << ' ';
1093  s << endl;
1094  s << setprecision(oldprec);
1095 }
1096 
1098  int format;
1099  s >> format;
1100  if (format != 1)
1102  " AstrometryTransformPolynomial::read : format is not 1 ");
1103 
1104  string order;
1105  s >> order >> _order;
1106  if (order != "order")
1108  " AstrometryTransformPolynomial::read : expecting \"order\" and found " + order);
1109  setOrder(_order);
1110  for (std::size_t k = 0; k < 2 * _nterms; ++k) s >> _coeffs[k];
1111 }
1112 
1113 ndarray::Array<double, 2, 2> AstrometryTransformPolynomial::toAstPolyMapCoefficients() const {
1114  std::size_t nCoeffs = _coeffs.size();
1115  ndarray::Array<double, 2, 2> result = ndarray::allocate(ndarray::makeVector(nCoeffs, std::size_t(4)));
1116 
1117  ndarray::Size k = 0;
1118  for (std::size_t iCoord = 0; iCoord < 2; ++iCoord) {
1119  for (std::size_t p = 0; p <= _order; ++p) {
1120  for (std::size_t py = 0; py <= p; ++py, ++k) {
1121  result[k][0] = coeff(p - py, py, iCoord);
1122  result[k][1] = iCoord + 1;
1123  result[k][2] = p - py;
1124  result[k][3] = py;
1125  }
1126  }
1127  }
1128 
1129  return result;
1130 }
1131 
1133  Frame const &domain,
1134  double const precision,
1135  std::size_t maxOrder,
1136  std::size_t nSteps) {
1137  StarMatchList sm;
1138  double xStart = domain.xMin;
1139  double yStart = domain.yMin;
1140  double xStep = domain.getWidth() / (nSteps - 1);
1141  double yStep = domain.getHeight() / (nSteps - 1);
1142  for (std::size_t i = 0; i < nSteps; ++i) {
1143  for (std::size_t j = 0; j < nSteps; ++j) {
1144  Point in(xStart + i * xStep, yStart + j * yStep);
1145  Point out(forward.apply(in));
1146  sm.push_back(StarMatch(out, in, nullptr, nullptr));
1147  }
1148  }
1149  std::size_t npairs = sm.size();
1150  std::size_t order;
1153  double chi2 = 0;
1154  double oldChi2 = std::numeric_limits<double>::infinity();
1155  for (order = 1; order <= maxOrder; ++order) {
1156  poly.reset(new AstrometryTransformPolynomial(order));
1157  auto success = poly->fit(sm);
1158  if (success == -1) {
1159  std::stringstream errMsg;
1160  errMsg << "Cannot fit a polynomial of order " << order << " with " << nSteps << "^2 points";
1161  throw pexExcept::RuntimeError(errMsg.str());
1162  }
1163  // compute the chi2 ignoring errors:
1164  chi2 = 0;
1165  for (auto const &i : sm) chi2 += i.point2.computeDist2(poly->apply((i.point1)));
1166  LOGLS_TRACE(_log, "inversePoly order " << order << ": " << chi2 << " / " << npairs << " = "
1167  << chi2 / npairs << " < " << precision * precision);
1168 
1169  if (chi2 / npairs < precision * precision) break;
1170 
1171  // If this triggers, we know we did not reach the required precision.
1172  if (chi2 > oldChi2) {
1173  LOGLS_WARN(_log, "inversePolyTransform: chi2 increases ("
1174  << chi2 << " > " << oldChi2 << "); ending fit with order: " << order);
1175  LOGLS_WARN(_log, "inversePolyTransform: requested precision not reached: "
1176  << chi2 << " / " << npairs << " = " << chi2 / npairs << " < "
1177  << precision * precision);
1178  poly = std::move(oldPoly);
1179  order--;
1180  break;
1181  } else {
1182  oldChi2 = chi2;
1183  // Clone it so we don't lose it in the next iteration.
1185  std::shared_ptr<AstrometryTransform>(poly->clone()));
1186  }
1187  }
1188  if (order > maxOrder)
1189  LOGLS_WARN(_log, "inversePolyTransform: Reached max order without reaching requested precision: "
1190  << chi2 << " / " << npairs << " = " << chi2 / npairs << " < "
1191  << precision * precision);
1192  return poly;
1193 }
1194 
1195 /**************** AstrometryTransformLinear ***************************************/
1196 /* AstrometryTransformLinear is a specialized constructor of AstrometryTransformPolynomial
1197  May be it could just disappear ??
1198 */
1199 
1200 AstrometryTransformLinear::AstrometryTransformLinear(const double Dx, const double Dy, const double A11,
1201  const double A12, const double A21, const double A22)
1203  dx() = Dx;
1204  a11() = A11;
1205  a12() = A12;
1206  dy() = Dy;
1207  a21() = A21;
1208  a22() = A22;
1209 }
1210 
1213  if (transform.getOrder() != 1)
1215  "Trying to build a AstrometryTransformLinear from a higher order transform. Aborting. ");
1217 }
1218 
1220  // There is a general routine in AstrometryTransformPolynomial that would do the job:
1221  // return AstrometryTransformLinear(AstrometryTransformPolynomial::operator*(right));
1222  // however, we are using this composition of linear stuff heavily in
1223  // AstrometryTransform::linearApproximation, itself used in inverseTransform::apply.
1224  // So, for sake of efficiency, and since it is easy, we take a shortcut:
1226  apply(right.Dx(), right.Dy(), result.dx(), result.dy());
1227  result.a11() = A11() * right.A11() + A12() * right.A21();
1228  result.a12() = A11() * right.A12() + A12() * right.A22();
1229  result.a21() = A21() * right.A11() + A22() * right.A21();
1230  result.a22() = A21() * right.A12() + A22() * right.A22();
1231  return result;
1232 }
1233 
1235  const double) const {
1236  derivative = *this;
1237  derivative.coeff(0, 0, 0) = 0;
1238  derivative.coeff(0, 0, 1) = 0;
1239 }
1240 
1242  return *this;
1243 }
1244 
1246  //
1247  // (T1,M1) * (T2,M2) = 1 i.e (0,1) implies
1248  // T1 = -M1 * T2
1249  // M1 = M2^(-1)
1250  //
1251 
1252  double a11 = A11();
1253  double a12 = A12();
1254  double a21 = A21();
1255  double a22 = A22();
1256  double d = (a11 * a22 - a12 * a21);
1257  if (d == 0) {
1258  LOGL_FATAL(_log,
1259  "AstrometryTransformLinear::invert singular transformation: transform contents will be "
1260  "dumped to stderr.");
1261  dump(cerr);
1262  }
1263 
1264  AstrometryTransformLinear result(0, 0, a22 / d, -a12 / d, -a21 / d, a11 / d);
1265  double rdx, rdy;
1266  result.apply(Dx(), Dy(), rdx, rdy);
1267  result.dx() = -rdx;
1268  result.dy() = -rdy;
1269  return result;
1270 }
1271 
1273  const Frame &) const {
1275 }
1276 
1278  throw pexExcept::NotFoundError("AstrometryTransformLinearRot::fit not implemented! aborting");
1279 }
1280 
1282  std::size_t npairs = starMatchList.size();
1283  if (npairs < 3) {
1284  LOGLS_FATAL(_log, "AstrometryTransformLinearShift::fit trying to fit a linear transform with only "
1285  << npairs << " matches.");
1286  return -1;
1287  }
1288 
1289  double sumr2 = 0; /* used to compute chi2 without relooping */
1290  /* loop on pairs and fill */
1291  Eigen::VectorXd B(2);
1292  B.setZero();
1293  Eigen::MatrixXd A(2, 2);
1294  A.setZero();
1295 
1296  for (auto const &it : starMatchList) {
1297  FatPoint const &point1 = it.point1;
1298  FatPoint const &point2 = it.point2;
1299  double deltax = point2.x - point1.x;
1300  double deltay = point2.y - point1.y;
1301  double vxx = point1.vx + point2.vx;
1302  double vyy = point1.vy + point2.vy;
1303  double vxy = point1.vxy + point2.vxy;
1304  double det = vxx * vyy - vxy * vxy;
1305  double wxx = vyy / det;
1306  double wyy = vxx / det;
1307  double wxy = -vxy / det;
1308  B(0) += deltax * wxx + wxy * deltay;
1309  B(1) += deltay * wyy + wxy * deltax;
1310  A(0, 0) += wxx;
1311  A(1, 1) += wyy;
1312  A(0, 1) += wxy;
1313  sumr2 += deltax * deltax * wxx + deltay * deltay * wyy + 2. * wxy * deltax * deltay;
1314  }
1315  double det = A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
1316  if (det <= 0) return -1;
1317  double tmp = A(0, 0);
1318  A(0, 0) = A(1, 1) / det;
1319  A(1, 1) = tmp / det;
1320  A(0, 1) = A(1, 0) = -A(0, 1) / det;
1321  Eigen::VectorXd sol = A * B;
1322  (*this) = AstrometryTransformLinearShift(sol(0), sol(1));
1323  return (sumr2 - sol.dot(B)); // chi2 compact form
1324 }
1325 
1327  const double scaleFactor) {
1328  double c = scaleFactor * std::cos(angleRad);
1329  double s = scaleFactor * std::sin(angleRad);
1330  a11() = a22() = c;
1331  a21() = s;
1332  a12() = -s;
1333 
1334  // we want that the center does not move : transform+M*C = C ==> transform = C - M*C
1335  Point a_point(0., 0.);
1336  if (center) a_point = *center;
1337 
1338  dx() = dy() = 0;
1339  AstrometryTransformPolynomial::apply(a_point.x, a_point.y, dx(), dy()); // compute M*C
1340  dx() = a_point.x - Dx();
1341  dy() = a_point.y - dy();
1342 }
1343 
1344 static double deg2rad(double degree) { return degree * M_PI / 180.; }
1345 
1346 static double rad2deg(double rad) { return rad * 180. / M_PI; }
1347 
1348 /************* WCS transform ******************/
1349 /************** LinPixelToTan *******************/
1350 
1351 /* Implementation note : it seemed wise to incorporate
1352  the radians to degreess convertion into the linPixelToTan
1353  part (right in the constructor), and to do the
1354  opposite operation in the LinPart routine.
1355  When I was coding the fit, I realized that it was a
1356  bad idea. then I realized that the fitting routine
1357  itself was probably useless. Finaly, for a persistor,
1358  it seems bizarre that the stored data is different
1359  from what convention (angles in degrees for WCS's)
1360  would expect.
1361  So, no more "automatic" degrees to radians and
1362  radians to degrees conversion. They are explicitely
1363  done in apply (for TanPixelToRaDec and TanRaDecToPixel).
1364  This is a minor concern though....
1365 */
1366 BaseTanWcs::BaseTanWcs(AstrometryTransformLinear const &pixToTan, Point const &tangentPoint,
1367  const AstrometryTransformPolynomial *corrections) {
1368  // the angles returned by linPixelToTan should be in degrees.
1369  linPixelToTan = pixToTan;
1370  ra0 = deg2rad(tangentPoint.x);
1371  dec0 = deg2rad(tangentPoint.y);
1372  cos0 = std::cos(dec0);
1373  sin0 = std::sin(dec0);
1374  corr = nullptr;
1375  if (corrections) corr.reset(new AstrometryTransformPolynomial(*corrections));
1376 }
1377 
1378 /* with some sort of smart pointer ro handle "corr", we could remove the
1379  copy constructor, the operator = and the destructor */
1380 
1381 // ": AstrometryTransform" suppresses a warning
1383  corr = nullptr;
1384  *this = original;
1385 }
1386 
1387 void BaseTanWcs::operator=(const BaseTanWcs &original) {
1388  linPixelToTan = original.linPixelToTan;
1389  ra0 = original.ra0;
1390  dec0 = original.dec0;
1391  cos0 = std::cos(dec0);
1392  sin0 = std::sin(dec0);
1393  corr = nullptr;
1394  if (original.corr) corr.reset(new AstrometryTransformPolynomial(*original.corr));
1395 }
1396 
1397 void BaseTanWcs::apply(const double xIn, const double yIn, double &xOut, double &yOut) const {
1398  double l, m; // radians in the tangent plane
1399  pixToTangentPlane(xIn, yIn, l, m); // l, m in degrees.
1400  l = deg2rad(l);
1401  m = deg2rad(m); // now in radians
1402  // Code inspired from worldpos.c in wcssubs (ancestor of the wcslib)
1403  /* At variance with wcslib, it collapses the projection to a plane
1404  and expression of sidereal cooordinates into a single set of
1405  operations. */
1406  double dect = cos0 - m * sin0;
1407  if (dect == 0) {
1408  LOGL_WARN(_log, "No sidereal coordinates at pole!");
1409  xOut = 0;
1410  yOut = 0;
1411  return;
1412  }
1413  double rat = ra0 + atan2(l, dect);
1414  dect = atan(std::cos(rat - ra0) * (m * cos0 + sin0) / dect);
1415  if (rat - ra0 > M_PI) rat -= (2. * M_PI);
1416  if (rat - ra0 < -M_PI) rat += (2. * M_PI);
1417  if (rat < 0.0) rat += (2. * M_PI);
1418  // convert to degree
1419  xOut = rad2deg(rat);
1420  yOut = rad2deg(dect);
1421 }
1422 
1423 Point BaseTanWcs::getTangentPoint() const { return Point(rad2deg(ra0), rad2deg(dec0)); }
1424 
1426 
1428  corr = std::move(corrections);
1429 }
1430 
1432  /* CRPIX's are defined by:
1433  ( CD1_1 CD1_2 ) (x - crpix1)
1434  transformed = ( ) * ( )
1435  ( CD2_1 CD2_2 ) (y - crpix2)
1436 
1437  so that CrPix is the point which transforms to (0,0)
1438  */
1440  return Point(inverse.Dx(), inverse.Dy());
1441 }
1442 
1444 
1446  : _skyWcs(skyWcs) {}
1447 
1448 void AstrometryTransformSkyWcs::apply(const double xIn, const double yIn, double &xOut, double &yOut) const {
1449  auto const outCoord = _skyWcs->pixelToSky(geom::Point2D(xIn, yIn));
1450  xOut = outCoord[0].asDegrees();
1451  yOut = outCoord[1].asDegrees();
1452 }
1453 
1455  stream << "AstrometryTransformSkyWcs(" << *_skyWcs << ")";
1456 }
1457 
1458 double AstrometryTransformSkyWcs::fit(const StarMatchList &starMatchList) {
1459  throw LSST_EXCEPT(pex::exceptions::LogicError, "Not implemented");
1460 }
1461 
1464 }
1465 
1466 /*************************** TanPixelToRaDec ***************/
1467 
1469  const AstrometryTransformPolynomial *corrections)
1470  : BaseTanWcs(pixToTan, tangentPoint, corrections) {}
1471 
1472 // ": AstrometryTransform" suppresses a warning
1474 
1476  AstrometryTransformLinear const &right) const {
1477  if (right.getOrder() == 1) {
1478  return std::make_unique<TanPixelToRaDec>((*this) * (right));
1479  } else {
1480  return std::unique_ptr<AstrometryTransform>(nullptr);
1481  }
1482 }
1483 
1485  TanPixelToRaDec result(*this);
1486  result.linPixelToTan = result.linPixelToTan * right;
1487  return result;
1488 }
1489 
1491  if (corr != nullptr) {
1492  LOGL_WARN(_log, "You are inverting a TanPixelToRaDec with corrections.");
1493  LOGL_WARN(_log, "The inverse you get ignores the corrections!");
1494  }
1496 }
1497 
1501 }
1502 
1504  const Frame &region) const {
1505  if (!corr)
1508  else
1509  return std::unique_ptr<AstrometryTransform>(new AstrometryTransformInverse(this, precision, region));
1510 }
1511 
1513  if (corr)
1514  return (*corr) * linPixelToTan;
1515  else
1516  return linPixelToTan;
1517 }
1518 
1519 void TanPixelToRaDec::pixToTangentPlane(double xPixel, double yPixel, double &xTangentPlane,
1520  double &yTangentPlane) const {
1521  // xTangentPlane, yTangentPlane in degrees.
1522  linPixelToTan.apply(xPixel, yPixel, xTangentPlane, yTangentPlane);
1523  if (corr) {
1524  double xtmp = xTangentPlane;
1525  double ytmp = yTangentPlane;
1526  corr->apply(xtmp, ytmp, xTangentPlane, yTangentPlane); // still in degrees.
1527  }
1528 }
1529 
1532  new TanPixelToRaDec(getLinPart(), getTangentPoint(), corr.get()));
1533 }
1534 
1536  stream << " TanPixelToRaDec, lin part :" << endl << linPixelToTan;
1537  Point tp = getTangentPoint();
1538  stream << " tangent point " << tp.x << ' ' << tp.y << endl;
1539  Point crpix = getCrPix();
1540  stream << " crpix " << crpix.x << ' ' << crpix.y << endl;
1541  if (corr) stream << "PV correction: " << endl << *corr;
1542 }
1543 
1545  /* OK we could implement this routine, but it is
1546  probably useless since to do the match, we have to
1547  project from sky to tangent plane. When a match is
1548  found, it is easier to carry out the fit in the
1549  tangent plane, rather than going back to the celestial
1550  sphere (and reproject to fit...). Anyway if this
1551  message shows up, we'll think about it.
1552  */
1554  "TanPixelToRaDec::fit is NOT implemented (although it is doable)) ");
1555  return -1;
1556 }
1557 
1558 /*************************** TanSipPixelToRaDec ***************/
1559 
1561  const AstrometryTransformPolynomial *corrections)
1562  : BaseTanWcs(pixToTan, tangentPoint, corrections) {}
1563 
1564 // ": AstrometryTransform" suppresses a warning
1566 
1567 /* Would require some checks before cooking up something more efficient
1568  than just a linear approximation */
1569 #if 0
1571 {
1572  if (&region) {}
1574 }
1575 #endif
1576 
1578  const Frame &region)
1579  const { /* We have not implemented (yet) the reverse corrections available in SIP */
1580  return std::unique_ptr<AstrometryTransform>(new AstrometryTransformInverse(this, precision, region));
1581 }
1582 
1584  if (corr)
1585  return AstrometryTransformPolynomial(linPixelToTan) * (*corr);
1586  else
1587  return linPixelToTan;
1588 }
1589 
1590 void TanSipPixelToRaDec::pixToTangentPlane(double xPixel, double yPixel, double &xTangentPlane,
1591  double &yTangentPlane) const {
1592  // xTangentPlane, yTangentPlane returned in degrees
1593  if (corr) {
1594  double xtmp, ytmp;
1595  corr->apply(xPixel, yPixel, xtmp, ytmp);
1596  linPixelToTan.apply(xtmp, ytmp, xTangentPlane, yTangentPlane);
1597  } else
1598  linPixelToTan.apply(xPixel, yPixel, xTangentPlane, yTangentPlane);
1599 }
1600 
1604 }
1605 
1607  stream << " TanSipPixelToRaDec, lin part :" << endl << linPixelToTan;
1608  Point tp = getTangentPoint();
1609  stream << " tangent point " << tp.x << ' ' << tp.y << endl;
1610  Point crpix = getCrPix();
1611  stream << " crpix " << crpix.x << ' ' << crpix.y << endl;
1612  if (corr) stream << "PV correction: " << endl << *corr;
1613 }
1614 
1616  /* OK we could implement this routine, but it is
1617  probably useless since to do the match, we have to
1618  project from sky to tangent plane. When a match is
1619  found, it is easier to carry out the fit in the
1620  tangent plane, rather than going back to the celestial
1621  sphere (and reproject to fit...). Anyway if this
1622  message shows up, we'll think about it.
1623  */
1625  "TanSipPixelToRaDec::fit is NOT implemented (although it is doable)) ");
1626  return -1;
1627 }
1628 
1629 /*************** reverse transform of TanPixelToRaDec: TanRaDecToPixel ********/
1630 
1632  : linTan2Pix(tan2Pix) {
1633  setTangentPoint(tangentPoint);
1634 }
1635 
1636 void TanRaDecToPixel::setTangentPoint(Point const &tangentPoint) {
1637  /* the radian to degrees conversion after projection
1638  is handled in apply */
1639  ra0 = deg2rad(tangentPoint.x);
1640  dec0 = deg2rad(tangentPoint.y);
1641  cos0 = std::cos(dec0);
1642  sin0 = std::sin(dec0);
1643 }
1644 
1646  ra0 = dec0 = 0;
1647  cos0 = 1;
1648  sin0 = 0;
1649 }
1650 
1651 Point TanRaDecToPixel::getTangentPoint() const { return Point(rad2deg(ra0), rad2deg(dec0)); }
1652 
1654 
1655 // Use analytic derivatives, computed at the same time as the transform itself
1657  /* this routine is very similar to apply, but also propagates errors.
1658  The deg2rad and rad2deg are ignored for errors because they act as
1659  2 global scalings that cancel each other.
1660  Derivatives were computed using maple:
1661 
1662  l1 := sin(a - a0)*cos(d);
1663  m1 := sin(d)*sin(d0)+cos(d)*cos(d0)*cos(a-a0);
1664  l2 := sin(d)*cos(d0)-cos(d)*sin(d0)*cos(a-a0);
1665  simplify(diff(l1/m1,a));
1666  simplify(diff(l1/m1,d));
1667  simplify(diff(l2/m1,a));
1668  simplify(diff(l2/m1,d));
1669 
1670  Checked against AstrometryTransform::transformPosAndErrors (dec 09)
1671  */
1672  double ra = deg2rad(in.x);
1673  double dec = deg2rad(in.y);
1674  if (ra - ra0 > M_PI) ra -= (2. * M_PI);
1675  if (ra - ra0 < -M_PI) ra += (2. * M_PI);
1676  // Code inspired from worldpos.c in wcssubs (ancestor of the wcslib)
1677  // The same code is copied in ::apply()
1678 
1679  double coss = std::cos(dec);
1680  double sins = std::sin(dec);
1681  double sinda = std::sin(ra - ra0);
1682  double cosda = std::cos(ra - ra0);
1683  double l = sinda * coss;
1684  double m = sins * sin0 + coss * cos0 * cosda;
1685  l = l / m;
1686  m = (sins * cos0 - coss * sin0 * cosda) / m;
1687 
1688  // derivatives
1689  double deno =
1690  sq(sin0) - sq(coss) + sq(coss * cos0) * (1 + sq(cosda)) + 2 * sins * sin0 * coss * cos0 * cosda;
1691  double a11 = coss * (cosda * sins * sin0 + coss * cos0) / deno;
1692  double a12 = -sinda * sin0 / deno;
1693  double a21 = coss * sinda * sins / deno;
1694  double a22 = cosda / deno;
1695 
1696  FatPoint tmp;
1697  tmp.vx = a11 * (a11 * in.vx + 2 * a12 * in.vxy) + a12 * a12 * in.vy;
1698  tmp.vy = a21 * a21 * in.vx + a22 * a22 * in.vy + 2. * a21 * a22 * in.vxy;
1699  tmp.vxy = a21 * a11 * in.vx + a22 * a12 * in.vy + (a21 * a12 + a11 * a22) * in.vxy;
1700 
1701  // l and m are now coordinates in the tangent plane, in radians.
1702  tmp.x = rad2deg(l);
1703  tmp.y = rad2deg(m);
1704 
1705  linTan2Pix.transformPosAndErrors(tmp, out);
1706 }
1707 
1708 void TanRaDecToPixel::apply(const double xIn, const double yIn, double &xOut, double &yOut) const {
1709  double ra = deg2rad(xIn);
1710  double dec = deg2rad(yIn);
1711  if (ra - ra0 > M_PI) ra -= (2. * M_PI);
1712  if (ra - ra0 < -M_PI) ra += (2. * M_PI);
1713  // Code inspired from worldpos.c in wcssubs (ancestor of the wcslib)
1714  // The same code is copied in ::transformPosAndErrors()
1715  double coss = std::cos(dec);
1716  double sins = std::sin(dec);
1717  double l = std::sin(ra - ra0) * coss;
1718  double m = sins * sin0 + coss * cos0 * std::cos(ra - ra0);
1719  l = l / m;
1720  m = (sins * cos0 - coss * sin0 * std::cos(ra - ra0)) / m;
1721  // l and m are now coordinates in the tangent plane, in radians.
1722  l = rad2deg(l);
1723  m = rad2deg(m);
1724  linTan2Pix.apply(l, m, xOut, yOut);
1725 }
1726 
1729 }
1730 
1732  Point tp = getTangentPoint();
1733  stream << " tan2pix " << linTan2Pix << " tangent point " << tp.x << ' ' << tp.y << endl;
1734 }
1735 
1739 }
1740 
1744 }
1745 
1748 }
1749 
1752  "TanRaDecToPixel::fit is NOT implemented (although it is doable)) ");
1753  return -1;
1754 }
1755 
1756 /************* a "run-time" transform, that does not require to modify this file */
1757 
1759  : _userFun(userFun), _userData(userData) {}
1760 
1761 void UserTransform::apply(const double xIn, const double yIn, double &xOut, double &yOut) const {
1762  _userFun(xIn, yIn, xOut, yOut, _userData);
1763 }
1764 
1766  stream << "UserTransform with user function @ " << _userFun << "and userData@ " << _userData << endl;
1767 }
1768 
1771  "UserTransform::fit is NOT implemented (and will never be)) ");
1772  return -1;
1773 }
1774 
1777 }
1778 
1779 /*************************************************************/
1780 
1782  ifstream s(fileName.c_str());
1783  if (!s)
1785  " astrometryTransformRead : cannot open " + fileName);
1786  try {
1788  s.close();
1789  return res;
1792  std::string(e.what()) + " in file " + fileName);
1793  }
1794 }
1795 
1797  std::string type;
1798  s >> type;
1799  if (s.fail())
1801  "astrometryTransformRead : could not find a AstrometryTransformtype");
1802  if (type == "AstrometryTransformIdentity") {
1804  res->read(s);
1805  return std::move(res);
1806  } else if (type == "AstrometryTransformPolynomial") {
1808  res->read(s);
1809  return std::move(res);
1810  } else
1812  " astrometryTransformRead : No reader for AstrometryTransform type " + type);
1813 }
1814 } // namespace jointcal
1815 } // namespace lsst
#define LOGLS_WARN(logger, message)
Log a warn-level message using an iostream-based interface.
Definition: Log.h:648
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 ...
def write(self, patchRef, catalog)
Write the output.
AstrometryTransformPolynomial operator+(AstrometryTransformPolynomial const &right) const
Addition.
std::unique_ptr< AstrometryTransform > inverseTransform(const double precision, const Frame &region) const
returns an inverse transform. Numerical if not overloaded.
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 ...
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
void() AstrometryTransformFun(const double, const double, double &, double &, const void *)
signature of the user-provided routine that actually does the coordinate transform for UserTransform...
The transformation that handles pixels to sideral transformations (Gnomonic, possibly with polynomial...
T atan(T... args)
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 apply(const double xIn, const double yIn, double &xOut, double &yOut) const
implements an iterative (Gauss-Newton) solver.
double dec
Definition: Match.cc:41
void paramDerivatives(Point const &where, double *dx, double *dy) const override
Derivative w.r.t parameters.
double coeff(std::size_t powX, std::size_t powY, std::size_t whichCoord) const
access to coefficients (read only)
#define LOG_LOGGER
Definition: Log.h:703
void computeDerivative(Point const &where, AstrometryTransformLinear &derivative, const double step=0.01) const override
specialised analytic routine
A hanger for star associations.
Definition: StarMatch.h:54
implements the linear transformations (6 real coefficients).
A point in a plane.
Definition: Point.h:36
void apply(const double xIn, const double yIn, double &xOut, double &yOut) const
return second(first(xIn,yIn))
double getArea() const
Definition: Frame.cc:118
AstrometryTransformLinear linPixelToTan
Low-level polynomials (including special polynomials) in C++.
Definition: Basis1d.h:26
void setOrder(std::size_t order)
Sets the polynomial order (the highest sum of exponents of the largest monomial). ...
std::unique_ptr< AstrometryTransform > composeAndReduce(AstrometryTransformLinear const &right) const
Return a reduced composition of newTransform = this(right()), or nullptr if it cannot be reduced...
std::unique_ptr< AstrometryTransform > clone() const
returns a copy (allocated by new) of the transformation.
bool isIntegerShift(const AstrometryTransform *transform)
Shorthand test to tell if a transform is a simple integer shift.
double fit(StarMatchList const &starMatchList)
fits a transform to a std::list of Point pairs (p1,p2, the Point fields in StarMatch).
table::PointKey< double > crpix
Definition: OldWcs.cc:131
T fail(T... args)
void apply(const double xIn, const double yIn, double &xOut, double &yOut) const override
AstrometryTransformPolynomial getPixelToTangentPlane() const
the transformation from pixels to tangent plane (degrees)
FatPoint point2
2 points
Definition: StarMatch.h:60
AstrometryTransformPolynomial(std::size_t order=1)
Default transform : identity for all orders (>=1 ).
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 =0
Transform from pixels to tangent plane (degrees)
T endl(T... args)
T left(T... args)
double fit(StarMatchList const &starMatchList)
fits a transform to a std::list of Point pairs (p1,p2, the Point fields in StarMatch).
int y
Definition: SpanSet.cc:49
table::Key< int > a
std::unique_ptr< AstrometryTransform > clone() const
returns a copy (allocated by new) of the transformation.
void apply(const double xIn, const double yIn, double &xOut, double &yOut) const
Transform pixels to ICRS RA, Dec in degrees.
STL namespace.
T precision(T... args)
double xMin
coordinate of boundary.
Definition: Frame.h:41
TanPixelToRaDec operator*(AstrometryTransformLinear const &right) const
composition with AstrometryTransformLinear
T end(T... args)
#define LOGL_ERROR(logger, message...)
Log a error-level message using a varargs/printf style interface.
Definition: Log.h:552
void apply(const double xIn, const double yIn, double &xOut, double &yOut) const
T floor(T... args)
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.
void dump(std::ostream &stream=std::cout) const override
print out of coefficients in a readable form.
AstrometryTransformInverse(const AstrometryTransform *direct, const double precision, const Frame &region)
long double & coeff(std::size_t powX, std::size_t powY)
STL class.
double paramRef(Eigen::Index const i) const override
int min
A Point with uncertainties.
Definition: FatPoint.h:34
#define M_PI
Definition: ListMatch.cc:31
void write(std::ostream &s) const override
void setCorrections(std::unique_ptr< AstrometryTransformPolynomial > corrections)
Assign the correction polynomial (what it means is left to derived classes)
std::ostream & operator<<(std::ostream &os, CameraSysPrefix const &detSysPrefix)
Definition: CameraSys.cc:47
std::unique_ptr< AstrometryTransform > roughInverse(const Frame &region) const
Overload the "generic routine" (available for all AstrometryTransform types.
std::unique_ptr< AstrometryTransformPolynomial > corr
T resize(T... args)
T atan2(T... args)
STL class.
UserTransform(AstrometryTransformFun &fun, const void *userData)
the transform routine and extra data that it may need.
LSST DM logging module built on log4cxx.
double x
coordinate
Definition: Point.h:41
T min(T... args)
pairs of points
T sin(T... args)
void computeDerivative(Point const &where, AstrometryTransformLinear &derivative, const double step=0.01) const
specialised analytic routine
T at(T... args)
AstrometryTransformLinear getLinPart() const
The Linear part (corresponding to CD&#39;s and CRPIX&#39;s)
double getHeight() const
size along y axis
Definition: Frame.h:57
std::size_t getOrder() const
Returns the polynomial order.
TanPixelToRaDec inverted() const
exact typed inverse:
Reports attempts to access elements using an invalid key.
Definition: Runtime.h:151
T push_back(T... args)
virtual char const * what(void) const noexcept
Return a character string summarizing this exception.
Definition: Exception.cc:99
void setTangentPoint(Point const &tangentPoint)
Resets the projection (or tangent) point.
STL class.
std::unique_ptr< AstrometryTransform > roughInverse(const Frame &region) const
Overload the "generic routine" (available for all AstrometryTransform types.
AstrometryTransformLinear getLinPart() const
The Linear part (corresponding to CD&#39;s and CRPIX&#39;s)
AstrometryTransformPolynomial operator*(AstrometryTransformPolynomial const &right) const
Composition (internal stuff in quadruple precision)
std::unique_ptr< AstrometryTransform > clone() const
returns a copy (allocated by new) of the transformation.
#define LOGL_FATAL(logger, message...)
Log a fatal-level message using a varargs/printf style interface.
Definition: Log.h:568
virtual void pixToTangentPlane(double xPixel, double yPixel, double &xTangentPlane, double &yTangentPlane) const
transforms from pixel space to tangent plane (degrees)
rectangle with sides parallel to axes.
Definition: Frame.h:38
AstrometryTransformLinear normalizeCoordinatesTransform(const Frame &frame)
Returns the transformation that maps the input frame along both axes to [-1,1].
void dump(std::ostream &stream) const
dumps the transform coefficients to stream.
A base class for image defects.
AstrometryTransformComposition(AstrometryTransform const &second, AstrometryTransform const &first)
#define LOGLS_TRACE(logger, message)
Log a trace-level message using an iostream-based interface.
Definition: Log.h:588
AstrometryTransformPolynomial getPixelToTangentPlane() const
the transformation from pixels to tangent plane (degrees)
table::Key< int > type
Definition: Detector.cc:163
std::unique_ptr< AstrometryTransform > clone() const
returns a copy (allocated by new) of the transformation.
Point getCrPix() const
Get the pixel origin of the WCS (CRPIX in FITS WCS terminology, but zero-based)
void dump(std::ostream &stream) const
dumps the transform coefficients to stream.
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...
std::unique_ptr< AstrometryTransform > clone() const override
returns a copy (allocated by new) of the transformation.
int const step
AstrometryTransformLinear()
the default constructor constructs the do-nothing transformation.
void apply(const double xIn, const double yIn, double &xOut, double &yOut) const override
std::unique_ptr< AstrometryTransform > compose(AstrometryTransform const &left, AstrometryTransform const &right)
Returns a pointer to a composition of transforms, representing left(right()).
PolyXY(AstrometryTransformPolynomial const &transform, std::size_t whichCoord)
def getParams(record, galType='sersic')
Definition: compareModel.py:17
T str(T... args)
table::Key< table::Array< int > > degree
Definition: PsfexPsf.cc:360
This one is the Tangent Plane (called gnomonic) projection (from celestial sphere to tangent plane) ...
int max
just here to provide a specialized constructor, and fit.
T reset(T... args)
T cos(T... args)
T dynamic_pointer_cast(T... args)
void apply(const double xIn, const double yIn, double &xOut, double &yOut) const
#define LOGL_WARN(logger, message...)
Log a warn-level message using a varargs/printf style interface.
Definition: Log.h:536
Reports errors in the logical structure of the program.
Definition: Runtime.h:46
T infinity(T... args)
T fabs(T... args)
T max(T... args)
T move(T... args)
std::unique_ptr< AstrometryTransform > clone() const
returns a copy (allocated by new) of the transformation.
double fit(StarMatchList const &starMatchList)
Not implemented yet, because we do it otherwise.
void dump(ostream &stream) const
dumps the transform coefficients to stream.
double x
double getWidth() const
size along x axis
Definition: Frame.h:54
T count(T... args)
std::unique_ptr< AstrometryTransform > inverseTransform(const double precision, const Frame &region) const
Inverse transform: returns a TanPixelToRaDec.
void dump(std::ostream &stream=std::cout) const
dumps the transform coefficients to stream.
BaseTanWcs(AstrometryTransformLinear const &pixToTan, Point const &tangentPoint, const AstrometryTransformPolynomial *corrections=nullptr)
void transformPosAndErrors(const FatPoint &in, FatPoint &out) const
transform with analytical derivatives
A do-nothing transformation. It anyway has dummy routines to mimick a AstrometryTransform.
std::unique_ptr< AstrometryTransform > astrometryTransformRead(const std::string &fileName)
The virtual constructor from a file.
T insert(T... args)
Point getCenter() const
Center of the frame.
Definition: Frame.h:60
T size(T... args)
std::size_t getOrder() const
double fit(StarMatchList const &starMatchList)
guess what
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
STL class.
AstrometryTransformLinear inverted() const
returns the inverse: T1 = T2.inverted();
void dump(std::ostream &stream) const
dumps the transform coefficients to stream.
double fit(StarMatchList const &starMatchList)
guess what
a virtual (interface) class for geometric transformations.
T begin(T... args)
std::unique_ptr< AstrometryTransform > inverseTransform(double, const Frame &) const
Inverse transform: returns the direct one!
T pow(T... args)
double fit(StarMatchList const &starMatchList)
fits a transform to a std::list of Point pairs (p1,p2, the Point fields in StarMatch).
T c_str(T... args)
virtual void pixToTangentPlane(double xPixel, double yPixel, double &xTangentPlane, double &yTangentPlane) const
transforms from pixel space to tangent plane (degrees)
AstrometryTransformLinear operator*(AstrometryTransformLinear const &right) const
enables to combine linear tranformations: T1=T2*T3 is legal.
virtual std::unique_ptr< AstrometryTransform > clone() const
returns a copy (allocated by new) of the transformation.
virtual void transformPosAndErrors(const FatPoint &in, FatPoint &out) const override
a mix of apply and Derivative
std::shared_ptr< afw::geom::SkyWcs > getSkyWcs() const
Point getTangentPoint() const
Get the sky origin (CRVAL in FITS WCS terminology) in degrees.
AstrometryTransformPolynomial operator-(AstrometryTransformPolynomial const &right) const
Subtraction.
void write(std::ostream &s) const override
Reports invalid arguments.
Definition: Runtime.h:66
AstrometryTransformSkyWcs(std::shared_ptr< afw::geom::SkyWcs > skyWcs)
std::unique_ptr< AstrometryTransform > composeAndReduce(AstrometryTransformPolynomial const &right) const
Return a reduced composition of newTransform = this(right()), or nullptr if it cannot be reduced...
int m
Definition: SpanSet.cc:49
std::unique_ptr< AstrometryTransform > roughInverse(const Frame &) const
Overload the "generic routine".
T transform(T... args)
T sqrt(T... args)
table::Key< int > transform
double fit(StarMatchList const &starMatchList) override
guess what
Private class to handle AstrometryTransform compositions (i.e.
#define LOGLS_FATAL(logger, message)
Log a fatal-level message using an iostream-based interface.
Definition: Log.h:688
TanRaDecToPixel inverted() const
approximate inverse : it ignores corrections;
double fit(const StarMatchList &starMatchList) override
Not implemented; throws pex::exceptions::LogicError.
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).
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...
#define LOG_GET(logger)
Returns a Log object associated with logger.
Definition: Log.h:75
T setprecision(T... args)
AstrometryTransformLinear linearApproximation(Point const &where, const double step=0.01) const
linear (local) approximation.
STL class.
virtual AstrometryTransformLinear linearApproximation(Point const &where, const double step=0.01) const override
linear approximation.
py::object result
Definition: _schema.cc:429
void dump(std::ostream &stream=std::cout) const override
dumps the transform coefficients to stream.
virtual std::unique_ptr< AstrometryTransform > clone() const =0
returns a copy (allocated by new) of the transformation.
virtual void dump(std::ostream &stream=std::cout) const =0
dumps the transform coefficients to stream.
just here to provide specialized constructors. AstrometryTransformLinear fit routine.
STL class.
void operator=(const BaseTanWcs &original)
double coeffOrZero(std::size_t powX, std::size_t powY, std::size_t whichCoord) const
read access, zero if beyond order
virtual void apply(const double xIn, const double yIn, double &xOut, double &yOut) const =0
T reserve(T... args)
Reports errors that are due to events beyond the control of the program.
Definition: Runtime.h:104
void dump(ostream &stream=cout) const
dumps the transform coefficients to stream.
Point getTangentPoint() const
tangent point coordinates (degrees)
T emplace_back(T... args)
long double coeff(std::size_t powX, std::size_t powY) const