LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
SipTransform.cc
Go to the documentation of this file.
1 // -*- LSST-C++ -*-
2 
3 /*
4  * LSST Data Management System
5  * Copyright 2016 LSST/AURA
6  *
7  * This product includes software developed by the
8  * LSST Project (http://www.lsst.org/).
9  *
10  * This program is free software: you can redistribute it and/or modify
11  * it under the terms of the GNU General Public License as published by
12  * the Free Software Foundation, either version 3 of the License, or
13  * (at your option) any later version.
14  *
15  * This program is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the LSST License Statement and
21  * the GNU General Public License along with this program. If not,
22  * see <http://www.lsstcorp.org/LegalNotices/>.
23  */
24 
25 #include <sstream>
26 
27 #include "lsst/geom/Point.h"
28 #include "lsst/geom/Angle.h"
29 #include "lsst/geom/SpherePoint.h"
32 #include "lsst/afw/geom/SkyWcs.h"
35 
36 namespace lsst {
37 namespace meas {
38 namespace astrom {
39 
41  // The implementation for transformPixels is identical for
42  // SipForwardTransform and SipReverseTransform. That's pretty obvious for
43  // the pixel origin and CD matrix, which are the same in both cases, but
44  // it wasn't obvious to me until I did the math that the polynomial
45  // transforms are composed with the affine transform the same way.
46  auto sInv = s.inverted();
47  _pixelOrigin = s.getLinear()(_pixelOrigin - sInv.getTranslation());
48  _cdMatrix = _cdMatrix * sInv.getLinear();
49  _poly = compose(s.getLinear(), compose(getPoly(), sInv.getLinear()));
50 }
51 
53  geom::Point2D const& pixelOrigin,
54  geom::LinearTransform const& cdMatrix) {
55  auto forwardSipPoly = compose(geom::AffineTransform(cdMatrix.inverted()),
57  // Subtracting 1 here accounts for the extra terms outside the sum in the
58  // transform definition (see class docs) - note that you can fold those
59  // terms into the sum by adding 1 from the A_10 and B_01 terms.
60  forwardSipPoly._xCoeffs(1, 0) -= 1;
61  forwardSipPoly._yCoeffs(0, 1) -= 1;
62  return SipForwardTransform(pixelOrigin, cdMatrix, forwardSipPoly);
63 }
64 
66  geom::Point2D const& pixelOrigin,
67  geom::LinearTransform const& cdMatrix) {
68  auto forwardSipPoly =
70  compose(scaled.getPoly(),
71  scaled.getInputScaling() * geom::AffineTransform(geom::Extent2D(pixelOrigin))));
72  // Account for the terms outside the sum in the definition (see comment
73  // earlier in the file for more explanation).
74  forwardSipPoly._xCoeffs(1, 0) -= 1;
75  forwardSipPoly._yCoeffs(0, 1) -= 1;
76  return SipForwardTransform(pixelOrigin, cdMatrix, forwardSipPoly);
77 }
78 
80  geom::Point2D pixelOrigin(-scaled.getOutputScalingInverse().getTranslation());
82  return convert(scaled, pixelOrigin, cdMatrix);
83 }
84 
87  return geom::AffineTransform(_cdMatrix) * (geom::AffineTransform() + _poly.linearize(tail(in))) * tail;
88 }
89 
92  return getCdMatrix()(geom::Extent2D(duv) + getPoly()(duv));
93 }
94 
97  result.transformPixelsInPlace(s);
98  return result;
99 }
100 
102  geom::Point2D const& pixelOrigin,
103  geom::LinearTransform const& cdMatrix) {
104  auto reverseSipPoly = compose(geom::AffineTransform(-geom::Extent2D(pixelOrigin)),
105  compose(poly, geom::AffineTransform(cdMatrix)));
106  // Account for the terms outside the sum in the definition (see comment
107  // earlier in the file for more explanation).
108  reverseSipPoly._xCoeffs(1, 0) -= 1;
109  reverseSipPoly._yCoeffs(0, 1) -= 1;
110  return SipReverseTransform(pixelOrigin, cdMatrix, reverseSipPoly);
111 }
112 
114  geom::Point2D const& pixelOrigin,
115  geom::LinearTransform const& cdMatrix) {
116  auto reverseSipPoly =
118  compose(scaled.getPoly(), scaled.getInputScaling() * geom::AffineTransform(cdMatrix)));
119  // Account for the terms outside the sum in the definition (see comment
120  // earlier in the file for more explanation).
121  reverseSipPoly._xCoeffs(1, 0) -= 1;
122  reverseSipPoly._yCoeffs(0, 1) -= 1;
123  return SipReverseTransform(pixelOrigin, cdMatrix, reverseSipPoly);
124 }
125 
127  return convert(scaled, geom::Point2D(scaled.getOutputScalingInverse().getTranslation()),
128  scaled.getInputScaling().getLinear());
129 }
130 
133  result.transformPixelsInPlace(s);
134  result._cdInverse = result._cdMatrix.inverted();
135  return result;
136 }
137 
140  (geom::AffineTransform() + _poly.linearize(_cdInverse(in))) * _cdInverse;
141 }
142 
144  geom::Point2D UV = _cdInverse(xy);
145  return geom::Extent2D(UV) + geom::Extent2D(getPixelOrigin()) + getPoly()(UV);
146 }
147 
149  SipReverseTransform const& sipReverse,
150  geom::SpherePoint const& skyOrigin) {
151  if (!sipForward.getPixelOrigin().asEigen().isApprox(sipReverse.getPixelOrigin().asEigen())) {
152  std::ostringstream oss;
153  oss << "SIP forward and reverse transforms have inconsistent CRPIX: " << sipForward.getPixelOrigin()
154  << " != " << sipReverse.getPixelOrigin();
156  }
157  if (!sipForward.getCdMatrix().getMatrix().isApprox(sipReverse.getCdMatrix().getMatrix())) {
158  std::ostringstream oss;
159  oss << "SIP forward and reverse transforms have inconsistent CD matrix: " << sipForward.getCdMatrix()
160  << "\n!=\n"
161  << sipReverse.getCdMatrix();
163  }
164  Eigen::MatrixXd sipA(ndarray::asEigenMatrix(sipForward.getPoly().getXCoeffs()));
165  Eigen::MatrixXd sipB(ndarray::asEigenMatrix(sipForward.getPoly().getYCoeffs()));
166  Eigen::MatrixXd sipAP(ndarray::asEigenMatrix(sipReverse.getPoly().getXCoeffs()));
167  Eigen::MatrixXd sipBP(ndarray::asEigenMatrix(sipReverse.getPoly().getYCoeffs()));
168 
169  return afw::geom::makeTanSipWcs(sipForward.getPixelOrigin(), skyOrigin,
170  sipForward.getCdMatrix().getMatrix(), sipA, sipB, sipAP, sipBP);
171 }
172 
174  geom::AffineTransform const& s) {
175  auto affineTransform22 = afw::geom::makeTransform(s);
176  return afw::geom::makeModifiedWcs(*affineTransform22->inverted(), wcs, true);
177 }
178 
180  geom::Extent2I const& dimensions) {
181  geom::Extent2D offset;
182  switch (nQuarter % 4) {
183  case 0:
184  offset = geom::Extent2D(0, 0);
185  break;
186  case 1:
187  offset = geom::Extent2D(dimensions.getY() - 1, 0);
188  break;
189  case 2:
190  offset = geom::Extent2D(dimensions - geom::Extent2I(1, 1));
191  break;
192  case 3:
193  offset = geom::Extent2D(0, dimensions.getX() - 1);
194  break;
195  }
196  auto rot = geom::LinearTransform::makeRotation(nQuarter * 90.0 * geom::degrees);
197  return transformWcsPixels(wcs, geom::AffineTransform(rot, offset));
198 }
199 
200 } // namespace astrom
201 } // namespace meas
202 } // namespace lsst
py::object result
Definition: _schema.cc:429
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
afw::table::PointKey< int > dimensions
Definition: GaussianPsf.cc:48
table::Key< table::Array< std::uint8_t > > wcs
Definition: SkyWcs.cc:66
A 2-dimensional celestial WCS that transform pixels to ICRS RA/Dec, using the LSST standard for pixel...
Definition: SkyWcs.h:117
An affine coordinate transformation consisting of a linear transformation and an offset.
AffineTransform const inverted() const
Return the inverse transform.
Extent2D const & getTranslation() const noexcept
LinearTransform const & getLinear() const noexcept
EigenVector const & asEigen() const noexcept(IS_ELEMENT_NOTHROW_COPYABLE)
Return a fixed-size Eigen representation of the coordinate object.
A 2D linear coordinate transformation.
static LinearTransform makeRotation(Angle t) noexcept
LinearTransform const inverted() const
Return the inverse transform.
Matrix const & getMatrix() const noexcept
Point in an unspecified spherical coordinate system.
Definition: SpherePoint.h:57
A 2-d coordinate transform represented by a pair of standard polynomials (one for each coordinate).
geom::AffineTransform linearize(geom::Point2D const &in) const
Return an approximate affine transform at the given point.
ndarray::Array< double const, 2, 2 > getYCoeffs() const
2-D polynomial coefficients that compute the output x coordinate.
ndarray::Array< double const, 2, 2 > getXCoeffs() const
2-D polynomial coefficients that compute the output x coordinate.
A 2-d coordinate transform represented by a lazy composition of an AffineTransform,...
geom::AffineTransform const & getOutputScalingInverse() const
Return the affine transform applied to points after the polynomial transform.
geom::AffineTransform const & getInputScaling() const
Return the first affine transform applied to input points.
PolynomialTransform const & getPoly() const
Return the polynomial transform applied after the input scaling.
A transform that maps pixel coordinates to intermediate world coordinates according to the SIP conven...
Definition: SipTransform.h:136
SipForwardTransform transformPixels(geom::AffineTransform const &s) const
Return a new forward SIP transform that includes a transformation of the pixel coordinate system by t...
Definition: SipTransform.cc:95
geom::AffineTransform linearize(geom::Point2D const &in) const
Return an approximate affine transform at the given point.
Definition: SipTransform.cc:85
geom::Point2D operator()(geom::Point2D const &uv) const
Apply the transform to a point.
Definition: SipTransform.cc:90
static SipForwardTransform convert(PolynomialTransform const &poly, geom::Point2D const &pixelOrigin, geom::LinearTransform const &cdMatrix)
Convert a PolynomialTransform to an equivalent SipForwardTransform.
Definition: SipTransform.cc:52
SipForwardTransform(geom::Point2D const &pixelOrigin, geom::LinearTransform const &cdMatrix, PolynomialTransform const &forwardSipPoly)
Construct a SipForwardTransform from its components.
Definition: SipTransform.h:175
A transform that maps intermediate world coordinates to pixel coordinates according to the SIP conven...
Definition: SipTransform.h:246
SipReverseTransform(geom::Point2D const &pixelOrigin, geom::LinearTransform const &cdMatrix, PolynomialTransform const &reverseSipPoly)
Construct a SipReverseTransform from its components.
Definition: SipTransform.h:285
SipReverseTransform transformPixels(geom::AffineTransform const &s) const
Return a new reverse SIP transform that includes a transformation of the pixel coordinate system by t...
geom::AffineTransform linearize(geom::Point2D const &in) const
Return an approximate affine transform at the given point.
static SipReverseTransform convert(PolynomialTransform const &poly, geom::Point2D const &pixelOrigin, geom::LinearTransform const &cdMatrix)
Convert a PolynomialTransform to an equivalent SipReverseTransform.
geom::Point2D operator()(geom::Point2D const &xy) const
Apply the transform to a point.
void transformPixelsInPlace(geom::AffineTransform const &s)
Definition: SipTransform.cc:40
geom::LinearTransform const & getCdMatrix() const
Return the CD matrix of the transform.
Definition: SipTransform.h:56
geom::Point2D const & getPixelOrigin() const
Return the pixel origin (CRPIX, but zero-indexed) of the transform.
Definition: SipTransform.h:51
geom::LinearTransform _cdMatrix
Definition: SipTransform.h:94
PolynomialTransform const & getPoly() const
Return the polynomial component of the transform (A,B) or (AP,BP).
Definition: SipTransform.h:61
Reports invalid arguments.
Definition: Runtime.h:66
std::shared_ptr< TransformPoint2ToPoint2 > makeTransform(lsst::geom::AffineTransform const &affine)
Wrap an lsst::geom::AffineTransform as a Transform.
std::shared_ptr< SkyWcs > makeTanSipWcs(lsst::geom::Point2D const &crpix, lsst::geom::SpherePoint const &crval, Eigen::Matrix2d const &cdMatrix, Eigen::MatrixXd const &sipA, Eigen::MatrixXd const &sipB)
Construct a TAN-SIP SkyWcs with forward SIP distortion terms and an iterative inverse.
Definition: SkyWcs.cc:538
std::shared_ptr< SkyWcs > makeModifiedWcs(TransformPoint2ToPoint2 const &pixelTransform, SkyWcs const &wcs, bool modifyActualPixels)
Create a new SkyWcs whose pixels are transformed by pixelTransform, as described below.
Definition: SkyWcs.cc:486
Low-level polynomials (including special polynomials) in C++.
Definition: Basis1d.h:26
Extent< double, 2 > Extent2D
Definition: Extent.h:400
constexpr AngleUnit degrees
constant with units of degrees
Definition: Angle.h:109
std::shared_ptr< afw::geom::SkyWcs > transformWcsPixels(afw::geom::SkyWcs const &wcs, geom::AffineTransform const &s)
Create a new SkyWcs whose pixel coordinate system has been transformed via an affine transform.
PolynomialTransform compose(geom::AffineTransform const &t1, PolynomialTransform const &t2)
Return a PolynomialTransform that is equivalent to the composition t1(t2())
std::shared_ptr< afw::geom::SkyWcs > rotateWcsPixelsBy90(afw::geom::SkyWcs const &wcs, int nQuarter, geom::Extent2I const &dimensions)
Return a new SkyWcs that represents a rotation of the image it corresponds to about the image's cente...
std::shared_ptr< afw::geom::SkyWcs > makeWcs(SipForwardTransform const &sipForward, SipReverseTransform const &sipReverse, geom::SpherePoint const &skyOrigin)
Create a new TAN SIP Wcs from a pair of SIP transforms and the sky origin.
A base class for image defects.
T str(T... args)