Coverage for python / lsst / meas / algorithms / subtractBackground.py: 21%
147 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-18 08:51 +0000
1#
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__ = (
24 "SubtractBackgroundConfig",
25 "SubtractBackgroundTask",
26 "backgroundFlatContext",
27 "TooManyMaskedPixelsError",
28)
30import itertools
31import math
33from contextlib import contextmanager
34import numpy
35from scipy.ndimage import median_filter
37from lsstDebug import getDebugFrame
38from lsst.utils import suppress_deprecations
40import lsst.afw.display as afwDisplay
41import lsst.afw.image as afwImage
42import lsst.afw.math as afwMath
43import lsst.pex.config as pexConfig
44import lsst.pipe.base as pipeBase
47@contextmanager
48def backgroundFlatContext(maskedImage, doApply, backgroundToPhotometricRatio=None):
49 """Context manager to convert from photometric-flattened to background-
50 flattened image.
52 Parameters
53 ----------
54 maskedImage : `lsst.afw.image.MaskedImage`
55 Masked image (image + mask + variance) to convert from a
56 photometrically flat image to an image suitable for background
57 subtraction.
58 doApply : `bool`
59 Apply the conversion? If False, this context manager will not
60 do anything.
61 backgroundToPhotometricRatio : `lsst.afw.image.Image`, optional
62 Image to multiply a photometrically-flattened image by to obtain a
63 background-flattened image.
64 Only used if ``doApply`` is ``True``.
66 Yields
67 ------
68 maskedImage : `lsst.afw.image.MaskedImage`
69 Masked image converted into an image suitable for background
70 subtraction.
72 Raises
73 ------
74 RuntimeError if doApply is True and no ratio is supplied.
75 ValueError if the ratio is not an `lsst.afw.image.Image`.
76 """
77 if doApply:
78 if backgroundToPhotometricRatio is None:
79 raise RuntimeError("backgroundFlatContext called with doApply=True, "
80 "but without a backgroundToPhotometricRatio")
81 if not isinstance(backgroundToPhotometricRatio, afwImage.Image):
82 raise ValueError("The backgroundToPhotometricRatio must be an lsst.afw.image.Image")
84 maskedImage *= backgroundToPhotometricRatio
86 try:
87 yield maskedImage
88 finally:
89 if doApply:
90 maskedImage /= backgroundToPhotometricRatio
93class TooManyMaskedPixelsError(pipeBase.AlgorithmError):
94 """Raised when all pixels in the image are masked and no background
95 can be estimated.
96 """
98 @property
99 def metadata(self) -> dict:
100 """There is no metadata associated with this error.
101 """
102 return {}
105class SubtractBackgroundConfig(pexConfig.Config):
106 """Config for SubtractBackgroundTask
108 Many of these fields match fields in `lsst.afw.math.BackgroundControl`,
109 the control class for `lsst.afw.math.makeBackground`
110 """
111 statisticsProperty = pexConfig.ChoiceField(
112 doc="type of statistic to use for grid points",
113 dtype=str, default="MEANCLIP",
114 allowed={
115 "MEANCLIP": "clipped mean",
116 "MEAN": "unclipped mean",
117 "MEDIAN": "median",
118 }
119 )
120 undersampleStyle = pexConfig.ChoiceField(
121 doc="behaviour if there are too few points in grid for requested interpolation style",
122 dtype=str, default="REDUCE_INTERP_ORDER",
123 allowed={
124 "THROW_EXCEPTION": "throw an exception if there are too few points",
125 "REDUCE_INTERP_ORDER": "use an interpolation style with a lower order.",
126 "INCREASE_NXNYSAMPLE": "Increase the number of samples used to make the interpolation grid.",
127 },
128 )
129 binSize = pexConfig.RangeField(
130 doc="how large a region of the sky should be used for each background point",
131 dtype=int, default=128, min=1,
132 )
133 binSizeX = pexConfig.RangeField(
134 doc=("Sky region size to be used for each background point in X direction. "
135 "If 0, the binSize config is used."),
136 dtype=int, default=0, min=0,
137 )
138 binSizeY = pexConfig.RangeField(
139 doc=("Sky region size to be used for each background point in Y direction. "
140 "If 0, the binSize config is used."),
141 dtype=int, default=0, min=0,
142 )
143 algorithm = pexConfig.ChoiceField(
144 doc="how to interpolate the background values. This maps to an enum; see afw::math::Background",
145 dtype=str, default="AKIMA_SPLINE", optional=True,
146 allowed={
147 "CONSTANT": "Use a single constant value",
148 "LINEAR": "Use linear interpolation",
149 "NATURAL_SPLINE": "cubic spline with zero second derivative at endpoints",
150 "AKIMA_SPLINE": "higher-level nonlinear spline that is more robust to outliers",
151 "NONE": "No background estimation is to be attempted",
152 },
153 )
154 ignoredPixelMask = pexConfig.ListField(
155 doc="Names of mask planes to ignore while estimating the background",
156 dtype=str, default=["BAD", "EDGE", "DETECTED", "DETECTED_NEGATIVE", "NO_DATA", ],
157 )
158 isNanSafe = pexConfig.Field(
159 doc="Ignore NaNs when estimating the background",
160 dtype=bool, default=False,
161 )
163 useApprox = pexConfig.Field(
164 doc="Use Approximate (Chebyshev) to model background.",
165 dtype=bool, default=True,
166 )
167 approxOrderX = pexConfig.Field(
168 doc="Approximation order in X for background Chebyshev (valid only with useApprox=True)",
169 dtype=int, default=6,
170 )
171 # Note: Currently X- and Y-orders must be equal due to a limitation in math::Chebyshev1Function2
172 # The following is being added so that the weighting attribute can also be configurable for the
173 # call to afwMath.ApproximateControl
174 approxOrderY = pexConfig.Field(
175 doc="Approximation order in Y for background Chebyshev (valid only with useApprox=True)",
176 dtype=int, default=-1,
177 )
178 weighting = pexConfig.Field(
179 doc="Use inverse variance weighting in calculation (valid only with useApprox=True)",
180 dtype=bool, default=True,
181 )
182 doApplyFlatBackgroundRatio = pexConfig.Field(
183 doc="Convert from a photometrically flat image to one suitable to background subtraction? "
184 "If True, then a backgroundToPhotometricRatio must be supplied to the task run method.",
185 dtype=bool,
186 default=False,
187 )
188 doFilterSuperPixels = pexConfig.Field(
189 doc="Remove outliers from the binned image.",
190 dtype=bool,
191 default=False,
192 )
193 superPixelFilterSize = pexConfig.Field(
194 doc="Size of the median filter to use to remove outliers from the binned image.",
195 dtype=int,
196 default=3,
197 )
200class SubtractBackgroundTask(pipeBase.Task):
201 """Subtract the background from an exposure
202 """
203 ConfigClass = SubtractBackgroundConfig
204 _DefaultName = "subtractBackground"
206 def __init__(self, config=None, *, name=None, parentTask=None, log=None):
207 super().__init__(config, name=name, parentTask=parentTask, log=log)
208 self.binSizeX = self.config.binSize if self.config.binSizeX == 0 else self.config.binSizeX
209 self.binSizeY = self.config.binSize if self.config.binSizeY == 0 else self.config.binSizeY
211 def run(self, exposure, background=None, stats=True, statsKeys=None, backgroundToPhotometricRatio=None):
212 """Fit and subtract the background of an exposure.
214 Parameters
215 ----------
216 exposure : `lsst.afw.image.Exposure`
217 Exposure whose background is to be subtracted.
218 background : `lsst.afw.math.BackgroundList`
219 Initial background model already subtracted. May be None if no background
220 has been subtracted.
221 stats : `bool`
222 If True then measure the mean and variance of the full background model and
223 record the results in the exposure's metadata.
224 statsKeys : `tuple`
225 Key names used to store the mean and variance of the background in the
226 exposure's metadata (another tuple); if None then use ("BGMEAN", "BGVAR");
227 ignored if stats is false.
228 backgroundToPhotometricRatio : `lsst.afw.image.Image`, optional
229 Image to multiply a photometrically-flattened image by to obtain a
230 background-flattened image.
231 Only used if config.doApplyFlatBackgroundRatio = True.
233 Returns
234 -------
235 background : `lsst.afw.math.BackgroundList`
236 Full background model (initial model with changes), contained in an
237 `lsst.pipe.base.Struct`.
238 """
239 if background is None:
240 background = afwMath.BackgroundList()
242 maskedImage = exposure.maskedImage
244 with backgroundFlatContext(
245 maskedImage,
246 self.config.doApplyFlatBackgroundRatio,
247 backgroundToPhotometricRatio=backgroundToPhotometricRatio,
248 ):
249 fitBg = self.fitBackground(maskedImage)
250 if self.config.doFilterSuperPixels:
251 filterSuperPixels(maskedImage.getBBox(), fitBg,
252 superPixelFilterSize=self.config.superPixelFilterSize)
254 maskedImage -= fitBg.getImageF(self.config.algorithm, self.config.undersampleStyle)
256 actrl = fitBg.getBackgroundControl().getApproximateControl()
257 background.append((fitBg, getattr(afwMath.Interpolate, self.config.algorithm),
258 fitBg.getAsUsedUndersampleStyle(), actrl.getStyle(),
259 actrl.getOrderX(), actrl.getOrderY(), actrl.getWeighting()))
261 if stats:
262 self._addStats(exposure, background, statsKeys=statsKeys)
264 subFrame = getDebugFrame(self._display, "subtracted")
265 if subFrame:
266 subDisp = afwDisplay.getDisplay(frame=subFrame)
267 subDisp.mtv(exposure, title="subtracted")
269 bgFrame = getDebugFrame(self._display, "background")
270 if bgFrame:
271 bgDisp = afwDisplay.getDisplay(frame=bgFrame)
272 bgImage = background.getImage()
273 bgDisp.mtv(bgImage, title="background")
275 return pipeBase.Struct(
276 background=background,
277 )
279 def _addStats(self, exposure, background, statsKeys=None):
280 """Add statistics about the background to the exposure's metadata
282 Parameters
283 ----------
284 exposure : `lsst.afw.image.Exposure`
285 Exposure whose background was subtracted.
286 background : `lsst.afw.math.BackgroundList`
287 Background model
288 statsKeys : `tuple`
289 Key names used to store the mean and variance of the background in
290 the exposure's metadata (a tuple); if None then use
291 ("BGMEAN", "BGVAR"); ignored if stats is false.
292 """
293 netBgImg = background.getImage()
294 if statsKeys is None:
295 statsKeys = ("BGMEAN", "BGVAR")
296 mnkey, varkey = statsKeys
297 meta = exposure.getMetadata()
298 s = afwMath.makeStatistics(netBgImg, afwMath.MEAN | afwMath.VARIANCE)
299 bgmean = s.getValue(afwMath.MEAN)
300 bgvar = s.getValue(afwMath.VARIANCE)
301 meta.addDouble(mnkey, bgmean)
302 meta.addDouble(varkey, bgvar)
304 def fitBackground(self, maskedImage, nx=0, ny=0, algorithm=None):
305 """Estimate the background of a masked image
307 Parameters
308 ----------
309 maskedImage : `lsst.afw.image.maskedImage`
310 Masked image whose background is to be computed
311 nx : 'int`
312 Number of x bands; if 0 compute from width and `self.config.binSizeX`
313 ny : `int`
314 Number of y bands; if 0 compute from height and `self.config.binSizeY`
315 algorithm : `str`
316 Name of interpolation algorithm; if None use `self.config.algorithm`
318 Returns
319 -------
320 bg : `lsst.afw.math.Background`
321 A fit background
323 Raises
324 ------
325 RuntimeError
326 Raised if lsst.afw.math.makeBackground returns None, an indicator
327 of failure.
328 """
330 if not nx:
331 nx = math.ceil(maskedImage.getWidth() / self.binSizeX)
332 if not ny:
333 ny = math.ceil(maskedImage.getHeight() / self.binSizeY)
335 unsubFrame = getDebugFrame(self._display, "unsubtracted")
336 if unsubFrame:
337 unsubDisp = afwDisplay.getDisplay(frame=unsubFrame)
338 unsubDisp.mtv(maskedImage, title="unsubtracted")
339 xPosts = numpy.rint(numpy.linspace(0, maskedImage.getWidth() + 1, num=nx, endpoint=True))
340 yPosts = numpy.rint(numpy.linspace(0, maskedImage.getHeight() + 1, num=ny, endpoint=True))
341 with unsubDisp.Buffering():
342 for (xMin, xMax), (yMin, yMax) in itertools.product(zip(xPosts[:-1], xPosts[1:]),
343 zip(yPosts[:-1], yPosts[1:])):
344 unsubDisp.line([(xMin, yMin), (xMin, yMax), (xMax, yMax), (xMax, yMin), (xMin, yMin)])
346 sctrl = afwMath.StatisticsControl()
347 badMask = maskedImage.mask.getPlaneBitMask(self.config.ignoredPixelMask)
349 sctrl.setAndMask(badMask)
350 sctrl.setNanSafe(self.config.isNanSafe)
352 self.log.debug("Ignoring mask planes: %s", ", ".join(self.config.ignoredPixelMask))
353 if (maskedImage.mask.getArray() & badMask).all():
354 raise TooManyMaskedPixelsError("All pixels masked. Cannot estimate background.")
356 if algorithm is None:
357 algorithm = self.config.algorithm
359 # TODO: DM-22814. This call to a deprecated BackgroundControl constructor
360 # is necessary to support the algorithm parameter; it # should be replaced with
361 #
362 # afwMath.BackgroundControl(nx, ny, sctrl, self.config.statisticsProperty)
363 #
364 # when algorithm has been deprecated and removed.
365 with suppress_deprecations():
366 bctrl = afwMath.BackgroundControl(algorithm, nx, ny,
367 self.config.undersampleStyle, sctrl,
368 self.config.statisticsProperty)
370 # TODO: The following check should really be done within lsst.afw.math.
371 # With the current code structure, it would need to be accounted for in the doGetImage()
372 # function in BackgroundMI.cc (which currently only checks against the interpolation settings,
373 # which is not appropriate when useApprox=True)
374 # and/or the makeApproximate() function in afw/Approximate.cc.
375 # See ticket DM-2920: "Clean up code in afw for Approximate background
376 # estimation" (which includes a note to remove the following and the
377 # similar checks in pipe_tasks/matchBackgrounds.py once implemented)
378 #
379 # Check that config setting of approxOrder/binSize make sense
380 # (i.e. ngrid (= shortDimension/binSize) > approxOrderX) and perform
381 # appropriate undersampleStlye behavior.
382 if self.config.useApprox:
383 if self.config.approxOrderY not in (self.config.approxOrderX, -1):
384 raise ValueError("Error: approxOrderY not in (approxOrderX, -1)")
385 order = self.config.approxOrderX
386 minNumberGridPoints = order + 1
387 if min(nx, ny) <= order:
388 self.log.warning("Too few points in grid to constrain fit: min(nx, ny) < approxOrder) "
389 "[min(%d, %d) < %d]", nx, ny, order)
390 if self.config.undersampleStyle == "THROW_EXCEPTION":
391 raise ValueError("Too few points in grid (%d, %d) for order (%d) and binSize (%d, %d)" %
392 (nx, ny, order, self.binSizeX, self.binSizeY))
393 elif self.config.undersampleStyle == "REDUCE_INTERP_ORDER":
394 if order < 1:
395 raise ValueError("Cannot reduce approxOrder below 0. "
396 "Try using undersampleStyle = \"INCREASE_NXNYSAMPLE\" instead?")
397 order = min(nx, ny) - 1
398 self.log.warning("Reducing approxOrder to %d", order)
399 elif self.config.undersampleStyle == "INCREASE_NXNYSAMPLE":
400 # Reduce bin size to the largest acceptable square bins
401 newBinSize = min(maskedImage.getWidth(), maskedImage.getHeight())//(minNumberGridPoints-1)
402 if newBinSize < 1:
403 raise ValueError("Binsize must be greater than 0")
404 newNx = maskedImage.getWidth()//newBinSize + 1
405 newNy = maskedImage.getHeight()//newBinSize + 1
406 bctrl.setNxSample(newNx)
407 bctrl.setNySample(newNy)
408 self.log.warning("Decreasing binSize from (%d, %d) to %d for a grid of (%d, %d)",
409 self.binSizeX, self.binSizeY, newBinSize, newNx, newNy)
411 actrl = afwMath.ApproximateControl(afwMath.ApproximateControl.CHEBYSHEV, order, order,
412 self.config.weighting)
413 bctrl.setApproximateControl(actrl)
415 bg = afwMath.makeBackground(maskedImage, bctrl)
416 if bg is None:
417 raise RuntimeError("lsst.afw.math.makeBackground failed to fit a background model")
418 return bg
421def filterSuperPixels(bbox, background, superPixelFilterSize=3):
422 """Remove outliers from the binned background model.
424 Parameters
425 ----------
426 bbox : `lsst.geom.Box2I`
427 Bounding box of the original image.
428 background : `lsst.afw.math.BackgroundMI`
429 Fit and binned background image, which will be modified in place.
430 superPixelFilterSize : `int`, optional
431 Size of the median filter to use, in pixels.
432 """
433 statsImg = background.getStatsImage()
434 # scipy's median_filter can't handle NaN values
435 bad = numpy.isnan(statsImg.image.array)
436 if numpy.count_nonzero(bad) > 0:
437 statsImg.image.array[bad] = numpy.nanmedian(statsImg.image.array)
438 statsImg.image.array = median_filter(statsImg.image.array, mode="reflect", size=superPixelFilterSize)
439 background = afwMath.BackgroundMI(bbox, statsImg)