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
BackgroundMI.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 
32 #include <iostream>
33 #include <limits>
34 #include <vector>
35 #include <cmath>
36 #include "lsst/utils/ieee.h"
42 
43 namespace lsst {
44 namespace ex = pex::exceptions;
45 
46 namespace afw {
47 namespace math {
48 
49 namespace {
50 
51  // Given two vectors x and y, with some nans in y we want vectors x' and y' that correspond to the data
52  // without the nans basic idea is that 'x' is the values, and 'y' is the ref (where nan checking happens)
53  // cullNan(x, y, x', y')
54  void cullNan(std::vector<double> const &values, std::vector<double> const &refs,
55  std::vector<double> &culledValues, std::vector<double> &culledRefs,
56  double const defaultValue=std::numeric_limits<double>::quiet_NaN()
57  ) {
58  if (culledValues.capacity() == 0) {
59  culledValues.reserve(refs.size());
60  } else {
61  culledValues.clear();
62  }
63  if (culledRefs.capacity() == 0) {
64  culledRefs.reserve(refs.size());
65  } else {
66  culledRefs.clear();
67  }
68 
69  bool const haveDefault = !lsst::utils::isnan(defaultValue);
70 
71  for (std::vector<double>::const_iterator pVal = values.begin(), pRef = refs.begin();
72  pRef != refs.end(); ++pRef, ++pVal) {
73  if (!lsst::utils::isnan(*pRef)) {
74  culledValues.push_back(*pVal);
75  culledRefs.push_back(*pRef);
76  } else if(haveDefault) {
77  culledValues.push_back(*pVal);
78  culledRefs.push_back(defaultValue);
79  } else {
80  ; // drop a NaN
81  }
82  }
83  }
84 }
85 
108 template<typename ImageT>
109 BackgroundMI::BackgroundMI(ImageT const& img,
110  BackgroundControl const& bgCtrl
111  ) :
112  Background(img, bgCtrl), _statsImage(image::MaskedImage<InternalPixelT>())
113 {
114  // =============================================================
115  // Loop over the cells in the image, computing statistical properties
116  // of each cell in turn and using them to set _statsImage
117  int const nxSample = bgCtrl.getNxSample();
118  int const nySample = bgCtrl.getNySample();
119  _statsImage = image::MaskedImage<InternalPixelT>(nxSample, nySample);
120 
123 
124  for (int iX = 0; iX < nxSample; ++iX) {
125  for (int iY = 0; iY < nySample; ++iY) {
126  ImageT subimg = ImageT(img, geom::Box2I(geom::Point2I(_xorig[iX], _yorig[iY]),
128 
129  std::pair<double, double> res = makeStatistics(subimg, bgCtrl.getStatisticsProperty() | ERRORS,
130  *bgCtrl.getStatisticsControl()).getResult();
131  im(iX, iY) = res.first;
132  var(iX, iY) = res.second;
133  }
134  }
135 }
140  image::MaskedImage<InternalPixelT> const& statsImage
141  ) :
142  Background(imageBBox, statsImage.getWidth(), statsImage.getHeight()),
143  _statsImage(statsImage)
144 {
145 }
146 
148  UndersampleStyle const undersampleStyle,
149  int const iX, std::vector<int> const& ypix) const
150 {
152 
153  int const height = _imgBBox.getHeight();
154  _gridColumns[iX].resize(height);
155 
156  // Set _grid as a transitional measure
157  std::vector<double> _grid(_statsImage.getHeight());
158  std::copy(im.col_begin(iX), im.col_end(iX), _grid.begin());
159 
160  // remove nan from the grid values before computing columns
161  // if we do it here (ie. in _setGridColumns), it should
162  // take care of all future occurrences, so we don't need to do this elsewhere
163  std::vector<double> ycenTmp, gridTmp;
164  cullNan(_ycen, _grid, ycenTmp, gridTmp);
165 
166  PTR(Interpolate) intobj;
167  try {
168  intobj = makeInterpolate(ycenTmp, gridTmp, interpStyle);
169  } catch(pex::exceptions::OutOfRangeError &e) {
170  switch (undersampleStyle) {
171  case THROW_EXCEPTION:
172  LSST_EXCEPT_ADD(e, "setting _gridcolumns");
173  throw;
174  case REDUCE_INTERP_ORDER:
175  {
176  if (gridTmp.empty()) {
177  // Set the column to NaN. We'll deal with this properly when interpolating in x
178  ycenTmp.push_back(0);
179  gridTmp.push_back(std::numeric_limits<double>::quiet_NaN());
180 
181  intobj = makeInterpolate(ycenTmp, gridTmp, Interpolate::CONSTANT);
182  break;
183  } else {
184  return _setGridColumns(lookupMaxInterpStyle(gridTmp.size()), undersampleStyle, iX, ypix);
185  }
186  }
187  case INCREASE_NXNYSAMPLE:
188  LSST_EXCEPT_ADD(e, "The BackgroundControl UndersampleStyle INCREASE_NXNYSAMPLE is not supported.");
189  throw;
190  default:
191  LSST_EXCEPT_ADD(e, str(boost::format("The selected BackgroundControl "
192  "UndersampleStyle %d is not defined.") % undersampleStyle));
193  throw;
194  }
195  } catch(ex::Exception &e) {
196  LSST_EXCEPT_ADD(e, "setting _gridcolumns");
197  throw;
198  }
199 
200  for (int iY = 0; iY < height; ++iY) {
201  _gridColumns[iX][iY] = intobj->interpolate(ypix[iY]);
202  }
203 }
204 
208 void BackgroundMI::operator+=(float const delta
209  )
210 {
211  _statsImage += delta;
212 }
213 
217 void BackgroundMI::operator-=(float const delta
218  )
219 {
220  _statsImage -= delta;
221 }
222 
231 double BackgroundMI::getPixel(Interpolate::Style const interpStyle,
232  int const x,
233  int const y
234  ) const
235 {
236  (void)getImage<InternalPixelT>(interpStyle); // setup the interpolation
237 
238  // build an interpobj along the row y and get the x'th value
239  int const nxSample = _statsImage.getWidth();
240  std::vector<double> bg_x(nxSample);
241  for (int iX = 0; iX < nxSample; iX++) {
242  bg_x[iX] = _gridColumns[iX][y];
243  }
244  std::vector<double> xcenTmp, bgTmp;
245  cullNan(_xcen, bg_x, xcenTmp, bgTmp);
246 
247  try {
248  PTR(Interpolate) intobj = makeInterpolate(xcenTmp, bgTmp, interpStyle);
249  return static_cast<double>(intobj->interpolate(x));
250  } catch(ex::Exception &e) {
251  LSST_EXCEPT_ADD(e, "in getPixel()");
252  throw;
253  }
254 }
255 /*
256  * Worker routine for getImage
257  */
258 template<typename PixelT>
260  geom::Box2I const& bbox,
261  Interpolate::Style const interpStyle_, // Style of the interpolation
262  UndersampleStyle const undersampleStyle // Behaviour if there are too few points
263  ) const
264 {
265  if (!_imgBBox.contains(bbox)) {
266  throw LSST_EXCEPT(ex::LengthError,
267  str(boost::format("BBox (%d:%d,%d:%d) out of range (%d:%d,%d:%d)") %
268  bbox.getMinX() % bbox.getMaxX() % bbox.getMinY() % bbox.getMaxY() %
269  _imgBBox.getMinX() % _imgBBox.getMaxX() %
270  _imgBBox.getMinY() % _imgBBox.getMaxY()));
271  }
272  int const nxSample = _statsImage.getWidth();
273  int const nySample = _statsImage.getHeight();
274  Interpolate::Style interpStyle = interpStyle_; // not const -- may be modified if REDUCE_INTERP_ORDER
275 
276  /*
277  * Save the as-used interpStyle and undersampleStyle.
278  *
279  * N.b. The undersampleStyle may actually be overridden for some columns of the statsImage if they
280  * have too few good values. This doesn't prevent you reproducing the results of getImage() by
281  * calling getImage(getInterpStyle(), getUndersampleStyle())
282  */
283  _asUsedInterpStyle = interpStyle;
284  _asUsedUndersampleStyle = undersampleStyle;
285  /*
286  * Check if the requested nx,ny are sufficient for the requested interpolation style,
287  * making suitable adjustments
288  */
289  bool const isXundersampled = (nxSample < lookupMinInterpPoints(interpStyle));
290  bool const isYundersampled = (nySample < lookupMinInterpPoints(interpStyle));
291 
292  switch (undersampleStyle) {
293  case THROW_EXCEPTION:
294  if (isXundersampled && isYundersampled) {
295  throw LSST_EXCEPT(ex::InvalidParameterError,
296  "nxSample and nySample have too few points for requested interpolation style.");
297  } else if (isXundersampled) {
298  throw LSST_EXCEPT(ex::InvalidParameterError,
299  "nxSample has too few points for requested interpolation style.");
300  } else if (isYundersampled) {
301  throw LSST_EXCEPT(ex::InvalidParameterError,
302  "nySample has too few points for requested interpolation style.");
303  }
304  break;
305  case REDUCE_INTERP_ORDER:
306  if (isXundersampled || isYundersampled) {
307  Interpolate::Style const xStyle = lookupMaxInterpStyle(nxSample);
308  Interpolate::Style const yStyle = lookupMaxInterpStyle(nySample);
309  interpStyle = (nxSample < nySample) ? xStyle : yStyle;
310  _asUsedInterpStyle = interpStyle;
311  }
312  break;
313  case INCREASE_NXNYSAMPLE:
314  if (isXundersampled || isYundersampled) {
315  throw LSST_EXCEPT(ex::InvalidParameterError,
316  "The BackgroundControl UndersampleStyle INCREASE_NXNYSAMPLE is not supported.");
317  }
318  break;
319  default:
320  throw LSST_EXCEPT(ex::InvalidParameterError,
321  str(boost::format("The selected BackgroundControl "
322  "UndersampleStyle %d is not defined.") % undersampleStyle));
323  }
324 
325  // if we're approximating, don't bother with the rest of the interp-related work. Return from here.
326  if (_bctrl->getApproximateControl()->getStyle() != ApproximateControl::UNKNOWN) {
327  return doGetApproximate<PixelT>(*_bctrl->getApproximateControl(), _asUsedUndersampleStyle)->getImage();
328  }
329 
330  // =============================================================
331  // --> We'll store nxSample fully-interpolated columns to interpolate the rows over
332  // make a vector containing the y pixel coords for the column
333  int const width = _imgBBox.getWidth();
334  int const height = _imgBBox.getHeight();
335  auto const bboxOff = bbox.getMin() - _imgBBox.getMin();
336 
337  std::vector<int> ypix(height);
338  for (int iY = 0; iY < height; ++iY) {
339  ypix[iY] = iY;
340  }
341 
342  _gridColumns.resize(width);
343  for (int iX = 0; iX < nxSample; ++iX) {
344  _setGridColumns(interpStyle, undersampleStyle, iX, ypix);
345  }
346 
347  // create a shared_ptr to put the background image in and return to caller
348  // start with xy0 = 0 and set final xy0 later
351 
352  // go through row by row
353  // - interpolate on the gridcolumns that were pre-computed by the constructor
354  // - copy the values to an ImageT to return to the caller.
355  std::vector<double> xcenTmp, bgTmp;
356 
357  // N.b. There's no API to set defaultValue to other than NaN (due to issues with persistence
358  // that I don't feel like fixing; #2825). If we want to address this, this is the place
359  // to start, but note that NaN is treated specially -- it means, "Interpolate" so to allow
360  // us to put a NaN into the outputs some changes will be needed
361  double defaultValue = std::numeric_limits<double>::quiet_NaN();
362 
363  for (int y = 0, iY = bboxOff.getY(); y < bbox.getHeight(); ++y, ++iY) {
364  // build an interp object for this row
365  std::vector<double> bg_x(nxSample);
366  for (int iX = 0; iX < nxSample; iX++) {
367  bg_x[iX] = static_cast<double>(_gridColumns[iX][iY]);
368  }
369  cullNan(_xcen, bg_x, xcenTmp, bgTmp, defaultValue);
370 
371  PTR(Interpolate) intobj;
372  try {
373  intobj = makeInterpolate(xcenTmp, bgTmp, interpStyle);
374  } catch(pex::exceptions::OutOfRangeError &e) {
375  switch (undersampleStyle) {
376  case THROW_EXCEPTION:
377  LSST_EXCEPT_ADD(e, str(boost::format("Interpolating in y (iY = %d)") % iY));
378  throw;
379  case REDUCE_INTERP_ORDER:
380  {
381  if (bgTmp.empty()) {
382  xcenTmp.push_back(0);
383  bgTmp.push_back(defaultValue);
384 
385  intobj = makeInterpolate(xcenTmp, bgTmp, Interpolate::CONSTANT);
386  break;
387  } else {
388  intobj = makeInterpolate(xcenTmp, bgTmp, lookupMaxInterpStyle(bgTmp.size()));
389  }
390  }
391  break;
392  case INCREASE_NXNYSAMPLE:
393  LSST_EXCEPT_ADD(e, "The BackgroundControl UndersampleStyle INCREASE_NXNYSAMPLE is not supported.");
394  throw;
395  default:
396  LSST_EXCEPT_ADD(e, str(boost::format("The selected BackgroundControl "
397  "UndersampleStyle %d is not defined.") % undersampleStyle));
398  throw;
399  }
400  } catch(ex::Exception &e) {
401  LSST_EXCEPT_ADD(e, str(boost::format("Interpolating in y (iY = %d)") % iY));
402  throw;
403  }
404 
405  // fill the image with interpolated values
406  for (int iX = bboxOff.getX(), x = 0; x < bbox.getWidth(); ++iX, ++x) {
407  (*bg)(x, y) = static_cast<PixelT>(intobj->interpolate(iX));
408  }
409  }
410  bg->setXY0(bbox.getMin());
411 
412  return bg;
413 }
414 
415 /************************************************************************************************************/
416 
417 template<typename PixelT>
419  ApproximateControl const& actrl, /* Approximation style */
420  UndersampleStyle const undersampleStyle /* Behaviour if there are too few points */
421  ) const
422 {
423  auto const localBBox = afw::geom::Box2I(afw::geom::Point2I(0, 0), _imgBBox.getDimensions());
424  return makeApproximate(_xcen, _ycen, _statsImage, localBBox, actrl);
425 }
426 
428 /*
429  * Create the versions we need of _get{Approximate,Image} and Explicit instantiations
430  *
431  */
432 #define CREATE_BACKGROUND(m, v, TYPE) \
433  template BackgroundMI::BackgroundMI(image::Image<TYPE> const& img, \
434  BackgroundControl const& bgCtrl); \
435  template BackgroundMI::BackgroundMI(image::MaskedImage<TYPE> const& img, \
436  BackgroundControl const& bgCtrl); \
437  PTR(image::Image<TYPE>) \
438  BackgroundMI::_getImage( \
439  geom::Box2I const& bbox, \
440  Interpolate::Style const interpStyle, /* Style of the interpolation */ \
441  UndersampleStyle const undersampleStyle, /* Behaviour if there are too few points */ \
442  TYPE /* disambiguate */ \
443  ) const \
444  { \
445  return BackgroundMI::doGetImage<TYPE>(bbox, interpStyle, undersampleStyle); \
446  }
447 
448 #define CREATE_getApproximate(m, v, TYPE) \
449 PTR(Approximate<TYPE>) BackgroundMI::_getApproximate( \
450  ApproximateControl const& actrl, /* Approximation style */ \
451  UndersampleStyle const undersampleStyle, /* Behaviour if there are too few points */ \
452  TYPE /* disambiguate */ \
453  ) const \
454  { \
455  return BackgroundMI::doGetApproximate<TYPE>(actrl, undersampleStyle); \
456  }
457 
458 BOOST_PP_SEQ_FOR_EACH(CREATE_BACKGROUND, , LSST_makeBackground_getImage_types)
459 BOOST_PP_SEQ_FOR_EACH(CREATE_getApproximate, , LSST_makeBackground_getApproximate_types)
460 
462 }}} // lsst::afw::math
int y
BackgroundMI(ImageT const &img, BackgroundControl const &bgCtrl)
Constructor for BackgroundMI.
std::vector< int > _ysize
y size ...
Definition: Background.h:340
int getWidth() const
Return the number of columns in the image.
Definition: MaskedImage.h:901
ImagePtr getImage(bool const noThrow=false) const
Return a (Ptr to) the MaskedImage&#39;s image.
Definition: MaskedImage.h:869
Interpolate values for a set of x,y vector&lt;&gt;s.
Definition: Interpolate.h:37
std::vector< int > _xorig
x origin pix coords of sub images
Definition: Background.h:337
#define PTR(...)
Definition: base.h:41
boost::shared_ptr< Interpolate > makeInterpolate(std::vector< double > const &x, std::vector< double > const &y, Interpolate::Style const style=Interpolate::AKIMA_SPLINE)
Definition: Interpolate.cc:353
SelectEigenView< T >::Type copy(Eigen::EigenBase< T > const &other)
Copy an arbitrary Eigen expression into a new EigenView.
Definition: eigen.h:390
#define LSST_makeBackground_getApproximate_types
Definition: Background.h:353
boost::shared_ptr< StatisticsControl > getStatisticsControl()
Definition: Background.h:206
std::vector< double > _ycen
y center ...
Definition: Background.h:336
int isnan(T t)
Definition: ieee.h:110
Pass parameters to a Background object.
Definition: Background.h:61
Approximate values for a MaskedImage.
Definition: Approximate.h:38
A virtual base class to evaluate image background levels.
Definition: Background.h:231
An integer coordinate rectangle.
Definition: Box.h:53
table::Key< table::Array< Kernel::Pixel > > image
Definition: FixedKernel.cc:117
Include errors of requested quantities.
Definition: Statistics.h:65
y_iterator col_begin(int x) const
Definition: Image.h:334
VariancePtr getVariance(bool const noThrow=false) const
Return a (Ptr to) the MaskedImage&#39;s variance.
Definition: MaskedImage.h:890
double getPixel(Interpolate::Style const style, int const x, int const y) const
Method to retrieve the background level at a pixel coord.
int getWidth() const
Return the number of columns in the image.
Definition: Image.h:237
Control how to make an approximation.
Definition: Approximate.h:47
int getHeight() const
Definition: Box.h:155
std::vector< int > _xsize
x size of sub images
Definition: Background.h:339
virtual void operator-=(float const delta)
Subtract a scalar from the Background (equivalent to subtracting a constant from the original image) ...
y_iterator col_end(int x) const
Return an y_iterator to the end of the y&#39;th row.
Definition: Image.h:339
Property getStatisticsProperty() const
Definition: Background.h:209
geom::Box2I _imgBBox
size and origin of input image
Definition: Background.h:330
lsst::afw::image::MaskedImage< InternalPixelT > _statsImage
Definition: Background.h:434
#define LSST_makeBackground_getImage_types
Definition: Background.h:352
Interpolate::Style lookupMaxInterpStyle(int const n)
Get the highest order Interpolation::Style available for &#39;n&#39; points.
Definition: Interpolate.cc:285
int x
int lookupMinInterpPoints(Interpolate::Style const style)
Get the minimum number of points needed to use the requested interpolation style. ...
Definition: Interpolate.cc:309
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
std::vector< double > _xcen
x center pix coords of sub images
Definition: Background.h:335
void _setGridColumns(Interpolate::Style const interpStyle, UndersampleStyle const undersampleStyle, int const iX, std::vector< int > const &ypix) const
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
virtual void operator+=(float const delta)
Add a scalar to the Background (equivalent to adding a constant to the original image) ...
Statistics makeStatistics(afwImage::Mask< afwImage::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl)
Specialization to handle Masks.
Definition: Statistics.cc:1082
Compute Image Statistics.
Implementation of the Class MaskedImage.
int getHeight() const
Return the number of rows in the image.
Definition: Image.h:239
A class to evaluate image background levels.
Definition: Background.h:403
std::vector< std::vector< double > > _gridColumns
Definition: Background.h:435
int getHeight() const
Return the number of rows in the image.
Definition: MaskedImage.h:903
std::vector< int > _yorig
y origin ...
Definition: Background.h:338
#define LSST_EXCEPT_ADD(e, m)
Definition: Exception.h:51
geom::Extent2I getDimensions() const
Return the image&#39;s size; useful for passing to constructors.
Definition: Image.h:298