LSST Applications  21.0.0-147-g0e635eb1+1acddb5be5,22.0.0+052faf71bd,22.0.0+1ea9a8b2b2,22.0.0+6312710a6c,22.0.0+729191ecac,22.0.0+7589c3a021,22.0.0+9f079a9461,22.0.1-1-g7d6de66+b8044ec9de,22.0.1-1-g87000a6+536b1ee016,22.0.1-1-g8e32f31+6312710a6c,22.0.1-10-gd060f87+016f7cdc03,22.0.1-12-g9c3108e+df145f6f68,22.0.1-16-g314fa6d+c825727ab8,22.0.1-19-g93a5c75+d23f2fb6d8,22.0.1-19-gb93eaa13+aab3ef7709,22.0.1-2-g8ef0a89+b8044ec9de,22.0.1-2-g92698f7+9f079a9461,22.0.1-2-ga9b0f51+052faf71bd,22.0.1-2-gac51dbf+052faf71bd,22.0.1-2-gb66926d+6312710a6c,22.0.1-2-gcb770ba+09e3807989,22.0.1-20-g32debb5+b8044ec9de,22.0.1-23-gc2439a9a+fb0756638e,22.0.1-3-g496fd5d+09117f784f,22.0.1-3-g59f966b+1e6ba2c031,22.0.1-3-g849a1b8+f8b568069f,22.0.1-3-gaaec9c0+c5c846a8b1,22.0.1-32-g5ddfab5d3+60ce4897b0,22.0.1-4-g037fbe1+64e601228d,22.0.1-4-g8623105+b8044ec9de,22.0.1-5-g096abc9+d18c45d440,22.0.1-5-g15c806e+57f5c03693,22.0.1-7-gba73697+57f5c03693,master-g6e05de7fdc+c1283a92b8,master-g72cdda8301+729191ecac,w.2021.39
LSST Data Management Base Package
transformFactory.cc
Go to the documentation of this file.
1 // -*- LSST-C++ -*-
2 /*
3  * LSST Data Management System
4  * See COPYRIGHT file at the top of the source tree.
5  *
6  * This product includes software developed by the
7  * LSST Project (http://www.lsst.org/).
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 LSST License Statement and
20  * the GNU General Public License along with this program. If not,
21  * see <http://www.lsstcorp.org/LegalNotices/>.
22  */
23 
24 #include <sstream>
25 #include <cmath>
26 
27 #include "astshim.h"
28 
29 #include "Eigen/Core"
30 
31 #include "lsst/afw/geom/Endpoint.h"
33 #include "lsst/pex/exceptions.h"
34 
35 #include "ndarray.h"
36 #include "ndarray/eigen.h"
37 
38 namespace lsst {
39 namespace afw {
40 namespace geom {
41 
42 namespace {
43 /*
44  * Print a vector to a stream.
45  *
46  * The exact details of the representation are unspecified and subject to
47  * change, but the following may be regarded as typical:
48  *
49  * [1.0, -3.560, 42.0]
50  *
51  * @tparam T the element type. Must support stream output.
52  */
53 template <typename T>
55  os << '[';
56  bool first = true;
57  for (T element : v) {
58  if (first) {
59  first = false;
60  } else {
61  os << ", ";
62  }
63  os << element;
64  }
65  os << ']';
66  return os;
67 }
68 
69 /*
70  * Convert a Matrix to the equivalent ndarray.
71  *
72  * @param matrix The matrix to convert.
73  * @returns an ndarray containing a copy of the data in `matrix`
74  *
75  * @note Returning a copy is safer that returning a view because ndarray cannot share
76  * management of the memory with Eigen. So as long as the matrix is small,
77  * making a copy is preferred even if a view would suffice.
78  */
79 template <typename Derived>
81  Eigen::MatrixBase<Derived> const &matrix) {
83  ndarray::allocate(ndarray::makeVector(matrix.rows(), matrix.cols()));
84  ndarray::asEigenMatrix(array) = matrix;
85  return array;
86 }
87 
88 /*
89  * Tests whether polynomial coefficients match the expected format.
90  *
91  * @param coeffs radial polynomial coefficients.
92  * @returns `true` if either `coeffs.size()` = 0, or `coeffs.size()` > 1,
93  * `coeffs[0]` = 0, and `coeffs[1]` &ne; 0. `false` otherwise.
94  */
95 bool areRadialCoefficients(std::vector<double> const &coeffs) noexcept {
96  if (coeffs.empty()) {
97  return true;
98  } else {
99  return coeffs.size() > 1 && coeffs[0] == 0.0 && coeffs[1] != 0.0;
100  }
101 }
102 
103 /*
104  * Make a one-dimensional polynomial distortion.
105  *
106  * The Mapping computes a scalar function
107  * @f[ f(x) = \sum_{i=1}^{N} \mathrm{coeffs[i]} \ x^i @f]
108  *
109  * @param coeffs radial polynomial coefficients. Must have `size` > 1,
110  * `coeffs[0]` = 0, and `coeffs[1]` &ne; 0.
111  * @returns the function represented by `coeffs`. The Mapping shall have an
112  * inverse, which may be approximate.
113  *
114  * @exceptsafe Provides basic exception safety.
115  *
116  * @warning Input to this function is not validated.
117  */
118 ast::PolyMap makeOneDDistortion(std::vector<double> const &coeffs) {
119  int const nCoeffs = coeffs.size() - 1; // ignore coeffs[0]
120  ndarray::Array<double, 2, 2> const polyCoeffs = ndarray::allocate(ndarray::makeVector(nCoeffs, 3));
121  for (size_t i = 1; i < coeffs.size(); ++i) {
122  polyCoeffs[i - 1][0] = coeffs[i];
123  polyCoeffs[i - 1][1] = 1;
124  polyCoeffs[i - 1][2] = i;
125  }
126 
127  return ast::PolyMap(polyCoeffs, 1, "IterInverse=1, TolInverse=1e-8, NIterInverse=30");
128 }
129 
130 } // namespace
131 
133  lsst::geom::Point2D const &inPoint) {
134  auto outPoint = original.applyForward(inPoint);
135  Eigen::Matrix2d jacobian = original.getJacobian(inPoint);
136  for (int i = 0; i < 2; ++i) {
137  if (!std::isfinite(outPoint[i])) {
138  std::ostringstream buffer;
139  buffer << "Transform ill-defined: " << inPoint << " -> " << outPoint;
141  }
142  }
143  if (!jacobian.allFinite()) {
144  std::ostringstream buffer;
145  buffer << "Transform not continuous at " << inPoint << ": J = " << jacobian;
147  }
148 
149  // y(x) = J (x - x0) + y0 = J x + (y0 - J x0)
150  auto offset = outPoint.asEigen() - jacobian * inPoint.asEigen();
151  return lsst::geom::AffineTransform(jacobian, offset);
152 }
153 
155  auto const offset = lsst::geom::Point2D(affine.getTranslation());
156  auto const jacobian = affine.getLinear().getMatrix();
157 
159  auto const map = ast::MatrixMap(toNdArray(jacobian))
160  .then(ast::ShiftMap(toEndpoint.dataFromPoint(offset)))
161  .simplified();
162  return std::make_shared<TransformPoint2ToPoint2>(*map);
163 }
164 
166  if (!areRadialCoefficients(coeffs)) {
167  std::ostringstream buffer;
168  buffer << "Invalid coefficient vector: " << coeffs;
170  }
171 
172  if (coeffs.empty()) {
173  return std::make_shared<TransformPoint2ToPoint2>(ast::UnitMap(2));
174  } else {
175  // distortion is a radial polynomial with center at focal plane center;
176  // the polynomial has an iterative inverse
177  std::vector<double> center = {0.0, 0.0};
178  ast::PolyMap const distortion = makeOneDDistortion(coeffs);
179  return std::make_shared<TransformPoint2ToPoint2>(*ast::makeRadialMapping(center, distortion));
180  }
181 }
182 
184  std::vector<double> const &inverseCoeffs) {
185  if (forwardCoeffs.empty() != inverseCoeffs.empty()) {
186  throw LSST_EXCEPT(
188  "makeRadialTransform must have either both empty or both non-empty coefficient vectors.");
189  }
190  if (forwardCoeffs.empty()) {
191  // no forward or inverse coefficients, so no distortion
192  return std::make_shared<TransformPoint2ToPoint2>(ast::UnitMap(2));
193  }
194 
195  if (!areRadialCoefficients(forwardCoeffs)) {
196  std::ostringstream buffer;
197  buffer << "Invalid forward coefficient vector: " << forwardCoeffs;
199  }
200  if (!areRadialCoefficients(inverseCoeffs)) {
201  std::ostringstream buffer;
202  buffer << "Invalid inverse coefficient vector: " << inverseCoeffs;
204  }
205  // distortion is a 1-d radial polynomial centered at focal plane center;
206  // the polynomial has coefficients specified for both the forward and inverse directions
207  std::vector<double> center = {0.0, 0.0};
208  ast::PolyMap const forward = makeOneDDistortion(forwardCoeffs);
209  auto inverse = makeOneDDistortion(inverseCoeffs).inverted();
210  auto distortion = ast::TranMap(forward, *inverse);
211  return std::make_shared<TransformPoint2ToPoint2>(*ast::makeRadialMapping(center, distortion));
212 }
213 
215  return std::make_shared<TransformPoint2ToPoint2>(ast::UnitMap(2));
216 }
217 
218 } // namespace geom
219 } // namespace afw
220 } // namespace lsst
double element[2]
Definition: BaseTable.cc:90
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
std::ostream * os
Definition: Schema.cc:557
SeriesMap then(Mapping const &next) const
Return a series compound mapping this(first(input)) containing shallow copies of the original.
Definition: Mapping.cc:37
std::shared_ptr< Mapping > simplified() const
Return a simplied version of the mapping (which may be a compound Mapping such as a CmpMap).
Definition: Mapping.h:248
MatrixMap is a form of Mapping which performs a general linear transformation.
Definition: MatrixMap.h:42
PolyMap is a Mapping which performs a general polynomial transformation.
Definition: PolyMap.h:49
ShiftMap is a linear Mapping which shifts each axis by a specified constant value.
Definition: ShiftMap.h:40
TranMap is a Mapping which combines the forward transformation of a supplied Mapping with the inverse...
Definition: TranMap.h:49
A UnitMap is a unit (null) Mapping that has no effect on the coordinates supplied to it.
Definition: UnitMap.h:44
An endpoint for lsst::geom::Point2D.
Definition: Endpoint.h:261
Transform LSST spatial data, such as lsst::geom::Point2D and lsst::geom::SpherePoint,...
Definition: Transform.h:68
Eigen::MatrixXd getJacobian(FromPoint const &x) const
The Jacobian matrix of this Transform.
Definition: Transform.cc:122
ToPoint applyForward(FromPoint const &point) const
Transform one point in the forward direction ("from" to "to")
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.
Matrix const & getMatrix() const noexcept
Reports invalid arguments.
Definition: Runtime.h:66
T empty(T... args)
T isfinite(T... args)
std::shared_ptr< Mapping > makeRadialMapping(std::vector< double > const &center, Mapping const &mapping1d)
Construct a radially symmetric mapping from a 1-dimensional mapping.
Definition: functional.cc:48
std::shared_ptr< TransformPoint2ToPoint2 > makeTransform(lsst::geom::AffineTransform const &affine)
Wrap an lsst::geom::AffineTransform as a Transform.
std::shared_ptr< TransformPoint2ToPoint2 > makeRadialTransform(std::vector< double > const &coeffs)
A purely radial polynomial distortion.
std::ostream & operator<<(std::ostream &os, GenericEndpoint const &endpoint)
Print "GenericEndpoint(_n_)" to the ostream where _n_ is the number of axes, e.g. "GenericAxes(4)".
Definition: Endpoint.cc:239
std::shared_ptr< TransformPoint2ToPoint2 > makeIdentityTransform()
Trivial Transform x → x.
lsst::geom::AffineTransform linearizeTransform(TransformPoint2ToPoint2 const &original, lsst::geom::Point2D const &inPoint)
Approximate a Transform by its local linearization.
Point< double, 2 > Point2D
Definition: Point.h:324
double Scalar
Typedefs to be used for probability and parameter values.
Definition: common.h:44
A base class for image defects.
T size(T... args)
T str(T... args)