PolynomialTransform.cc
Go to the documentation of this file.
1 // -*- LSST-C++ -*-
2
3 /*
4  * LSST Data Management System
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
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
41  PolynomialTransform poly = other.getPoly();
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(),
48  compose(poly, geom::AffineTransform(geom::Point2D() - other.getPixelOrigin())));
49 }
50
52  PolynomialTransform poly = other.getPoly();
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;
57  return compose(geom::AffineTransform(geom::Extent2D(other.getPixelOrigin())),
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) {
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
y
int y
Definition: SpanSet.cc:49
lsst::meas::astrom::PolynomialTransform::convert
static PolynomialTransform convert(ScaledPolynomialTransform const &other)
Convert a ScaledPolynomialTransform to an equivalent PolynomialTransform.
Definition: PolynomialTransform.cc:36
lsst::meas::astrom::ScaledPolynomialTransform::linearize
geom::AffineTransform linearize(geom::Point2D const &in) const
Return an approximate affine transform at the given point.
Definition: PolynomialTransform.cc:206
SipTransform.h
lsst::meas::astrom::ScaledPolynomialTransform::swap
void swap(ScaledPolynomialTransform &other)
Definition: PolynomialTransform.cc:200
polynomialUtils.h
lsst::meas::astrom::PolynomialTransform::getOrder
int getOrder() const
Return the order of the polynomials.
Definition: PolynomialTransform.h:107
lsst::geom::LinearTransform
A 2D linear coordinate transformation.
Definition: LinearTransform.h:69
PolynomialTransform.h
AffineTransform.h
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
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.pex.config.history.format
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
lsst::meas::astrom::ScaledPolynomialTransform::convert
static ScaledPolynomialTransform convert(PolynomialTransform const &poly)
Convert a PolynomialTransform to an equivalent ScaledPolynomialTransform.
Definition: PolynomialTransform.cc:172
lsst::geom::Point2D
Point< double, 2 > Point2D
Definition: Point.h:324
lsst::afw::image::X
@ X
Definition: ImageUtils.h:36
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
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
Extent.h
ndarray
Definition: SpanSetFunctorGetters.h:303
z
double z
Definition: Match.cc:44
x
double x
Definition: ChebyshevBoundedField.cc:277
lsst::meas::astrom::PolynomialTransform::operator=
PolynomialTransform & operator=(PolynomialTransform const &other)
Copy assignment.
Definition: PolynomialTransform.cc:107
lsst.pex::exceptions::LengthError
Reports attempts to exceed implementation-defined length limits for some classes.
Definition: Runtime.h:76
lsst::meas::astrom::PolynomialTransform::operator()
geom::Point2D operator()(geom::Point2D const &in) const
Apply the transform to a point.
Definition: PolynomialTransform.cc:157
other
ItemVariant const * other
Definition: Schema.cc:56
lsst::afw::image::Y
@ Y
Definition: ImageUtils.h:36
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
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::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::swap
T swap(T... args)
lsst::meas::astrom::PolynomialTransform::PolynomialTransform
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.
Definition: PolynomialTransform.cc:74
lsst::meas::astrom::SipTransformBase::getCdMatrix
geom::LinearTransform const & getCdMatrix() const
Return the CD matrix of the transform.
Definition: SipTransform.h:56
lsst::geom::polynomials
Low-level polynomials (including special polynomials) in C++.
Definition: Basis1d.h:26
lsst::meas::astrom::ScaledPolynomialTransform::operator()
geom::Point2D operator()(geom::Point2D const &in) const
Apply the transform to a point.
Definition: PolynomialTransform.cc:210
lsst::meas::astrom::ScaledPolynomialTransform::ScaledPolynomialTransform
ScaledPolynomialTransform(PolynomialTransform const &poly, geom::AffineTransform const &inputScaling, geom::AffineTransform const &outputScalingInverse)
Construct a new ScaledPolynomialTransform from its constituents.
Definition: PolynomialTransform.cc:195
lsst::geom::Point< double, 2 >
lsst::meas::astrom::detail::computePowers
void computePowers(Eigen::VectorXd &r, double x)
Fill an array with integer powers of x, so .
Definition: polynomialUtils.cc:33
Point.h
lsst::meas::astrom::PolynomialTransform::compose
friend PolynomialTransform compose(geom::AffineTransform const &t1, PolynomialTransform const &t2)
Return a PolynomialTransform that is equivalent to the composition t1(t2())
Definition: PolynomialTransform.cc:214
lsst::meas::astrom::detail::BinomialMatrix
A class that computes binomial coefficients up to a certain power.
Definition: polynomialUtils.h:84
lsst::meas::astrom::SipReverseTransform
A transform that maps intermediate world coordinates to pixel coordinates according to the SIP conven...
Definition: SipTransform.h:246
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.pipe.drivers.constructCalibs.getSize
def getSize(dimList)
Definition: constructCalibs.py:176
lsst::meas::astrom::PolynomialTransform::swap
void swap(PolynomialTransform &other)
Lightweight swap.
Definition: PolynomialTransform.cc:122
m
int m
Definition: SpanSet.cc:49
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