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
subtractBackground.py
Go to the documentation of this file.
2# LSST Data Management System
3#
4# Copyright 2008-2017 AURA/LSST.
5#
6# This product includes software developed by the
7# LSST Project (http://www.lsst.org/).
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the LSST License Statement and
20# the GNU General Public License along with this program. If not,
21# see <https://www.lsstcorp.org/LegalNotices/>.
22#
23__all__ = ("SubtractBackgroundConfig", "SubtractBackgroundTask")
24
25import itertools
26
27import numpy
28
29from lsstDebug import getDebugFrame
30from lsst.utils import suppress_deprecations
31
32import lsst.afw.display as afwDisplay
33import lsst.afw.image as afwImage
34import lsst.afw.math as afwMath
35import lsst.pex.config as pexConfig
36import lsst.pipe.base as pipeBase
37
38
39class SubtractBackgroundConfig(pexConfig.Config):
40 """Config for SubtractBackgroundTask
41
42 Many of these fields match fields in `lsst.afw.math.BackgroundControl`,
43 the control class for `lsst.afw.math.makeBackground`
44 """
45 statisticsProperty = pexConfig.ChoiceField(
46 doc="type of statistic to use for grid points",
47 dtype=str, default="MEANCLIP",
48 allowed={
49 "MEANCLIP": "clipped mean",
50 "MEAN": "unclipped mean",
51 "MEDIAN": "median",
52 }
53 )
54 undersampleStyle = pexConfig.ChoiceField(
55 doc="behaviour if there are too few points in grid for requested interpolation style",
56 dtype=str, default="REDUCE_INTERP_ORDER",
57 allowed={
58 "THROW_EXCEPTION": "throw an exception if there are too few points",
59 "REDUCE_INTERP_ORDER": "use an interpolation style with a lower order.",
60 "INCREASE_NXNYSAMPLE": "Increase the number of samples used to make the interpolation grid.",
61 },
62 )
63 binSize = pexConfig.RangeField(
64 doc="how large a region of the sky should be used for each background point",
65 dtype=int, default=128, min=1,
66 )
67 binSizeX = pexConfig.RangeField(
68 doc=("Sky region size to be used for each background point in X direction. "
69 "If 0, the binSize config is used."),
70 dtype=int, default=0, min=0,
71 )
72 binSizeY = pexConfig.RangeField(
73 doc=("Sky region size to be used for each background point in Y direction. "
74 "If 0, the binSize config is used."),
75 dtype=int, default=0, min=0,
76 )
77 algorithm = pexConfig.ChoiceField(
78 doc="how to interpolate the background values. This maps to an enum; see afw::math::Background",
79 dtype=str, default="AKIMA_SPLINE", optional=True,
80 allowed={
81 "CONSTANT": "Use a single constant value",
82 "LINEAR": "Use linear interpolation",
83 "NATURAL_SPLINE": "cubic spline with zero second derivative at endpoints",
84 "AKIMA_SPLINE": "higher-level nonlinear spline that is more robust to outliers",
85 "NONE": "No background estimation is to be attempted",
86 },
87 )
88 ignoredPixelMask = pexConfig.ListField(
89 doc="Names of mask planes to ignore while estimating the background",
90 dtype=str, default=["BAD", "EDGE", "DETECTED", "DETECTED_NEGATIVE", "NO_DATA", ],
91 itemCheck=lambda x: x in afwImage.Mask().getMaskPlaneDict().keys(),
92 )
93 isNanSafe = pexConfig.Field(
94 doc="Ignore NaNs when estimating the background",
95 dtype=bool, default=False,
96 )
97
98 useApprox = pexConfig.Field(
99 doc="Use Approximate (Chebyshev) to model background.",
100 dtype=bool, default=True,
101 )
102 approxOrderX = pexConfig.Field(
103 doc="Approximation order in X for background Chebyshev (valid only with useApprox=True)",
104 dtype=int, default=6,
105 )
106 # Note: Currently X- and Y-orders must be equal due to a limitation in math::Chebyshev1Function2
107 # The following is being added so that the weighting attribute can also be configurable for the
108 # call to afwMath.ApproximateControl
109 approxOrderY = pexConfig.Field(
110 doc="Approximation order in Y for background Chebyshev (valid only with useApprox=True)",
111 dtype=int, default=-1,
112 )
113 weighting = pexConfig.Field(
114 doc="Use inverse variance weighting in calculation (valid only with useApprox=True)",
115 dtype=bool, default=True,
116 )
117
118
119class SubtractBackgroundTask(pipeBase.Task):
120 """Subtract the background from an exposure
121 """
122 ConfigClass = SubtractBackgroundConfig
123 _DefaultName = "subtractBackground"
124
125 def run(self, exposure, background=None, stats=True, statsKeys=None):
126 """Fit and subtract the background of an exposure.
127
128 Parameters
129 ----------
130 exposure : `lsst.afw.image.Exposure`
131 Exposure whose background is to be subtracted.
132 background : `lsst.afw.math.BackgroundList`
133 Initial background model already subtracted. May be None if no background
134 has been subtracted.
135 stats : `bool`
136 If True then measure the mean and variance of the full background model and
137 record the results in the exposure's metadata.
138 statsKeys : `tuple`
139 Key names used to store the mean and variance of the background in the
140 exposure's metadata (another tuple); if None then use ("BGMEAN", "BGVAR");
141 ignored if stats is false.
142
143 Returns
144 -------
145 background : `lsst.afw.math.BackgroundLst`
146 Full background model (initial model with changes), contained in an
147 `lsst.pipe.base.Struct`.
148 """
149 if background is None:
150 background = afwMath.BackgroundList()
151
152 maskedImage = exposure.getMaskedImage()
153 fitBg = self.fitBackgroundfitBackground(maskedImage)
154 maskedImage -= fitBg.getImageF(self.config.algorithm, self.config.undersampleStyle)
155
156 actrl = fitBg.getBackgroundControl().getApproximateControl()
157 background.append((fitBg, getattr(afwMath.Interpolate, self.config.algorithm),
158 fitBg.getAsUsedUndersampleStyle(), actrl.getStyle(),
159 actrl.getOrderX(), actrl.getOrderY(), actrl.getWeighting()))
160
161 if stats:
162 self._addStats_addStats(exposure, background, statsKeys=statsKeys)
163
164 subFrame = getDebugFrame(self._display, "subtracted")
165 if subFrame:
166 subDisp = afwDisplay.getDisplay(frame=subFrame)
167 subDisp.mtv(exposure, title="subtracted")
168
169 bgFrame = getDebugFrame(self._display, "background")
170 if bgFrame:
171 bgDisp = afwDisplay.getDisplay(frame=bgFrame)
172 bgImage = background.getImage()
173 bgDisp.mtv(bgImage, title="background")
174
175 return pipeBase.Struct(
176 background=background,
177 )
178
179 def _addStats(self, exposure, background, statsKeys=None):
180 """Add statistics about the background to the exposure's metadata
181
182 Parameters
183 ----------
184 exposure : `lsst.afw.image.Exposure`
185 Exposure whose background was subtracted.
186 background : `lsst.afw.math.BackgroundList`
187 Background model
188 statsKeys : `tuple`
189 Key names used to store the mean and variance of the background in
190 the exposure's metadata (a tuple); if None then use
191 ("BGMEAN", "BGVAR"); ignored if stats is false.
192 """
193 netBgImg = background.getImage()
194 if statsKeys is None:
195 statsKeys = ("BGMEAN", "BGVAR")
196 mnkey, varkey = statsKeys
197 meta = exposure.getMetadata()
198 s = afwMath.makeStatistics(netBgImg, afwMath.MEAN | afwMath.VARIANCE)
199 bgmean = s.getValue(afwMath.MEAN)
200 bgvar = s.getValue(afwMath.VARIANCE)
201 meta.addDouble(mnkey, bgmean)
202 meta.addDouble(varkey, bgvar)
203
204 def fitBackground(self, maskedImage, nx=0, ny=0, algorithm=None):
205 """Estimate the background of a masked image
206
207 Parameters
208 ----------
209 maskedImage : `lsst.afw.image.maskedImage`
210 Masked image whose background is to be computed
211 nx : 'int`
212 Number of x bands; if 0 compute from width and `self.config.binSizeX`
213 ny : `int`
214 Number of y bands; if 0 compute from height and `self.config.binSizeY`
215 algorithm : `str`
216 Name of interpolation algorithm; if None use `self.config.algorithm`
217
218 Returns
219 -------
221 A fit background
222
223 Raises
224 ------
225 RuntimeError
226 Raised if lsst.afw.math.makeBackground returns None, an indicator
227 of failure.
228 """
229
230 binSizeX = self.config.binSize if self.config.binSizeX == 0 else self.config.binSizeX
231 binSizeY = self.config.binSize if self.config.binSizeY == 0 else self.config.binSizeY
232
233 if not nx:
234 nx = maskedImage.getWidth()//binSizeX + 1
235 if not ny:
236 ny = maskedImage.getHeight()//binSizeY + 1
237
238 unsubFrame = getDebugFrame(self._display, "unsubtracted")
239 if unsubFrame:
240 unsubDisp = afwDisplay.getDisplay(frame=unsubFrame)
241 unsubDisp.mtv(maskedImage, title="unsubtracted")
242 xPosts = numpy.rint(numpy.linspace(0, maskedImage.getWidth() + 1, num=nx, endpoint=True))
243 yPosts = numpy.rint(numpy.linspace(0, maskedImage.getHeight() + 1, num=ny, endpoint=True))
244 with unsubDisp.Buffering():
245 for (xMin, xMax), (yMin, yMax) in itertools.product(zip(xPosts[:-1], xPosts[1:]),
246 zip(yPosts[:-1], yPosts[1:])):
247 unsubDisp.line([(xMin, yMin), (xMin, yMax), (xMax, yMax), (xMax, yMin), (xMin, yMin)])
248
250 badMask = maskedImage.mask.getPlaneBitMask(self.config.ignoredPixelMask)
251
252 sctrl.setAndMask(badMask)
253 sctrl.setNanSafe(self.config.isNanSafe)
254
255 self.log.debug("Ignoring mask planes: %s", ", ".join(self.config.ignoredPixelMask))
256 if (maskedImage.mask.getArray() & badMask).all():
257 raise pipeBase.TaskError("All pixels masked. Cannot estimate background")
258
259 if algorithm is None:
260 algorithm = self.config.algorithm
261
262 # TODO: DM-22814. This call to a deprecated BackgroundControl constructor
263 # is necessary to support the algorithm parameter; it # should be replaced with
264 #
265 # afwMath.BackgroundControl(nx, ny, sctrl, self.config.statisticsProperty)
266 #
267 # when algorithm has been deprecated and removed.
268 with suppress_deprecations():
269 bctrl = afwMath.BackgroundControl(algorithm, nx, ny,
270 self.config.undersampleStyle, sctrl,
271 self.config.statisticsProperty)
272
273 # TODO: The following check should really be done within lsst.afw.math.
274 # With the current code structure, it would need to be accounted for in the doGetImage()
275 # function in BackgroundMI.cc (which currently only checks against the interpolation settings,
276 # which is not appropriate when useApprox=True)
277 # and/or the makeApproximate() function in afw/Approximate.cc.
278 # See ticket DM-2920: "Clean up code in afw for Approximate background
279 # estimation" (which includes a note to remove the following and the
280 # similar checks in pipe_tasks/matchBackgrounds.py once implemented)
281 #
282 # Check that config setting of approxOrder/binSize make sense
283 # (i.e. ngrid (= shortDimension/binSize) > approxOrderX) and perform
284 # appropriate undersampleStlye behavior.
285 if self.config.useApprox:
286 if self.config.approxOrderY not in (self.config.approxOrderX, -1):
287 raise ValueError("Error: approxOrderY not in (approxOrderX, -1)")
288 order = self.config.approxOrderX
289 minNumberGridPoints = order + 1
290 if min(nx, ny) <= order:
291 self.log.warning("Too few points in grid to constrain fit: min(nx, ny) < approxOrder) "
292 "[min(%d, %d) < %d]", nx, ny, order)
293 if self.config.undersampleStyle == "THROW_EXCEPTION":
294 raise ValueError("Too few points in grid (%d, %d) for order (%d) and binSize (%d, %d)" %
295 (nx, ny, order, binSizeX, binSizeY))
296 elif self.config.undersampleStyle == "REDUCE_INTERP_ORDER":
297 if order < 1:
298 raise ValueError("Cannot reduce approxOrder below 0. "
299 "Try using undersampleStyle = \"INCREASE_NXNYSAMPLE\" instead?")
300 order = min(nx, ny) - 1
301 self.log.warning("Reducing approxOrder to %d", order)
302 elif self.config.undersampleStyle == "INCREASE_NXNYSAMPLE":
303 # Reduce bin size to the largest acceptable square bins
304 newBinSize = min(maskedImage.getWidth(), maskedImage.getHeight())//(minNumberGridPoints-1)
305 if newBinSize < 1:
306 raise ValueError("Binsize must be greater than 0")
307 newNx = maskedImage.getWidth()//newBinSize + 1
308 newNy = maskedImage.getHeight()//newBinSize + 1
309 bctrl.setNxSample(newNx)
310 bctrl.setNySample(newNy)
311 self.log.warning("Decreasing binSize from (%d, %d) to %d for a grid of (%d, %d)",
312 binSizeX, binSizeY, newBinSize, newNx, newNy)
313
314 actrl = afwMath.ApproximateControl(afwMath.ApproximateControl.CHEBYSHEV, order, order,
315 self.config.weighting)
316 bctrl.setApproximateControl(actrl)
317
318 bg = afwMath.makeBackground(maskedImage, bctrl)
319 if bg is None:
320 raise RuntimeError("lsst.afw.math.makeBackground failed to fit a background model")
321 return bg
int min
A class to contain the data, WCS, and other information needed to describe an image of the sky.
Definition: Exposure.h:72
Represent a 2-dimensional array of bitmask pixels.
Definition: Mask.h:77
Control how to make an approximation.
Definition: Approximate.h:48
Pass parameters to a Background object.
Definition: Background.h:56
A virtual base class to evaluate image background levels.
Definition: Background.h:235
Pass parameters to a Statistics object.
Definition: Statistics.h:92
def fitBackground(self, maskedImage, nx=0, ny=0, algorithm=None)
def run(self, exposure, background=None, stats=True, statsKeys=None)
def _addStats(self, exposure, background, statsKeys=None)
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
std::shared_ptr< Background > makeBackground(ImageT const &img, BackgroundControl const &bgCtrl)
A convenience function that uses function overloading to make the correct type of Background.
Definition: Background.h:526
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
bool all(CoordinateExpr< N > const &expr) noexcept
Return true if all elements are true.
def getDebugFrame(debugDisplay, name)
Definition: lsstDebug.py:95