LSST Applications g0b6bd0c080+a72a5dd7e6,g1182afd7b4+2a019aa3bb,g17e5ecfddb+2b8207f7de,g1d67935e3f+06cf436103,g38293774b4+ac198e9f13,g396055baef+6a2097e274,g3b44f30a73+6611e0205b,g480783c3b1+98f8679e14,g48ccf36440+89c08d0516,g4b93dc025c+98f8679e14,g5c4744a4d9+a302e8c7f0,g613e996a0d+e1c447f2e0,g6c8d09e9e7+25247a063c,g7271f0639c+98f8679e14,g7a9cd813b8+124095ede6,g9d27549199+a302e8c7f0,ga1cf026fa3+ac198e9f13,ga32aa97882+7403ac30ac,ga786bb30fb+7a139211af,gaa63f70f4e+9994eb9896,gabf319e997+ade567573c,gba47b54d5d+94dc90c3ea,gbec6a3398f+06cf436103,gc6308e37c7+07dd123edb,gc655b1545f+ade567573c,gcc9029db3c+ab229f5caf,gd01420fc67+06cf436103,gd877ba84e5+06cf436103,gdb4cecd868+6f279b5b48,ge2d134c3d5+cc4dbb2e3f,ge448b5faa6+86d1ceac1d,gecc7e12556+98f8679e14,gf3ee170dca+25247a063c,gf4ac96e456+ade567573c,gf9f5ea5b4d+ac198e9f13,gff490e6085+8c2580be5c,w.2022.27
LSST Data Management Base Package
Approximate.cc
Go to the documentation of this file.
1// -*- LSST-C++ -*-
2
3/*
4 * LSST Data Management System
5 * Copyright 2008-2015 AURA/LSST.
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
11 * it under the terms of the GNU General Public License as published by
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 <https://www.lsstcorp.org/LegalNotices/>.
23 */
24
25/*
26 * Approximate values for a set of x,y vector<>s
27 */
28#include <limits>
29#include <algorithm>
30#include <cmath>
31#include <memory>
32#include <numeric>
33#include "Eigen/Core"
34#include "Eigen/LU"
35#include "boost/format.hpp"
36#include "lsst/pex/exceptions.h"
40
41namespace lsst {
42namespace ex = pex::exceptions;
43namespace afw {
44namespace math {
45
46ApproximateControl::ApproximateControl(Style style, int orderX, int orderY, bool weighting)
47 : _style(style), _orderX(orderX), _orderY(orderY < 0 ? orderX : orderY), _weighting(weighting) {
48 if (_orderX != _orderY) {
50 str(boost::format("X- and Y-orders must be equal (%d != %d) "
51 "due to a limitation in math::Chebyshev1Function2") %
52 _orderX % _orderY));
53 }
54}
55
56namespace {
60template <typename PixelT>
61class ApproximateChebyshev : public Approximate<PixelT> {
62 template <typename T>
64 std::vector<double> const& yVec,
65 image::MaskedImage<T> const& im,
67 ApproximateControl const& ctrl);
68
69public:
70 ~ApproximateChebyshev() override;
71
72private:
74
75 ApproximateChebyshev(std::vector<double> const& xVec, std::vector<double> const& yVec,
77 ApproximateControl const& ctrl);
79 int orderX, int orderY) const override;
81 int orderX, int orderY) const override;
82};
83
84namespace {
85// N.b. physically inlining these routines into ApproximateChebyshev
86// causes clang++ 3.1 problems; I suspect a clang bug (see http://llvm.org/bugs/show_bug.cgi?id=14162)
87inline void solveMatrix_Eigen(Eigen::MatrixXd& a, Eigen::VectorXd& b, Eigen::Map<Eigen::VectorXd>& c) {
88 Eigen::PartialPivLU<Eigen::MatrixXd> lu(a);
89 c = lu.solve(b);
90}
91} // namespace
92
96template <typename PixelT>
97ApproximateChebyshev<PixelT>::ApproximateChebyshev(
98 std::vector<double> const& xVec,
99 std::vector<double> const& yVec,
101 lsst::geom::Box2I const& bbox,
102 ApproximateControl const& ctrl
103 )
104 : Approximate<PixelT>(xVec, yVec, bbox, ctrl),
105 _poly(math::Chebyshev1Function2<double>(ctrl.getOrderX(), lsst::geom::Box2D(bbox))) {
106#if !defined(NDEBUG)
107 {
108 std::vector<double> const& coeffs = _poly.getParameters();
109 assert(std::accumulate(coeffs.begin(), coeffs.end(), 0.0) ==
110 0.0); // i.e. coeffs is initialised to 0.0
111 }
112#endif
113 int const nTerm = _poly.getNParameters(); // number of terms in polynomial
114 int const nData = im.getWidth() * im.getHeight(); // number of data points
115 /*
116 * N.b. in the comments,
117 * i runs over the 0..nTerm-1 coefficients
118 * alpha runs over the 0..nData-1 data points
119 */
120
121 /*
122 * We need the value of the polynomials evaluated at every data point, so it's more
123 * efficient to pre-calculate the values: termCoeffs[i][alpha]
124 */
125 std::vector<std::vector<double>> termCoeffs(nTerm);
126
127 for (int i = 0; i != nTerm; ++i) {
128 termCoeffs[i].reserve(nData);
129 }
130
131 for (int iy = 0; iy != im.getHeight(); ++iy) {
132 double const y = yVec[iy];
133
134 for (int ix = 0; ix != im.getWidth(); ++ix) {
135 double const x = xVec[ix];
136
137 for (int i = 0; i != nTerm; ++i) {
138 _poly.setParameter(i, 1.0);
139 termCoeffs[i].push_back(_poly(x, y));
140 _poly.setParameter(i, 0.0);
141 }
142 }
143 }
144 // We'll solve A*c = b
145 Eigen::MatrixXd A;
146 A.setZero(nTerm, nTerm); // We'll solve A*c = b
147 Eigen::VectorXd b;
148 b.setZero(nTerm);
149 /*
150 * Go through the data accumulating the values of the A and b matrix/vector
151 */
152 int alpha = 0;
153 for (int iy = 0; iy != im.getHeight(); ++iy) {
155 end = im.row_end(iy);
156 ptr != end; ++ptr, ++alpha) {
157 double const val = ptr.image();
158 double const ivar = ctrl.getWeighting() ? 1 / ptr.variance() : 1.0;
159 if (!std::isfinite(val + ivar)) {
160 continue;
161 }
162
163 for (int i = 0; i != nTerm; ++i) {
164 double const c_i = termCoeffs[i][alpha];
165 double const tmp = c_i * ivar;
166 b(i) += val * tmp;
167 A(i, i) += c_i * tmp;
168 for (int j = 0; j < i; ++j) {
169 double const c_j = termCoeffs[j][alpha];
170 A(i, j) += c_j * tmp;
171 }
172 }
173 }
174 }
175 if (A.isZero()) {
176 if (ctrl.getWeighting()) {
177 throw LSST_EXCEPT(pex::exceptions::RuntimeError,
178 "No valid points to fit. Variance is likely zero. Try weighting=False");
179 } else {
180 throw LSST_EXCEPT(pex::exceptions::RuntimeError,
181 "No valid points to fit (even though weighting is False). "
182 "Check that binSize & approxOrderX settings are appropriate for image size.");
183 }
184 }
185 // We only filled out the lower triangular part of A
186 for (int j = 0; j != nTerm; ++j) {
187 for (int i = j + 1; i != nTerm; ++i) {
188 A(j, i) = A(i, j);
189 }
190 }
191 /*
192 * OK, now all we ned do is solve that...
193 */
194 std::vector<double> cvec(nTerm);
195 Eigen::Map<Eigen::VectorXd> c(&cvec[0], nTerm); // N.b. c shares memory with cvec
196
197 solveMatrix_Eigen(A, b, c);
198
199 _poly.setParameters(cvec);
200}
201
203template <typename PixelT>
204ApproximateChebyshev<PixelT>::~ApproximateChebyshev() = default;
205
216template <typename PixelT>
218ApproximateChebyshev<PixelT>::doGetImage(int orderX, int orderY) const {
220 if (orderY < 0) orderY = Approximate<PixelT>::_ctrl.getOrderY();
221
222 math::Chebyshev1Function2<double> poly = (orderX == Approximate<PixelT>::_ctrl.getOrderX() &&
224 ? _poly
225 : _poly.truncate(orderX);
226
228
230 for (int iy = 0; iy != im->getHeight(); ++iy) {
231 double const y = iy;
232
233 int ix = 0;
234 for (typename ImageT::x_iterator ptr = im->row_begin(iy), end = im->row_end(iy); ptr != end;
235 ++ptr, ++ix) {
236 double const x = ix;
237
238 *ptr = poly(x, y);
239 }
240 }
241
242 return im;
243}
254template <typename PixelT>
256ApproximateChebyshev<PixelT>::doGetMaskedImage(int orderX, int orderY) const {
258
261
262 for (int iy = 0; iy != im->getHeight(); ++iy) {
263 double const y = iy;
264
265 int ix = 0;
266 for (typename MImageT::Image::x_iterator ptr = im->row_begin(iy), end = im->row_end(iy); ptr != end;
267 ++ptr, ++ix) {
268 double const x = ix;
269
270 *ptr = _poly(x, y);
271 }
272 }
273
274 return mi;
275}
276} // namespace
277
278template <typename PixelT>
280 std::vector<double> const& y,
282 lsst::geom::Box2I const& bbox,
283 ApproximateControl const& ctrl) {
284 switch (ctrl.getStyle()) {
287 new ApproximateChebyshev<PixelT>(x, y, im, bbox, ctrl));
288 default:
290 str(boost::format("Unknown ApproximationStyle: %d") % ctrl.getStyle()));
291 }
292}
294/*
295 * Explicit instantiations
296 *
297 */
298#define INSTANTIATE(PIXEL_T) \
299 template std::shared_ptr<Approximate<PIXEL_T>> makeApproximate( \
300 std::vector<double> const& x, std::vector<double> const& y, \
301 image::MaskedImage<PIXEL_T> const& im, lsst::geom::Box2I const& bbox, \
302 ApproximateControl const& ctrl)
303
304INSTANTIATE(float);
305// INSTANTIATE(int);
306
308} // namespace math
309} // namespace afw
310} // namespace lsst
AmpInfoBoxKey bbox
Definition: Amplifier.cc:117
int end
table::Key< int > orderX
double x
#define INSTANTIATE(FROMSYS, TOSYS)
Definition: Detector.cc:484
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
uint64_t * ptr
Definition: RangeSet.cc:88
int y
Definition: SpanSet.cc:48
table::Key< int > b
table::Key< int > a
T accumulate(T... args)
T begin(T... args)
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:51
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:73
int getHeight() const
Return the number of rows in the image.
Definition: MaskedImage.h:1056
int getWidth() const
Return the number of columns in the image.
Definition: MaskedImage.h:1054
x_iterator row_end(int y) const
Return an x_iterator to the end of the image.
Definition: MaskedImage.cc:631
x_iterator row_begin(int y) const
Return an x_iterator to the start of the image.
Definition: MaskedImage.cc:621
ImagePtr getImage() const
Return a (shared_ptr to) the MaskedImage's image.
Definition: MaskedImage.h:1018
Control how to make an approximation.
Definition: Approximate.h:48
int getOrderX() const
Return the order of approximation to use in the x-direction.
Definition: Approximate.h:70
ApproximateControl(Style style, int orderX, int orderY=-1, bool weighting=true)
ctor
Definition: Approximate.cc:46
Style getStyle() const
Return the Style.
Definition: Approximate.h:66
Style
Choose the type of approximation to use.
Definition: Approximate.h:51
@ CHEBYSHEV
Use a 2-D Chebyshev polynomial.
Definition: Approximate.h:53
int getOrderY() const
Return the order of approximation to use in the y-direction.
Definition: Approximate.h:74
Approximate values for a MaskedImage.
Definition: Approximate.h:109
ApproximateControl const _ctrl
desired approximation algorithm
Definition: Approximate.h:148
lsst::geom::Box2I const _bbox
Domain for approximation.
Definition: Approximate.h:147
2-dimensional weighted sum of Chebyshev polynomials of the first kind.
An integer coordinate rectangle.
Definition: Box.h:55
Reports invalid arguments.
Definition: Runtime.h:66
T end(T... args)
T isfinite(T... args)
std::shared_ptr< Approximate< PixelT > > makeApproximate(std::vector< double > const &x, std::vector< double > const &y, image::MaskedImage< PixelT > const &im, lsst::geom::Box2I const &bbox, ApproximateControl const &ctrl)
Construct a new Approximate object, inferring the type from the type of the given MaskedImage.
Definition: Approximate.cc:279
Low-level polynomials (including special polynomials) in C++.
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
A base class for image defects.
ImageT val
Definition: CR.cc:146