LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
ChebyshevBoundedField.cc
Go to the documentation of this file.
1// -*- LSST-C++ -*-
2/*
3 * LSST Data Management System
4 * Copyright 2008-2014 LSST Corporation.
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 <memory>
25
26#include "ndarray/eigen.h"
34
35namespace lsst {
36namespace afw {
37
38template std::shared_ptr<math::ChebyshevBoundedField> table::io::PersistableFacade<
39 math::ChebyshevBoundedField>::dynamicCast(std::shared_ptr<table::io::Persistable> const&);
40
41namespace math {
42
44
45// ------------------ Constructors and helpers ---------------------------------------------------------------
46
47namespace {
48
49// Compute an affine transform that maps an arbitrary box to [-1,1]x[-1,1]
50lsst::geom::AffineTransform makeChebyshevRangeTransform(lsst::geom::Box2D const bbox) {
52 lsst::geom::LinearTransform::makeScaling(2.0 / bbox.getWidth(), 2.0 / bbox.getHeight()),
53 lsst::geom::Extent2D(-(2.0 * bbox.getCenterX()) / bbox.getWidth(),
54 -(2.0 * bbox.getCenterY()) / bbox.getHeight()));
55}
56
57} // namespace
58
60 ndarray::Array<double const, 2, 2> const& coefficients)
62 _toChebyshevRange(makeChebyshevRangeTransform(lsst::geom::Box2D(bbox))),
63 _coefficients(coefficients) {}
64
66 : BoundedField(bbox), _toChebyshevRange(makeChebyshevRangeTransform(lsst::geom::Box2D(bbox))) {}
67
68ChebyshevBoundedField::ChebyshevBoundedField(ChebyshevBoundedField const&) = default;
69ChebyshevBoundedField::ChebyshevBoundedField(ChebyshevBoundedField&&) = default;
71
72// ------------------ fit() and helpers ---------------------------------------------------------------------
73
74namespace {
75
76using Control = ChebyshevBoundedField::Control;
77using Packer = detail::TrapezoidalPacker;
78
79// fill an array with 1-d Chebyshev functions of the 1st kind T(x), evaluated at the given point x
80void evaluateBasis1d(ndarray::Array<double, 1, 1> const& t, double x) {
81 int const n = t.getSize<0>();
82 if (n > 0) {
83 t[0] = 1.0;
84 }
85 if (n > 1) {
86 t[1] = x;
87 }
88 for (int i = 2; i < n; ++i) {
89 t[i] = 2.0 * x * t[i - 1] - t[i - 2];
90 }
91}
92
93// Create a matrix of 2-d Chebyshev functions evaluated at a set of positions, with
94// Chebyshev order along columns and evaluation positions along rows. We pack the
95// 2-d functions using the TrapezoidalPacker class, because we don't want any columns
96// that correspond to coefficients that should be set to zero.
97ndarray::Array<double, 2, 2> makeMatrix(ndarray::Array<double const, 1> const& x,
98 ndarray::Array<double const, 1> const& y,
99 lsst::geom::AffineTransform const& toChebyshevRange,
100 Packer const& packer, Control const& ctrl) {
101 int const nPoints = x.getSize<0>();
102 ndarray::Array<double, 1, 1> tx = ndarray::allocate(packer.nx);
103 ndarray::Array<double, 1, 1> ty = ndarray::allocate(packer.ny);
104 ndarray::Array<double, 2, 2> out = ndarray::allocate(nPoints, packer.size);
105 // Loop over x and y together, computing T_i(x) and T_j(y) arrays for each point,
106 // then packing them together.
107 for (int p = 0; p < nPoints; ++p) {
108 lsst::geom::Point2D sxy = toChebyshevRange(lsst::geom::Point2D(x[p], y[p]));
109 evaluateBasis1d(tx, sxy.getX());
110 evaluateBasis1d(ty, sxy.getY());
111 packer.pack(out[p], tx, ty); // this sets a row of out to the packed outer product of tx and ty
112 }
113 return out;
114}
115
116// Create a matrix of 2-d Chebyshev functions evaluated on a grid of positions, with
117// Chebyshev order along columns and evaluation positions along rows. We pack the
118// 2-d functions using the TrapezoidalPacker class, because we don't want any columns
119// that correspond to coefficients that should be set to zero.
120ndarray::Array<double, 2, 2> makeMatrix(lsst::geom::Box2I const& bbox,
121 lsst::geom::AffineTransform const& toChebyshevRange,
122 Packer const& packer, Control const& ctrl) {
123 // Create a 2-d array that contains T_j(x) for each x value, with x values in rows and j in columns
124 ndarray::Array<double, 2, 2> tx = ndarray::allocate(bbox.getWidth(), packer.nx);
125 for (int x = bbox.getBeginX(), p = 0; p < bbox.getWidth(); ++p, ++x) {
126 evaluateBasis1d(tx[p], toChebyshevRange[lsst::geom::AffineTransform::XX] * x +
127 toChebyshevRange[lsst::geom::AffineTransform::X]);
128 }
129
130 // Loop over y values, and at each point, compute T_i(y), then loop over x and multiply by the T_j(x)
131 // we already computed and stored above.
132 ndarray::Array<double, 2, 2> out = ndarray::allocate(bbox.getArea(), packer.size);
133 ndarray::Array<double, 2, 2>::Iterator outIter = out.begin();
134 ndarray::Array<double, 1, 1> ty = ndarray::allocate(ctrl.orderY + 1);
135 for (int y = bbox.getBeginY(), i = 0; i < bbox.getHeight(); ++i, ++y) {
136 evaluateBasis1d(ty, toChebyshevRange[lsst::geom::AffineTransform::YY] * y +
137 toChebyshevRange[lsst::geom::AffineTransform::Y]);
138 for (int j = 0; j < bbox.getWidth(); ++j, ++outIter) {
139 // this sets a row of out to the packed outer product of tx and ty
140 packer.pack(*outIter, tx[j], ty);
141 }
142 }
143 return out;
144}
145
146} // namespace
147
149 ndarray::Array<double const, 1> const& x,
150 ndarray::Array<double const, 1> const& y,
151 ndarray::Array<double const, 1> const& z,
152 Control const& ctrl) {
153 // Initialize the result object, so we can make use of the AffineTransform it builds
155 // This packer object knows how to map the 2-d Chebyshev functions onto a 1-d array,
156 // using only those that the control says should have nonzero coefficients.
157 Packer const packer(ctrl);
158 // Create a "design matrix" for the linear least squares problem (A in min||Ax-b||)
159 ndarray::Array<double, 2, 2> matrix = makeMatrix(x, y, result->_toChebyshevRange, packer, ctrl);
160 // Solve the linear least squares problem.
162 // Unpack the solution into a 2-d matrix, with zeros for values we didn't fit.
163 result->_coefficients = packer.unpack(lstsq.getSolution());
164 return result;
165}
166
168 ndarray::Array<double const, 1> const& x,
169 ndarray::Array<double const, 1> const& y,
170 ndarray::Array<double const, 1> const& z,
171 ndarray::Array<double const, 1> const& w,
172 Control const& ctrl) {
173 // Initialize the result object, so we can make use of the AffineTransform it builds
175 // This packer object knows how to map the 2-d Chebyshev functions onto a 1-d array,
176 // using only those that the control says should have nonzero coefficients.
177 Packer const packer(ctrl);
178 // Create a "design matrix" for the linear least squares problem ('A' in min||Ax-b||)
179 ndarray::Array<double, 2, 2> matrix = makeMatrix(x, y, result->_toChebyshevRange, packer, ctrl);
180 // We want to do weighted least squares, so we multiply both the data vector 'b' and the
181 // matrix 'A' by the weights.
182 ndarray::asEigenArray(matrix).colwise() *= ndarray::asEigenArray(w);
183 ndarray::Array<double, 1, 1> wz = ndarray::copy(z);
184 ndarray::asEigenArray(wz) *= ndarray::asEigenArray(w);
185 // Solve the linear least squares problem.
187 // Unpack the solution into a 2-d matrix, with zeros for values we didn't fit.
188 result->_coefficients = packer.unpack(lstsq.getSolution());
189 return result;
190}
191
192template <typename T>
194 Control const& ctrl) {
195 // Initialize the result object, so we can make use of the AffineTransform it builds
196 lsst::geom::Box2I bbox = img.getBBox(image::PARENT);
198 // This packer object knows how to map the 2-d Chebyshev functions onto a 1-d array,
199 // using only those that the control says should have nonzero coefficients.
200 Packer const packer(ctrl);
201 ndarray::Array<double, 2, 2> matrix = makeMatrix(bbox, result->_toChebyshevRange, packer, ctrl);
202 // Flatten the data image into a 1-d vector.
203 ndarray::Array<double, 2, 2> imgCopy = ndarray::allocate(img.getArray().getShape());
204 imgCopy.deep() = img.getArray();
205 ndarray::Array<double const, 1, 1> z = ndarray::flatten<1>(imgCopy);
206 // Solve the linear least squares problem.
208 // Unpack the solution into a 2-d matrix, with zeros for values we didn't fit.
209 result->_coefficients = packer.unpack(lstsq.getSolution());
210 return result;
211}
212
213// ------------------ modifier factories ---------------------------------------------------------------
214
216 if (static_cast<std::size_t>(ctrl.orderX) >= _coefficients.getSize<1>()) {
218 (boost::format("New x order (%d) exceeds old x order (%d)") % ctrl.orderX %
219 (_coefficients.getSize<1>() - 1))
220 .str());
221 }
222 if (static_cast<std::size_t>(ctrl.orderY) >= _coefficients.getSize<0>()) {
224 (boost::format("New y order (%d) exceeds old y order (%d)") % ctrl.orderY %
225 (_coefficients.getSize<0>() - 1))
226 .str());
227 }
228 ndarray::Array<double, 2, 2> coefficients = ndarray::allocate(ctrl.orderY + 1, ctrl.orderX + 1);
229 coefficients.deep() = _coefficients[ndarray::view(0, ctrl.orderY + 1)(0, ctrl.orderX + 1)];
230 if (ctrl.triangular) {
231 Packer packer(ctrl);
232 ndarray::Array<double, 1, 1> packed = ndarray::allocate(packer.size);
233 packer.pack(packed, coefficients);
234 packer.unpack(coefficients, packed);
235 }
236 return std::make_shared<ChebyshevBoundedField>(getBBox(), coefficients);
237}
238
240 return std::make_shared<ChebyshevBoundedField>(bbox, _coefficients);
241}
242
243// ------------------ evaluate() and helpers ---------------------------------------------------------------
244
245namespace {
246
247// To evaluate a 1-d Chebyshev function without needing to have workspace, we use the
248// Clenshaw algorith, which is like going through the recurrence relation in reverse.
249// The CoeffGetter argument g is something that behaves like an array, providing access
250// to the coefficients.
251template <typename CoeffGetter>
252double evaluateFunction1d(CoeffGetter g, double x, int size) {
253 double b_kp2 = 0.0, b_kp1 = 0.0;
254 for (int k = (size - 1); k > 0; --k) {
255 double b_k = g[k] + 2 * x * b_kp1 - b_kp2;
256 b_kp2 = b_kp1;
257 b_kp1 = b_k;
258 }
259 return g[0] + x * b_kp1 - b_kp2;
260}
261
262// This class imitates a 1-d array, by running evaluateFunction1d on a nested dimension;
263// this lets us reuse the logic in evaluateFunction1d for both dimensions. Essentially,
264// we run evaluateFunction1d on a column of coefficients to evaluate T_i(x), then pass
265// the result of that to evaluateFunction1d with the results as the "coefficients" associated
266// with the T_j(y) functions.
267struct RecursionArrayImitator {
268 double operator[](int i) const {
269 return evaluateFunction1d(coefficients[i], x, coefficients.getSize<1>());
270 }
271
272 RecursionArrayImitator(ndarray::Array<double const, 2, 2> const& coefficients_, double x_)
273 : coefficients(coefficients_), x(x_) {}
274
275 ndarray::Array<double const, 2, 2> coefficients;
276 double x;
277};
278
279} // namespace
280
282 lsst::geom::Point2D p = _toChebyshevRange(position);
283 return evaluateFunction1d(RecursionArrayImitator(_coefficients, p.getX()), p.getY(),
284 _coefficients.getSize<0>());
285}
286
287// The integral of T_n(x) over [-1,1]:
288// https://en.wikipedia.org/wiki/Chebyshev_polynomials#Differentiation_and_integration
289double integrateTn(int n) {
290 if (n % 2 == 1)
291 return 0;
292 else
293 return 2.0 / (1.0 - double(n * n));
294}
295
297 double result = 0;
298 double determinant = getBBox().getArea() / 4.0;
299 for (ndarray::Size j = 0; j < _coefficients.getSize<0>(); j++) {
300 for (ndarray::Size i = 0; i < _coefficients.getSize<1>(); i++) {
301 result += _coefficients[j][i] * integrateTn(i) * integrateTn(j);
302 }
303 }
304 return result * determinant;
305}
306
307double ChebyshevBoundedField::mean() const { return integrate() / getBBox().getArea(); }
308
309// ------------------ persistence ---------------------------------------------------------------------------
310
311namespace {
312
313struct PersistenceHelper {
314 table::Schema schema;
315 table::Key<int> orderX;
317 table::Key<table::Array<double> > coefficients;
318
319 PersistenceHelper(int nx, int ny)
320 : schema(),
321 orderX(schema.addField<int>("order_x", "maximum Chebyshev function order in x")),
322 bbox(table::Box2IKey::addFields(schema, "bbox", "bounding box", "pixel")),
323 coefficients(schema.addField<table::Array<double> >(
324 "coefficients", "Chebyshev function coefficients, ordered by y then x", nx * ny)) {}
325
326 PersistenceHelper(table::Schema const& s)
327 : schema(s), orderX(s["order_x"]), bbox(s["bbox"]), coefficients(s["coefficients"]) {}
328};
329
330class ChebyshevBoundedFieldFactory : public table::io::PersistableFactory {
331public:
332 explicit ChebyshevBoundedFieldFactory(std::string const& name)
333 : afw::table::io::PersistableFactory(name) {}
334
335 std::shared_ptr<table::io::Persistable> read(InputArchive const& archive,
336 CatalogVector const& catalogs) const override {
337 LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
338 LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
339 table::BaseRecord const& record = catalogs.front().front();
340 PersistenceHelper const keys(record.getSchema());
341 lsst::geom::Box2I bbox(record.get(keys.bbox));
342 std::size_t nx = record.get(keys.orderX) + 1;
343 std::size_t ny = keys.coefficients.getSize() / nx;
344 LSST_ARCHIVE_ASSERT(nx * ny == keys.coefficients.getSize());
345 ndarray::Array<double, 2, 2> coefficients = ndarray::allocate(ny, nx);
346 ndarray::flatten<1>(coefficients) = record.get(keys.coefficients);
347 return std::make_shared<ChebyshevBoundedField>(bbox, coefficients);
348 }
349};
350
351std::string getChebyshevBoundedFieldPersistenceName() { return "ChebyshevBoundedField"; }
352
353ChebyshevBoundedFieldFactory registration(getChebyshevBoundedFieldPersistenceName());
354
355} // namespace
356
358 return getChebyshevBoundedFieldPersistenceName();
359}
360
361std::string ChebyshevBoundedField::getPythonModule() const { return "lsst.afw.math"; }
362
364 PersistenceHelper const keys(_coefficients.getSize<1>(), _coefficients.getSize<0>());
365 table::BaseCatalog catalog = handle.makeCatalog(keys.schema);
366 std::shared_ptr<table::BaseRecord> record = catalog.addNew();
367 record->set(keys.orderX, _coefficients.getSize<1>() - 1);
368 record->set(keys.bbox, getBBox());
369 (*record)[keys.coefficients].deep() = ndarray::flatten<1>(_coefficients);
370 handle.saveCatalog(catalog);
371}
372
373// ------------------ operators -----------------------------------------------------------------------------
374
376 return std::make_shared<ChebyshevBoundedField>(getBBox(), ndarray::copy(getCoefficients() * scale));
377}
378
380 auto rhsCasted = dynamic_cast<ChebyshevBoundedField const*>(&rhs);
381 if (!rhsCasted) return false;
382
383 return (getBBox() == rhsCasted->getBBox()) &&
384 (_coefficients.getShape() == rhsCasted->_coefficients.getShape()) &&
385 all(equal(_coefficients, rhsCasted->_coefficients));
386}
387
390 os << "ChebyshevBoundedField (" << _coefficients.getShape() << " coefficients in y,x)";
391 return os.str();
392}
393
394// ------------------ explicit instantiation ----------------------------------------------------------------
395
396#ifndef DOXYGEN
397
398#define INSTANTIATE(T) \
399 template std::shared_ptr<ChebyshevBoundedField> ChebyshevBoundedField::fit(image::Image<T> const& image, \
400 Control const& ctrl)
401
402INSTANTIATE(float);
403INSTANTIATE(double);
404
405#endif
406} // namespace math
407} // namespace afw
408} // namespace lsst
py::object result
Definition _schema.cc:429
AmpInfoBoxKey bbox
Definition Amplifier.cc:117
table::Key< int > orderX
ndarray::Array< double const, 2, 2 > coefficients
#define INSTANTIATE(FROMSYS, TOSYS)
Definition Detector.cc:509
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
double z
Definition Match.cc:44
std::ostream * os
Definition Schema.cc:557
int y
Definition SpanSet.cc:48
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition Persistable.h:48
lsst::geom::Box2I getBBox(ImageOrigin origin=PARENT) const
Definition ImageBase.h:445
A class to represent a 2-dimensional array of pixels.
Definition Image.h:51
An abstract base class for 2-d functions defined on an integer bounding boxes.
lsst::geom::Box2I getBBox() const
Return the bounding box that defines the region where the field is valid.
A control object used when fitting ChebyshevBoundedField to data (see ChebyshevBoundedField::fit)
int computeSize() const
Return the number of nonzero coefficients in the Chebyshev function defined by this object.
bool triangular
"if true, only include terms where the sum of the x and y order " "is less than or equal to max(order...
int orderY
"maximum Chebyshev function order in y" ;
int orderX
"maximum Chebyshev function order in x" ;
A BoundedField based on 2-d Chebyshev polynomials of the first kind.
static std::shared_ptr< ChebyshevBoundedField > fit(lsst::geom::Box2I const &bbox, ndarray::Array< double const, 1 > const &x, ndarray::Array< double const, 1 > const &y, ndarray::Array< double const, 1 > const &z, Control const &ctrl)
Fit a Chebyshev approximation to non-gridded data with equal weights.
std::shared_ptr< ChebyshevBoundedField > relocate(lsst::geom::Box2I const &bbox) const
Return a new ChebyshevBoundedField with domain set to the given bounding box.
bool operator==(BoundedField const &rhs) const override
BoundedFields (of the same sublcass) are equal if their bounding boxes and parameters are equal.
ndarray::Array< double const, 2, 2 > getCoefficients() const
Return the coefficient matrix.
std::shared_ptr< BoundedField > operator*(double const scale) const override
Return a scaled BoundedField.
void write(OutputArchiveHandle &handle) const override
Write the object to one or more catalogs.
std::string getPythonModule() const override
Return the fully-qualified Python module that should be imported to guarantee that its factory is reg...
double evaluate(lsst::geom::Point2D const &position) const override
Evaluate the field at the given point.
double mean() const override
Compute the mean of this function over its bounding-box.
std::shared_ptr< ChebyshevBoundedField > truncate(Control const &ctrl) const
Return a new ChebyshevBoudedField with maximum orders set by the given control object.
std::string getPersistenceName() const override
Return the unique name used to persist this object and look up its factory.
ChebyshevBoundedField(lsst::geom::Box2I const &bbox, ndarray::Array< double const, 2, 2 > const &coefficients)
Initialize the field from its bounding box an coefficients.
double integrate() const override
Compute the integral of this function over its bounding-box.
Solver for linear least-squares problems.
static LeastSquares fromDesignMatrix(ndarray::Array< T1, 2, C1 > const &design, ndarray::Array< T2, 1, C2 > const &data, Factorization factorization=NORMAL_EIGENSYSTEM)
Initialize from the design matrix and data vector given as ndarrays.
@ NORMAL_EIGENSYSTEM
Use the normal equations with a symmetric Eigensystem decomposition.
An object passed to Persistable::write to allow it to persist itself.
void saveCatalog(BaseCatalog const &catalog)
Save a catalog in the archive.
BaseCatalog makeCatalog(Schema const &schema)
Return a new, empty catalog with the given schema.
An affine coordinate transformation consisting of a linear transformation and an offset.
A floating-point coordinate rectangle geometry.
Definition Box.h:413
An integer coordinate rectangle.
Definition Box.h:55
int getArea() const
Definition Box.h:189
static LinearTransform makeScaling(double s) noexcept
Reports attempts to exceed implementation-defined length limits for some classes.
Definition Runtime.h:76
BoxKey< lsst::geom::Box2I > Box2IKey
Definition aggregates.h:283
double w
Definition CoaddPsf.cc:69
A helper class ChebyshevBoundedField, for mapping trapezoidal matrices to 1-d arrays.
std::shared_ptr< table::io::Persistable > read(table::io::InputArchive const &archive, table::io::CatalogVector const &catalogs) const override