LSST Applications g0f08755f38+82efc23009,g12f32b3c4e+e7bdf1200e,g1653933729+a8ce1bb630,g1a0ca8cf93+50eff2b06f,g28da252d5a+52db39f6a5,g2bbee38e9b+37c5a29d61,g2bc492864f+37c5a29d61,g2cdde0e794+c05ff076ad,g3156d2b45e+41e33cbcdc,g347aa1857d+37c5a29d61,g35bb328faa+a8ce1bb630,g3a166c0a6a+37c5a29d61,g3e281a1b8c+fb992f5633,g414038480c+7f03dfc1b0,g41af890bb2+11b950c980,g5fbc88fb19+17cd334064,g6b1c1869cb+12dd639c9a,g781aacb6e4+a8ce1bb630,g80478fca09+72e9651da0,g82479be7b0+04c31367b4,g858d7b2824+82efc23009,g9125e01d80+a8ce1bb630,g9726552aa6+8047e3811d,ga5288a1d22+e532dc0a0b,gae0086650b+a8ce1bb630,gb58c049af0+d64f4d3760,gc28159a63d+37c5a29d61,gcf0d15dbbd+2acd6d4d48,gd7358e8bfb+778a810b6e,gda3e153d99+82efc23009,gda6a2b7d83+2acd6d4d48,gdaeeff99f8+1711a396fd,ge2409df99d+6b12de1076,ge79ae78c31+37c5a29d61,gf0baf85859+d0a5978c5a,gf3967379c6+4954f8c433,gfb92a5be7c+82efc23009,gfec2e1e490+2aaed99252,w.2024.46
LSST Data Management Base Package
Loading...
Searching...
No Matches
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"
35
36namespace lsst {
37namespace meas {
38namespace 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
82 return convert(scaled, pixelOrigin, cdMatrix);
83}
84
89
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)),
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
130
133 result.transformPixelsInPlace(s);
134 result._cdInverse = result._cdMatrix.inverted();
135 return result;
136}
137
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())) {
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())) {
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
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.
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.
Matrix const & getMatrix() const noexcept
static LinearTransform makeRotation(Angle t) noexcept
LinearTransform const inverted() const
Return the inverse transform.
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 > getXCoeffs() const
2-D polynomial coefficients that compute the output x coordinate.
ndarray::Array< double const, 2, 2 > getYCoeffs() 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.
PolynomialTransform const & getPoly() const
Return the polynomial transform applied after the input scaling.
geom::AffineTransform const & getInputScaling() const
Return the first affine transform applied to input points.
A transform that maps pixel coordinates to intermediate world coordinates according to the SIP conven...
SipForwardTransform transformPixels(geom::AffineTransform const &s) const
Return a new forward 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.
geom::Point2D operator()(geom::Point2D const &uv) const
Apply the transform to a point.
static SipForwardTransform convert(PolynomialTransform const &poly, geom::Point2D const &pixelOrigin, geom::LinearTransform const &cdMatrix)
Convert a PolynomialTransform to an equivalent SipForwardTransform.
SipForwardTransform(geom::Point2D const &pixelOrigin, geom::LinearTransform const &cdMatrix, PolynomialTransform const &forwardSipPoly)
Construct a SipForwardTransform from its components.
A transform that maps intermediate world coordinates to pixel coordinates according to the SIP conven...
SipReverseTransform(geom::Point2D const &pixelOrigin, geom::LinearTransform const &cdMatrix, PolynomialTransform const &reverseSipPoly)
Construct a SipReverseTransform from its components.
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.
geom::Point2D const & getPixelOrigin() const
Return the pixel origin (CRPIX, but zero-indexed) of the transform.
PolynomialTransform const & getPoly() const
Return the polynomial component of the transform (A,B) or (AP,BP).
void transformPixelsInPlace(geom::AffineTransform const &s)
geom::LinearTransform const & getCdMatrix() const
Return the CD matrix of the transform.
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++.
AngleUnit constexpr degrees
constant with units of degrees
Definition Angle.h:110
Extent< double, 2 > Extent2D
Definition Extent.h:400
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.
T str(T... args)