LSST Applications  21.0.0+c4f5df5339,21.0.0+e70536a077,21.0.0-1-ga51b5d4+7c60f8a6ea,21.0.0-10-gcf60f90+74aa0801fd,21.0.0-12-g63909ac9+643a1044a5,21.0.0-15-gedb9d5423+1041c3824f,21.0.0-2-g103fe59+a356b2badb,21.0.0-2-g1367e85+6d3f3f41db,21.0.0-2-g45278ab+e70536a077,21.0.0-2-g5242d73+6d3f3f41db,21.0.0-2-g7f82c8f+8d7c04eab9,21.0.0-2-g8f08a60+9c9a9cfcc8,21.0.0-2-ga326454+8d7c04eab9,21.0.0-2-gde069b7+bedfc5e1fb,21.0.0-2-gecfae73+6cb6558258,21.0.0-2-gfc62afb+6d3f3f41db,21.0.0-3-g21c7a62+f6e98b25aa,21.0.0-3-g357aad2+bd62456bea,21.0.0-3-g4be5c26+6d3f3f41db,21.0.0-3-g65f322c+03a4076c01,21.0.0-3-g7d9da8d+c4f5df5339,21.0.0-3-gaa929c8+c6b98066dc,21.0.0-3-gc44e71e+a26d5c1aea,21.0.0-3-ge02ed75+04b527a9d5,21.0.0-35-g64f566ff+b27e5ef93e,21.0.0-4-g591bb35+04b527a9d5,21.0.0-4-g88306b8+8773047b2e,21.0.0-4-gc004bbf+80a0b7acb7,21.0.0-4-gccdca77+a5c54364a0,21.0.0-4-ge8fba5a+ccfc328107,21.0.0-5-gdf36809+87b8d260e6,21.0.0-6-g00874e7+7eeda2b6ba,21.0.0-6-g2d4f3f3+e70536a077,21.0.0-6-g4e60332+04b527a9d5,21.0.0-6-g5ef7dad+f53629abd8,21.0.0-7-gc8ca178+b63e69433b,21.0.0-8-gfbe0b4b+c6b98066dc,w.2021.06
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 
41 namespace lsst {
42 namespace ex = pex::exceptions;
43 namespace afw {
44 namespace math {
45 
46 ApproximateControl::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 
56 namespace {
60 template <typename PixelT>
61 class ApproximateChebyshev : public Approximate<PixelT> {
62  template <typename T>
64  std::vector<double> const& yVec,
65  image::MaskedImage<T> const& im,
66  lsst::geom::Box2I const& bbox,
67  ApproximateControl const& ctrl);
68 
69 public:
70  ~ApproximateChebyshev() override;
71 
72 private:
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 
84 namespace {
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)
87 inline 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 
96 template <typename PixelT>
97 ApproximateChebyshev<PixelT>::ApproximateChebyshev(
98  std::vector<double> const& xVec,
99  std::vector<double> const& yVec,
100  image::MaskedImage<PixelT> const& im,
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 
203 template <typename PixelT>
204 ApproximateChebyshev<PixelT>::~ApproximateChebyshev() {}
205 
216 template <typename PixelT>
218 ApproximateChebyshev<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 }
254 template <typename PixelT>
256 ApproximateChebyshev<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 
278 template <typename PixelT>
280  std::vector<double> const& y,
281  image::MaskedImage<PixelT> const& im,
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 
304 INSTANTIATE(float);
305 // INSTANTIATE(int);
306 
308 } // namespace math
309 } // namespace afw
310 } // namespace lsst
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
uint64_t * ptr
Definition: RangeSet.cc:88
T accumulate(T... args)
T begin(T... args)
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:58
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:1084
int getWidth() const
Return the number of columns in the image.
Definition: MaskedImage.h:1082
x_iterator row_end(int y) const
Return an x_iterator to the end of the image.
Definition: MaskedImage.cc:615
x_iterator row_begin(int y) const
Return an x_iterator to the start of the image.
Definition: MaskedImage.cc:605
ImagePtr getImage() const
Return a (shared_ptr to) the MaskedImage's image.
Definition: MaskedImage.h:1046
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++.
Definition: Basis1d.h:26
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
A base class for image defects.
AmpInfoBoxKey bbox
Definition: Amplifier.cc:117
#define INSTANTIATE(FROMSYS, TOSYS)
Definition: Detector.cc:484
ImageT val
Definition: CR.cc:146
int y
Definition: SpanSet.cc:49
table::Key< int > b
table::Key< int > a
int end
table::Key< int > orderX
double x