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
dynamicDetection.py
Go to the documentation of this file.
1 
2 __all__ = ["DynamicDetectionConfig", "DynamicDetectionTask"]
3 
4 import numpy as np
5 
6 from lsst.pex.config import Field, ConfigurableField
7 from lsst.pipe.base import Struct
8 
9 from .detection import SourceDetectionConfig, SourceDetectionTask
10 from .skyObjects import SkyObjectsTask
11 
12 from lsst.afw.detection import FootprintSet
13 from lsst.afw.table import SourceCatalog, SourceTable
14 from lsst.meas.base import ForcedMeasurementTask
15 
16 import lsst.afw.image
17 import lsst.afw.math
18 
19 
21  """Configuration for DynamicDetectionTask
22  """
23  prelimThresholdFactor = Field(dtype=float, default=0.5,
24  doc="Fraction of the threshold to use for first pass (to find sky objects)")
25  skyObjects = ConfigurableField(target=SkyObjectsTask, doc="Generate sky objects")
26  doBackgroundTweak = Field(dtype=bool, default=True,
27  doc="Tweak background level so median PSF flux of sky objects is zero?")
28  minNumSources = Field(dtype=int, default=10,
29  doc="Minimum number of sky sources in statistical sample; "
30  "if below this number, we refuse to modify the threshold.")
31 
32  def setDefaults(self):
33  SourceDetectionConfig.setDefaults(self)
34  self.skyObjectsskyObjects.nSources = 1000 # For good statistics
35 
36 
38  """Detection of sources on an image with a dynamic threshold
39 
40  We first detect sources using a lower threshold than normal (see config
41  parameter ``prelimThresholdFactor``) in order to identify good sky regions
42  (configurable ``skyObjects``). Then we perform forced PSF photometry on
43  those sky regions. Using those PSF flux measurements and estimated errors,
44  we set the threshold so that the stdev of the measurements matches the
45  median estimated error.
46 
47  Besides the usual initialisation of configurables, we also set up
48  the forced measurement which is deliberately not represented in
49  this Task's configuration parameters because we're using it as
50  part of the algorithm and we don't want to allow it to be modified.
51  """
52  ConfigClass = DynamicDetectionConfig
53  _DefaultName = "dynamicDetection"
54 
55  def __init__(self, *args, **kwargs):
56 
57  SourceDetectionTask.__init__(self, *args, **kwargs)
58  self.makeSubtask("skyObjects")
59 
60  # Set up forced measurement.
61  config = ForcedMeasurementTask.ConfigClass()
62  config.plugins.names = ['base_TransformedCentroid', 'base_PsfFlux', 'base_LocalBackground']
63  # We'll need the "centroid" and "psfFlux" slots
64  for slot in ("shape", "psfShape", "apFlux", "modelFlux", "gaussianFlux", "calibFlux"):
65  setattr(config.slots, slot, None)
66  config.copyColumns = {}
67  self.skySchemaskySchema = SourceTable.makeMinimalSchema()
68  self.skyMeasurementskyMeasurement = ForcedMeasurementTask(config=config, name="skyMeasurement", parentTask=self,
69  refSchema=self.skySchemaskySchema)
70 
71  def calculateThreshold(self, exposure, seed, sigma=None):
72  """Calculate new threshold
73 
74  This is the main functional addition to the vanilla
75  `SourceDetectionTask`.
76 
77  We identify sky objects and perform forced PSF photometry on
78  them. Using those PSF flux measurements and estimated errors,
79  we set the threshold so that the stdev of the measurements
80  matches the median estimated error.
81 
82  Parameters
83  ----------
84  exposure : `lsst.afw.image.Exposure`
85  Exposure on which we're detecting sources.
86  seed : `int`
87  RNG seed to use for finding sky objects.
88  sigma : `float`, optional
89  Gaussian sigma of smoothing kernel; if not provided,
90  will be deduced from the exposure's PSF.
91 
92  Returns
93  -------
94  result : `lsst.pipe.base.Struct`
95  Result struct with components:
96 
97  - ``multiplicative``: multiplicative factor to be applied to the
98  configured detection threshold (`float`).
99  - ``additive``: additive factor to be applied to the background
100  level (`float`).
101  """
102  # Make a catalog of sky objects
103  fp = self.skyObjects.run(exposure.maskedImage.mask, seed)
104  skyFootprints = FootprintSet(exposure.getBBox())
105  skyFootprints.setFootprints(fp)
106  table = SourceTable.make(self.skyMeasurementskyMeasurement.schema)
107  catalog = SourceCatalog(table)
108  catalog.reserve(len(skyFootprints.getFootprints()))
109  skyFootprints.makeSources(catalog)
110  key = catalog.getCentroidSlot().getMeasKey()
111  for source in catalog:
112  peaks = source.getFootprint().getPeaks()
113  assert len(peaks) == 1
114  source.set(key, peaks[0].getF())
115  source.updateCoord(exposure.getWcs())
116 
117  # Forced photometry on sky objects
118  self.skyMeasurementskyMeasurement.run(catalog, exposure, catalog, exposure.getWcs())
119 
120  # Calculate new threshold
121  fluxes = catalog["base_PsfFlux_instFlux"]
122  area = catalog["base_PsfFlux_area"]
123  bg = catalog["base_LocalBackground_instFlux"]
124 
125  good = (~catalog["base_PsfFlux_flag"] & ~catalog["base_LocalBackground_flag"]
126  & np.isfinite(fluxes) & np.isfinite(area) & np.isfinite(bg))
127 
128  if good.sum() < self.config.minNumSources:
129  self.log.warn("Insufficient good flux measurements (%d < %d) for dynamic threshold calculation",
130  good.sum(), self.config.minNumSources)
131  return Struct(multiplicative=1.0, additive=0.0)
132 
133  bgMedian = np.median((fluxes/area)[good])
134 
135  lq, uq = np.percentile((fluxes - bg*area)[good], [25.0, 75.0])
136  stdevMeas = 0.741*(uq - lq)
137  medianError = np.median(catalog["base_PsfFlux_instFluxErr"][good])
138  return Struct(multiplicative=medianError/stdevMeas, additive=bgMedian)
139 
140  def detectFootprints(self, exposure, doSmooth=True, sigma=None, clearMask=True, expId=None):
141  """Detect footprints with a dynamic threshold
142 
143  This varies from the vanilla ``detectFootprints`` method because we
144  do detection twice: one with a low threshold so that we can find
145  sky uncontaminated by objects, then one more with the new calculated
146  threshold.
147 
148  Parameters
149  ----------
150  exposure : `lsst.afw.image.Exposure`
151  Exposure to process; DETECTED{,_NEGATIVE} mask plane will be
152  set in-place.
153  doSmooth : `bool`, optional
154  If True, smooth the image before detection using a Gaussian
155  of width ``sigma``.
156  sigma : `float`, optional
157  Gaussian Sigma of PSF (pixels); used for smoothing and to grow
158  detections; if `None` then measure the sigma of the PSF of the
159  ``exposure``.
160  clearMask : `bool`, optional
161  Clear both DETECTED and DETECTED_NEGATIVE planes before running
162  detection.
163  expId : `int`, optional
164  Exposure identifier, used as a seed for the random number
165  generator. If absent, the seed will be the sum of the image.
166 
167  Return Struct contents
168  ----------------------
169  positive : `lsst.afw.detection.FootprintSet`
170  Positive polarity footprints (may be `None`)
171  negative : `lsst.afw.detection.FootprintSet`
172  Negative polarity footprints (may be `None`)
173  numPos : `int`
174  Number of footprints in positive or 0 if detection polarity was
175  negative.
176  numNeg : `int`
177  Number of footprints in negative or 0 if detection polarity was
178  positive.
179  background : `lsst.afw.math.BackgroundList`
180  Re-estimated background. `None` if
181  ``reEstimateBackground==False``.
182  factor : `float`
183  Multiplication factor applied to the configured detection
184  threshold.
185  prelim : `lsst.pipe.base.Struct`
186  Results from preliminary detection pass.
187  """
188  maskedImage = exposure.maskedImage
189 
190  if clearMask:
191  self.clearMaskclearMask(maskedImage.mask)
192  else:
193  oldDetected = maskedImage.mask.array & maskedImage.mask.getPlaneBitMask(["DETECTED",
194  "DETECTED_NEGATIVE"])
195 
196  with self.tempWideBackgroundContexttempWideBackgroundContext(exposure):
197  # Could potentially smooth with a wider kernel than the PSF in order to better pick up the
198  # wings of stars and galaxies, but for now sticking with the PSF as that's more simple.
199  psf = self.getPsfgetPsf(exposure, sigma=sigma)
200  convolveResults = self.convolveImageconvolveImage(maskedImage, psf, doSmooth=doSmooth)
201  middle = convolveResults.middle
202  sigma = convolveResults.sigma
203  prelim = self.applyThresholdapplyThreshold(middle, maskedImage.getBBox(), self.config.prelimThresholdFactor)
204  self.finalizeFootprintsfinalizeFootprints(maskedImage.mask, prelim, sigma, self.config.prelimThresholdFactor)
205 
206  # Calculate the proper threshold
207  # seed needs to fit in a C++ 'int' so pybind doesn't choke on it
208  seed = (expId if expId is not None else int(maskedImage.image.array.sum())) % (2**31 - 1)
209  threshResults = self.calculateThresholdcalculateThreshold(exposure, seed, sigma=sigma)
210  factor = threshResults.multiplicative
211  self.log.info("Modifying configured detection threshold by factor %f to %f",
212  factor, factor*self.config.thresholdValue)
213 
214  # Blow away preliminary (low threshold) detection mask
215  self.clearMaskclearMask(maskedImage.mask)
216  if not clearMask:
217  maskedImage.mask.array |= oldDetected
218 
219  # Rinse and repeat thresholding with new calculated threshold
220  results = self.applyThresholdapplyThreshold(middle, maskedImage.getBBox(), factor)
221  results.prelim = prelim
222  results.background = lsst.afw.math.BackgroundList()
223  if self.config.doTempLocalBackground:
224  self.applyTempLocalBackgroundapplyTempLocalBackground(exposure, middle, results)
225  self.finalizeFootprintsfinalizeFootprints(maskedImage.mask, results, sigma, factor)
226 
227  self.clearUnwantedResultsclearUnwantedResults(maskedImage.mask, results)
228 
229  if self.config.reEstimateBackground:
230  self.reEstimateBackgroundreEstimateBackground(maskedImage, results.background)
231 
232  self.displaydisplay(exposure, results, middle)
233 
234  if self.config.doBackgroundTweak:
235  # Re-do the background tweak after any temporary backgrounds have been restored
236  #
237  # But we want to keep any large-scale background (e.g., scattered light from bright stars)
238  # from being selected for sky objects in the calculation, so do another detection pass without
239  # either the local or wide temporary background subtraction; the DETECTED pixels will mark
240  # the area to ignore.
241  originalMask = maskedImage.mask.array.copy()
242  try:
243  self.clearMaskclearMask(exposure.mask)
244  convolveResults = self.convolveImageconvolveImage(maskedImage, psf, doSmooth=doSmooth)
245  tweakDetResults = self.applyThresholdapplyThreshold(convolveResults.middle, maskedImage.getBBox(), factor)
246  self.finalizeFootprintsfinalizeFootprints(maskedImage.mask, tweakDetResults, sigma, factor)
247  bgLevel = self.calculateThresholdcalculateThreshold(exposure, seed, sigma=sigma).additive
248  finally:
249  maskedImage.mask.array[:] = originalMask
250  self.tweakBackgroundtweakBackground(exposure, bgLevel, results.background)
251 
252  return results
253 
254  def tweakBackground(self, exposure, bgLevel, bgList=None):
255  """Modify the background by a constant value
256 
257  Parameters
258  ----------
259  exposure : `lsst.afw.image.Exposure`
260  Exposure for which to tweak background.
261  bgLevel : `float`
262  Background level to remove
263  bgList : `lsst.afw.math.BackgroundList`, optional
264  List of backgrounds to append to.
265 
266  Returns
267  -------
268  bg : `lsst.afw.math.BackgroundMI`
269  Constant background model.
270  """
271  self.log.info("Tweaking background by %f to match sky photometry", bgLevel)
272  exposure.image -= bgLevel
273  bgStats = lsst.afw.image.MaskedImageF(1, 1)
274  bgStats.set(bgLevel, 0, bgLevel)
275  bg = lsst.afw.math.BackgroundMI(exposure.getBBox(), bgStats)
276  bgData = (bg, lsst.afw.math.Interpolate.LINEAR, lsst.afw.math.REDUCE_INTERP_ORDER,
277  lsst.afw.math.ApproximateControl.UNKNOWN, 0, 0, False)
278  if bgList is not None:
279  bgList.append(bgData)
280  return bg
A set of Footprints, associated with a MaskedImage.
Definition: FootprintSet.h:53
A class to evaluate image background levels.
Definition: Background.h:434
def getPsf(self, exposure, sigma=None)
Definition: detection.py:387
def run(self, table, exposure, doSmooth=True, sigma=None, clearMask=True, expId=None)
Definition: detection.py:204
def convolveImage(self, maskedImage, psf, doSmooth=True)
Definition: detection.py:414
def applyTempLocalBackground(self, exposure, middle, results)
Definition: detection.py:318
def reEstimateBackground(self, maskedImage, backgrounds)
Definition: detection.py:586
def finalizeFootprints(self, mask, results, sigma, factor=1.0)
Definition: detection.py:524
def display(self, exposure, results, convolvedImage=None)
Definition: detection.py:264
def applyThreshold(self, middle, bbox, factor=1.0)
Definition: detection.py:475
def tweakBackground(self, exposure, bgLevel, bgList=None)
def detectFootprints(self, exposure, doSmooth=True, sigma=None, clearMask=True, expId=None)
def calculateThreshold(self, exposure, seed, sigma=None)
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.