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