LSSTApplications  20.0.0
LSSTDataManagementBasePackage
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
lsst::geom::degrees
constexpr AngleUnit degrees
constant with units of degrees
Definition: Angle.h:109
lsst::geom::LinearTransform::makeRotation
static LinearTransform makeRotation(Angle t) noexcept
Definition: LinearTransform.h:102
SipTransform.h
std::shared_ptr
STL class.
wcs
table::Key< table::Array< std::uint8_t > > wcs
Definition: SkyWcs.cc:71
lsst::geom::LinearTransform
A 2D linear coordinate transformation.
Definition: LinearTransform.h:69
AffineTransform.h
lsst::meas::astrom::transformWcsPixels
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.
Definition: SipTransform.cc:173
lsst::meas::astrom::PolynomialTransform::getYCoeffs
ndarray::Array< double const, 2, 2 > getYCoeffs() const
2-D polynomial coefficients that compute the output x coordinate.
Definition: PolynomialTransform.h:123
lsst::meas::astrom::SipForwardTransform::operator()
geom::Point2D operator()(geom::Point2D const &uv) const
Apply the transform to a point.
Definition: SipTransform.cc:90
lsst::meas::astrom::SipReverseTransform::transformPixels
SipReverseTransform transformPixels(geom::AffineTransform const &s) const
Return a new reverse SIP transform that includes a transformation of the pixel coordinate system by t...
Definition: SipTransform.cc:131
lsst::meas::astrom::ScaledPolynomialTransform::getOutputScalingInverse
geom::AffineTransform const & getOutputScalingInverse() const
Return the affine transform applied to points after the polynomial transform.
Definition: PolynomialTransform.h:215
SkyWcs.h
lsst::geom::AffineTransform::inverted
AffineTransform const inverted() const
Return the inverse transform.
Definition: AffineTransform.cc:53
lsst::meas::astrom::SipTransformBase::getPoly
PolynomialTransform const & getPoly() const
Return the polynomial component of the transform (A,B) or (AP,BP).
Definition: SipTransform.h:61
lsst::afw::geom::makeModifiedWcs
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:491
lsst::geom::AffineTransform::getTranslation
Extent2D const & getTranslation() const noexcept
Definition: AffineTransform.h:152
lsst::afw::geom::SkyWcs
A 2-dimensional celestial WCS that transform pixels to ICRS RA/Dec, using the LSST standard for pixel...
Definition: SkyWcs.h:117
lsst::meas::astrom::ScaledPolynomialTransform
A 2-d coordinate transform represented by a lazy composition of an AffineTransform,...
Definition: PolynomialTransform.h:157
lsst::meas::astrom::SipForwardTransform::SipForwardTransform
SipForwardTransform(geom::Point2D const &pixelOrigin, geom::LinearTransform const &cdMatrix, PolynomialTransform const &forwardSipPoly)
Construct a SipForwardTransform from its components.
Definition: SipTransform.h:175
lsst::meas::astrom::SipReverseTransform::convert
static SipReverseTransform convert(PolynomialTransform const &poly, geom::Point2D const &pixelOrigin, geom::LinearTransform const &cdMatrix)
Convert a PolynomialTransform to an equivalent SipReverseTransform.
Definition: SipTransform.cc:101
lsst::meas::astrom::SipForwardTransform
A transform that maps pixel coordinates to intermediate world coordinates according to the SIP conven...
Definition: SipTransform.h:136
lsst::geom::AffineTransform
An affine coordinate transformation consisting of a linear transformation and an offset.
Definition: AffineTransform.h:75
Angle.h
lsst::afw::geom::makeTanSipWcs
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:543
lsst::meas::astrom::SipTransformBase::_pixelOrigin
geom::Point2D _pixelOrigin
Definition: SipTransform.h:93
lsst::meas::astrom::SipTransformBase::_poly
PolynomialTransform _poly
Definition: SipTransform.h:95
transformFactory.h
lsst::geom::AffineTransform::getLinear
LinearTransform const & getLinear() const noexcept
Definition: AffineTransform.h:155
lsst::meas::astrom::SipReverseTransform::operator()
geom::Point2D operator()(geom::Point2D const &xy) const
Apply the transform to a point.
Definition: SipTransform.cc:143
LinearTransform.h
lsst::geom::LinearTransform::inverted
LinearTransform const inverted() const
Return the inverse transform.
Definition: LinearTransform.cc:45
lsst::meas::astrom::SipReverseTransform::SipReverseTransform
SipReverseTransform(geom::Point2D const &pixelOrigin, geom::LinearTransform const &cdMatrix, PolynomialTransform const &reverseSipPoly)
Construct a SipReverseTransform from its components.
Definition: SipTransform.h:285
lsst::meas::astrom::PolynomialTransform::linearize
geom::AffineTransform linearize(geom::Point2D const &in) const
Return an approximate affine transform at the given point.
Definition: PolynomialTransform.cc:129
dimensions
afw::table::PointKey< int > dimensions
Definition: GaussianPsf.cc:49
lsst::meas::astrom::SipTransformBase::transformPixelsInPlace
void transformPixelsInPlace(geom::AffineTransform const &s)
Definition: SipTransform.cc:40
lsst::meas::astrom::PolynomialTransform
A 2-d coordinate transform represented by a pair of standard polynomials (one for each coordinate).
Definition: PolynomialTransform.h:45
result
py::object result
Definition: _schema.cc:429
lsst::geom::LinearTransform::getMatrix
Matrix const & getMatrix() const noexcept
Definition: LinearTransform.h:151
lsst::meas::astrom::SipForwardTransform::linearize
geom::AffineTransform linearize(geom::Point2D const &in) const
Return an approximate affine transform at the given point.
Definition: SipTransform.cc:85
lsst::meas::astrom::SipTransformBase::getPixelOrigin
geom::Point2D const & getPixelOrigin() const
Return the pixel origin (CRPIX, but zero-indexed) of the transform.
Definition: SipTransform.h:51
lsst
A base class for image defects.
Definition: imageAlgorithm.dox:1
lsst::meas::astrom::ScaledPolynomialTransform::getPoly
PolynomialTransform const & getPoly() const
Return the polynomial transform applied after the input scaling.
Definition: PolynomialTransform.h:209
LSST_EXCEPT
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
std::ostringstream
STL class.
lsst::meas::astrom::makeWcs
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.
Definition: SipTransform.cc:148
lsst::meas::astrom::SipTransformBase::getCdMatrix
geom::LinearTransform const & getCdMatrix() const
Return the CD matrix of the transform.
Definition: SipTransform.h:56
lsst::geom::polynomials
Definition: Basis1d.h:26
lsst::pex::exceptions::InvalidParameterError
Reports invalid arguments.
Definition: Runtime.h:66
lsst::geom::Point< double, 2 >
lsst::meas::astrom::SipForwardTransform::convert
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
lsst::meas::astrom::PolynomialTransform::getXCoeffs
ndarray::Array< double const, 2, 2 > getXCoeffs() const
2-D polynomial coefficients that compute the output x coordinate.
Definition: PolynomialTransform.h:115
lsst::geom::Extent2D
Extent< double, 2 > Extent2D
Definition: Extent.h:400
std::ostringstream::str
T str(T... args)
SpherePoint.h
Point.h
lsst::geom::SpherePoint
Point in an unspecified spherical coordinate system.
Definition: SpherePoint.h:57
lsst::meas::astrom::rotateWcsPixelsBy90
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...
Definition: SipTransform.cc:179
lsst::meas::astrom::SipReverseTransform
A transform that maps intermediate world coordinates to pixel coordinates according to the SIP conven...
Definition: SipTransform.h:246
lsst::meas::astrom::SipReverseTransform::linearize
geom::AffineTransform linearize(geom::Point2D const &in) const
Return an approximate affine transform at the given point.
Definition: SipTransform.cc:138
lsst::meas::astrom::SipTransformBase::_cdMatrix
geom::LinearTransform _cdMatrix
Definition: SipTransform.h:94
lsst::afw::geom::makeTransform
std::shared_ptr< TransformPoint2ToPoint2 > makeTransform(lsst::geom::AffineTransform const &affine)
Wrap an lsst::geom::AffineTransform as a Transform.
Definition: transformFactory.cc:154
lsst::geom::Extent< double, 2 >
lsst::meas::astrom::ScaledPolynomialTransform::getInputScaling
geom::AffineTransform const & getInputScaling() const
Return the first affine transform applied to input points.
Definition: PolynomialTransform.h:212
lsst::meas::astrom::SipForwardTransform::transformPixels
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
lsst::meas::astrom::compose
PolynomialTransform compose(geom::AffineTransform const &t1, PolynomialTransform const &t2)
Return a PolynomialTransform that is equivalent to the composition t1(t2())
Definition: PolynomialTransform.cc:214