LSST Applications g0b6bd0c080+a72a5dd7e6,g1182afd7b4+2a019aa3bb,g17e5ecfddb+2b8207f7de,g1d67935e3f+06cf436103,g38293774b4+ac198e9f13,g396055baef+6a2097e274,g3b44f30a73+6611e0205b,g480783c3b1+98f8679e14,g48ccf36440+89c08d0516,g4b93dc025c+98f8679e14,g5c4744a4d9+a302e8c7f0,g613e996a0d+e1c447f2e0,g6c8d09e9e7+25247a063c,g7271f0639c+98f8679e14,g7a9cd813b8+124095ede6,g9d27549199+a302e8c7f0,ga1cf026fa3+ac198e9f13,ga32aa97882+7403ac30ac,ga786bb30fb+7a139211af,gaa63f70f4e+9994eb9896,gabf319e997+ade567573c,gba47b54d5d+94dc90c3ea,gbec6a3398f+06cf436103,gc6308e37c7+07dd123edb,gc655b1545f+ade567573c,gcc9029db3c+ab229f5caf,gd01420fc67+06cf436103,gd877ba84e5+06cf436103,gdb4cecd868+6f279b5b48,ge2d134c3d5+cc4dbb2e3f,ge448b5faa6+86d1ceac1d,gecc7e12556+98f8679e14,gf3ee170dca+25247a063c,gf4ac96e456+ade567573c,gf9f5ea5b4d+ac198e9f13,gff490e6085+8c2580be5c,w.2022.27
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
37namespace lsst {
38namespace ex = pex::exceptions;
39
40namespace afw {
41namespace math {
42
43namespace {
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')
48void 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
79template <typename ImageT>
80BackgroundMI::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
89 image::MaskedImage<InternalPixelT>::Image& im = *_statsImage.getImage();
90 image::MaskedImage<InternalPixelT>::Variance& var = *_statsImage.getVariance();
91
92 for (int iX = 0; iX < nxSample; ++iX) {
93 for (int iY = 0; iY < nySample; ++iY) {
94 ImageT subimg = ImageT(img,
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
111void BackgroundMI::_setGridColumns(Interpolate::Style const interpStyle,
112 UndersampleStyle const undersampleStyle, int const iX,
113 std::vector<int> const& ypix) const {
114 image::MaskedImage<InternalPixelT>::Image& im = *_statsImage.getImage();
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);
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 }
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
179template <typename PixelT>
180std::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;
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;
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;
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
340template <typename PixelT>
341std::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
380} // 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
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:73
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 size(T... args)