LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
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/SVD"
27 #include "Eigen/QR"
30 
31 namespace lsst { namespace afw { namespace geom {
32 
33 namespace poly = lsst::geom::polynomials;
34 
35 namespace {
36 
38  int order,
39  lsst::geom::Box2D const & box,
40  double svdThreshold,
43 ) {
44  // The scaled polynomial basis evaluates polynomials after mapping the
45  // input coordinates from the given box to [-1, 1]x[-1, 1] (for numerical
46  // stability).
47  auto basis = poly::ScaledPolynomialBasis2dYX(order, box);
48  auto workspace = basis.makeWorkspace();
49  Eigen::MatrixXd matrix = Eigen::MatrixXd::Zero(input.size(), basis.size());
50  Eigen::VectorXd xRhs(input.size());
51  Eigen::VectorXd yRhs(input.size());
52  for (int i = 0; i < matrix.rows(); ++i) {
53  basis.fill(input[i], matrix.row(i), workspace);
54  auto rhs = output[i] - input[i];
55  xRhs[i] = rhs.getX();
56  yRhs[i] = rhs.getY();
57  }
58  // Since we're not trying to null the zeroth- and first-order terms, the
59  // solution is just linear least squares, and we can do that with SVD.
60  auto decomp = matrix.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV);
61  if (svdThreshold >= 0) {
62  decomp.setThreshold(svdThreshold);
63  }
64  auto scaledX = makeFunction2d(basis, decomp.solve(xRhs));
65  auto scaledY = makeFunction2d(basis, decomp.solve(yRhs));
66  // On return, we simplify the polynomials by moving the remapping transform
67  // into the coefficients themselves.
68  return std::make_pair(simplified(scaledX), simplified(scaledY));
69 }
70 
71 // Return a vector of points on a grid, covering the given bounding box.
73  lsst::geom::Extent2I const & shape) {
74  if (shape.getX() <= 0 || shape.getY() <= 0) {
75  throw LSST_EXCEPT(
76  pex::exceptions::InvalidParameterError,
77  "Grid shape must be positive."
78  );
79  }
81  points.reserve(shape.getX()*shape.getY());
82  double const dx = bbox.getWidth()/shape.getX();
83  double const dy = bbox.getHeight()/shape.getY();
84  for (int iy = 0; iy < shape.getY(); ++iy) {
85  double const y = bbox.getMinY() + iy*dy;
86  for (int ix = 0; ix < shape.getX(); ++ix) {
87  points.emplace_back(bbox.getMinX() + ix*dx, y);
88  }
89  }
90  return points;
91 }
92 
93 // Make a polynomial object (with packed coefficients) from a square coefficients matrix.
94 poly::PolynomialFunction2dYX makePolynomialFromCoeffMatrix(ndarray::Array<double const, 2> const & coeffs) {
95  LSST_THROW_IF_NE(coeffs.getSize<0>(), coeffs.getSize<1>(), pex::exceptions::InvalidParameterError,
96  "Coefficient matrix must be square (%d != %d).");
97  poly::PolynomialBasis2dYX basis(coeffs.getSize<0>() - 1);
98  Eigen::VectorXd packed(basis.size());
99  for (auto const & i : basis.getIndices()) {
100  packed[i.flat] = coeffs[i.nx][i.ny];
101  }
102  return poly::PolynomialFunction2dYX(basis, packed);
103 }
104 
105 } // anonymous
106 
107 // Private implementation object for SipApproximation that manages the grid of points on which
108 // we evaluate the exact transform.
110 
111  // Set up the grid.
112  Grid(lsst::geom::Extent2I const & shape_, SipApproximation const & parent);
113 
114  lsst::geom::Extent2I const shape; // number of grid points in each dimension
115  std::vector<lsst::geom::Point2D> dpix1; // [pixel coords] - CRPIX
116  std::vector<lsst::geom::Point2D> siwc; // CD^{-1}([intermediate world coords])
117  std::vector<lsst::geom::Point2D> dpix2; // round-tripped version of dpix1 if useInverse, or exactly dpix1
118 };
119 
120 // Private implementation object for SipApproximation that manages the solution
122 
123  static std::unique_ptr<Solution> fit(int order_, double svdThreshold, SipApproximation const & parent);
124 
126  poly::PolynomialFunction2dYX const & b_,
127  poly::PolynomialFunction2dYX const & ap_,
128  poly::PolynomialFunction2dYX const & bp_) :
129  a(a_), b(b_), ap(ap_), bp(bp_)
130  {
131  LSST_THROW_IF_NE(a.getBasis().getOrder(), b.getBasis().getOrder(),
133  "A and B polynomials must have the same order (%d != %d).");
134  LSST_THROW_IF_NE(a.getBasis().getOrder(), ap.getBasis().getOrder(),
136  "A and AP polynomials must have the same order (%d != %d).");
137  LSST_THROW_IF_NE(a.getBasis().getOrder(), bp.getBasis().getOrder(),
139  "A and BP polynomials must have the same order (%d != %d).");
140  }
141 
143 
144  Workspace makeWorkspace() const { return a.makeWorkspace(); }
145 
147  return dpix + lsst::geom::Extent2D(a(dpix, ws), b(dpix, ws));
148  }
149 
151  return siwc + lsst::geom::Extent2D(ap(siwc, ws), bp(siwc, ws));
152  }
153 
158 };
159 
161  shape(shape_),
162  dpix1(makeGrid(parent._bbox, shape)),
163  siwc(parent._pixelToIwc->applyForward(dpix1))
164 {
165  // Apply the CRPIX offset to make pix1 into dpix1 (in-place)
166  std::for_each(dpix1.begin(), dpix1.end(), [&parent](lsst::geom::Point2D & p){ p -= parent._crpix; });
167 
168  if (parent._useInverse) {
169  // Set from the given inverse of the given pixels-to-iwc transform
170  // Note that at this point, siwc is still just iwc, because the scaling by cdInv is later.
171  dpix2 = parent._pixelToIwc->applyInverse(siwc);
172  // Apply the CRPIX offset to make pix1 into dpix2 (in-place)
173  std::for_each(dpix2.begin(), dpix2.end(), [&parent](lsst::geom::Point2D & p){ p -= parent._crpix; });
174  } else {
175  // Just make dpix2 = dpix1, and hence fit to the true inverse of pixels-to-iwc.
176  dpix2 = dpix1;
177  }
178 
179  // Apply the CD^{-1} transform to siwc
180  std::for_each(siwc.begin(), siwc.end(), [&parent](lsst::geom::Point2D & p){ p = parent._cdInv(p); });
181 }
182 
184  int order,
185  double svdThreshold,
186  SipApproximation const & parent
187 ) {
189  if (basis.size() > parent._grid->dpix1.size()) {
190  throw LSST_EXCEPT(
192  (boost::format("Number of parameters (%d) is larger than number of data points (%d)")
193  % (2*basis.size()) % (2*parent._grid->dpix1.size())).str()
194  );
195  }
196 
197  lsst::geom::Box2D boxFwd(parent._bbox);
198  boxFwd.shift(-parent._crpix);
199  auto fwd = fitSipOneDirection(order, boxFwd, svdThreshold, parent._grid->dpix1, parent._grid->siwc);
200 
201  lsst::geom::Box2D boxInv;
202  for (auto const & point : parent._grid->siwc) {
203  boxInv.include(point);
204  }
205  auto inv = fitSipOneDirection(order, boxInv, svdThreshold, parent._grid->siwc, parent._grid->dpix2);
206 
207  return std::make_unique<Solution>(fwd.first, fwd.second, inv.first, inv.second);
208 }
209 
212  lsst::geom::Point2D const & crpix,
213  Eigen::Matrix2d const & cd,
214  lsst::geom::Box2D const & bbox,
215  lsst::geom::Extent2I const & gridShape,
216  int order,
217  bool useInverse,
218  double svdThreshold
219 ) :
220  _useInverse(useInverse),
221  _pixelToIwc(std::move(pixelToIwc)),
222  _bbox(bbox),
223  _crpix(crpix),
224  _cdInv(lsst::geom::LinearTransform(cd).inverted()),
225  _grid(new Grid(gridShape, *this)),
226  _solution(Solution::fit(order, svdThreshold, *this))
227 {}
228 
231  lsst::geom::Point2D const & crpix,
232  Eigen::Matrix2d const & cd,
233  lsst::geom::Box2D const & bbox,
234  lsst::geom::Extent2I const & gridShape,
235  ndarray::Array<double const, 2> const & a,
236  ndarray::Array<double const, 2> const & b,
237  ndarray::Array<double const, 2> const & ap,
238  ndarray::Array<double const, 2> const & bp,
239  bool useInverse
240 ) :
241  _useInverse(useInverse),
242  _pixelToIwc(std::move(pixelToIwc)),
243  _bbox(bbox),
244  _crpix(crpix),
245  _cdInv(lsst::geom::LinearTransform(cd).inverted()),
246  _grid(new Grid(gridShape, *this)),
247  _solution(
248  new Solution(
249  makePolynomialFromCoeffMatrix(a),
250  makePolynomialFromCoeffMatrix(b),
251  makePolynomialFromCoeffMatrix(ap),
252  makePolynomialFromCoeffMatrix(bp)
253  )
254  )
255 {}
256 
258 
259 int SipApproximation::getOrder() const noexcept {
260  return _solution->a.getBasis().getOrder();
261 }
262 
263 double SipApproximation::getA(int p, int q) const {
264  return _solution->a[_solution->a.getBasis().index(p, q)];
265 }
266 
267 double SipApproximation::getB(int p, int q) const {
268  return _solution->b[_solution->b.getBasis().index(p, q)];
269 }
270 
271 double SipApproximation::getAP(int p, int q) const {
272  return _solution->ap[_solution->ap.getBasis().index(p, q)];
273 }
274 
275 double SipApproximation::getBP(int p, int q) const {
276  return _solution->bp[_solution->bp.getBasis().index(p, q)];
277 }
278 
279 namespace {
280 
281 template <typename F>
282 Eigen::MatrixXd makeCoefficientMatrix(std::size_t order, F getter) {
283  Eigen::MatrixXd result = Eigen::MatrixXd::Zero(order + 1, order + 1);
284  for (std::size_t p = 0; p <= order; ++p) {
285  for (std::size_t q = 0; q <= order - p; ++q) {
286  result(p, q) = getter(p, q);
287  }
288  }
289  return result;
290 }
291 
292 } // anonymous
293 
294 Eigen::MatrixXd SipApproximation::getA() const noexcept {
295  return makeCoefficientMatrix(
296  getOrder(),
297  [this](int p, int q) { return getA(p, q); }
298  );
299 }
300 
301 Eigen::MatrixXd SipApproximation::getB() const noexcept {
302  return makeCoefficientMatrix(
303  getOrder(),
304  [this](int p, int q) { return getB(p, q); }
305  );
306 }
307 
308 Eigen::MatrixXd SipApproximation::getAP() const noexcept {
309  return makeCoefficientMatrix(
310  getOrder(),
311  [this](int p, int q) { return getAP(p, q); }
312  );
313 }
314 
315 Eigen::MatrixXd SipApproximation::getBP() const noexcept {
316  return makeCoefficientMatrix(
317  getOrder(),
318  [this](int p, int q) { return getBP(p, q); }
319  );
320 }
321 
323  auto cd = _cdInv.inverted();
324  auto ws = _solution->makeWorkspace();
325  return cd(_solution->applyForward(pix - _crpix, ws));
326 }
327 
329  std::vector<lsst::geom::Point2D> const & pix) const {
330  auto ws = _solution->makeWorkspace();
332  iwc.reserve(pix.size());
333  auto cd = _cdInv.inverted();
334  for (auto const & point : pix) {
335  iwc.push_back(cd(_solution->applyForward(point - _crpix, ws)));
336  }
337  return iwc;
338 }
339 
341  auto ws = _solution->makeWorkspace();
342  return _solution->applyInverse(_cdInv(iwc), ws) + _crpix;
343 }
344 
346  std::vector<lsst::geom::Point2D> const & iwc) const {
347  auto ws = _solution->makeWorkspace();
349  pix.reserve(iwc.size());
350  for (auto const & point : iwc) {
351  pix.push_back(_solution->applyInverse(_cdInv(point), ws) + _crpix);
352  }
353  return pix;
354 }
355 
357  return lsst::geom::Extent2D(_bbox.getWidth()/_grid->shape.getX(),
358  _bbox.getHeight()/_grid->shape.getY());
359 }
360 
362  return _grid->shape;
363 }
364 
366  _grid = std::make_unique<Grid>(shape, *this);
367 }
368 
370  // We shrink the grid spacing by the given factor, which is not the same
371  // as increasing the number of grid points by that factor, because there
372  // is one more grid point that step in each dimension.
373  lsst::geom::Extent2I unit(1);
374  updateGrid((_grid->shape - unit)*f + unit);
375 }
376 
377 void SipApproximation::fit(int order, double svdThreshold) {
378  _solution = Solution::fit(order, svdThreshold, *this);
379 }
380 
382  std::pair<double, double> maxDiff(0.0, 0.0);
383  auto ws = _solution->makeWorkspace();
384  for (std::size_t i = 0; i < _grid->dpix1.size(); ++i) {
385  auto siwc2 = _solution->applyForward(_grid->dpix1[i], ws);
386  auto dpix2 = _solution->applyInverse(_grid->siwc[i], ws);
387  maxDiff.first = std::max(maxDiff.first, (_grid->siwc[i] - siwc2).computeNorm());
388  maxDiff.second = std::max(maxDiff.second, (_grid->dpix2[i] - dpix2).computeNorm());
389  }
390  return maxDiff;
391 }
392 
393 }}} // namespace lsst::afw::geom
py::object result
Definition: _schema.cc:430
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:131
table::Key< table::Array< double > > cd
Definition: OldWcs.cc:132
int y
Definition: SpanSet.cc:49
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++.
Definition: Basis1d.h:26
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.
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
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