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
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
#define INSTANTIATE(FROMSYS, TOSYS)
Definition Detector.cc:509
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
uint64_t * ptr
Definition RangeSet.cc:95
int y
Definition SpanSet.cc:48
table::Key< int > b
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:74
int getHeight() const
Return the number of rows in the image.
int getWidth() const
Return the number of columns in the image.
x_iterator row_end(int y) const
Return an x_iterator to the end of the image.
x_iterator row_begin(int y) const
Return an x_iterator to the start of the image.
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
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.
ApproximateControl const _ctrl
desired approximation algorithm
lsst::geom::Box2I const _bbox
Domain for approximation.
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.
Low-level polynomials (including special polynomials) in C++.
Extent< int, N > truncate(Extent< double, N > const &input) noexcept
Return the component-wise truncation (round towards zero).
Definition Extent.cc:100
ImageT val
Definition CR.cc:146