LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
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 <iostream>
29 #include <limits>
30 #include <vector>
31 #include <cmath>
37 
38 namespace lsst {
39 namespace ex = pex::exceptions;
40 
41 namespace afw {
42 namespace math {
43 
44 namespace {
45 
46 // Given two vectors x and y, with some nans in y we want vectors x' and y' that correspond to the data
47 // without the nans basic idea is that 'x' is the values, and 'y' is the ref (where nan checking happens)
48 // cullNan(x, y, x', y')
49 void cullNan(std::vector<double> const& values, std::vector<double> const& refs,
50  std::vector<double>& culledValues, std::vector<double>& culledRefs,
51  double const defaultValue = std::numeric_limits<double>::quiet_NaN()) {
52  if (culledValues.capacity() == 0) {
53  culledValues.reserve(refs.size());
54  } else {
55  culledValues.clear();
56  }
57  if (culledRefs.capacity() == 0) {
58  culledRefs.reserve(refs.size());
59  } else {
60  culledRefs.clear();
61  }
62 
63  bool const haveDefault = !std::isnan(defaultValue);
64 
65  for (std::vector<double>::const_iterator pVal = values.begin(), pRef = refs.begin(); pRef != refs.end();
66  ++pRef, ++pVal) {
67  if (!std::isnan(*pRef)) {
68  culledValues.push_back(*pVal);
69  culledRefs.push_back(*pRef);
70  } else if (haveDefault) {
71  culledValues.push_back(*pVal);
72  culledRefs.push_back(defaultValue);
73  } else {
74  ; // drop a NaN
75  }
76  }
77 }
78 } // namespace
79 
80 template <typename ImageT>
81 BackgroundMI::BackgroundMI(ImageT const& img, BackgroundControl const& bgCtrl)
82  : Background(img, bgCtrl), _statsImage(image::MaskedImage<InternalPixelT>()) {
83  // =============================================================
84  // Loop over the cells in the image, computing statistical properties
85  // of each cell in turn and using them to set _statsImage
86  int const nxSample = bgCtrl.getNxSample();
87  int const nySample = bgCtrl.getNySample();
88  _statsImage = image::MaskedImage<InternalPixelT>(nxSample, nySample);
89 
92 
93  for (int iX = 0; iX < nxSample; ++iX) {
94  for (int iY = 0; iY < nySample; ++iY) {
95  ImageT subimg = ImageT(img,
98  image::LOCAL);
99 
101  *bgCtrl.getStatisticsControl())
102  .getResult();
103  im(iX, iY) = res.first;
104  var(iX, iY) = res.second;
105  }
106  }
107 }
109  image::MaskedImage<InternalPixelT> const& statsImage)
110  : Background(imageBBox, statsImage.getWidth(), statsImage.getHeight()), _statsImage(statsImage) {}
111 
112 void BackgroundMI::_setGridColumns(Interpolate::Style const interpStyle,
113  UndersampleStyle const undersampleStyle, int const iX,
114  std::vector<int> const& ypix) const {
116 
117  int const height = _imgBBox.getHeight();
118  _gridColumns[iX].resize(height);
119 
120  // Set _grid as a transitional measure
121  std::vector<double> _grid(_statsImage.getHeight());
122  std::copy(im.col_begin(iX), im.col_end(iX), _grid.begin());
123 
124  // remove nan from the grid values before computing columns
125  // if we do it here (ie. in _setGridColumns), it should
126  // take care of all future occurrences, so we don't need to do this elsewhere
127  std::vector<double> ycenTmp, gridTmp;
128  cullNan(_ycen, _grid, ycenTmp, gridTmp);
129 
131  try {
132  intobj = makeInterpolate(ycenTmp, gridTmp, interpStyle);
133  } catch (pex::exceptions::OutOfRangeError& e) {
134  switch (undersampleStyle) {
135  case THROW_EXCEPTION:
136  LSST_EXCEPT_ADD(e, "setting _gridcolumns");
137  throw;
138  case REDUCE_INTERP_ORDER: {
139  if (gridTmp.empty()) {
140  // Set the column to NaN. We'll deal with this properly when interpolating in x
141  ycenTmp.push_back(0);
143 
144  intobj = makeInterpolate(ycenTmp, gridTmp, Interpolate::CONSTANT);
145  break;
146  } else {
147  return _setGridColumns(lookupMaxInterpStyle(gridTmp.size()), undersampleStyle, iX, ypix);
148  }
149  }
150  case INCREASE_NXNYSAMPLE:
152  e, "The BackgroundControl UndersampleStyle INCREASE_NXNYSAMPLE is not supported.");
153  throw;
154  default:
155  LSST_EXCEPT_ADD(e, str(boost::format("The selected BackgroundControl "
156  "UndersampleStyle %d is not defined.") %
157  undersampleStyle));
158  throw;
159  }
160  } catch (ex::Exception& e) {
161  LSST_EXCEPT_ADD(e, "setting _gridcolumns");
162  throw;
163  }
164 
165  for (int iY = 0; iY < height; ++iY) {
166  _gridColumns[iX][iY] = intobj->interpolate(ypix[iY]);
167  }
168 }
169 
171  _statsImage += delta;
172  return *this;
173 }
174 
176  _statsImage -= delta;
177  return *this;
178 }
179 
180 template <typename PixelT>
181 std::shared_ptr<image::Image<PixelT>> BackgroundMI::doGetImage(
182  lsst::geom::Box2I const& bbox,
183  Interpolate::Style const interpStyle_, // Style of the interpolation
184  UndersampleStyle const undersampleStyle // Behaviour if there are too few points
185  ) const {
186  if (!_imgBBox.contains(bbox)) {
187  throw LSST_EXCEPT(
189  str(boost::format("BBox (%d:%d,%d:%d) out of range (%d:%d,%d:%d)") % bbox.getMinX() %
190  bbox.getMaxX() % bbox.getMinY() % bbox.getMaxY() % _imgBBox.getMinX() %
192  }
193  int const nxSample = _statsImage.getWidth();
194  int const nySample = _statsImage.getHeight();
195  Interpolate::Style interpStyle = interpStyle_; // not const -- may be modified if REDUCE_INTERP_ORDER
196 
197  /*
198  * Save the as-used interpStyle and undersampleStyle.
199  *
200  * N.b. The undersampleStyle may actually be overridden for some columns of the statsImage if they
201  * have too few good values. This doesn't prevent you reproducing the results of getImage() by
202  * calling getImage(getInterpStyle(), getUndersampleStyle())
203  */
204  _asUsedInterpStyle = interpStyle;
205  _asUsedUndersampleStyle = undersampleStyle;
206  /*
207  * Check if the requested nx,ny are sufficient for the requested interpolation style,
208  * making suitable adjustments
209  */
210  bool const isXundersampled = (nxSample < lookupMinInterpPoints(interpStyle));
211  bool const isYundersampled = (nySample < lookupMinInterpPoints(interpStyle));
212 
213  switch (undersampleStyle) {
214  case THROW_EXCEPTION:
215  if (isXundersampled && isYundersampled) {
216  throw LSST_EXCEPT(
218  "nxSample and nySample have too few points for requested interpolation style.");
219  } else if (isXundersampled) {
221  "nxSample has too few points for requested interpolation style.");
222  } else if (isYundersampled) {
224  "nySample has too few points for requested interpolation style.");
225  }
226  break;
227  case REDUCE_INTERP_ORDER:
228  if (isXundersampled || isYundersampled) {
229  Interpolate::Style const xStyle = lookupMaxInterpStyle(nxSample);
230  Interpolate::Style const yStyle = lookupMaxInterpStyle(nySample);
231  interpStyle = (nxSample < nySample) ? xStyle : yStyle;
232  _asUsedInterpStyle = interpStyle;
233  }
234  break;
235  case INCREASE_NXNYSAMPLE:
236  if (isXundersampled || isYundersampled) {
237  throw LSST_EXCEPT(
239  "The BackgroundControl UndersampleStyle INCREASE_NXNYSAMPLE is not supported.");
240  }
241  break;
242  default:
244  str(boost::format("The selected BackgroundControl "
245  "UndersampleStyle %d is not defined.") %
246  undersampleStyle));
247  }
248 
249  // if we're approximating, don't bother with the rest of the interp-related work. Return from here.
250  if (_bctrl->getApproximateControl()->getStyle() != ApproximateControl::UNKNOWN) {
251  return doGetApproximate<PixelT>(*_bctrl->getApproximateControl(), _asUsedUndersampleStyle)
252  ->getImage();
253  }
254 
255  // =============================================================
256  // --> We'll store nxSample fully-interpolated columns to interpolate the rows over
257  // make a vector containing the y pixel coords for the column
258  int const width = _imgBBox.getWidth();
259  int const height = _imgBBox.getHeight();
260  auto const bboxOff = bbox.getMin() - _imgBBox.getMin();
261 
262  std::vector<int> ypix(height);
263  for (int iY = 0; iY < height; ++iY) {
264  ypix[iY] = iY;
265  }
266 
267  _gridColumns.resize(width);
268  for (int iX = 0; iX < nxSample; ++iX) {
269  _setGridColumns(interpStyle, undersampleStyle, iX, ypix);
270  }
271 
272  // create a shared_ptr to put the background image in and return to caller
273  // start with xy0 = 0 and set final xy0 later
276 
277  // go through row by row
278  // - interpolate on the gridcolumns that were pre-computed by the constructor
279  // - copy the values to an ImageT to return to the caller.
280  std::vector<double> xcenTmp, bgTmp;
281 
282  // N.b. There's no API to set defaultValue to other than NaN (due to issues with persistence
283  // that I don't feel like fixing; #2825). If we want to address this, this is the place
284  // to start, but note that NaN is treated specially -- it means, "Interpolate" so to allow
285  // us to put a NaN into the outputs some changes will be needed
286  double defaultValue = std::numeric_limits<double>::quiet_NaN();
287 
288  for (int y = 0, iY = bboxOff.getY(); y < bbox.getHeight(); ++y, ++iY) {
289  // build an interp object for this row
290  std::vector<double> bg_x(nxSample);
291  for (int iX = 0; iX < nxSample; iX++) {
292  bg_x[iX] = static_cast<double>(_gridColumns[iX][iY]);
293  }
294  cullNan(_xcen, bg_x, xcenTmp, bgTmp, defaultValue);
295 
297  try {
298  intobj = makeInterpolate(xcenTmp, bgTmp, interpStyle);
299  } catch (pex::exceptions::OutOfRangeError& e) {
300  switch (undersampleStyle) {
301  case THROW_EXCEPTION:
302  LSST_EXCEPT_ADD(e, str(boost::format("Interpolating in y (iY = %d)") % iY));
303  throw;
304  case REDUCE_INTERP_ORDER: {
305  if (bgTmp.empty()) {
306  xcenTmp.push_back(0);
307  bgTmp.push_back(defaultValue);
308 
309  intobj = makeInterpolate(xcenTmp, bgTmp, Interpolate::CONSTANT);
310  break;
311  } else {
312  intobj = makeInterpolate(xcenTmp, bgTmp, lookupMaxInterpStyle(bgTmp.size()));
313  }
314  } break;
315  case INCREASE_NXNYSAMPLE:
317  e,
318  "The BackgroundControl UndersampleStyle INCREASE_NXNYSAMPLE is not supported.");
319  throw;
320  default:
321  LSST_EXCEPT_ADD(e, str(boost::format("The selected BackgroundControl "
322  "UndersampleStyle %d is not defined.") %
323  undersampleStyle));
324  throw;
325  }
326  } catch (ex::Exception& e) {
327  LSST_EXCEPT_ADD(e, str(boost::format("Interpolating in y (iY = %d)") % iY));
328  throw;
329  }
330 
331  // fill the image with interpolated values
332  for (int iX = bboxOff.getX(), x = 0; x < bbox.getWidth(); ++iX, ++x) {
333  (*bg)(x, y) = static_cast<PixelT>(intobj->interpolate(iX));
334  }
335  }
336  bg->setXY0(bbox.getMin());
337 
338  return bg;
339 }
340 
341 template <typename PixelT>
342 std::shared_ptr<Approximate<PixelT>> BackgroundMI::doGetApproximate(
343  ApproximateControl const& actrl, /* Approximation style */
344  UndersampleStyle const undersampleStyle /* Behaviour if there are too few points */
345  ) const {
346  auto const localBBox = lsst::geom::Box2I(lsst::geom::Point2I(0, 0), _imgBBox.getDimensions());
347  return makeApproximate(_xcen, _ycen, _statsImage, localBBox, actrl);
348 }
349 
351 /*
352  * Create the versions we need of _get{Approximate,Image} and Explicit instantiations
353  *
354  */
355 #define CREATE_BACKGROUND(m, v, TYPE) \
356  template BackgroundMI::BackgroundMI(image::Image<TYPE> const& img, BackgroundControl const& bgCtrl); \
357  template BackgroundMI::BackgroundMI(image::MaskedImage<TYPE> const& img, \
358  BackgroundControl const& bgCtrl); \
359  std::shared_ptr<image::Image<TYPE>> BackgroundMI::_getImage( \
360  lsst::geom::Box2I const& bbox, \
361  Interpolate::Style const interpStyle, /* Style of the interpolation */ \
362  UndersampleStyle const undersampleStyle, /* Behaviour if there are too few points */ \
363  TYPE /* disambiguate */ \
364  ) const { \
365  return BackgroundMI::doGetImage<TYPE>(bbox, interpStyle, undersampleStyle); \
366  }
367 
368 #define CREATE_getApproximate(m, v, TYPE) \
369  std::shared_ptr<Approximate<TYPE>> BackgroundMI::_getApproximate( \
370  ApproximateControl const& actrl, /* Approximation style */ \
371  UndersampleStyle const undersampleStyle, /* Behaviour if there are too few points */ \
372  TYPE /* disambiguate */ \
373  ) const { \
374  return BackgroundMI::doGetApproximate<TYPE>(actrl, undersampleStyle); \
375  }
376 
379 
380 } // namespace math
382 } // namespace afw
383 } // 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:49
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:58
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
VariancePtr getVariance() const
Return a (shared_ptr to) the MaskedImage's variance.
Definition: MaskedImage.h:1079
ImagePtr getImage() const
Return a (shared_ptr to) the MaskedImage's image.
Definition: MaskedImage.h:1046
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:81
Value getResult(Property const prop=NOTHING) const
Return the value and error in the specified statistic (e.g.
Definition: Statistics.cc:931
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:275
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:354
@ ERRORS
Include errors of requested quantities.
Definition: Statistics.h:65
int lookupMinInterpPoints(Interpolate::Style const style)
Get the minimum number of points needed to use the requested interpolation style.
Definition: Interpolate.cc:314
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:343
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)