LSST Applications g0b6bd0c080+a72a5dd7e6,g1182afd7b4+2a019aa3bb,g17e5ecfddb+2b8207f7de,g1d67935e3f+06cf436103,g38293774b4+ac198e9f13,g396055baef+6a2097e274,g3b44f30a73+6611e0205b,g480783c3b1+98f8679e14,g48ccf36440+89c08d0516,g4b93dc025c+98f8679e14,g5c4744a4d9+a302e8c7f0,g613e996a0d+e1c447f2e0,g6c8d09e9e7+25247a063c,g7271f0639c+98f8679e14,g7a9cd813b8+124095ede6,g9d27549199+a302e8c7f0,ga1cf026fa3+ac198e9f13,ga32aa97882+7403ac30ac,ga786bb30fb+7a139211af,gaa63f70f4e+9994eb9896,gabf319e997+ade567573c,gba47b54d5d+94dc90c3ea,gbec6a3398f+06cf436103,gc6308e37c7+07dd123edb,gc655b1545f+ade567573c,gcc9029db3c+ab229f5caf,gd01420fc67+06cf436103,gd877ba84e5+06cf436103,gdb4cecd868+6f279b5b48,ge2d134c3d5+cc4dbb2e3f,ge448b5faa6+86d1ceac1d,gecc7e12556+98f8679e14,gf3ee170dca+25247a063c,gf4ac96e456+ade567573c,gf9f5ea5b4d+ac198e9f13,gff490e6085+8c2580be5c,w.2022.27
LSST Data Management Base Package
SipApproximation.cc
Go to the documentation of this file.
1// -*- lsst-c++ -*-
2/*
3 * Developed for the LSST Data Management System.
4 * This product includes software developed by the LSST Project
5 * (https://www.lsst.org).
6 * See the COPYRIGHT file at the top-level directory of this distribution
7 * for details of code ownership.
8 *
9 * This program is free software: you can redistribute it and/or modify
10 * it under the terms of the GNU General Public License as published by
11 * the Free Software Foundation, either version 3 of the License, or
12 * (at your option) any later version.
13 *
14 * This program is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the GNU General Public License
20 * along with this program. If not, see <https://www.gnu.org/licenses/>.
21 */
22
23#include <algorithm>
24#include <vector>
25
26#include "Eigen/QR"
29
30namespace lsst { namespace afw { namespace geom {
31
33
34namespace {
35
37 int order,
38 lsst::geom::Box2D const & box,
39 double svdThreshold,
42) {
43 // The scaled polynomial basis evaluates polynomials after mapping the
44 // input coordinates from the given box to [-1, 1]x[-1, 1] (for numerical
45 // stability).
47 auto workspace = basis.makeWorkspace();
48 Eigen::MatrixXd matrix = Eigen::MatrixXd::Zero(input.size(), basis.size());
49 Eigen::VectorXd xRhs(input.size());
50 Eigen::VectorXd yRhs(input.size());
51 for (int i = 0; i < matrix.rows(); ++i) {
52 basis.fill(input[i], matrix.row(i), workspace);
53 auto rhs = output[i] - input[i];
54 xRhs[i] = rhs.getX();
55 yRhs[i] = rhs.getY();
56 }
57 // Since we're not trying to null the zeroth- and first-order terms, the
58 // solution is just linear least squares, and we can do that with SVD.
59 auto decomp = matrix.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV);
60 if (svdThreshold >= 0) {
61 decomp.setThreshold(svdThreshold);
62 }
63 auto scaledX = makeFunction2d(basis, decomp.solve(xRhs));
64 auto scaledY = makeFunction2d(basis, decomp.solve(yRhs));
65 // On return, we simplify the polynomials by moving the remapping transform
66 // into the coefficients themselves.
67 return std::make_pair(simplified(scaledX), simplified(scaledY));
68}
69
70// Return a vector of points on a grid, covering the given bounding box.
72 lsst::geom::Extent2I const & shape) {
73 if (shape.getX() <= 0 || shape.getY() <= 0) {
74 throw LSST_EXCEPT(
75 pex::exceptions::InvalidParameterError,
76 "Grid shape must be positive."
77 );
78 }
80 points.reserve(shape.getX()*shape.getY());
81 double const dx = bbox.getWidth()/shape.getX();
82 double const dy = bbox.getHeight()/shape.getY();
83 for (int iy = 0; iy < shape.getY(); ++iy) {
84 double const y = bbox.getMinY() + iy*dy;
85 for (int ix = 0; ix < shape.getX(); ++ix) {
86 points.emplace_back(bbox.getMinX() + ix*dx, y);
87 }
88 }
89 return points;
90}
91
92// Make a polynomial object (with packed coefficients) from a square coefficients matrix.
93poly::PolynomialFunction2dYX makePolynomialFromCoeffMatrix(ndarray::Array<double const, 2> const & coeffs) {
94 LSST_THROW_IF_NE(coeffs.getSize<0>(), coeffs.getSize<1>(), pex::exceptions::InvalidParameterError,
95 "Coefficient matrix must be square (%d != %d).");
96 poly::PolynomialBasis2dYX basis(coeffs.getSize<0>() - 1);
97 Eigen::VectorXd packed(basis.size());
98 for (auto const & i : basis.getIndices()) {
99 packed[i.flat] = coeffs[i.nx][i.ny];
100 }
101 return poly::PolynomialFunction2dYX(basis, packed);
102}
103
104} // anonymous
105
106// Private implementation object for SipApproximation that manages the grid of points on which
107// we evaluate the exact transform.
109
110 // Set up the grid.
111 Grid(lsst::geom::Extent2I const & shape_, SipApproximation const & parent);
112
113 lsst::geom::Extent2I const shape; // number of grid points in each dimension
114 std::vector<lsst::geom::Point2D> dpix1; // [pixel coords] - CRPIX
115 std::vector<lsst::geom::Point2D> siwc; // CD^{-1}([intermediate world coords])
116 std::vector<lsst::geom::Point2D> dpix2; // round-tripped version of dpix1 if useInverse, or exactly dpix1
117};
118
119// Private implementation object for SipApproximation that manages the solution
121
122 static std::unique_ptr<Solution> fit(int order_, double svdThreshold, SipApproximation const & parent);
123
127 poly::PolynomialFunction2dYX const & bp_) :
128 a(a_), b(b_), ap(ap_), bp(bp_)
129 {
130 LSST_THROW_IF_NE(a.getBasis().getOrder(), b.getBasis().getOrder(),
132 "A and B polynomials must have the same order (%d != %d).");
133 LSST_THROW_IF_NE(a.getBasis().getOrder(), ap.getBasis().getOrder(),
135 "A and AP polynomials must have the same order (%d != %d).");
136 LSST_THROW_IF_NE(a.getBasis().getOrder(), bp.getBasis().getOrder(),
138 "A and BP polynomials must have the same order (%d != %d).");
139 }
140
142
144
146 return dpix + lsst::geom::Extent2D(a(dpix, ws), b(dpix, ws));
147 }
148
150 return siwc + lsst::geom::Extent2D(ap(siwc, ws), bp(siwc, ws));
151 }
152
157};
158
160 shape(shape_),
161 dpix1(makeGrid(parent._bbox, shape)),
162 siwc(parent._pixelToIwc->applyForward(dpix1))
163{
164 // Apply the CRPIX offset to make pix1 into dpix1 (in-place)
165 std::for_each(dpix1.begin(), dpix1.end(), [&parent](lsst::geom::Point2D & p){ p -= parent._crpix; });
166
167 if (parent._useInverse) {
168 // Set from the given inverse of the given pixels-to-iwc transform
169 // Note that at this point, siwc is still just iwc, because the scaling by cdInv is later.
170 dpix2 = parent._pixelToIwc->applyInverse(siwc);
171 // Apply the CRPIX offset to make pix1 into dpix2 (in-place)
172 std::for_each(dpix2.begin(), dpix2.end(), [&parent](lsst::geom::Point2D & p){ p -= parent._crpix; });
173 } else {
174 // Just make dpix2 = dpix1, and hence fit to the true inverse of pixels-to-iwc.
175 dpix2 = dpix1;
176 }
177
178 // Apply the CD^{-1} transform to siwc
179 std::for_each(siwc.begin(), siwc.end(), [&parent](lsst::geom::Point2D & p){ p = parent._cdInv(p); });
180}
181
183 int order,
184 double svdThreshold,
185 SipApproximation const & parent
186) {
188 if (basis.size() > parent._grid->dpix1.size()) {
189 throw LSST_EXCEPT(
191 (boost::format("Number of parameters (%d) is larger than number of data points (%d)")
192 % (2*basis.size()) % (2*parent._grid->dpix1.size())).str()
193 );
194 }
195
196 lsst::geom::Box2D boxFwd(parent._bbox);
197 boxFwd.shift(-parent._crpix);
198 auto fwd = fitSipOneDirection(order, boxFwd, svdThreshold, parent._grid->dpix1, parent._grid->siwc);
199
200 lsst::geom::Box2D boxInv;
201 for (auto const & point : parent._grid->siwc) {
202 boxInv.include(point);
203 }
204 auto inv = fitSipOneDirection(order, boxInv, svdThreshold, parent._grid->siwc, parent._grid->dpix2);
205
206 return std::make_unique<Solution>(fwd.first, fwd.second, inv.first, inv.second);
207}
208
212 Eigen::Matrix2d const & cd,
213 lsst::geom::Box2D const & bbox,
214 lsst::geom::Extent2I const & gridShape,
215 int order,
216 bool useInverse,
217 double svdThreshold
218) :
219 _useInverse(useInverse),
220 _pixelToIwc(std::move(pixelToIwc)),
221 _bbox(bbox),
222 _crpix(crpix),
223 _cdInv(lsst::geom::LinearTransform(cd).inverted()),
224 _grid(new Grid(gridShape, *this)),
225 _solution(Solution::fit(order, svdThreshold, *this))
226{}
227
231 Eigen::Matrix2d const & cd,
232 lsst::geom::Box2D const & bbox,
233 lsst::geom::Extent2I const & gridShape,
234 ndarray::Array<double const, 2> const & a,
235 ndarray::Array<double const, 2> const & b,
236 ndarray::Array<double const, 2> const & ap,
237 ndarray::Array<double const, 2> const & bp,
238 bool useInverse
239) :
240 _useInverse(useInverse),
241 _pixelToIwc(std::move(pixelToIwc)),
242 _bbox(bbox),
243 _crpix(crpix),
244 _cdInv(lsst::geom::LinearTransform(cd).inverted()),
245 _grid(new Grid(gridShape, *this)),
246 _solution(
247 new Solution(
248 makePolynomialFromCoeffMatrix(a),
249 makePolynomialFromCoeffMatrix(b),
250 makePolynomialFromCoeffMatrix(ap),
251 makePolynomialFromCoeffMatrix(bp)
252 )
253 )
254{}
255
256SipApproximation::~SipApproximation() noexcept = default;
257
258int SipApproximation::getOrder() const noexcept {
259 return _solution->a.getBasis().getOrder();
260}
261
262double SipApproximation::getA(int p, int q) const {
263 return _solution->a[_solution->a.getBasis().index(p, q)];
264}
265
266double SipApproximation::getB(int p, int q) const {
267 return _solution->b[_solution->b.getBasis().index(p, q)];
268}
269
270double SipApproximation::getAP(int p, int q) const {
271 return _solution->ap[_solution->ap.getBasis().index(p, q)];
272}
273
274double SipApproximation::getBP(int p, int q) const {
275 return _solution->bp[_solution->bp.getBasis().index(p, q)];
276}
277
278namespace {
279
280template <typename F>
281Eigen::MatrixXd makeCoefficientMatrix(std::size_t order, F getter) {
282 Eigen::MatrixXd result = Eigen::MatrixXd::Zero(order + 1, order + 1);
283 for (std::size_t p = 0; p <= order; ++p) {
284 for (std::size_t q = 0; q <= order - p; ++q) {
285 result(p, q) = getter(p, q);
286 }
287 }
288 return result;
289}
290
291} // anonymous
292
293Eigen::MatrixXd SipApproximation::getA() const noexcept {
294 return makeCoefficientMatrix(
295 getOrder(),
296 [this](int p, int q) { return getA(p, q); }
297 );
298}
299
300Eigen::MatrixXd SipApproximation::getB() const noexcept {
301 return makeCoefficientMatrix(
302 getOrder(),
303 [this](int p, int q) { return getB(p, q); }
304 );
305}
306
307Eigen::MatrixXd SipApproximation::getAP() const noexcept {
308 return makeCoefficientMatrix(
309 getOrder(),
310 [this](int p, int q) { return getAP(p, q); }
311 );
312}
313
314Eigen::MatrixXd SipApproximation::getBP() const noexcept {
315 return makeCoefficientMatrix(
316 getOrder(),
317 [this](int p, int q) { return getBP(p, q); }
318 );
319}
320
322 auto cd = _cdInv.inverted();
323 auto ws = _solution->makeWorkspace();
324 return cd(_solution->applyForward(pix - _crpix, ws));
325}
326
328 std::vector<lsst::geom::Point2D> const & pix) const {
329 auto ws = _solution->makeWorkspace();
331 iwc.reserve(pix.size());
332 auto cd = _cdInv.inverted();
333 for (auto const & point : pix) {
334 iwc.push_back(cd(_solution->applyForward(point - _crpix, ws)));
335 }
336 return iwc;
337}
338
340 auto ws = _solution->makeWorkspace();
341 return _solution->applyInverse(_cdInv(iwc), ws) + _crpix;
342}
343
345 std::vector<lsst::geom::Point2D> const & iwc) const {
346 auto ws = _solution->makeWorkspace();
348 pix.reserve(iwc.size());
349 for (auto const & point : iwc) {
350 pix.push_back(_solution->applyInverse(_cdInv(point), ws) + _crpix);
351 }
352 return pix;
353}
354
356 return lsst::geom::Extent2D(_bbox.getWidth()/_grid->shape.getX(),
357 _bbox.getHeight()/_grid->shape.getY());
358}
359
361 return _grid->shape;
362}
363
365 _grid = std::make_unique<Grid>(shape, *this);
366}
367
369 // We shrink the grid spacing by the given factor, which is not the same
370 // as increasing the number of grid points by that factor, because there
371 // is one more grid point that step in each dimension.
372 lsst::geom::Extent2I unit(1);
373 updateGrid((_grid->shape - unit)*f + unit);
374}
375
376void SipApproximation::fit(int order, double svdThreshold) {
377 _solution = Solution::fit(order, svdThreshold, *this);
378}
379
381 std::pair<double, double> maxDiff(0.0, 0.0);
382 auto ws = _solution->makeWorkspace();
383 for (std::size_t i = 0; i < _grid->dpix1.size(); ++i) {
384 auto siwc2 = _solution->applyForward(_grid->dpix1[i], ws);
385 auto dpix2 = _solution->applyInverse(_grid->siwc[i], ws);
386 maxDiff.first = std::max(maxDiff.first, (_grid->siwc[i] - siwc2).computeNorm());
387 maxDiff.second = std::max(maxDiff.second, (_grid->dpix2[i] - dpix2).computeNorm());
388 }
389 return maxDiff;
390}
391
392}}} // namespace lsst::afw::geom
py::object result
Definition: _schema.cc:429
AmpInfoBoxKey bbox
Definition: Amplifier.cc:117
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
table::PointKey< double > crpix
Definition: OldWcs.cc:129
table::Key< table::Array< double > > cd
Definition: OldWcs.cc:130
int y
Definition: SpanSet.cc:48
table::Key< int > b
table::Key< int > a
#define LSST_THROW_IF_NE(N1, N2, EXC_CLASS, MSG)
Check whether the given values are equal, and throw an LSST Exception if they are not.
Definition: asserts.h:38
A fitter and results class for approximating a general Transform in a form compatible with FITS WCS p...
void fit(int order, double svdThreshold=-1)
Obtain a new solution at the given order with the current grid.
lsst::geom::Extent2I getGridShape() const noexcept
Return the number of grid points in x and y.
Eigen::MatrixXd getBP() const noexcept
Return the coefficients of the reverse transform polynomial.
Eigen::MatrixXd getAP() const noexcept
Return the coefficients of the reverse transform polynomial.
Eigen::MatrixXd getA() const noexcept
Return the coefficients of the forward transform polynomial.
int getOrder() const noexcept
Return the polynomial order of the current solution (same for forward and reverse).
lsst::geom::Extent2D getGridStep() const noexcept
Return the distance between grid points in pixels.
void updateGrid(lsst::geom::Extent2I const &shape)
Update the grid to the given number of points in x and y.
void refineGrid(int factor=2)
Update the grid by making it finer by a given integer factor.
lsst::geom::Point2D applyInverse(lsst::geom::Point2D const &iwcs) const
Convert a point from intermediate world coordinates to pixels.
Eigen::MatrixXd getB() const noexcept
Return the coefficients of the forward transform polynomial.
std::pair< double, double > computeMaxDeviation() const noexcept
Return the maximum deviation of the solution from the exact transform on the current grid.
lsst::geom::Point2D applyForward(lsst::geom::Point2D const &pix) const
Convert a point from pixels to intermediate world coordinates.
SipApproximation(std::shared_ptr< TransformPoint2ToPoint2 > pixelToIwc, lsst::geom::Point2D const &crpix, Eigen::Matrix2d const &cd, lsst::geom::Box2D const &bbox, lsst::geom::Extent2I const &gridShape, int order, bool useInverse=true, double svdThreshold=-1)
Construct a new approximation by fitting on a grid of points.
A floating-point coordinate rectangle geometry.
Definition: Box.h:413
void shift(Extent2D const &offset)
Shift the position of the box by the given offset.
Definition: Box.cc:350
double getWidth() const noexcept
1-d interval accessors
Definition: Box.h:529
void include(Point2D const &point) noexcept
Expand this to ensure that this->contains(point).
Definition: Box.cc:380
double getHeight() const noexcept
1-d interval accessors
Definition: Box.h:530
LinearTransform const inverted() const
Return the inverse transform.
A 2-d function defined by a series expansion and its coefficients.
Definition: Function2d.h:42
Workspace makeWorkspace() const
Allocate workspace that can be passed to operator() to avoid repeated memory allocations.
Definition: Function2d.h:107
Basis const & getBasis() const
Return the associated Basis2d object.
Definition: Function2d.h:101
typename Basis::Workspace Workspace
Type returned by makeWorkspace().
Definition: Function2d.h:52
A Basis2d formed from the product of a Basis1d for each of x and y, truncated at the sum of their ord...
Definition: PackedBasis2d.h:75
Reports invalid arguments.
Definition: Runtime.h:66
Reports errors in the logical structure of the program.
Definition: Runtime.h:46
T emplace_back(T... args)
T for_each(T... args)
T make_pair(T... args)
T max(T... args)
Low-level polynomials (including special polynomials) in C++.
Function2d< Basis > makeFunction2d(Basis const &basis, Eigen::VectorXd const &coefficients)
Create a Function2d of the appropriate type from a Basis2d and an Eigen object containing coefficient...
Definition: Function2d.h:155
PolynomialFunction1d simplified(ScaledPolynomialFunction1d const &f)
Calculate the standard polynomial function that is equivalent to a scaled standard polynomial functio...
PolynomialFunction2d< PackingOrder::YX > PolynomialFunction2dYX
A Function2d for standard polynomials, ordered via PackingOrder::YX.
ScaledPolynomialBasis2d< PackingOrder::YX > ScaledPolynomialBasis2dYX
A Basis2d for scaled standard polynomials, ordered via PackingOrder::YX.
Extent< double, 2 > Extent2D
Definition: Extent.h:400
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
A base class for image defects.
STL namespace.
T push_back(T... args)
T reserve(T... args)
T size(T... args)
table::Key< table::Array< double > > basis
Definition: PsfexPsf.cc:361
Grid(lsst::geom::Extent2I const &shape_, SipApproximation const &parent)
std::vector< lsst::geom::Point2D > dpix1
std::vector< lsst::geom::Point2D > dpix2
std::vector< lsst::geom::Point2D > siwc
static std::unique_ptr< Solution > fit(int order_, double svdThreshold, SipApproximation const &parent)
lsst::geom::Point2D applyInverse(lsst::geom::Point2D const &siwc, Workspace &ws) const
Solution(poly::PolynomialFunction2dYX const &a_, poly::PolynomialFunction2dYX const &b_, poly::PolynomialFunction2dYX const &ap_, poly::PolynomialFunction2dYX const &bp_)
lsst::geom::Point2D applyForward(lsst::geom::Point2D const &dpix, Workspace &ws) const
poly::PolynomialFunction2dYX::Workspace Workspace
table::Key< int > order