LSST Applications  21.0.0-147-g0e635eb1+1acddb5be5,22.0.0+052faf71bd,22.0.0+1ea9a8b2b2,22.0.0+6312710a6c,22.0.0+729191ecac,22.0.0+7589c3a021,22.0.0+9f079a9461,22.0.1-1-g7d6de66+b8044ec9de,22.0.1-1-g87000a6+536b1ee016,22.0.1-1-g8e32f31+6312710a6c,22.0.1-10-gd060f87+016f7cdc03,22.0.1-12-g9c3108e+df145f6f68,22.0.1-16-g314fa6d+c825727ab8,22.0.1-19-g93a5c75+d23f2fb6d8,22.0.1-19-gb93eaa13+aab3ef7709,22.0.1-2-g8ef0a89+b8044ec9de,22.0.1-2-g92698f7+9f079a9461,22.0.1-2-ga9b0f51+052faf71bd,22.0.1-2-gac51dbf+052faf71bd,22.0.1-2-gb66926d+6312710a6c,22.0.1-2-gcb770ba+09e3807989,22.0.1-20-g32debb5+b8044ec9de,22.0.1-23-gc2439a9a+fb0756638e,22.0.1-3-g496fd5d+09117f784f,22.0.1-3-g59f966b+1e6ba2c031,22.0.1-3-g849a1b8+f8b568069f,22.0.1-3-gaaec9c0+c5c846a8b1,22.0.1-32-g5ddfab5d3+60ce4897b0,22.0.1-4-g037fbe1+64e601228d,22.0.1-4-g8623105+b8044ec9de,22.0.1-5-g096abc9+d18c45d440,22.0.1-5-g15c806e+57f5c03693,22.0.1-7-gba73697+57f5c03693,master-g6e05de7fdc+c1283a92b8,master-g72cdda8301+729191ecac,w.2021.39
LSST Data Management Base Package
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 
25 /*
26  * Background estimation class code
27  */
28 #include <limits>
29 #include <vector>
30 #include <cmath>
36 
37 namespace lsst {
38 namespace ex = pex::exceptions;
39 
40 namespace afw {
41 namespace math {
42 
43 namespace {
44 
45 // Given two vectors x and y, with some nans in y we want vectors x' and y' that correspond to the data
46 // without the nans basic idea is that 'x' is the values, and 'y' is the ref (where nan checking happens)
47 // cullNan(x, y, x', y')
48 void cullNan(std::vector<double> const& values, std::vector<double> const& refs,
49  std::vector<double>& culledValues, std::vector<double>& culledRefs,
50  double const defaultValue = std::numeric_limits<double>::quiet_NaN()) {
51  if (culledValues.capacity() == 0) {
52  culledValues.reserve(refs.size());
53  } else {
54  culledValues.clear();
55  }
56  if (culledRefs.capacity() == 0) {
57  culledRefs.reserve(refs.size());
58  } else {
59  culledRefs.clear();
60  }
61 
62  bool const haveDefault = !std::isnan(defaultValue);
63 
64  for (std::vector<double>::const_iterator pVal = values.begin(), pRef = refs.begin(); pRef != refs.end();
65  ++pRef, ++pVal) {
66  if (!std::isnan(*pRef)) {
67  culledValues.push_back(*pVal);
68  culledRefs.push_back(*pRef);
69  } else if (haveDefault) {
70  culledValues.push_back(*pVal);
71  culledRefs.push_back(defaultValue);
72  } else {
73  ; // drop a NaN
74  }
75  }
76 }
77 } // namespace
78 
79 template <typename ImageT>
80 BackgroundMI::BackgroundMI(ImageT const& img, BackgroundControl const& bgCtrl)
81  : Background(img, bgCtrl), _statsImage(image::MaskedImage<InternalPixelT>()) {
82  // =============================================================
83  // Loop over the cells in the image, computing statistical properties
84  // of each cell in turn and using them to set _statsImage
85  int const nxSample = bgCtrl.getNxSample();
86  int const nySample = bgCtrl.getNySample();
87  _statsImage = image::MaskedImage<InternalPixelT>(nxSample, nySample);
88 
91 
92  for (int iX = 0; iX < nxSample; ++iX) {
93  for (int iY = 0; iY < nySample; ++iY) {
94  ImageT subimg = ImageT(img,
97  image::LOCAL);
98 
100  *bgCtrl.getStatisticsControl())
101  .getResult();
102  im(iX, iY) = res.first;
103  var(iX, iY) = res.second;
104  }
105  }
106 }
108  image::MaskedImage<InternalPixelT> const& statsImage)
109  : Background(imageBBox, statsImage.getWidth(), statsImage.getHeight()), _statsImage(statsImage) {}
110 
111 void BackgroundMI::_setGridColumns(Interpolate::Style const interpStyle,
112  UndersampleStyle const undersampleStyle, int const iX,
113  std::vector<int> const& ypix) const {
115 
116  int const height = _imgBBox.getHeight();
117  _gridColumns[iX].resize(height);
118 
119  // Set _grid as a transitional measure
120  std::vector<double> _grid(_statsImage.getHeight());
121  std::copy(im.col_begin(iX), im.col_end(iX), _grid.begin());
122 
123  // remove nan from the grid values before computing columns
124  // if we do it here (ie. in _setGridColumns), it should
125  // take care of all future occurrences, so we don't need to do this elsewhere
126  std::vector<double> ycenTmp, gridTmp;
127  cullNan(_ycen, _grid, ycenTmp, gridTmp);
128 
130  try {
131  intobj = makeInterpolate(ycenTmp, gridTmp, interpStyle);
132  } catch (pex::exceptions::OutOfRangeError& e) {
133  switch (undersampleStyle) {
134  case THROW_EXCEPTION:
135  LSST_EXCEPT_ADD(e, "setting _gridcolumns");
136  throw;
137  case REDUCE_INTERP_ORDER: {
138  if (gridTmp.empty()) {
139  // Set the column to NaN. We'll deal with this properly when interpolating in x
140  ycenTmp.push_back(0);
142 
143  intobj = makeInterpolate(ycenTmp, gridTmp, Interpolate::CONSTANT);
144  break;
145  } else {
146  return _setGridColumns(lookupMaxInterpStyle(gridTmp.size()), undersampleStyle, iX, ypix);
147  }
148  }
149  case INCREASE_NXNYSAMPLE:
151  e, "The BackgroundControl UndersampleStyle INCREASE_NXNYSAMPLE is not supported.");
152  throw;
153  default:
154  LSST_EXCEPT_ADD(e, str(boost::format("The selected BackgroundControl "
155  "UndersampleStyle %d is not defined.") %
156  undersampleStyle));
157  throw;
158  }
159  } catch (ex::Exception& e) {
160  LSST_EXCEPT_ADD(e, "setting _gridcolumns");
161  throw;
162  }
163 
164  for (int iY = 0; iY < height; ++iY) {
165  _gridColumns[iX][iY] = intobj->interpolate(ypix[iY]);
166  }
167 }
168 
170  _statsImage += delta;
171  return *this;
172 }
173 
175  _statsImage -= delta;
176  return *this;
177 }
178 
179 template <typename PixelT>
180 std::shared_ptr<image::Image<PixelT>> BackgroundMI::doGetImage(
181  lsst::geom::Box2I const& bbox,
182  Interpolate::Style const interpStyle_, // Style of the interpolation
183  UndersampleStyle const undersampleStyle // Behaviour if there are too few points
184  ) const {
185  if (!_imgBBox.contains(bbox)) {
186  throw LSST_EXCEPT(
188  str(boost::format("BBox (%d:%d,%d:%d) out of range (%d:%d,%d:%d)") % bbox.getMinX() %
189  bbox.getMaxX() % bbox.getMinY() % bbox.getMaxY() % _imgBBox.getMinX() %
191  }
192  int const nxSample = _statsImage.getWidth();
193  int const nySample = _statsImage.getHeight();
194  Interpolate::Style interpStyle = interpStyle_; // not const -- may be modified if REDUCE_INTERP_ORDER
195 
196  /*
197  * Save the as-used interpStyle and undersampleStyle.
198  *
199  * N.b. The undersampleStyle may actually be overridden for some columns of the statsImage if they
200  * have too few good values. This doesn't prevent you reproducing the results of getImage() by
201  * calling getImage(getInterpStyle(), getUndersampleStyle())
202  */
203  _asUsedInterpStyle = interpStyle;
204  _asUsedUndersampleStyle = undersampleStyle;
205  /*
206  * Check if the requested nx,ny are sufficient for the requested interpolation style,
207  * making suitable adjustments
208  */
209  bool const isXundersampled = (nxSample < lookupMinInterpPoints(interpStyle));
210  bool const isYundersampled = (nySample < lookupMinInterpPoints(interpStyle));
211 
212  switch (undersampleStyle) {
213  case THROW_EXCEPTION:
214  if (isXundersampled && isYundersampled) {
215  throw LSST_EXCEPT(
217  "nxSample and nySample have too few points for requested interpolation style.");
218  } else if (isXundersampled) {
220  "nxSample has too few points for requested interpolation style.");
221  } else if (isYundersampled) {
223  "nySample has too few points for requested interpolation style.");
224  }
225  break;
226  case REDUCE_INTERP_ORDER:
227  if (isXundersampled || isYundersampled) {
228  Interpolate::Style const xStyle = lookupMaxInterpStyle(nxSample);
229  Interpolate::Style const yStyle = lookupMaxInterpStyle(nySample);
230  interpStyle = (nxSample < nySample) ? xStyle : yStyle;
231  _asUsedInterpStyle = interpStyle;
232  }
233  break;
234  case INCREASE_NXNYSAMPLE:
235  if (isXundersampled || isYundersampled) {
236  throw LSST_EXCEPT(
238  "The BackgroundControl UndersampleStyle INCREASE_NXNYSAMPLE is not supported.");
239  }
240  break;
241  default:
243  str(boost::format("The selected BackgroundControl "
244  "UndersampleStyle %d is not defined.") %
245  undersampleStyle));
246  }
247 
248  // if we're approximating, don't bother with the rest of the interp-related work. Return from here.
249  if (_bctrl->getApproximateControl()->getStyle() != ApproximateControl::UNKNOWN) {
250  return doGetApproximate<PixelT>(*_bctrl->getApproximateControl(), _asUsedUndersampleStyle)
251  ->getImage();
252  }
253 
254  // =============================================================
255  // --> We'll store nxSample fully-interpolated columns to interpolate the rows over
256  // make a vector containing the y pixel coords for the column
257  int const width = _imgBBox.getWidth();
258  int const height = _imgBBox.getHeight();
259  auto const bboxOff = bbox.getMin() - _imgBBox.getMin();
260 
261  std::vector<int> ypix(height);
262  for (int iY = 0; iY < height; ++iY) {
263  ypix[iY] = iY;
264  }
265 
266  _gridColumns.resize(width);
267  for (int iX = 0; iX < nxSample; ++iX) {
268  _setGridColumns(interpStyle, undersampleStyle, iX, ypix);
269  }
270 
271  // create a shared_ptr to put the background image in and return to caller
272  // start with xy0 = 0 and set final xy0 later
275 
276  // go through row by row
277  // - interpolate on the gridcolumns that were pre-computed by the constructor
278  // - copy the values to an ImageT to return to the caller.
279  std::vector<double> xcenTmp, bgTmp;
280 
281  // N.b. There's no API to set defaultValue to other than NaN (due to issues with persistence
282  // that I don't feel like fixing; #2825). If we want to address this, this is the place
283  // to start, but note that NaN is treated specially -- it means, "Interpolate" so to allow
284  // us to put a NaN into the outputs some changes will be needed
285  double defaultValue = std::numeric_limits<double>::quiet_NaN();
286 
287  for (int y = 0, iY = bboxOff.getY(); y < bbox.getHeight(); ++y, ++iY) {
288  // build an interp object for this row
289  std::vector<double> bg_x(nxSample);
290  for (int iX = 0; iX < nxSample; iX++) {
291  bg_x[iX] = static_cast<double>(_gridColumns[iX][iY]);
292  }
293  cullNan(_xcen, bg_x, xcenTmp, bgTmp, defaultValue);
294 
296  try {
297  intobj = makeInterpolate(xcenTmp, bgTmp, interpStyle);
298  } catch (pex::exceptions::OutOfRangeError& e) {
299  switch (undersampleStyle) {
300  case THROW_EXCEPTION:
301  LSST_EXCEPT_ADD(e, str(boost::format("Interpolating in y (iY = %d)") % iY));
302  throw;
303  case REDUCE_INTERP_ORDER: {
304  if (bgTmp.empty()) {
305  xcenTmp.push_back(0);
306  bgTmp.push_back(defaultValue);
307 
308  intobj = makeInterpolate(xcenTmp, bgTmp, Interpolate::CONSTANT);
309  break;
310  } else {
311  intobj = makeInterpolate(xcenTmp, bgTmp, lookupMaxInterpStyle(bgTmp.size()));
312  }
313  } break;
314  case INCREASE_NXNYSAMPLE:
316  e,
317  "The BackgroundControl UndersampleStyle INCREASE_NXNYSAMPLE is not supported.");
318  throw;
319  default:
320  LSST_EXCEPT_ADD(e, str(boost::format("The selected BackgroundControl "
321  "UndersampleStyle %d is not defined.") %
322  undersampleStyle));
323  throw;
324  }
325  } catch (ex::Exception& e) {
326  LSST_EXCEPT_ADD(e, str(boost::format("Interpolating in y (iY = %d)") % iY));
327  throw;
328  }
329 
330  // fill the image with interpolated values
331  for (int iX = bboxOff.getX(), x = 0; x < bbox.getWidth(); ++iX, ++x) {
332  (*bg)(x, y) = static_cast<PixelT>(intobj->interpolate(iX));
333  }
334  }
335  bg->setXY0(bbox.getMin());
336 
337  return bg;
338 }
339 
340 template <typename PixelT>
341 std::shared_ptr<Approximate<PixelT>> BackgroundMI::doGetApproximate(
342  ApproximateControl const& actrl, /* Approximation style */
343  UndersampleStyle const undersampleStyle /* Behaviour if there are too few points */
344  ) const {
345  auto const localBBox = lsst::geom::Box2I(lsst::geom::Point2I(0, 0), _imgBBox.getDimensions());
346  return makeApproximate(_xcen, _ycen, _statsImage, localBBox, actrl);
347 }
348 
350 /*
351  * Create the versions we need of _get{Approximate,Image} and Explicit instantiations
352  *
353  */
354 #define CREATE_BACKGROUND(m, v, TYPE) \
355  template BackgroundMI::BackgroundMI(image::Image<TYPE> const& img, BackgroundControl const& bgCtrl); \
356  template BackgroundMI::BackgroundMI(image::MaskedImage<TYPE> const& img, \
357  BackgroundControl const& bgCtrl); \
358  std::shared_ptr<image::Image<TYPE>> BackgroundMI::_getImage( \
359  lsst::geom::Box2I const& bbox, \
360  Interpolate::Style const interpStyle, /* Style of the interpolation */ \
361  UndersampleStyle const undersampleStyle, /* Behaviour if there are too few points */ \
362  TYPE /* disambiguate */ \
363  ) const { \
364  return BackgroundMI::doGetImage<TYPE>(bbox, interpStyle, undersampleStyle); \
365  }
366 
367 #define CREATE_getApproximate(m, v, TYPE) \
368  std::shared_ptr<Approximate<TYPE>> BackgroundMI::_getApproximate( \
369  ApproximateControl const& actrl, /* Approximation style */ \
370  UndersampleStyle const undersampleStyle, /* Behaviour if there are too few points */ \
371  TYPE /* disambiguate */ \
372  ) const { \
373  return BackgroundMI::doGetApproximate<TYPE>(actrl, undersampleStyle); \
374  }
375 
378 
379 } // namespace math
381 } // namespace afw
382 } // namespace lsst
AmpInfoBoxKey bbox
Definition: Amplifier.cc:117
#define LSST_makeBackground_getApproximate_types
Definition: Background.h:385
#define LSST_makeBackground_getImage_types
Definition: Background.h:384
double x
#define LSST_EXCEPT_ADD(e, m)
Add the current location and a message to an existing exception before rethrowing it.
Definition: Exception.h:54
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
int y
Definition: SpanSet.cc:48
T begin(T... args)
T capacity(T... args)
y_iterator col_begin(int x) const
Return an y_iterator to the start of the y'th row.
Definition: ImageBase.h:413
y_iterator col_end(int x) const
Return an y_iterator to the start of the y'th row.
Definition: ImageBase.h:416
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:51
int getHeight() const
Return the number of rows in the image.
Definition: MaskedImage.h:1056
int getWidth() const
Return the number of columns in the image.
Definition: MaskedImage.h:1054
VariancePtr getVariance() const
Return a (shared_ptr to) the MaskedImage's variance.
Definition: MaskedImage.h:1051
ImagePtr getImage() const
Return a (shared_ptr to) the MaskedImage's image.
Definition: MaskedImage.h:1018
Pass parameters to a Background object.
Definition: Background.h:56
std::shared_ptr< StatisticsControl > getStatisticsControl()
Definition: Background.h:211
Property getStatisticsProperty() const
Definition: Background.h:214
A virtual base class to evaluate image background levels.
Definition: Background.h:235
std::vector< double > _ycen
y center ...
Definition: Background.h:369
std::vector< double > _xcen
x center pix coords of sub images
Definition: Background.h:368
std::vector< int > _xsize
x size of sub images
Definition: Background.h:372
lsst::geom::Box2I _imgBBox
size and origin of input image
Definition: Background.h:363
std::vector< int > _xorig
x origin pix coords of sub images
Definition: Background.h:370
UndersampleStyle _asUsedUndersampleStyle
the undersampleStyle we actually used
Definition: Background.h:366
std::vector< int > _ysize
y size ...
Definition: Background.h:373
float InternalPixelT
type used for any internal images, and returned by getApproximate
Definition: Background.h:262
std::vector< int > _yorig
y origin ...
Definition: Background.h:371
Interpolate::Style _asUsedInterpStyle
the style we actually used
Definition: Background.h:365
std::shared_ptr< BackgroundControl > _bctrl
control info set by user.
Definition: Background.h:364
A class to evaluate image background levels.
Definition: Background.h:434
BackgroundMI & operator-=(float const delta) override
Subtract a scalar from the Background (equivalent to subtracting a constant from the original image)
BackgroundMI & operator+=(float const delta) override
Add a scalar to the Background (equivalent to adding a constant to the original image)
BackgroundMI(ImageT const &img, BackgroundControl const &bgCtrl)
Constructor for BackgroundMI.
Definition: BackgroundMI.cc:80
Value getResult(Property const prop=NOTHING) const
Return the value and error in the specified statistic (e.g.
Definition: Statistics.cc:922
An integer coordinate rectangle.
Definition: Box.h:55
int getMinY() const noexcept
Definition: Box.h:158
int getHeight() const noexcept
Definition: Box.h:188
Point2I const getMin() const noexcept
Definition: Box.h:156
int getMinX() const noexcept
Definition: Box.h:157
int getWidth() const noexcept
Definition: Box.h:187
int getMaxX() const noexcept
Definition: Box.h:161
bool contains(Point2I const &point) const noexcept
Return true if the box contains the point.
Definition: Box.cc:114
int getMaxY() const noexcept
Definition: Box.h:162
Extent2I const getDimensions() const noexcept
Definition: Box.h:186
Provides consistent interface for LSST exceptions.
Definition: Exception.h:107
Reports invalid arguments.
Definition: Runtime.h:66
Reports attempts to exceed implementation-defined length limits for some classes.
Definition: Runtime.h:76
Reports attempts to access elements outside a valid range of indices.
Definition: Runtime.h:89
T clear(T... args)
T copy(T... args)
T empty(T... args)
T end(T... args)
T isnan(T... args)
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
Interpolate::Style lookupMaxInterpStyle(int const n)
Get the highest order Interpolation::Style available for 'n' points.
Definition: Interpolate.cc:274
Statistics makeStatistics(lsst::afw::image::Image< Pixel > const &img, lsst::afw::image::Mask< image::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Handle a watered-down front-end to the constructor (no variance)
Definition: Statistics.h:359
@ ERRORS
Include errors of requested quantities.
Definition: Statistics.h:64
int lookupMinInterpPoints(Interpolate::Style const style)
Get the minimum number of points needed to use the requested interpolation style.
Definition: Interpolate.cc:313
std::shared_ptr< Interpolate > makeInterpolate(std::vector< double > const &x, std::vector< double > const &y, Interpolate::Style const style=Interpolate::AKIMA_SPLINE)
A factory function to make Interpolate objects.
Definition: Interpolate.cc:342
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
BOOST_PP_SEQ_FOR_EACH(INSTANTIATE_COLUMNVIEW_SCALAR, _, BOOST_PP_TUPLE_TO_SEQ(AFW_TABLE_SCALAR_FIELD_TYPE_N, AFW_TABLE_SCALAR_FIELD_TYPE_TUPLE)) BOOST_PP_SEQ_FOR_EACH(INSTANTIATE_COLUMNVIEW_ARRAY
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
A base class for image defects.
T push_back(T... args)
T quiet_NaN(T... args)
T reserve(T... args)
T resize(T... args)
T size(T... args)