Coverage for python / lsst / meas / algorithms / dynamicDetection.py: 14%
295 statements
« prev ^ index » next coverage.py v7.12.0, created at 2025-12-12 01:22 -0800
2__all__ = [
3 "DynamicDetectionConfig",
4 "DynamicDetectionTask",
5 "InsufficientSourcesError",
6 "ZeroFootprintError",
7]
9import numpy as np
11from lsst.pex.config import Field, ConfigurableField, FieldValidationError
13from .detection import SourceDetectionConfig, SourceDetectionTask
14from .skyObjects import SkyObjectsTask
15from .subtractBackground import TooManyMaskedPixelsError
17from lsst.afw.detection import FootprintSet
18from lsst.afw.geom import makeCdMatrix, makeSkyWcs, SpanSet
19from lsst.afw.table import SourceCatalog, SourceTable
20from lsst.meas.base import ForcedMeasurementTask
21from lsst.pipe.base import AlgorithmError, Struct
23import lsst.afw.image
24import lsst.afw.math
25import lsst.geom as geom
28class InsufficientSourcesError(AlgorithmError):
29 """Raised if an insufficient number of sky sources are found for
30 dynamic detection.
32 Parameters
33 ----------
34 msg : `str`
35 Error message.
36 nGoodPix : `int`
37 Number of good pixels (i.e. not NO_DATA or BAD).
38 nPix : `int`
39 Total number of pixels.
40 **kwargs : `dict`, optional
41 Additional keyword arguments to initialize the Exception base class.
42 """
43 def __init__(self, msg, nGoodPix, nPix, **kwargs):
44 self.msg = msg
45 self._metadata = kwargs
46 super().__init__(msg, **kwargs)
47 self._metadata["nGoodPix"] = int(nGoodPix)
48 self._metadata["nPix"] = int(nPix)
50 def __str__(self):
51 # Exception doesn't handle **kwargs, so we need a custom str.
52 return f"{self.msg}: {self.metadata}"
54 @property
55 def metadata(self):
56 for key, value in self._metadata.items():
57 if not isinstance(value, (int, float, str)):
58 raise TypeError(f"{key} is of type {type(value)}, but only (int, float, str) are allowed.")
59 return self._metadata
62class ZeroFootprintError(AlgorithmError):
63 """Raised if no footprints are detected in the image.
65 Parameters
66 ----------
67 msg : `str`
68 Error message.
69 **kwargs : `dict`, optional
70 Additional keyword arguments to initialize the Exception base class.
71 """
72 def __init__(self, msg, **kwargs):
73 self.msg = msg
74 self._metadata = kwargs
75 super().__init__(msg, **kwargs)
77 def __str__(self):
78 # Exception doesn't handle **kwargs, so we need a custom str.
79 return f"{self.msg}: {self.metadata}"
81 @property
82 def metadata(self):
83 for key, value in self._metadata.items():
84 if not isinstance(value, (int, float, str)):
85 raise TypeError(f"{key} is of type {type(value)}, but only (int, float, str) are allowed.")
86 return self._metadata
89class DynamicDetectionConfig(SourceDetectionConfig):
90 """Configuration for DynamicDetectionTask
91 """
92 prelimThresholdFactor = Field(dtype=float, default=0.5,
93 doc="Factor by which to multiply the main detection threshold "
94 "(thresholdValue) to use for first pass (to find sky objects).")
95 prelimNegMultiplier = Field(dtype=float, default=2.5,
96 doc="Multiplier for the negative (relative to positive) polarity "
97 "detections threshold to use for first pass (to find sky objects).")
98 skyObjects = ConfigurableField(target=SkyObjectsTask, doc="Generate sky objects.")
99 minGoodPixelFraction = Field(dtype=float, default=0.005,
100 doc="Minimum fraction of 'good' pixels required to be deemed "
101 "worthwhile for an attempt at further processing.")
102 doThresholdScaling = Field(dtype=bool, default=True,
103 doc="Scale the threshold level to get empirically measured S/N requested?")
104 minThresholdScaleFactor = Field(dtype=float, default=0.1, optional=True,
105 doc="Mininum threshold scaling allowed (i.e. it will be set to this "
106 "if the computed value is smaller than it). Set to None for no limit.")
107 maxThresholdScaleFactor = Field(dtype=float, default=10.0, optional=True,
108 doc="Maximum threshold scaling allowed (i.e. it will be set to this "
109 "if the computed value is greater than it). Set to None for no limit.")
110 doBackgroundTweak = Field(dtype=bool, default=True,
111 doc="Tweak background level so median PSF flux of sky objects is zero?")
112 minBackgroundTweak = Field(dtype=float, default=-100.0, optional=True,
113 doc="Mininum background tweak allowed (i.e. it will be set to this "
114 "if the computed value is smaller than it). Set to None for no limit.")
115 maxBackgroundTweak = Field(dtype=float, default=100.0, optional=True,
116 doc="Maximum background tweak allowed (i.e. it will be set to this "
117 "if the computed value is greater than it). Set to None for no limit.")
118 minFractionSources = Field(dtype=float, default=0.02,
119 doc="Minimum fraction of the requested number of sky sources for dynamic "
120 "detection to be considered a success. If the number of good sky sources "
121 "identified falls below this threshold, an InsufficientSourcesError error "
122 "is raised so that this dataId is no longer considered in downstream "
123 "processing.")
124 doBrightPrelimDetection = Field(dtype=bool, default=True,
125 doc="Do initial bright detection pass where footprints are grown "
126 "by brightGrowFactor?")
127 brightMultiplier = Field(dtype=float, default=2000.0,
128 doc="Multiplier to apply to the prelimThresholdFactor for the "
129 "\"bright\" detections stage (want this to be large to only "
130 "detect the brightest sources).")
131 brightNegFactor = Field(dtype=float, default=2.2,
132 doc="Factor by which to multiply the threshold for the negative polatiry "
133 "detections for the \"bright\" detections stage (this needs to be fairly "
134 "low given the nature of the negative polarity detections in the very "
135 "large positive polarity threshold).")
136 brightGrowFactor = Field(dtype=int, default=40,
137 doc="Factor by which to grow the footprints of sources detected in the "
138 "\"bright\" detections stage (want this to be large to mask wings of "
139 "bright sources).")
140 brightMaskFractionMax = Field(dtype=float, default=0.95,
141 doc="Maximum allowed fraction of masked pixes from the \"bright\" "
142 "detection stage (to mask regions unsuitable for sky sourcess). "
143 "If this fraction is exeeded, the detection threshold for this stage "
144 "will be increased by bisectFactor until the fraction of masked "
145 "pixels drops below this threshold.")
146 bisectFactor = Field(dtype=float, default=1.2,
147 doc="Factor by which to increase thresholds in brightMaskFractionMax loop.")
148 allowMaskErode = Field(dtype=bool, default=True,
149 doc="Crowded/large fill-factor fields make it difficult to find "
150 "suitable locations to lay down sky objects. To allow for best effort "
151 "sky source placement, if True, this allows for a slight erosion of "
152 "the detection masks.")
153 maxPeakToFootRatio = Field(dtype=float, default=150.0,
154 doc="Maximum ratio of peak per footprint in the detection mask. "
155 "This is to help prevent single contiguous footprints that nothing "
156 "can be done with (i.e. deblending will be skipped). If the current "
157 "detection plane does not satisfy this constraint, the detection "
158 "threshold is increased iteratively until it is. This behaviour is "
159 "intended to be an effective no-op for most \"typical\" scenes/standard "
160 "quality observations, but can avoid total meltdown in, e.g. very "
161 "crowded regions.")
163 def setDefaults(self):
164 SourceDetectionConfig.setDefaults(self)
165 self.skyObjects.nSources = 1000 # For good statistics
166 for maskStr in ["SAT"]:
167 if maskStr not in self.skyObjects.avoidMask:
168 self.skyObjects.avoidMask.append(maskStr)
170 def validate(self):
171 super().validate()
173 if self.doApplyFlatBackgroundRatio:
174 raise ValueError("DynamicDetectionTask does not support doApplyFlatBackgroundRatio.")
176 if self.doThresholdScaling:
177 if self.minThresholdScaleFactor and self.maxThresholdScaleFactor:
178 if self.minThresholdScaleFactor > self.maxThresholdScaleFactor:
179 msg = "minThresholdScaleFactor must be <= maxThresholdScaleFactor"
180 raise FieldValidationError(
181 DynamicDetectionConfig.doThresholdScaling, self, msg,
182 )
184 if self.doBackgroundTweak:
185 if self.minBackgroundTweak and self.maxBackgroundTweak:
186 if self.minBackgroundTweak > self.maxBackgroundTweak:
187 msg = "minBackgroundTweak must be <= maxBackgroundTweak"
188 raise FieldValidationError(
189 DynamicDetectionConfig.doBackgroundTweak, self, msg,
190 )
193class DynamicDetectionTask(SourceDetectionTask):
194 """Detection of sources on an image with a dynamic threshold
196 We first detect sources using a lower threshold than normal (see config
197 parameter ``prelimThresholdFactor``) in order to identify good sky regions
198 (configurable ``skyObjects``). Then we perform forced PSF photometry on
199 those sky regions. Using those PSF flux measurements and estimated errors,
200 we set the threshold so that the stdev of the measurements matches the
201 median estimated error.
203 Besides the usual initialisation of configurables, we also set up
204 the forced measurement which is deliberately not represented in
205 this Task's configuration parameters because we're using it as
206 part of the algorithm and we don't want to allow it to be modified.
207 """
208 ConfigClass = DynamicDetectionConfig
209 _DefaultName = "dynamicDetection"
211 def __init__(self, *args, **kwargs):
213 SourceDetectionTask.__init__(self, *args, **kwargs)
214 self.makeSubtask("skyObjects")
216 # Set up forced measurement.
217 config = ForcedMeasurementTask.ConfigClass()
218 config.plugins.names = ["base_TransformedCentroid", "base_PsfFlux"]
219 # We'll need the "centroid" and "psfFlux" slots
220 for slot in ("shape", "psfShape", "apFlux", "modelFlux", "gaussianFlux", "calibFlux"):
221 setattr(config.slots, slot, None)
222 config.copyColumns = {}
223 self.skySchema = SourceTable.makeMinimalSchema()
224 self.skyMeasurement = ForcedMeasurementTask(config=config, name="skyMeasurement", parentTask=self,
225 refSchema=self.skySchema)
227 def calculateThreshold(self, exposure, seed, sigma=None, minFractionSourcesFactor=1.0,
228 isBgTweak=False, nPixMaskErode=None, maxMaskErodeIter=10):
229 """Calculate new threshold
231 This is the main functional addition to the vanilla
232 `SourceDetectionTask`.
234 We identify sky objects and perform forced PSF photometry on
235 them. Using those PSF flux measurements and estimated errors,
236 we set the threshold so that the stdev of the measurements
237 matches the median estimated error.
239 Parameters
240 ----------
241 exposure : `lsst.afw.image.Exposure`
242 Exposure on which we're detecting sources.
243 seed : `int`
244 RNG seed to use for finding sky objects.
245 sigma : `float`, optional
246 Gaussian sigma of smoothing kernel; if not provided,
247 will be deduced from the exposure's PSF.
248 minFractionSourcesFactor : `float`
249 Change the fraction of required sky sources from that set in
250 ``self.config.minFractionSources`` by this factor. NOTE: this
251 is intended for use in the background tweak pass (the detection
252 threshold is much lower there, so many more pixels end up marked
253 as DETECTED or DETECTED_NEGATIVE, leaving less room for sky
254 object placement).
255 isBgTweak : `bool`
256 Set to ``True`` for the background tweak pass (for more helpful
257 log messages).
258 nPixMaskErode : `int`, optional
259 Number of pixels by which to erode the detection masks on each
260 iteration of best-effort sky object placement.
261 maxMaskErodeIter : `int`, optional
262 Maximum number of iterations for the detection mask erosion.
264 Returns
265 -------
266 result : `lsst.pipe.base.Struct`
267 Result struct with components:
269 ``multiplicative``
270 Multiplicative factor to be applied to the
271 configured detection threshold (`float`).
272 ``additive``
273 Additive factor to be applied to the background
274 level (`float`).
276 Raises
277 ------
278 InsufficientSourcesError
279 Raised if the number of good sky sources found is less than the
280 minimum fraction
281 (``self.config.minFractionSources``*``minFractionSourcesFactor``)
282 of the number requested (``self.skyObjects.config.nSources``).
283 """
284 wcsIsNone = exposure.getWcs() is None
285 if wcsIsNone: # create a dummy WCS as needed by ForcedMeasurementTask
286 self.log.info("WCS for exposure is None. Setting a dummy WCS for dynamic detection.")
287 exposure.setWcs(makeSkyWcs(crpix=geom.Point2D(0, 0),
288 crval=geom.SpherePoint(0, 0, geom.degrees),
289 cdMatrix=makeCdMatrix(scale=1e-5*geom.degrees)))
290 minNumSources = int(self.config.minFractionSources*self.skyObjects.config.nSources)
291 # Reduce the number of sky sources required if requested, but ensure
292 # a minumum of 3.
293 if minFractionSourcesFactor != 1.0:
294 minNumSources = max(3, int(minNumSources*minFractionSourcesFactor))
295 fp = self.skyObjects.run(exposure.maskedImage.mask, seed)
297 if self.config.allowMaskErode:
298 detectedMaskPlanes = ["DETECTED", "DETECTED_NEGATIVE"]
299 mask = exposure.maskedImage.mask
300 for nIter in range(maxMaskErodeIter):
301 if nIter > 0:
302 fp = self.skyObjects.run(mask, seed)
303 if len(fp) < int(2*minNumSources): # Allow for measurement failures
304 self.log.info("Current number of sky sources is below 2*minimum required "
305 "(%d < %d, allowing for some subsequent measurement failures). "
306 "Allowing erosion of detected mask planes for sky placement "
307 "nIter: %d [of %d max]",
308 len(fp), 2*minNumSources, nIter, maxMaskErodeIter)
309 if nPixMaskErode is None:
310 if len(fp) == 0:
311 nPixMaskErode = 4
312 elif len(fp) < int(0.75*minNumSources):
313 nPixMaskErode = 2
314 else:
315 nPixMaskErode = 1
316 for maskName in detectedMaskPlanes:
317 # Compute the eroded detection mask plane using SpanSet
318 detectedMaskBit = mask.getPlaneBitMask(maskName)
319 detectedMaskSpanSet = SpanSet.fromMask(mask, detectedMaskBit)
320 detectedMaskSpanSet = detectedMaskSpanSet.eroded(nPixMaskErode)
321 # Clear the detected mask plane
322 detectedMask = mask.getMaskPlane(maskName)
323 mask.clearMaskPlane(detectedMask)
324 # Set the mask plane to the eroded one
325 detectedMaskSpanSet.setMask(mask, detectedMaskBit)
326 else:
327 break
329 skyFootprints = FootprintSet(exposure.getBBox())
330 skyFootprints.setFootprints(fp)
331 table = SourceTable.make(self.skyMeasurement.schema)
332 catalog = SourceCatalog(table)
333 catalog.reserve(len(skyFootprints.getFootprints()))
334 skyFootprints.makeSources(catalog)
335 key = catalog.getCentroidSlot().getMeasKey()
336 for source in catalog:
337 peaks = source.getFootprint().getPeaks()
338 assert len(peaks) == 1
339 source.set(key, peaks[0].getF())
340 # Coordinate covariance is not used, so don't bother calulating it.
341 source.updateCoord(exposure.getWcs(), include_covariance=False)
343 # Forced photometry on sky objects
344 self.skyMeasurement.run(catalog, exposure, catalog, exposure.getWcs())
346 # Calculate new threshold
347 fluxes = catalog["base_PsfFlux_instFlux"]
348 area = catalog["base_PsfFlux_area"]
349 good = (~catalog["base_PsfFlux_flag"] & np.isfinite(fluxes))
351 if good.sum() < minNumSources:
352 if not isBgTweak:
353 msg = (f"Insufficient good sky source flux measurements ({good.sum()} < "
354 f"{minNumSources}) for dynamic threshold calculation.")
355 else:
356 msg = (f"Insufficient good sky source flux measurements ({good.sum()} < "
357 f"{minNumSources}) for background tweak calculation.")
359 nPix = exposure.mask.array.size
360 badPixelMask = lsst.afw.image.Mask.getPlaneBitMask(["NO_DATA", "BAD"])
361 nGoodPix = np.sum(exposure.mask.array & badPixelMask == 0)
362 if nGoodPix/nPix > 0.2:
363 detectedPixelMask = lsst.afw.image.Mask.getPlaneBitMask(["DETECTED", "DETECTED_NEGATIVE"])
364 nDetectedPix = np.sum((exposure.mask.array & detectedPixelMask != 0)
365 & (exposure.mask.array & badPixelMask == 0))
366 msg += (" However, {} of {} ({:.3f}%) pixels are not marked NO_DATA or BAD, "
367 "so there should be sufficient area to locate suitable sky sources. "
368 "Note that {} of {} ({:.3f}%) \"good\" pixels were marked "
369 "as DETECTED or DETECTED_NEGATIVE.".format(
370 nGoodPix, nPix, 100.0*nGoodPix/nPix,
371 nDetectedPix, nGoodPix, 100.0*nDetectedPix/nGoodPix))
372 raise InsufficientSourcesError(msg, nGoodPix, nPix)
373 raise InsufficientSourcesError(msg, nGoodPix, nPix)
375 if not isBgTweak:
376 self.log.info("Number of good sky sources used for dynamic detection: %d (of %d requested).",
377 good.sum(), self.skyObjects.config.nSources)
378 else:
379 self.log.info("Number of good sky sources used for dynamic detection background tweak:"
380 " %d (of %d requested).", good.sum(), self.skyObjects.config.nSources)
382 bgMedian = np.median((fluxes/area)[good])
383 lq, uq = np.percentile(fluxes[good], [25.0, 75.0])
384 stdevMeas = 0.741*(uq - lq)
385 medianError = np.median(catalog["base_PsfFlux_instFluxErr"][good])
386 if wcsIsNone:
387 exposure.setWcs(None)
388 return Struct(multiplicative=medianError/stdevMeas, additive=bgMedian)
390 def detectFootprints(self, exposure, doSmooth=True, sigma=None, clearMask=True, expId=None,
391 background=None, backgroundToPhotometricRatio=None):
392 """Detect footprints with a dynamic threshold
394 This varies from the vanilla ``detectFootprints`` method because we
395 do detection three times: first with a high threshold to detect
396 "bright" (both positive and negative, the latter to identify very
397 over-subtracted regions) sources for which we grow the DETECTED and
398 DETECTED_NEGATIVE masks significantly to account for wings. Second,
399 with a low threshold to mask all non-empty regions of the image. These
400 two masks are combined and used to identify regions of sky
401 uncontaminated by objects. A final round of detection is then done
402 with the new calculated threshold.
404 Parameters
405 ----------
406 exposure : `lsst.afw.image.Exposure`
407 Exposure to process; DETECTED{,_NEGATIVE} mask plane will be
408 set in-place.
409 doSmooth : `bool`, optional
410 If True, smooth the image before detection using a Gaussian
411 of width ``sigma``.
412 sigma : `float`, optional
413 Gaussian Sigma of PSF (pixels); used for smoothing and to grow
414 detections; if `None` then measure the sigma of the PSF of the
415 ``exposure``.
416 clearMask : `bool`, optional
417 Clear both DETECTED and DETECTED_NEGATIVE planes before running
418 detection.
419 expId : `int`, optional
420 Exposure identifier, used as a seed for the random number
421 generator. If absent, the seed will be the sum of the image.
422 background : `lsst.afw.math.BackgroundList`, optional
423 Background that was already subtracted from the exposure; will be
424 modified in-place if ``reEstimateBackground=True``.
425 backgroundToPhotometricRatio : `lsst.afw.image.Image`, optional
426 Unused; if set will Raise.
428 Returns
429 -------
430 results : `lsst.pipe.base.Struct`
431 The results `~lsst.pipe.base.Struct` contains:
433 ``positive``
434 Positive polarity footprints.
435 (`lsst.afw.detection.FootprintSet` or `None`)
436 ``negative``
437 Negative polarity footprints.
438 (`lsst.afw.detection.FootprintSet` or `None`)
439 ``numPos``
440 Number of footprints in positive or 0 if detection polarity was
441 negative. (`int`)
442 ``numNeg``
443 Number of footprints in negative or 0 if detection polarity was
444 positive. (`int`)
445 ``background``
446 Re-estimated background. `None` or the input ``background``
447 if ``reEstimateBackground==False``.
448 (`lsst.afw.math.BackgroundList`)
449 ``factor``
450 Multiplication factor applied to the configured detection
451 threshold. (`float`)
452 ``prelim``
453 Results from preliminary detection pass.
454 (`lsst.pipe.base.Struct`)
455 """
456 if backgroundToPhotometricRatio is not None:
457 raise RuntimeError("DynamicDetectionTask does not support backgroundToPhotometricRatio.")
458 maskedImage = exposure.maskedImage
460 if clearMask:
461 self.clearMask(maskedImage.mask)
462 else:
463 oldDetected = maskedImage.mask.array & maskedImage.mask.getPlaneBitMask(["DETECTED",
464 "DETECTED_NEGATIVE"])
465 nPix = maskedImage.mask.array.size
466 badPixelMask = lsst.afw.image.Mask.getPlaneBitMask(["NO_DATA", "BAD"])
467 nGoodPix = np.sum(maskedImage.mask.array & badPixelMask == 0)
468 self.log.info("Number of good data pixels (i.e. not NO_DATA or BAD): {} ({:.2f}% of total)".
469 format(nGoodPix, 100*nGoodPix/nPix))
470 if nGoodPix/nPix < self.config.minGoodPixelFraction:
471 msg = (f"Image has a very low good pixel fraction ({nGoodPix} of {nPix}), so not worth further "
472 "consideration")
473 raise TooManyMaskedPixelsError(msg)
475 with self.tempWideBackgroundContext(exposure):
476 # Could potentially smooth with a wider kernel than the PSF in
477 # order to better pick up the wings of stars and galaxies, but for
478 # now sticking with the PSF as that's more simple.
479 psf = self.getPsf(exposure, sigma=sigma)
480 convolveResults = self.convolveImage(maskedImage, psf, doSmooth=doSmooth)
482 if self.config.doThresholdScaling:
483 if self.config.doBrightPrelimDetection:
484 brightDetectedMask = self._computeBrightDetectionMask(maskedImage, convolveResults)
485 else:
486 prelim = None
487 factor = 1.0
489 # seed needs to fit in a C++ 'int' so pybind doesn't choke on it
490 seed = (expId if expId is not None else int(maskedImage.image.array.sum())) % (2**31 - 1)
492 middle = convolveResults.middle
493 sigma = convolveResults.sigma
494 if self.config.doThresholdScaling:
495 prelim = self.applyThreshold(
496 middle, maskedImage.getBBox(), factor=self.config.prelimThresholdFactor,
497 factorNeg=self.config.prelimNegMultiplier*self.config.prelimThresholdFactor
498 )
499 self.finalizeFootprints(
500 maskedImage.mask, prelim, sigma, factor=self.config.prelimThresholdFactor,
501 factorNeg=self.config.prelimNegMultiplier*self.config.prelimThresholdFactor
502 )
503 if self.config.doBrightPrelimDetection:
504 # Combine prelim and bright detection masks for multiplier
505 # determination.
506 maskedImage.mask.array |= brightDetectedMask
508 # Calculate the proper threshold
509 threshResults = self.calculateThreshold(exposure, seed, sigma=sigma)
510 if (self.config.minThresholdScaleFactor
511 and threshResults.multiplicative < self.config.minThresholdScaleFactor):
512 self.log.warning("Measured threshold scaling factor (%.2f) is outside [min, max] "
513 "bounds [%.2f, %.2f]. Setting factor to lower limit: %.2f.",
514 threshResults.multiplicative, self.config.minThresholdScaleFactor,
515 self.config.maxThresholdScaleFactor, self.config.minThresholdScaleFactor)
516 factor = self.config.minThresholdScaleFactor
517 elif (self.config.maxThresholdScaleFactor
518 and threshResults.multiplicative > self.config.maxThresholdScaleFactor):
519 self.log.warning("Measured threshold scaling factor (%.2f) is outside [min, max] "
520 "bounds [%.2f, %.2f]. Setting factor to upper limit: %.2f.",
521 threshResults.multiplicative, self.config.minThresholdScaleFactor,
522 self.config.maxThresholdScaleFactor, self.config.maxThresholdScaleFactor)
523 factor = self.config.maxThresholdScaleFactor
524 else:
525 factor = threshResults.multiplicative
526 self.log.info("Modifying configured detection threshold by factor %.2f to %.2f",
527 factor, factor*self.config.thresholdValue)
529 growOverride = None
530 inFinalize = True
531 while inFinalize:
532 inFinalize = False
533 # Blow away preliminary (low threshold) detection mask
534 self.clearMask(maskedImage.mask)
535 if not clearMask:
536 maskedImage.mask.array |= oldDetected
538 # Rinse and repeat thresholding with new calculated threshold
539 results = self.applyThreshold(middle, maskedImage.getBBox(), factor)
540 results.prelim = prelim
541 results.background = background if background is not None else lsst.afw.math.BackgroundList()
542 if self.config.doTempLocalBackground:
543 self.applyTempLocalBackground(exposure, middle, results)
544 self.finalizeFootprints(maskedImage.mask, results, sigma, factor=factor,
545 growOverride=growOverride)
546 if results.numPos == 0:
547 msg = "No footprints were detected, so further processing would be moot"
548 raise ZeroFootprintError(msg)
549 else:
550 self.log.warning("nPeaks/nFootprint = %.2f (max is %.1f)",
551 results.numPosPeaks/results.numPos,
552 self.config.maxPeakToFootRatio)
553 if results.numPosPeaks/results.numPos > self.config.maxPeakToFootRatio:
554 if results.numPosPeaks/results.numPos > 3*self.config.maxPeakToFootRatio:
555 factor *= 1.4
556 else:
557 factor *= 1.2
558 if factor > 2.0:
559 if growOverride is None:
560 growOverride = 0.75*self.config.nSigmaToGrow
561 else:
562 growOverride *= 0.75
563 self.log.warning("Decreasing nSigmaToGrow to %.2f", growOverride)
564 if factor >= 5:
565 self.log.warning("New theshold value would be > 5 times the initially requested "
566 "one (%.2f > %.2f). Leaving inFinalize iteration without "
567 "getting the number of peaks per footprint below %.1f",
568 factor*self.config.thresholdValue, self.config.thresholdValue,
569 self.config.maxPeakToFootRatio)
570 inFinalize = False
571 else:
572 inFinalize = True
573 self.log.warning("numPosPeaks/numPos (%d) > maxPeakPerFootprint (%.1f). "
574 "Increasing threshold factor to %.2f and re-running,",
575 results.numPosPeaks/results.numPos,
576 self.config.maxPeakToFootRatio, factor)
578 self.clearUnwantedResults(maskedImage.mask, results)
580 if self.config.reEstimateBackground:
581 self.reEstimateBackground(maskedImage, results.background)
583 self.display(exposure, results, middle)
585 # Re-do the background tweak after any temporary backgrounds have
586 # been restored.
587 #
588 # But we want to keep any large-scale background (e.g., scattered
589 # light from bright stars) from being selected for sky objects in
590 # the calculation, so do another detection pass without either the
591 # local or wide temporary background subtraction; the DETECTED
592 # pixels will mark the area to ignore.
594 # The following if/else is to workaround the fact that it is
595 # currently not possible to persist an empty BackgroundList, so
596 # we instead set the value of the backround tweak to 0.0 if
597 # doBackgroundTweak is False and call the tweakBackground function
598 # regardless to get at least one background into the list (do we
599 # need a TODO here?).
600 if self.config.doBackgroundTweak:
601 originalMask = maskedImage.mask.array.copy()
602 try:
603 self.clearMask(exposure.mask)
604 convolveResults = self.convolveImage(maskedImage, psf, doSmooth=doSmooth)
605 tweakDetResults = self.applyThreshold(convolveResults.middle, maskedImage.getBBox(), factor)
606 self.finalizeFootprints(maskedImage.mask, tweakDetResults, sigma, factor=factor)
607 bgLevel = self.calculateThreshold(exposure, seed, sigma=sigma, minFractionSourcesFactor=0.5,
608 isBgTweak=True).additive
609 if self.config.minBackgroundTweak and bgLevel < self.config.minBackgroundTweak:
610 self.log.warning("Measured background tweak (%.2f) is outside [min, max] bounds "
611 "[%.2f, %.2f]. Setting tweak to lower limit: %.2f.", bgLevel,
612 self.config.minBackgroundTweak, self.config.maxBackgroundTweak,
613 self.config.minBackgroundTweak)
614 bgLevel = self.config.minBackgroundTweak
615 if self.config.maxBackgroundTweak and bgLevel > self.config.maxBackgroundTweak:
616 self.log.warning("Measured background tweak (%.2f) is outside [min, max] bounds "
617 "[%.2f, %.2f]. Setting tweak to upper limit: %.2f.", bgLevel,
618 self.config.minBackgroundTweak, self.config.maxBackgroundTweak,
619 self.config.maxBackgroundTweak)
620 bgLevel = self.config.maxBackgroundTweak
621 finally:
622 maskedImage.mask.array[:] = originalMask
623 else:
624 bgLevel = 0.0
625 self.tweakBackground(exposure, bgLevel, results.background)
627 return results
629 def tweakBackground(self, exposure, bgLevel, bgList=None):
630 """Modify the background by a constant value
632 Parameters
633 ----------
634 exposure : `lsst.afw.image.Exposure`
635 Exposure for which to tweak background.
636 bgLevel : `float`
637 Background level to remove
638 bgList : `lsst.afw.math.BackgroundList`, optional
639 List of backgrounds to append to.
641 Returns
642 -------
643 bg : `lsst.afw.math.BackgroundMI`
644 Constant background model.
645 """
646 if bgLevel != 0.0:
647 self.log.info("Tweaking background by %f to match sky photometry", bgLevel)
648 exposure.image -= bgLevel
649 bgStats = lsst.afw.image.MaskedImageF(1, 1)
650 bgStats.set(bgLevel, 0, bgLevel)
651 bg = lsst.afw.math.BackgroundMI(exposure.getBBox(), bgStats)
652 bgData = (bg, lsst.afw.math.Interpolate.LINEAR, lsst.afw.math.REDUCE_INTERP_ORDER,
653 lsst.afw.math.ApproximateControl.UNKNOWN, 0, 0, False)
654 if bgList is not None:
655 bgList.append(bgData)
656 return bg
658 def _computeBrightDetectionMask(self, maskedImage, convolveResults):
659 """Perform an initial bright source detection pass.
661 Perform an initial bright object detection pass using a high detection
662 threshold. The footprints in this pass are grown significantly more
663 than is typical to account for wings around bright sources. The
664 negative polarity detections in this pass help in masking severely
665 over-subtracted regions.
667 A maximum fraction of masked pixel from this pass is ensured via
668 the config ``brightMaskFractionMax``. If the masked pixel fraction is
669 above this value, the detection thresholds here are increased by
670 ``bisectFactor`` in a while loop until the detected masked fraction
671 falls below this value.
673 Parameters
674 ----------
675 maskedImage : `lsst.afw.image.MaskedImage`
676 Masked image on which to run the detection.
677 convolveResults : `lsst.pipe.base.Struct`
678 The results of the self.convolveImage function with attributes:
680 ``middle``
681 Convolved image, without the edges
682 (`lsst.afw.image.MaskedImage`).
683 ``sigma``
684 Gaussian sigma used for the convolution (`float`).
686 Returns
687 -------
688 brightDetectedMask : `numpy.ndarray`
689 Boolean array representing the union of the bright detection pass
690 DETECTED and DETECTED_NEGATIVE masks.
691 """
692 # Initialize some parameters.
693 brightPosFactor = (
694 self.config.prelimThresholdFactor*self.config.brightMultiplier/self.config.bisectFactor
695 )
696 brightNegFactor = self.config.brightNegFactor/self.config.bisectFactor
697 nPix = 1
698 nPixDet = 1
699 nPixDetNeg = 1
700 brightMaskFractionMax = self.config.brightMaskFractionMax
702 # Loop until masked fraction is smaller than
703 # brightMaskFractionMax, increasing the thresholds by
704 # config.bisectFactor on each iteration (rarely necessary
705 # for current defaults).
706 while nPixDetNeg/nPix > brightMaskFractionMax or nPixDet/nPix > brightMaskFractionMax:
707 self.clearMask(maskedImage.mask)
708 brightPosFactor *= self.config.bisectFactor
709 brightNegFactor *= self.config.bisectFactor
710 prelimBright = self.applyThreshold(convolveResults.middle, maskedImage.getBBox(),
711 factor=brightPosFactor, factorNeg=brightNegFactor)
712 self.finalizeFootprints(
713 maskedImage.mask, prelimBright, convolveResults.sigma*self.config.brightGrowFactor,
714 factor=brightPosFactor, factorNeg=brightNegFactor
715 )
716 # Check that not too many pixels got masked.
717 nPix = maskedImage.mask.array.size
718 nPixDet = countMaskedPixels(maskedImage, "DETECTED")
719 self.log.info("Number (%) of bright DETECTED pix: {} ({:.1f}%)".
720 format(nPixDet, 100*nPixDet/nPix))
721 nPixDetNeg = countMaskedPixels(maskedImage, "DETECTED_NEGATIVE")
722 self.log.info("Number (%) of bright DETECTED_NEGATIVE pix: {} ({:.1f}%)".
723 format(nPixDetNeg, 100*nPixDetNeg/nPix))
724 if nPixDetNeg/nPix > brightMaskFractionMax or nPixDet/nPix > brightMaskFractionMax:
725 self.log.warn("Too high a fraction (%.1f > %.1f) of pixels were masked with current "
726 "\"bright\" detection round thresholds. Increasing by a factor of %.2f "
727 "and trying again.", max(nPixDetNeg, nPixDet)/nPix,
728 brightMaskFractionMax, self.config.bisectFactor)
730 # Save the mask planes from the "bright" detection round, then
731 # clear them before moving on to the "prelim" detection phase.
732 brightDetectedMask = (maskedImage.mask.array
733 & maskedImage.mask.getPlaneBitMask(["DETECTED", "DETECTED_NEGATIVE"]))
734 self.clearMask(maskedImage.mask)
735 return brightDetectedMask
738def countMaskedPixels(maskedIm, maskPlane):
739 """Count the number of pixels in a given mask plane.
741 Parameters
742 ----------
743 maskedIm : `lsst.afw.image.MaskedImage`
744 Masked image to examine.
745 maskPlane : `str`
746 Name of the mask plane to examine.
748 Returns
749 -------
750 nPixMasked : `int`
751 Number of pixels with ``maskPlane`` bit set.
752 """
753 maskBit = maskedIm.mask.getPlaneBitMask(maskPlane)
754 nPixMasked = np.sum(np.bitwise_and(maskedIm.mask.array, maskBit))/maskBit
755 return nPixMasked