LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
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 
29 #include <limits>
30 #include <algorithm>
31 #include <numeric>
32 #include "Eigen/Core"
33 #include "Eigen/LU"
34 #include "boost/format.hpp"
35 #include "boost/shared_ptr.hpp"
36 #include "lsst/utils/ieee.h"
37 #include "lsst/pex/exceptions.h"
41 
42 namespace lsst {
43 namespace ex = pex::exceptions;
44 namespace afw {
45 namespace math {
46 
49  int orderX,
50  int orderY,
51  bool weighting
52  ) :
53  _style(style), _orderX(orderX), _orderY(orderY < 0 ? orderX : orderY), _weighting(weighting) {
54  if (_orderX != _orderY) {
55  throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
56  str(boost::format("X- and Y-orders must be equal (%d != %d) "
57  "due to a limitation in math::Chebyshev1Function2")
58  % _orderX % _orderY));
59  }
60 }
61 
62 /************************************************************************************************************/
63 
64 namespace {
68 template<typename PixelT>
69 class ApproximateChebyshev : public Approximate<PixelT> {
70  template<typename T>
71  friend PTR(Approximate<T>)
72  math::makeApproximate(std::vector<double> const &xVec, std::vector<double> const &yVec,
73  image::MaskedImage<T> const& im, geom::Box2I const& bbox,
74  ApproximateControl const& ctrl);
75 public:
76  virtual ~ApproximateChebyshev();
77 private:
78  math::Chebyshev1Function2<double> _poly;
79 
80  ApproximateChebyshev(std::vector<double> const &xVec, std::vector<double> const &yVec,
81  image::MaskedImage<PixelT> const& im, geom::Box2I const& bbox,
82  ApproximateControl const& ctrl);
83  virtual PTR(image::Image<typename Approximate<PixelT>::OutPixelT>)
84  doGetImage(int orderX, int orderY) const;
85  virtual PTR(image::MaskedImage<typename Approximate<PixelT>::OutPixelT>)
86  doGetMaskedImage(int orderX, int orderY) const;
87 };
88 
89 /************************************************************************************************************/
90 
91 namespace {
92  // N.b. physically inlining these routines into ApproximateChebyshev
93  // causes clang++ 3.1 problems; I suspect a clang bug (see http://llvm.org/bugs/show_bug.cgi?id=14162)
94  inline void
95  solveMatrix_Eigen(Eigen::MatrixXd &a,
96  Eigen::VectorXd &b,
97  Eigen::Map<Eigen::VectorXd> &c
98  ) {
99  Eigen::PartialPivLU<Eigen::MatrixXd> lu(a);
100  c = lu.solve(b);
101  }
102 }
103 
107 template<typename PixelT>
108 ApproximateChebyshev<PixelT>::ApproximateChebyshev(
109  std::vector<double> const &xVec,
110  std::vector<double> const &yVec,
111  image::MaskedImage<PixelT> const& im,
112  geom::Box2I const& bbox,
113  ApproximateControl const& ctrl
114  )
115  : Approximate<PixelT>(xVec, yVec, bbox, ctrl),
116  _poly(math::Chebyshev1Function2<double>(ctrl.getOrderX(), geom::Box2D(bbox)))
117 {
118 #if !defined(NDEBUG)
119  {
120  std::vector<double> const& coeffs = _poly.getParameters();
121  assert(std::accumulate(coeffs.begin(), coeffs.end(), 0.0) == 0.0); // i.e. coeffs is initialised to 0.0
122  }
123 #endif
124  int const nTerm = _poly.getNParameters(); // number of terms in polynomial
125  int const nData = im.getWidth()*im.getHeight(); // number of data points
126  /*
127  * N.b. in the comments,
128  * i runs over the 0..nTerm-1 coefficients
129  * alpha runs over the 0..nData-1 data points
130  */
131 
132  /*
133  * We need the value of the polynomials evaluated at every data point, so it's more
134  * efficient to pre-calculate the values: termCoeffs[i][alpha]
135  */
136  std::vector<std::vector<double> > termCoeffs(nTerm);
137 
138  for (int i = 0; i != nTerm; ++i) {
139  termCoeffs[i].reserve(nData);
140  }
141 
142  for (int iy = 0; iy != im.getHeight(); ++iy) {
143  double const y = yVec[iy];
144 
145  for (int ix = 0; ix != im.getWidth(); ++ix) {
146  double const x = xVec[ix];
147 
148  for (int i = 0; i != nTerm; ++i) {
149  _poly.setParameter(i, 1.0);
150  termCoeffs[i].push_back(_poly(x, y));
151  _poly.setParameter(i, 0.0);
152  }
153  }
154  }
155  // We'll solve A*c = b
156  Eigen::MatrixXd A; A.setZero(nTerm, nTerm); // We'll solve A*c = b
157  Eigen::VectorXd b; b.setZero(nTerm);
158  /*
159  * Go through the data accumulating the values of the A and b matrix/vector
160  */
161  int alpha = 0;
162  for (int iy = 0; iy != im.getHeight(); ++iy) {
163  for (typename image::MaskedImage<PixelT>::const_x_iterator ptr = im.row_begin(iy),
164  end = im.row_end(iy); ptr != end; ++ptr, ++alpha) {
165  double const val = ptr.image();
166  double const ivar = ctrl.getWeighting() ? 1/ptr.variance() : 1.0;
167  if (!lsst::utils::isfinite(val + ivar)) {
168  continue;
169  }
170 
171  for (int i = 0; i != nTerm; ++i) {
172  double const c_i = termCoeffs[i][alpha];
173  double const tmp = c_i*ivar;
174  b(i) += val*tmp;
175  A(i, i) += c_i*tmp;
176  for (int j = 0; j < i; ++j) {
177  double const c_j = termCoeffs[j][alpha];
178  A(i, j) += c_j*tmp;
179  }
180  }
181  }
182  }
183  if (A.isZero()) {
184  if (ctrl.getWeighting()) {
185  throw LSST_EXCEPT(pex::exceptions::RuntimeError,
186  "No valid points to fit. Variance is likely zero. Try weighting=False");
187  }
188  else {
189  throw LSST_EXCEPT(pex::exceptions::RuntimeError,
190  "No valid points to fit (even though weighting is False). "
191  "Check that binSize & approxOrderX settings are appropriate for image size.");
192  }
193  }
194  // We only filled out the lower triangular part of A
195  for (int j = 0; j != nTerm; ++j) {
196  for (int i = j + 1; i != nTerm; ++i) {
197  A(j, i) = A(i, j);
198  }
199  }
200  /*
201  * OK, now all we ned do is solve that...
202  */
203  std::vector<double> cvec(nTerm);
204  Eigen::Map<Eigen::VectorXd> c(&cvec[0], nTerm); // N.b. c shares memory with cvec
205 
206  solveMatrix_Eigen(A, b, c);
207 
208  _poly.setParameters(cvec);
209 }
210 
212 template<typename PixelT>
213 ApproximateChebyshev<PixelT>::~ApproximateChebyshev() {
214 }
215 
223 template<typename PixelT>
224 PTR(image::Image<typename Approximate<PixelT>::OutPixelT>)
225 ApproximateChebyshev<PixelT>::doGetImage(int orderX,
226  int orderY
227  ) const
228 {
229  if (orderX < 0) orderX = Approximate<PixelT>::_ctrl.getOrderX();
230  if (orderY < 0) orderY = Approximate<PixelT>::_ctrl.getOrderY();
231 
233  (orderX == Approximate<PixelT>::_ctrl.getOrderX() &&
234  orderY == Approximate<PixelT>::_ctrl.getOrderY()) ? _poly : _poly.truncate(orderX);
235 
237 
238  PTR(ImageT) im(new ImageT(Approximate<PixelT>::_bbox));
239  for (int iy = 0; iy != im->getHeight(); ++iy) {
240  double const y = iy;
241 
242  int ix = 0;
243  for (typename ImageT::x_iterator ptr = im->row_begin(iy),
244  end = im->row_end(iy); ptr != end; ++ptr, ++ix) {
245  double const x = ix;
246 
247  *ptr = poly(x, y);
248  }
249  }
250 
251  return im;
252 }
261 template<typename PixelT>
262 PTR(image::MaskedImage<typename Approximate<PixelT>::OutPixelT>)
263 ApproximateChebyshev<PixelT>::doGetMaskedImage(
264  int orderX,
265  int orderY
266  ) const
267 {
269 
270  PTR(MImageT) mi(new MImageT(Approximate<PixelT>::_bbox));
271  PTR(typename MImageT::Image) im = mi->getImage();
272 
273  for (int iy = 0; iy != im->getHeight(); ++iy) {
274  double const y = iy;
275 
276  int ix = 0;
277  for (typename MImageT::Image::x_iterator ptr = im->row_begin(iy),
278  end = im->row_end(iy); ptr != end; ++ptr, ++ix) {
279  double const x = ix;
280 
281  *ptr = _poly(x, y);
282  }
283  }
284 
285  return mi;
286 }
287 }
288 
289 /************************************************************************************************************/
293 template<typename PixelT>
294 PTR(Approximate<PixelT>)
295 makeApproximate(std::vector<double> const &x,
296  std::vector<double> const &y,
297  image::MaskedImage<PixelT> const& im,
298  geom::Box2I const& bbox,
299  ApproximateControl const& ctrl
300  )
301 {
302  switch (ctrl.getStyle()) {
303  case ApproximateControl::CHEBYSHEV:
304  return PTR(Approximate<PixelT>)(new ApproximateChebyshev<PixelT>(x, y, im, bbox, ctrl));
305  default:
306  throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
307  str(boost::format("Unknown ApproximationStyle: %d") % ctrl.getStyle()));
308  }
309 }
311 /*
312  * Explicit instantiations
313  *
314  */
315 #define INSTANTIATE(PIXEL_T) \
316  template \
317  PTR(Approximate<PIXEL_T>) makeApproximate( \
318  std::vector<double> const &x, std::vector<double> const &y, \
319  image::MaskedImage<PIXEL_T> const& im, \
320  geom::Box2I const& bbox, \
321  ApproximateControl const& ctrl)
322 
323 INSTANTIATE(float);
324 //INSTANTIATE(int);
325 
327 
328 }}} // lsst::afw::math
int y
2-dimensional weighted sum of Chebyshev polynomials of the first kind.
int getWidth() const
Return the number of columns in the image.
Definition: MaskedImage.h:901
#define PTR(...)
Definition: base.h:41
for(FootprintList::const_iterator ptr=feet->begin(), end=feet->end();ptr!=end;++ptr)
Definition: saturated.cc:82
Define a collection of useful Functions.
x_iterator row_end(int y) const
Return an x_iterator to the end of the image.
Definition: MaskedImage.cc:752
Approximate values for a MaskedImage.
Definition: Approximate.h:38
Matrix alpha
An integer coordinate rectangle.
Definition: Box.h:53
table::Key< table::Array< Kernel::Pixel > > image
Definition: FixedKernel.cc:117
ApproximateControl(Style style, int orderX, int orderY=-1, bool weighting=true)
ctor
Definition: Approximate.cc:48
math::Chebyshev1Function2< double > _poly
Definition: Approximate.cc:78
Control how to make an approximation.
Definition: Approximate.h:47
geom::Box2I const & _bbox
Definition: fits_io_mpl.h:80
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:77
Style
Choose the type of approximation to use.
Definition: Approximate.h:50
int isfinite(T t)
Definition: ieee.h:100
int x
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
#define INSTANTIATE(T)
boost::shared_ptr< Approximate< PixelT > > makeApproximate(std::vector< double > const &x, std::vector< double > const &y, image::MaskedImage< PixelT > const &im, 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:295
afw::table::Key< double > b
x_iterator row_begin(int y) const
Return an x_iterator to the start of the image.
Definition: MaskedImage.cc:742
Implementation of the Class MaskedImage.
table::Key< int > orderX
int getHeight() const
Return the number of rows in the image.
Definition: MaskedImage.h:903
ImageT val
Definition: CR.cc:154
A class to represent a 2-dimensional array of pixels.
Definition: PSF.h:43
Include files required for standard LSST Exception handling.