LSSTApplications  17.0+124,17.0+14,17.0+73,18.0.0+37,18.0.0+80,18.0.0-4-g68ffd23+4,18.1.0-1-g0001055+12,18.1.0-1-g03d53ef+5,18.1.0-1-g1349e88+55,18.1.0-1-g2505f39+44,18.1.0-1-g5315e5e+4,18.1.0-1-g5e4b7ea+14,18.1.0-1-g7e8fceb+4,18.1.0-1-g85f8cd4+48,18.1.0-1-g8ff0b9f+4,18.1.0-1-ga2c679d+1,18.1.0-1-gd55f500+35,18.1.0-10-gb58edde+2,18.1.0-11-g0997b02+4,18.1.0-13-gfe4edf0b+12,18.1.0-14-g259bd21+21,18.1.0-19-gdb69f3f+2,18.1.0-2-g5f9922c+24,18.1.0-2-gd3b74e5+11,18.1.0-2-gfbf3545+32,18.1.0-26-g728bddb4+5,18.1.0-27-g6ff7ca9+2,18.1.0-3-g52aa583+25,18.1.0-3-g8ea57af+9,18.1.0-3-gb69f684+42,18.1.0-3-gfcaddf3+6,18.1.0-32-gd8786685a,18.1.0-4-gf3f9b77+6,18.1.0-5-g1dd662b+2,18.1.0-5-g6dbcb01+41,18.1.0-6-gae77429+3,18.1.0-7-g9d75d83+9,18.1.0-7-gae09a6d+30,18.1.0-9-gc381ef5+4,w.2019.45
LSSTDataManagementBasePackage
PolynomialTransform.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 "lsst/geom/Extent.h"
26 #include "lsst/geom/Point.h"
31 
32 namespace lsst {
33 namespace meas {
34 namespace astrom {
35 
37  return compose(scaled.getOutputScalingInverse(), compose(scaled.getPoly(), scaled.getInputScaling()));
38 }
39 
42  // Adding 1 here accounts for the extra terms outside the sum in the SIP
43  // transform definition (see SipForwardTransform docs) - note that you can
44  // fold those terms into the sum by adding 1 from the A_10 and B_01 terms.
45  poly._xCoeffs(1, 0) += 1;
46  poly._yCoeffs(0, 1) += 1;
47  return compose(other.getCdMatrix(),
49 }
50 
53  // Account for the terms outside the sum in the SIP definition (see comment
54  // earlier in the file for more explanation).
55  poly._xCoeffs(1, 0) += 1;
56  poly._yCoeffs(0, 1) += 1;
58  compose(poly, other._cdInverse));
59 }
60 
61 PolynomialTransform::PolynomialTransform(int order) : _xCoeffs(), _yCoeffs(), _u(), _v() {
62  if (order < 0) {
63  throw LSST_EXCEPT(pex::exceptions::LengthError, "PolynomialTransform order must be >= 0");
64  }
65  // Delay allocation until after error checking.
66  _xCoeffs = ndarray::allocate(order + 1, order + 1);
67  _yCoeffs = ndarray::allocate(order + 1, order + 1);
68  _xCoeffs.deep() = 0;
69  _yCoeffs.deep() = 0;
70  _u = Eigen::VectorXd(order + 1);
71  _v = Eigen::VectorXd(order + 1);
72 }
73 
74 PolynomialTransform::PolynomialTransform(ndarray::Array<double const, 2, 0> const& xCoeffs,
75  ndarray::Array<double const, 2, 0> const& yCoeffs)
76  : _xCoeffs(ndarray::copy(xCoeffs)),
77  _yCoeffs(ndarray::copy(yCoeffs)),
78  _u(_xCoeffs.getSize<0>()),
79  _v(_xCoeffs.getSize<0>()) {
80  if (xCoeffs.getShape() != yCoeffs.getShape()) {
81  throw LSST_EXCEPT(
83  (boost::format("X and Y coefficient matrices must have the same shape: "
84  " (%d,%d) != (%d,%d)") %
85  xCoeffs.getSize<0>() % xCoeffs.getSize<1>() % yCoeffs.getSize<0>() % yCoeffs.getSize<1>())
86  .str());
87  }
88  if (_xCoeffs.getSize<1>() != _xCoeffs.getSize<0>()) {
90  (boost::format("Coefficient matrices must be triangular, not trapezoidal: "
91  " %d != %d ") %
92  _xCoeffs.getSize<0>() % _xCoeffs.getSize<1>())
93  .str());
94  }
95 }
96 
98  : _xCoeffs(ndarray::copy(other.getXCoeffs())),
99  _yCoeffs(ndarray::copy(other.getYCoeffs())),
100  _u(other._u.size()),
101  _v(other._v.size()) {}
102 
103 PolynomialTransform::PolynomialTransform(PolynomialTransform&& other) : _xCoeffs(), _yCoeffs(), _u(), _v() {
104  this->swap(other);
105 }
106 
108  if (&other != this) {
109  PolynomialTransform tmp(other);
110  tmp.swap(*this);
111  }
112  return *this;
113 }
114 
116  if (&other != this) {
117  other.swap(*this);
118  }
119  return *this;
120 }
121 
123  _xCoeffs.swap(other._xCoeffs);
124  _yCoeffs.swap(other._yCoeffs);
125  _u.swap(other._u);
126  _v.swap(other._v);
127 }
128 
130  double xu = 0.0, xv = 0.0, yu = 0.0, yv = 0.0, x = 0.0, y = 0.0;
131  int const order = getOrder();
132  detail::computePowers(_u, in.getX());
133  detail::computePowers(_v, in.getY());
134  for (int p = 0; p <= order; ++p) {
135  for (int q = 0; q <= order; ++q) {
136  if (p > 0) {
137  xu += _xCoeffs(p, q) * p * _u[p - 1] * _v[q];
138  yu += _yCoeffs(p, q) * p * _u[p - 1] * _v[q];
139  }
140  if (q > 0) {
141  xv += _xCoeffs(p, q) * q * _u[p] * _v[q - 1];
142  yv += _yCoeffs(p, q) * q * _u[p] * _v[q - 1];
143  }
144  x += _xCoeffs(p, q) * _u[p] * _v[q];
145  y += _yCoeffs(p, q) * _u[p] * _v[q];
146  }
147  }
148  geom::LinearTransform linear;
149  linear.getMatrix()(0, 0) = xu;
150  linear.getMatrix()(0, 1) = xv;
151  linear.getMatrix()(1, 0) = yu;
152  linear.getMatrix()(1, 1) = yv;
153  geom::Point2D origin(x, y);
154  return geom::AffineTransform(linear, origin - linear(in));
155 }
156 
158  int const order = getOrder();
159  detail::computePowers(_u, in.getX());
160  detail::computePowers(_v, in.getY());
161  double x = 0;
162  double y = 0;
163  for (int p = 0; p <= order; ++p) {
164  for (int q = 0; q <= order; ++q) {
165  x += _xCoeffs(p, q) * _u[p] * _v[q];
166  y += _yCoeffs(p, q) * _u[p] * _v[q];
167  }
168  }
169  return geom::Point2D(x, y);
170 }
171 
174 }
175 
178  geom::AffineTransform(geom::Point2D(0, 0) - sipForward.getPixelOrigin()),
179  geom::AffineTransform(sipForward.getCdMatrix()));
180  // Account for the terms outside the sum in the SIP definition (see comment
181  // earlier in the file for more explanation).
182  result._poly._xCoeffs(1, 0) += 1;
183  result._poly._yCoeffs(0, 1) += 1;
184  return result;
185 }
186 
188  ScaledPolynomialTransform result(sipReverse.getPoly(), geom::AffineTransform(sipReverse._cdInverse),
190  result._poly._xCoeffs(1, 0) += 1;
191  result._poly._yCoeffs(0, 1) += 1;
192  return result;
193 }
194 
196  geom::AffineTransform const& inputScaling,
197  geom::AffineTransform const& outputScalingInverse)
198  : _poly(poly), _inputScaling(inputScaling), _outputScalingInverse(outputScalingInverse) {}
199 
201  _poly.swap(other._poly);
202  std::swap(_inputScaling, other._inputScaling);
203  std::swap(_outputScalingInverse, other._outputScalingInverse);
204 }
205 
207  return _outputScalingInverse * _poly.linearize(_inputScaling(in)) * _inputScaling;
208 }
209 
211  return _outputScalingInverse(_poly(_inputScaling(in)));
212 }
213 
215  typedef geom::AffineTransform AT;
217 
218  result._xCoeffs.deep() = t2._xCoeffs * t1[AT::XX] + t2._yCoeffs * t1[AT::XY];
219  result._yCoeffs.deep() = t2._xCoeffs * t1[AT::YX] + t2._yCoeffs * t1[AT::YY];
220  result._xCoeffs(0, 0) += t1[AT::X];
221  result._yCoeffs(0, 0) += t1[AT::Y];
222  return result;
223 }
224 
226  typedef geom::AffineTransform AT;
227  int const order = t1.getOrder();
228  if (order < 1) {
229  PolynomialTransform t1a(1);
230  t1a._xCoeffs(0, 0) = t1._xCoeffs(0, 0);
231  t1a._yCoeffs(0, 0) = t1._yCoeffs(0, 0);
232  return compose(t1a, t2);
233  }
234  detail::BinomialMatrix binomial(order);
235  // For each of these, (e.g.) a[n] == pow(a, n)
236  auto const t2u = detail::computePowers(t2[AT::X], order);
237  auto const t2v = detail::computePowers(t2[AT::Y], order);
238  auto const t2uu = detail::computePowers(t2[AT::XX], order);
239  auto const t2uv = detail::computePowers(t2[AT::XY], order);
240  auto const t2vu = detail::computePowers(t2[AT::YX], order);
241  auto const t2vv = detail::computePowers(t2[AT::YY], order);
243  for (int p = 0; p <= order; ++p) {
244  for (int m = 0; m <= p; ++m) {
245  for (int j = 0; j <= m; ++j) {
246  for (int q = 0; p + q <= order; ++q) {
247  for (int n = 0; n <= q; ++n) {
248  for (int k = 0; k <= n; ++k) {
249  double z = binomial(p, m) * t2u[p - m] * binomial(m, j) * t2uu[j] * t2uv[m - j] *
250  binomial(q, n) * t2v[q - n] * binomial(n, k) * t2vu[k] * t2vv[n - k];
251  result._xCoeffs(j + k, m + n - j - k) += t1._xCoeffs(p, q) * z;
252  result._yCoeffs(j + k, m + n - j - k) += t1._yCoeffs(p, q) * z;
253  } // k
254  } // n
255  } // q
256  } // j
257  } // m
258  } // p
259  return result;
260 }
261 
262 } // namespace astrom
263 } // namespace meas
264 } // namespace lsst
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
PolynomialTransform const & getPoly() const
Return the polynomial component of the transform (A,B) or (AP,BP).
Definition: SipTransform.h:61
Matrix const & getMatrix() const noexcept
geom::Point2D operator()(geom::Point2D const &in) const
Apply the transform to a point.
Low-level polynomials (including special polynomials) in C++.
Definition: Basis1d.h:26
An affine coordinate transformation consisting of a linear transformation and an offset.
T swap(T... args)
A class that computes binomial coefficients up to a certain power.
Reports attempts to exceed implementation-defined length limits for some classes. ...
Definition: Runtime.h:76
int y
Definition: SpanSet.cc:49
PolynomialTransform const & getPoly() const
Return the polynomial transform applied after the input scaling.
PolynomialTransform compose(geom::AffineTransform const &t1, PolynomialTransform const &t2)
Return a PolynomialTransform that is equivalent to the composition t1(t2())
ItemVariant const * other
Definition: Schema.cc:56
double z
Definition: Match.cc:44
geom::AffineTransform const & getInputScaling() const
Return the first affine transform applied to input points.
Point< double, 2 > Point2D
Definition: Point.h:324
geom::Point2D const & getPixelOrigin() const
Return the pixel origin (CRPIX, but zero-indexed) of the transform.
Definition: SipTransform.h:51
A transform that maps pixel coordinates to intermediate world coordinates according to the SIP conven...
Definition: SipTransform.h:136
geom::AffineTransform const & getOutputScalingInverse() const
Return the affine transform applied to points after the polynomial transform.
PolynomialTransform(ndarray::Array< double const, 2, 0 > const &xCoeffs, ndarray::Array< double const, 2, 0 > const &yCoeffs)
Construct a new transform from existing coefficient arrays.
A base class for image defects.
geom::LinearTransform const & getCdMatrix() const
Return the CD matrix of the transform.
Definition: SipTransform.h:56
void swap(PolynomialTransform &other)
Lightweight swap.
PolynomialTransform & operator=(PolynomialTransform const &other)
Copy assignment.
ScaledPolynomialTransform(PolynomialTransform const &poly, geom::AffineTransform const &inputScaling, geom::AffineTransform const &outputScalingInverse)
Construct a new ScaledPolynomialTransform from its constituents.
static PolynomialTransform convert(ScaledPolynomialTransform const &other)
Convert a ScaledPolynomialTransform to an equivalent PolynomialTransform.
double x
A 2-d coordinate transform represented by a lazy composition of an AffineTransform, a PolynomialTransform, and another AffineTransform.
geom::AffineTransform linearize(geom::Point2D const &in) const
Return an approximate affine transform at the given point.
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
A transform that maps intermediate world coordinates to pixel coordinates according to the SIP conven...
Definition: SipTransform.h:246
void swap(ScaledPolynomialTransform &other)
geom::Point2D operator()(geom::Point2D const &in) const
Apply the transform to a 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.
int m
Definition: SpanSet.cc:49
static ScaledPolynomialTransform convert(PolynomialTransform const &poly)
Convert a PolynomialTransform to an equivalent ScaledPolynomialTransform.
geom::AffineTransform linearize(geom::Point2D const &in) const
Return an approximate affine transform at the given point.
int getOrder() const
Return the order of the polynomials.
py::object result
Definition: _schema.cc:429
friend PolynomialTransform compose(geom::AffineTransform const &t1, PolynomialTransform const &t2)
Return a PolynomialTransform that is equivalent to the composition t1(t2())
A 2D linear coordinate transformation.
void computePowers(Eigen::VectorXd &r, double x)
Fill an array with integer powers of x, so .
A 2-d coordinate transform represented by a pair of standard polynomials (one for each coordinate)...