LSSTApplications  19.0.0-10-g4a5fae6+3,19.0.0-10-g920eed2,19.0.0-11-g48a0200+2,19.0.0-18-gfc4e62b+16,19.0.0-2-g3b2f90d+2,19.0.0-2-gd671419+6,19.0.0-20-g5a5a17ab+14,19.0.0-21-g2644856+17,19.0.0-24-g0913cb1,19.0.0-24-g878c510+4,19.0.0-25-g6c8df7140+1,19.0.0-25-gb330496+4,19.0.0-3-g2b32d65+6,19.0.0-3-g8227491+15,19.0.0-3-g9c54d0d+15,19.0.0-3-gca68e65+11,19.0.0-3-gcfc5f51+6,19.0.0-3-ge110943+14,19.0.0-3-ge74d124,19.0.0-30-g9c3fd16+5,19.0.0-4-g06f5963+6,19.0.0-4-g10df615,19.0.0-4-g3d16501+17,19.0.0-4-g4a9c019+6,19.0.0-4-g5a8b323,19.0.0-4-g66397f0+1,19.0.0-4-g8557e14,19.0.0-4-g8964aba+16,19.0.0-4-ge404a01+15,19.0.0-5-g40f3a5a,19.0.0-5-g4db63b3,19.0.0-5-gb9eeb60,19.0.0-5-gfb03ce7+16,19.0.0-6-gbaebbfb+15,19.0.0-61-gec4c6e08+5,19.0.0-7-g039c0b5+15,19.0.0-7-gbea9075+4,19.0.0-7-gc567de5+16,19.0.0-72-g37abf38+2,19.0.0-9-g463f923+15,v20.0.0.rc1
LSSTDataManagementBasePackage
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
std::shared_ptr
STL class.
ast::Mapping::simplified
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
lsst::afw::geom::Transform
Transform LSST spatial data, such as lsst::geom::Point2D and lsst::geom::SpherePoint,...
Definition: Transform.h:68
std::vector
STL class.
std::vector::size
T size(T... args)
ast::Mapping::then
SeriesMap then(Mapping const &next) const
Return a series compound mapping this(first(input)) containing shallow copies of the original.
Definition: Mapping.cc:37
lsst::afw::table._match.first
first
Definition: _match.py:76
lsst::afw
Definition: imageAlgorithm.dox:1
lsst::meas::modelfit::Scalar
double Scalar
Typedefs to be used for probability and parameter values.
Definition: common.h:44
lsst::afw::geom::Transform::applyForward
ToPoint applyForward(FromPoint const &point) const
Transform one point in the forward direction ("from" to "to")
Definition: Transform.cc:79
lsst::afw::geom::makeRadialTransform
std::shared_ptr< TransformPoint2ToPoint2 > makeRadialTransform(std::vector< double > const &coeffs)
A purely radial polynomial distortion.
Definition: transformFactory.cc:165
lsst::afw::geom::linearizeTransform
lsst::geom::AffineTransform linearizeTransform(TransformPoint2ToPoint2 const &original, lsst::geom::Point2D const &inPoint)
Approximate a Transform by its local linearization.
Definition: transformFactory.cc:132
lsst::geom::AffineTransform::getTranslation
Extent2D const & getTranslation() const noexcept
Definition: AffineTransform.h:152
Endpoint.h
lsst::geom::Point2D
Point< double, 2 > Point2D
Definition: Point.h:324
lsst::afw::geom::Point2Endpoint
An endpoint for lsst::geom::Point2D.
Definition: Endpoint.h:261
ast::PolyMap
PolyMap is a Mapping which performs a general polynomial transformation.
Definition: PolyMap.h:49
astshim.h
ast::MatrixMap
MatrixMap is a form of Mapping which performs a general linear transformation.
Definition: MatrixMap.h:42
lsst::geom::AffineTransform
An affine coordinate transformation consisting of a linear transformation and an offset.
Definition: AffineTransform.h:75
lsst::afw::geom::operator<<
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:240
exceptions.h
std::isfinite
T isfinite(T... args)
std::ostream
STL class.
transformFactory.h
ast::ShiftMap
ShiftMap is a linear Mapping which shifts each axis by a specified constant value.
Definition: ShiftMap.h:40
lsst::geom::AffineTransform::getLinear
LinearTransform const & getLinear() const noexcept
Definition: AffineTransform.h:155
ast::makeRadialMapping
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
lsst::geom::LinearTransform::getMatrix
Matrix const & getMatrix() const noexcept
Definition: LinearTransform.h:151
lsst
A base class for image defects.
Definition: imageAlgorithm.dox:1
LSST_EXCEPT
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
std::ostringstream
STL class.
lsst::geom
Definition: AffineTransform.h:36
lsst::afw::geom.python.transform.toEndpoint
toEndpoint
Definition: transform.py:60
os
std::ostream * os
Definition: Schema.cc:746
lsst::pex::exceptions::InvalidParameterError
Reports invalid arguments.
Definition: Runtime.h:66
lsst::geom::Point
A coordinate class intended to represent absolute positions.
Definition: CoordinateBase.h:39
std::vector::empty
T empty(T... args)
std::ostringstream::str
T str(T... args)
element
double element[2]
Definition: BaseTable.cc:91
lsst::afw::geom::Transform::getJacobian
Eigen::MatrixXd getJacobian(FromPoint const &x) const
The Jacobian matrix of this Transform.
Definition: Transform.cc:123
ast::TranMap
TranMap is a Mapping which combines the forward transformation of a supplied Mapping with the inverse...
Definition: TranMap.h:49
lsst::afw::geom::makeIdentityTransform
std::shared_ptr< TransformPoint2ToPoint2 > makeIdentityTransform()
Trivial Transform x → x.
Definition: transformFactory.cc:214
lsst::afw::geom::makeTransform
std::shared_ptr< TransformPoint2ToPoint2 > makeTransform(lsst::geom::AffineTransform const &affine)
Wrap an lsst::geom::AffineTransform as a Transform.
Definition: transformFactory.cc:154
ast::UnitMap
A UnitMap is a unit (null) Mapping that has no effect on the coordinates supplied to it.
Definition: UnitMap.h:44