LSSTApplications  20.0.0
LSSTDataManagementBasePackage
PhotometryTransform.cc
Go to the documentation of this file.
1 // -*- LSST-C++ -*-
2 /*
3  * This file is part of jointcal.
4  *
5  * Developed for the LSST Data Management System.
6  * This product includes software developed by the LSST Project
7  * (https://www.lsst.org).
8  * See the COPYRIGHT file at the top-level directory of this distribution
9  * for details of code ownership.
10  *
11  * This program is free software: you can redistribute it and/or modify
12  * it under the terms of the GNU General Public License as published by
13  * the Free Software Foundation, either version 3 of the License, or
14  * (at your option) any later version.
15  *
16  * This program is distributed in the hope that it will be useful,
17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19  * GNU General Public License for more details.
20  *
21  * You should have received a copy of the GNU General Public License
22  * along with this program. If not, see <https://www.gnu.org/licenses/>.
23  */
24 
25 #include "ndarray.h"
26 #include "Eigen/Core"
27 
29 
30 #include "lsst/jointcal/Point.h"
32 
33 namespace lsst {
34 namespace jointcal {
35 
36 // ------------------ PhotometryTransformChebyshev helpers ---------------------------------------------------
37 
38 namespace {
39 
40 // To evaluate a 1-d Chebyshev function without needing to have workspace, we use the
41 // Clenshaw algorith, which is like going through the recurrence relation in reverse.
42 // The CoeffGetter argument g is something that behaves like an array, providing access
43 // to the coefficients.
44 template <typename CoeffGetter>
45 double evaluateFunction1d(CoeffGetter g, double x, std::size_t size) {
46  double b_kp2 = 0.0, b_kp1 = 0.0;
47  for (std::size_t k = (size - 1); k > 0; --k) {
48  double b_k = g[k] + 2 * x * b_kp1 - b_kp2;
49  b_kp2 = b_kp1;
50  b_kp1 = b_k;
51  }
52  return g[0] + x * b_kp1 - b_kp2;
53 }
54 
55 // This class imitates a 1-d array, by running evaluateFunction1d on a nested dimension;
56 // this lets us reuse the logic in evaluateFunction1d for both dimensions. Essentially,
57 // we run evaluateFunction1d on a column of coefficients to evaluate T_i(x), then pass
58 // the result of that to evaluateFunction1d with the results as the "coefficients" associated
59 // with the T_j(y) functions.
60 struct RecursionArrayImitator {
61  double operator[](Eigen::Index i) const {
62  return evaluateFunction1d(coefficients[i], x, coefficients.getSize<1>());
63  }
64 
65  RecursionArrayImitator(ndarray::Array<double const, 2, 2> const &coefficients_, double x_)
66  : coefficients(coefficients_), x(x_) {}
67 
68  ndarray::Array<double const, 2, 2> coefficients;
69  double x;
70 };
71 
72 // Compute an affine transform that maps an arbitrary box to [-1,1]x[-1,1]
73 geom::AffineTransform makeChebyshevRangeTransform(geom::Box2D const &bbox) {
74  return geom::AffineTransform(
75  geom::LinearTransform::makeScaling(2.0 / bbox.getWidth(), 2.0 / bbox.getHeight()),
76  geom::Extent2D(-(2.0 * bbox.getCenterX()) / bbox.getWidth(),
77  -(2.0 * bbox.getCenterY()) / bbox.getHeight()));
78 }
79 
80 // Initialize a "unit" Chebyshev
81 ndarray::Array<double, 2, 2> _initializeChebyshev(size_t order, bool identity) {
82  ndarray::Array<double, 2, 2> coeffs = ndarray::allocate(ndarray::makeVector(order + 1, order + 1));
83  coeffs.deep() = 0.0;
84  if (identity) {
85  coeffs[0][0] = 1;
86  }
87  return coeffs;
88 }
89 } // namespace
90 
92  bool identity)
93  : _bbox(bbox),
94  _toChebyshevRange(makeChebyshevRangeTransform(bbox)),
95  _coefficients(_initializeChebyshev(order, identity)),
96  _order(order),
97  _nParameters((order + 1) * (order + 2) / 2) {}
98 
100  geom::Box2D const &bbox)
101  : _bbox(bbox),
102  _toChebyshevRange(makeChebyshevRangeTransform(bbox)),
103  _coefficients(coefficients),
104  _order(coefficients.size() - 1),
105  _nParameters((_order + 1) * (_order + 2) / 2) {}
106 
107 void PhotometryTransformChebyshev::offsetParams(Eigen::VectorXd const &delta) {
108  // NOTE: the indexing in this method and computeParameterDerivatives must be kept consistent!
109  Eigen::VectorXd::Index k = 0;
110  for (ndarray::Size j = 0; j <= _order; ++j) {
111  ndarray::Size const iMax = _order - j; // to save re-computing `i+j <= order` every inner step.
112  for (ndarray::Size i = 0; i <= iMax; ++i, ++k) {
113  _coefficients[j][i] -= delta[k];
114  }
115  }
116 }
117 
118 namespace {
119 // The integral of T_n(x) over [-1,1]:
120 // https://en.wikipedia.org/wiki/Chebyshev_polynomials#Differentiation_and_integration
121 double integrateTn(int n) {
122  if (n % 2 == 1)
123  return 0;
124  else
125  return 2.0 / (1.0 - static_cast<double>(n * n));
126 }
127 
136 Eigen::VectorXd computeIntegralTn(ndarray::Size order, double x) {
137  // We need to compute through order+2, because the integral recurrance relation uses Tn+1(x).
138  Eigen::VectorXd Tn(order + 2);
139 
140  // The nth chebyshev polynomial evaluated at x.
141  Tn[0] = 1;
142  Tn[1] = x;
143  for (ndarray::Size i = 2; i <= order + 1; ++i) {
144  Tn[i] = 2 * x * Tn[i - 1] - Tn[i - 2];
145  }
146 
147  // The integral of the nth chebyshev polynomial at x.
148  Eigen::VectorXd integralTn(order + 1);
149  integralTn[0] = x;
150  if (order > 0) { // have to evaluate this separately from the recurrence relation below.
151  integralTn[1] = 0.5 * x * x;
152  }
153  for (ndarray::Size i = 2; i <= order; ++i) {
154  // recurrence relation: integral(Tn(x)) = n*T[n+1](x)/(n^2 - 1) - x*T[n](x)/(n-1)
155  integralTn[i] = i * Tn[i + 1] / (i * i - 1) - x * Tn[i] / (i - 1);
156  }
157 
158  return integralTn;
159 }
160 
161 } // namespace
162 
163 double PhotometryTransformChebyshev::oneIntegral(double x, double y) const {
164  geom::Point2D point = _toChebyshevRange(geom::Point2D(x, y));
165 
166  auto integralTnx = computeIntegralTn(_order, point.getX());
167  auto integralTmy = computeIntegralTn(_order, point.getY());
168 
169  // NOTE: the indexing in this method and offsetParams must be kept consistent!
170  // Roll up the x and y terms
171  double result = 0;
172  for (ndarray::Size j = 0; j <= _order; ++j) {
173  ndarray::Size const iMax = _order - j; // to save re-computing `i+j <= order` every inner step.
174  for (ndarray::Size i = 0; i <= iMax; ++i) {
175  result += _coefficients[j][i] * integralTnx[i] * integralTmy[j];
176  }
177  }
178  return result;
179 }
180 
182  double result = 0;
183 
184  result += oneIntegral(bbox.getMaxX(), bbox.getMaxY());
185  result += oneIntegral(bbox.getMinX(), bbox.getMinY());
186  result -= oneIntegral(bbox.getMaxX(), bbox.getMinY());
187  result -= oneIntegral(bbox.getMinX(), bbox.getMaxY());
188 
189  // scale factor due to the change of limits in the integral
190  return result / _toChebyshevRange.getLinear().computeDeterminant();
191 }
192 
194  double result = 0;
195  double determinant = _bbox.getArea() / 4.0;
196  for (ndarray::Size j = 0; j < _coefficients.getSize<0>(); j++) {
197  for (ndarray::Size i = 0; i < _coefficients.getSize<1>(); i++) {
198  result += _coefficients[j][i] * integrateTn(i) * integrateTn(j);
199  }
200  }
201  return result * determinant;
202 }
203 
205  return integrate(bbox) / bbox.getArea();
206 }
207 
208 double PhotometryTransformChebyshev::mean() const { return integrate() / _bbox.getArea(); }
209 
211  Eigen::VectorXd parameters(_nParameters);
212  // NOTE: the indexing in this method and offsetParams must be kept consistent!
213  Eigen::VectorXd::Index k = 0;
214  for (ndarray::Size j = 0; j <= _order; ++j) {
215  ndarray::Size const iMax = _order - j; // to save re-computing `i+j <= order` every inner step.
216  for (ndarray::Size i = 0; i <= iMax; ++i, ++k) {
217  parameters[k] = _coefficients[j][i];
218  }
219  }
220 
221  return parameters;
222 }
223 
224 double PhotometryTransformChebyshev::computeChebyshev(double x, double y) const {
225  geom::Point2D p = _toChebyshevRange(geom::Point2D(x, y));
226  return evaluateFunction1d(RecursionArrayImitator(_coefficients, p.getX()), p.getY(),
227  _coefficients.getSize<0>());
228 }
229 
231  double x, double y, Eigen::Ref<Eigen::VectorXd> derivatives) const {
232  geom::Point2D p = _toChebyshevRange(geom::Point2D(x, y));
233  // Algorithm: compute all the individual components recursively (since we'll need them anyway),
234  // then combine them into the final answer vectors.
235  Eigen::VectorXd Tnx(_order + 1);
236  Eigen::VectorXd Tmy(_order + 1);
237  Tnx[0] = 1;
238  Tmy[0] = 1;
239  if (_order >= 1) {
240  Tnx[1] = p.getX();
241  Tmy[1] = p.getY();
242  }
243  for (ndarray::Size i = 2; i <= _order; ++i) {
244  Tnx[i] = 2 * p.getX() * Tnx[i - 1] - Tnx[i - 2];
245  Tmy[i] = 2 * p.getY() * Tmy[i - 1] - Tmy[i - 2];
246  }
247 
248  // NOTE: the indexing in this method and offsetParams must be kept consistent!
249  Eigen::VectorXd::Index k = 0;
250  for (ndarray::Size j = 0; j <= _order; ++j) {
251  ndarray::Size const iMax = _order - j; // to save re-computing `i+j <= order` every inner step.
252  for (ndarray::Size i = 0; i <= iMax; ++i, ++k) {
253  derivatives[k] = Tmy[j] * Tnx[i];
254  }
255  }
256 }
257 
258 } // namespace jointcal
259 } // namespace lsst
y
int y
Definition: SpanSet.cc:49
lsst::jointcal::PhotometryTransformChebyshev::PhotometryTransformChebyshev
PhotometryTransformChebyshev(size_t order, geom::Box2D const &bbox, bool identity)
Create a Chebyshev transform with terms up to order in (x*y).
Definition: PhotometryTransform.cc:91
PhotometryTransform.h
lsst::afw::math::integrateTn
double integrateTn(int n)
Definition: ChebyshevBoundedField.cc:290
lsst::geom::LinearTransform::makeScaling
static LinearTransform makeScaling(double s) noexcept
Definition: LinearTransform.h:94
lsst::geom::AffineTransform
An affine coordinate transformation consisting of a linear transformation and an offset.
Definition: AffineTransform.h:75
lsst::jointcal::PhotometryTransformChebyshev::getParameters
Eigen::VectorXd getParameters() const override
Get a copy of the parameters of this model, in the same order as offsetParams.
Definition: PhotometryTransform.cc:210
Point.h
lsst::jointcal::PhotometryTransformChebyshev::offsetParams
void offsetParams(Eigen::VectorXd const &delta) override
Offset the parameters by some (negative) amount during fitting.
Definition: PhotometryTransform.cc:107
lsst::geom::AffineTransform::getLinear
LinearTransform const & getLinear() const noexcept
Definition: AffineTransform.h:155
TrapezoidalPacker.h
lsst::jointcal::PhotometryTransformChebyshev::mean
double mean() const
Definition: PhotometryTransform.cc:208
lsst::jointcal::PhotometryTransformChebyshev::computeChebyshev
double computeChebyshev(double x, double y) const
Return the value of this polynomial at x,y.
Definition: PhotometryTransform.cc:224
lsst::jointcal
Definition: Associations.h:49
result
py::object result
Definition: _schema.cc:429
lsst
A base class for image defects.
Definition: imageAlgorithm.dox:1
coefficients
ndarray::Array< double const, 2, 2 > coefficients
Definition: PhotometryTransform.cc:68
lsst::geom::LinearTransform::computeDeterminant
double computeDeterminant() const noexcept
Return the determinant of the 2x2 matrix.
Definition: LinearTransform.cc:54
lsst::jointcal::PhotometryTransformChebyshev::integrate
double integrate() const
Definition: PhotometryTransform.cc:193
lsst::geom::Point< double, 2 >
x
double x
Definition: PhotometryTransform.cc:69
std::size_t
lsst::geom::Box2D
A floating-point coordinate rectangle geometry.
Definition: Box.h:413
lsst::geom::Box2D::getArea
double getArea() const noexcept
1-d interval accessors
Definition: Box.h:531
lsst::geom::Extent< double, 2 >
bbox
AmpInfoBoxKey bbox
Definition: Amplifier.cc:117
lsst::jointcal::PhotometryTransformChebyshev::computeChebyshevDerivatives
void computeChebyshevDerivatives(double x, double y, Eigen::Ref< Eigen::VectorXd > derivatives) const
Set the derivatives of this polynomial at x,y.
Definition: PhotometryTransform.cc:230