LSSTApplications  18.0.0+106,18.0.0+50,19.0.0,19.0.0+1,19.0.0+10,19.0.0+11,19.0.0+13,19.0.0+17,19.0.0+2,19.0.0-1-g20d9b18+6,19.0.0-1-g425ff20,19.0.0-1-g5549ca4,19.0.0-1-g580fafe+6,19.0.0-1-g6fe20d0+1,19.0.0-1-g7011481+9,19.0.0-1-g8c57eb9+6,19.0.0-1-gb5175dc+11,19.0.0-1-gdc0e4a7+9,19.0.0-1-ge272bc4+6,19.0.0-1-ge3aa853,19.0.0-10-g448f008b,19.0.0-12-g6990b2c,19.0.0-2-g0d9f9cd+11,19.0.0-2-g3d9e4fb2+11,19.0.0-2-g5037de4,19.0.0-2-gb96a1c4+3,19.0.0-2-gd955cfd+15,19.0.0-3-g2d13df8,19.0.0-3-g6f3c7dc,19.0.0-4-g725f80e+11,19.0.0-4-ga671dab3b+1,19.0.0-4-gad373c5+3,19.0.0-5-ga2acb9c+2,19.0.0-5-gfe96e6c+2,w.2020.01
LSSTDataManagementBasePackage
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,
76  image::MaskedImage<PixelT> const& im, lsst::geom::Box2I const& bbox,
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),
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()) {
178  "No valid points to fit. Variance is likely zero. Try weighting=False");
179  } else {
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 {
219  if (orderX < 0) orderX = Approximate<PixelT>::_ctrl.getOrderX();
220  if (orderY < 0) orderY = Approximate<PixelT>::_ctrl.getOrderY();
221 
223  orderY == Approximate<PixelT>::_ctrl.getOrderY())
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 
260  std::shared_ptr<typename MImageT::Image> im = mi->getImage();
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
AmpInfoBoxKey bbox
Definition: Amplifier.cc:117
2-dimensional weighted sum of Chebyshev polynomials of the first kind.
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
uint64_t * ptr
Definition: RangeSet.cc:88
bool getWeighting() const
Return whether inverse variance weighting will be used in calculation.
Definition: Approximate.h:78
Use a 2-D Chebyshev polynomial.
Definition: Approximate.h:53
Low-level polynomials (including special polynomials) in C++.
Definition: Basis1d.h:26
int getOrderX() const
Return the order of approximation to use in the x-direction.
Definition: Approximate.h:70
table::Key< int > b
int y
Definition: SpanSet.cc:49
Style getStyle() const
Return the Style.
Definition: Approximate.h:66
table::Key< int > a
ImageT val
Definition: CR.cc:146
T end(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
Approximate values for a MaskedImage.
Definition: Approximate.h:39
T push_back(T... args)
Control how to make an approximation.
Definition: Approximate.h:48
int getHeight() const
Return the number of rows in the image.
Definition: MaskedImage.h:1095
A base class for image defects.
x_iterator row_begin(int y) const
Return an x_iterator to the start of the image.
Definition: MaskedImage.cc:612
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:73
Style
Choose the type of approximation to use.
Definition: Approximate.h:51
T isfinite(T... args)
double x
x_iterator row_end(int y) const
Return an x_iterator to the end of the image.
Definition: MaskedImage.cc:622
ApproximateControl(Style style, int orderX, int orderY=-1, bool weighting=true)
ctor
Definition: Approximate.cc:46
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
An const iterator to the MaskedImage.
Definition: MaskedImage.h:102
T begin(T... args)
Reports invalid arguments.
Definition: Runtime.h:66
int getWidth() const
Return the number of columns in the image.
Definition: MaskedImage.h:1093
T accumulate(T... args)
An integer coordinate rectangle.
Definition: Box.h:55
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:58
int end
table::Key< int > orderX
#define INSTANTIATE(FROMSYS, TOSYS)
Definition: Detector.cc:484
T reserve(T... args)
Reports errors that are due to events beyond the control of the program.
Definition: Runtime.h:104