43 #include "Eigen/Cholesky"
59 if (shift ==
nullptr)
return false;
61 static const double eps = 1e-5;
66 static Point dumb(4000, 4000);
68 fabs(dumb.
x + dx - shift->
apply(dumb).x) < eps &&
fabs(dumb.
y + dy - shift->
apply(dumb).y) < eps)
76 Frame AstrometryTransform::apply(
Frame const &inputframe,
bool inscribed)
const {
78 double xtmin1, xtmax1, ytmin1, ytmax1;
79 apply(inputframe.
xMin, inputframe.
yMin, xtmin1, ytmin1);
80 apply(inputframe.
xMax, inputframe.
yMax, xtmax1, ytmax1);
84 double xtmin2, xtmax2, ytmin2, ytmax2;
85 apply(inputframe.
xMin, inputframe.
yMax, xtmin2, ytmax2);
86 apply(inputframe.
xMax, inputframe.
yMin, xtmax2, ytmin2);
90 if (inscribed)
return fr1 * fr2;
101 double eps =
x * 0.01;
102 if (eps == 0) eps = 0.01;
105 apply(
x + eps,
y, dxdx, dydx);
109 apply(
x,
y + eps, dxdy, dydy);
112 return ((dxdx * dydy - dxdy * dydx) / (eps * eps));
119 const double step)
const {
123 apply(
x,
y, xp0, yp0);
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;
137 const double step)
const {
138 Point outwhere = apply(where);
140 computeDerivative(where, der,
step);
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;
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();
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];
190 vOut[xx] = b11 * a11 + b12 * a12;
191 vOut[xy] = b11 * a21 + b12 * a22;
192 vOut[yy] = b21 * a21 + b22 * a22;
198 Point centerIn = apply(centerOut);
216 void AstrometryTransform::getParams(
double *params)
const {
218 for (
std::size_t i = 0; i < npar; ++i) params[i] = paramRef(i);
221 void AstrometryTransform::offsetParams(Eigen::VectorXd
const &delta) {
223 for (
std::size_t i = 0; i < npar; ++i) paramRef(i) += delta[i];
226 double AstrometryTransform::paramRef(Eigen::Index
const)
const {
228 std::string(
"AstrometryTransform::paramRef should never be called "));
231 double &AstrometryTransform::paramRef(Eigen::Index
const) {
233 "AstrometryTransform::paramRef should never be called ");
236 void AstrometryTransform::paramDerivatives(
Point const &,
double *,
double *)
const {
238 "AstrometryTransform::paramDerivatives() should never be called ");
253 "AstrometryTransform::write, something went wrong for file " + fileName);
259 "AstrometryTransform::write(ostream), should never be called. MEans that it is missing in some "
275 const Frame ®ion);
279 void apply(
const double xIn,
const double yIn,
double &xOut,
double &yOut)
const;
281 void print(
ostream &out)
const;
294 return _direct->clone();
304 const Frame ®ion)
const {
310 _direct = direct->
clone();
311 _roughInverse = _direct->roughInverse(region);
317 _direct = model._direct->clone();
318 _roughInverse = model._roughInverse->clone();
319 precision2 = model.precision2;
325 _direct = model._direct->clone();
326 _roughInverse = model._roughInverse->clone();
327 precision2 = model.precision2;
332 Point outGuess = _roughInverse->apply(in);
339 Point inGuess = _direct->apply(outGuess);
340 _direct->computeDerivative(outGuess, directDer);
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);
354 stream <<
" AstrometryTransformInverse of :" <<
endl << *_direct <<
endl;
359 "Cannot fit a AstrometryTransformInverse. Use StarMatchList::inverseTransform instead.");
381 void apply(
const double xIn,
const double yIn,
double &xOut,
double &yOut)
const;
393 _first =
first.clone();
398 double &yOut)
const {
400 _first->apply(xIn, yIn, xout, yout);
401 _second->apply(xout, yout, xOut, yOut);
405 stream <<
"Composed AstrometryTransform consisting of:" <<
std::endl;
406 _first->print(stream);
407 _second->print(stream);
412 return _first->fit(starMatchList);
416 return std::make_unique<AstrometryTransformComposition>(*_second, *_first);
432 if (composition ==
nullptr)
433 return std::make_unique<AstrometryTransformComposition>(
left,
right);
440 const double)
const {
445 const double)
const {
451 return std::make_shared<ast::UnitMap>(2);
455 stream <<
"AstrometryTransformIdentity 1" <<
endl;
463 " AstrometryTransformIdentity::read : format is not 1 ");
474 _coeffs.
resize(2 * _nterms, 0.);
491 auto pix = std::make_shared<BaseStar>(
x,
y, 0, 0);
494 auto tp = std::make_shared<BaseStar>(xtr, ytr, 0, 0);
509 double xStart = domain.
xMin;
510 double yStart = domain.
yMin;
511 double xStep = domain.
getWidth() / (nSteps + 1);
512 double yStep = domain.
getHeight() / (nSteps + 1);
524 poly.fit(starMatchList);
528 void AstrometryTransformPolynomial::computeMonomials(
double xIn,
double yIn,
double *monomial)
const {
542 for (
std::size_t iy = 0; iy <= _order - ix; ++iy) {
543 monomial[k] = xx * yy;
554 _nterms = (_order + 1) * (_order + 2) / 2;
559 _coeffs.
resize(2 * _nterms);
562 for (
std::size_t k = 0; k < _nterms; ++k) _coeffs[k] = 0;
566 _coeffs[k] = old_coeffs[k];
567 _coeffs[k + _nterms] = old_coeffs[k + old_nterms];
573 double &yOut)
const {
583 double monomials[_nterms];
584 computeMonomials(xIn, yIn, monomials);
588 const double *c = &_coeffs[0];
589 const double *pm = &monomials[0];
591 for (
int k = _nterms; k--;) xOut += (*(pm++)) * (*(c++));
593 for (
int k = _nterms; k--;) yOut += (*(pm++)) * (*(c++));
602 derivative.
dx() = derivative.
dy() = 0;
606 double dermx[2 * _nterms];
607 double *dermy = dermx + _nterms;
608 double xin = where.
x;
609 double yin = where.
y;
616 dermx[k] = ix * xxm1;
620 for (
std::size_t iy = 1; iy <= _order - ix; ++iy) {
621 dermx[k] = ix * xxm1 * yym1 * yin;
622 dermy[k] = iy * xx * yym1;
627 if (ix >= 1) xxm1 *= xin;
633 const double *mx = &dermx[0];
634 const double *my = &dermy[0];
635 const double *c = &_coeffs[0];
637 double a11 = 0, a12 = 0;
638 for (
int k = _nterms; k--;) {
639 a11 += (*(mx++)) * (*c);
640 a12 += (*(my++)) * (*(c++));
642 derivative.
a11() = a11;
643 derivative.
a12() = a12;
645 double a21 = 0, a22 = 0;
648 for (
int k = _nterms; k--;) {
649 a21 += (*(mx++)) * (*c);
650 a22 += (*(my++)) * (*(c++));
652 derivative.
a21() = a21;
653 derivative.
a22() = a22;
669 double monomials[_nterms];
673 double dermx[2 * _nterms];
674 double *dermy = dermx + _nterms;
683 dermx[k] = ix * xxm1;
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;
698 if (ix >= 1) xxm1 *= xin;
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++));
707 for (
int k = _nterms; k--;) yout += (*(pm++)) * (*(c++));
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++));
721 double a21 = 0, a22 = 0;
724 for (
int k = _nterms; k--;) {
725 a21 += (*(mx++)) * (*c);
726 a22 += (*(my++)) * (*(c++));
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;
742 assert((degX + degY <= _order) && whichCoord < 2);
746 return _coeffs[(degX + degY) * (degX + degY + 1) / 2 + degY + whichCoord * _nterms];
751 assert((degX + degY <= _order) && whichCoord < 2);
752 return _coeffs[(degX + degY) * (degX + degY + 1) / 2 + degY + whichCoord * _nterms];
758 assert(whichCoord < 2);
759 if (degX + degY <= _order)
760 return _coeffs[(degX + degY) * (degX + degY + 1) / 2 + degY + whichCoord * _nterms];
766 assert(i < 2 * Eigen::Index(_nterms));
771 assert(i < 2 * Eigen::Index(_nterms));
777 computeMonomials(where.
x, where.
y, dx);
779 dy[_nterms + k] = dx[k];
780 dx[_nterms + k] = dy[k] = 0;
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;
804 bool printed =
false;
808 if (coefficient != 0 ||
isnan(coefficient)) {
810 if (printed && p +
py != 0) stream <<
" + ";
812 stream << coefficient << monomialString(p -
py,
py);
821 if (_order > 0) stream <<
"Linear determinant = " <<
determinant();
840 static AstrometryTransformLinear shiftAndNormalize(StarMatchList
const &starMatchList) {
846 for (
auto it = starMatchList.begin(); it != starMatchList.end(); ++it) {
847 const StarMatch &a_match = *it;
848 Point
const &point1 = a_match.point1;
855 if (count == 0)
return AstrometryTransformLinear();
861 return AstrometryTransformLinearScale(2. / xspan, 2. / yspan) *
862 AstrometryTransformLinearShift(-xav, -yav);
865 static double sq(
double x) {
return x *
x; }
867 double AstrometryTransformPolynomial::computeFit(StarMatchList
const &starMatchList,
868 AstrometryTransform
const &shiftToCenter,
869 const bool useErrors) {
870 Eigen::MatrixXd A(2 * _nterms, 2 * _nterms);
872 Eigen::VectorXd B(2 * _nterms);
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;
883 computeMonomials(point1.x, point1.y, 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;
896 apply(point1.x, point1.y, tr1.x, tr1.y);
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;
902 double bxcoeff = wxx * resx + wxy * resy;
903 double bycoeff = wyy * resy + wxy * resx;
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];
910 B(j) += bxcoeff * monomials[j];
911 B(j + _nterms) += bycoeff * monomials[j];
914 Eigen::LDLT<Eigen::MatrixXd, Eigen::Lower> factor(A);
916 if (factor.info() != Eigen::Success) {
917 LOGL_ERROR(_log,
"AstrometryTransformPolynomial::fit could not factorize");
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));
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.");
936 computeFit(starMatchList, conditionner,
false);
937 computeFit(starMatchList, conditionner,
true);
938 double chi2 = computeFit(starMatchList, conditionner,
true);
940 (*this) = (*this) * conditionner;
941 if (starMatchList.
size() == _nterms)
return 0;
948 return std::make_unique<AstrometryTransformLinear>((*
this) * (
right));
950 return std::make_unique<AstrometryTransformPolynomial>((*
this) * (
right));
984 assert(powX + powY <=
order);
985 return coeffs.
at((powX + powY) * (powX + powY + 1) / 2 + powY);
989 assert(powX + powY <=
order);
990 return coeffs.
at((powX + powY) * (powX + powY + 1) / 2 + powY);
996 static void operator+=(PolyXY &left,
const PolyXY &right) {
998 assert(
left.getOrder() >= rdeg);
1000 for (
std::size_t j = 0; j <= rdeg - i; ++j)
left.getCoefficient(i, j) +=
right.getCoefficient(i, j);
1004 static PolyXY operator*(
const long double &
a,
const PolyXY &polyXY) {
1014 static PolyXY product(
const PolyXY &p1,
const PolyXY &p2) {
1017 PolyXY
result(deg1 + deg2);
1022 result.getCoefficient(i1 + i2, j1 + j2) +=
1023 p1.getCoefficient(i1, j1) * p2.getCoefficient(i2, j2);
1031 powers[0].getCoefficient(0, 0) = 1L;
1036 static PolyXY composition(
const PolyXY &polyXY,
const PolyXY &polyX,
const PolyXY &polyY) {
1038 PolyXY
result(pdeg *
max(polyX.getOrder(), polyY.getOrder()));
1041 computePowers(polyX, pdeg, pXPowers);
1042 computePowers(polyY, pdeg, pYPowers);
1045 result += polyXY.getCoefficient(px,
py) * product(pXPowers.
at(px), pYPowers.
at(
py));
1062 PolyXY rx(composition(plx, prx, pry));
1063 PolyXY ry(composition(ply, prx, pry));
1077 if (_order >=
right._order) {
1086 return (
right + (*
this));
1093 for (
std::size_t j = 0; j <= res._order - i; ++j) {
1102 return std::make_shared<ast::PolyMap>(toAstPolyMapCoefficients(), inverse->toAstPolyMapCoefficients());
1106 s <<
" AstrometryTransformPolynomial 1" <<
endl;
1107 s <<
"order " << _order <<
endl;
1110 for (
std::size_t k = 0; k < 2 * _nterms; ++k) s << _coeffs[k] <<
' ';
1120 " AstrometryTransformPolynomial::read : format is not 1 ");
1123 s >>
order >> _order;
1124 if (
order !=
"order")
1126 " AstrometryTransformPolynomial::read : expecting \"order\" and found " +
order);
1128 for (
std::size_t k = 0; k < 2 * _nterms; ++k) s >> _coeffs[k];
1131 ndarray::Array<double, 2, 2> AstrometryTransformPolynomial::toAstPolyMapCoefficients()
const {
1133 ndarray::Array<double, 2, 2>
result = ndarray::allocate(ndarray::makeVector(nCoeffs,
std::size_t(4)));
1135 ndarray::Size k = 0;
1136 for (
std::size_t iCoord = 0; iCoord < 2; ++iCoord) {
1140 result[k][1] = iCoord + 1;
1151 Frame const &domain,
1156 double xStart = domain.
xMin;
1157 double yStart = domain.
yMin;
1158 double xStep = domain.
getWidth() / (nSteps - 1);
1159 double yStep = domain.
getHeight() / (nSteps - 1);
1162 Point in(xStart + i * xStep, yStart + j * yStep);
1175 auto success =
poly->fit(sm);
1176 if (success == -1) {
1178 errMsg <<
"Cannot fit a polynomial of order " <<
order <<
" with " << nSteps <<
"^2 points";
1183 for (
auto const &i : sm) chi2 += i.point2.computeDist2(
poly->apply((i.point1)));
1184 LOGLS_TRACE(_log,
"inversePoly order " <<
order <<
": " << chi2 <<
" / " << npairs <<
" = "
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 <<
" < "
1202 oldPoly = dynamic_pointer_cast<AstrometryTransformPolynomial>(
1206 if (
order > maxOrder)
1207 LOGLS_WARN(_log,
"inversePolyTransform: Reached max order without reaching requested precision: "
1208 << chi2 <<
" / " << npairs <<
" = " << chi2 / npairs <<
" < "
1219 const double A12,
const double A21,
const double A22)
1233 "Trying to build a AstrometryTransformLinear from a higher order transform. Aborting. ");
1253 const double)
const {
1277 "AstrometryTransformLinear::invert singular transformation: transform contents will be "
1278 "dumped to stderr.");
1291 const Frame &)
const {
1298 stream <<
"Linear AstrometryTransform:" <<
std::endl;
1313 LOGLS_FATAL(_log,
"AstrometryTransformLinearShift::fit trying to fit a linear transform with only "
1314 << npairs <<
" matches.");
1320 Eigen::VectorXd B(2);
1322 Eigen::MatrixXd A(2, 2);
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;
1342 sumr2 += deltax * deltax * wxx + deltay * deltay * wyy + 2. * wxy * deltax * deltay;
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;
1352 return (sumr2 - sol.dot(B));
1356 const double scaleFactor) {
1357 double c = scaleFactor *
std::cos(angleRad);
1358 double s = scaleFactor *
std::sin(angleRad);
1364 Point a_point(0., 0.);
1365 if (center) a_point = *center;
1369 dx() = a_point.
x -
Dx();
1370 dy() = a_point.
y -
dy();
1375 static double rad2deg(
double rad) {
return rad * 180. /
M_PI; }
1399 ra0 = deg2rad(tangentPoint.
x);
1400 dec0 = deg2rad(tangentPoint.
y);
1437 LOGL_WARN(_log,
"No sidereal coordinates at pole!");
1446 if (rat < 0.0) rat += (2. *
M_PI);
1448 xOut = rad2deg(rat);
1449 yOut = rad2deg(dect);
1469 return Point(inverse.
Dx(), inverse.
Dy());
1475 : _skyWcs(skyWcs) {}
1478 auto const outCoord = _skyWcs->pixelToSky(
geom::Point2D(xIn, yIn));
1479 xOut = outCoord[0].asDegrees();
1480 yOut = outCoord[1].asDegrees();
1484 out <<
"AstrometryTransformSkyWcs(" << *_skyWcs <<
")";
1499 :
BaseTanWcs(pixToTan, tangentPoint, corrections) {}
1506 if (
right.getOrder() == 1) {
1507 return std::make_unique<TanPixelToRaDec>((*
this) * (
right));
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!");
1533 const Frame ®ion)
const {
1549 double &yTangentPlane)
const {
1553 double xtmp = xTangentPlane;
1554 double ytmp = yTangentPlane;
1555 corr->apply(xtmp, ytmp, xTangentPlane, yTangentPlane);
1567 stream <<
" tangent point " << tp.
x <<
' ' << tp.
y <<
endl;
1570 if (
corr) stream <<
"PV correction: " <<
endl << *
corr;
1583 "TanPixelToRaDec::fit is NOT implemented (although it is doable)) ");
1591 :
BaseTanWcs(pixToTan, tangentPoint, corrections) {}
1607 const Frame ®ion)
1620 double &yTangentPlane)
const {
1624 corr->apply(xPixel, yPixel, xtmp, ytmp);
1638 stream <<
" tangent point " << tp.
x <<
' ' << tp.
y <<
endl;
1641 if (
corr) stream <<
"PV correction: " <<
endl << *
corr;
1654 "TanSipPixelToRaDec::fit is NOT implemented (although it is doable)) ");
1661 : linTan2Pix(tan2Pix) {
1668 ra0 = deg2rad(tangentPoint.
x);
1669 dec0 = deg2rad(tangentPoint.
y);
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);
1712 double l = sinda * coss;
1713 double m = sins * sin0 + coss * cos0 * cosda;
1715 m = (sins * cos0 - coss * sin0 * cosda) /
m;
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;
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;
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);
1746 double l =
std::sin(ra - ra0) * coss;
1747 double m = sins * sin0 + coss * cos0 *
std::cos(ra - ra0);
1749 m = (sins * cos0 - coss * sin0 *
std::cos(ra - ra0)) /
m;
1753 linTan2Pix.
apply(l,
m, xOut, yOut);
1762 stream <<
"tan2pix " << linTan2Pix <<
std::endl;
1763 stream <<
"tangent point " << tp.
x <<
' ' << tp.
y <<
endl;
1782 "TanRaDecToPixel::fit is NOT implemented (although it is doable)) ");
1789 : _userFun(userFun), _userData(userData) {}
1792 _userFun(xIn, yIn, xOut, yOut, _userData);
1796 stream <<
"UserTransform with user function @ " << _userFun <<
"and userData@ " << _userData <<
endl;
1801 "UserTransform::fit is NOT implemented (and will never be)) ");
1815 " astrometryTransformRead : cannot open " + fileName);
1831 "astrometryTransformRead : could not find a AstrometryTransformtype");
1832 if (
type ==
"AstrometryTransformIdentity") {
1836 }
else if (
type ==
"AstrometryTransformPolynomial") {
1842 " astrometryTransformRead : No reader for AstrometryTransform type " +
type);
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
LSST DM logging module built on log4cxx.
#define LOGLS_WARN(logger, message)
Log a warn-level message using an iostream-based interface.
#define LOGL_WARN(logger, message...)
Log a warn-level message using a varargs/printf style interface.
#define LOGLS_FATAL(logger, message)
Log a fatal-level message using an iostream-based interface.
#define LOG_GET(logger)
Returns a Log object associated with logger.
#define LOGL_ERROR(logger, message...)
Log a error-level message using a varargs/printf style interface.
#define LOGLS_TRACE(logger, message)
Log a trace-level message using an iostream-based interface.
#define LOGL_FATAL(logger, message...)
Log a fatal-level message using a varargs/printf style interface.
table::PointKey< double > crpix
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.
rectangle with sides parallel to axes.
Point getCenter() const
Center of the frame.
double getHeight() const
size along y axis
double xMin
coordinate of boundary.
double getWidth() const
size along x axis
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.
The transformation that handles pixels to sideral transformations (Gnomonic, possibly with polynomial...
std::unique_ptr< AstrometryTransform > inverseTransform(const double precision, const Frame ®ion) 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 ®ion) 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 ®ion) const
Overload the "generic routine" (available for all AstrometryTransform types.
std::unique_ptr< AstrometryTransform > inverseTransform(const double precision, const Frame ®ion) 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 ®ion) 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)
virtual char const * what(void) const noexcept
Return a character string summarizing this exception.
Reports invalid arguments.
Reports errors in the logical structure of the program.
Reports attempts to access elements using an invalid key.
Reports errors that are due to events beyond the control of the program.
T emplace_back(T... args)
std::ostream & operator<<(std::ostream &os, CameraSysPrefix const &detSysPrefix)
void write(OutputArchiveHandle &handle) const override
FilterProperty & operator=(FilterProperty const &)=default
Low-level polynomials (including special polynomials) in C++.
AstrometryTransformLinear normalizeCoordinatesTransform(const Frame &frame)
Returns the transformation that maps the input frame along both axes to [-1,1].
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)
A base class for image defects.
T setprecision(T... args)
table::Key< table::Array< int > > degree