LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
subtractBackground.py
Go to the documentation of this file.
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__ = ("SubtractBackgroundConfig", "SubtractBackgroundTask")
24 
25 import itertools
26 
27 import numpy
28 
29 from lsstDebug import getDebugFrame
30 from lsst.utils import suppress_deprecations
31 
32 import lsst.afw.display as afwDisplay
33 import lsst.afw.image as afwImage
34 import lsst.afw.math as afwMath
35 import lsst.pex.config as pexConfig
36 import lsst.pipe.base as pipeBase
37 
38 
39 class 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 
119 class 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  -------
220  bg : `lsst.afw.math.Background`
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 
249  sctrl = afwMath.StatisticsControl()
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
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
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.
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
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
bool all(CoordinateExpr< N > const &expr) noexcept
Return true if all elements are true.
def getDebugFrame(debugDisplay, name)
Definition: lsstDebug.py:95