Coverage for python/lsst/meas/algorithms/dynamicDetection.py: 12%
325 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-05-30 08:51 +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 detectedPosPixelMask = lsst.afw.image.Mask.getPlaneBitMask(["DETECTED"])
367 detectedNegPixelMask = lsst.afw.image.Mask.getPlaneBitMask(["DETECTED_NEGATIVE"])
368 nDetectedPix = np.sum(((exposure.mask.array & detectedPosPixelMask != 0)
369 | (exposure.mask.array & detectedNegPixelMask != 0))
370 & (exposure.mask.array & badPixelMask == 0))
371 msg += (" However, {} of {} ({:.3f}%) pixels are not marked NO_DATA or BAD, "
372 "so there should be sufficient area to locate suitable sky sources. "
373 "Note that {} of {} ({:.3f}%) \"good\" pixels were marked "
374 "as DETECTED or DETECTED_NEGATIVE.".format(
375 nGoodPix, nPix, 100.0*nGoodPix/nPix,
376 nDetectedPix, nGoodPix, 100.0*nDetectedPix/nGoodPix))
377 raise InsufficientSourcesError(msg, nGoodPix, nPix)
378 raise InsufficientSourcesError(msg, nGoodPix, nPix)
380 if not isBgTweak:
381 self.log.info("Number of good sky sources used for dynamic detection: %d (of %d requested).",
382 good.sum(), self.skyObjects.config.nSources)
383 else:
384 self.log.info("Number of good sky sources used for dynamic detection background tweak:"
385 " %d (of %d requested).", good.sum(), self.skyObjects.config.nSources)
387 bgMedian = np.median((fluxes/area)[good])
388 lq, uq = np.percentile(fluxes[good], [25.0, 75.0])
389 stdevMeas = 0.741*(uq - lq)
390 medianError = np.median(catalog["base_PsfFlux_instFluxErr"][good])
391 if wcsIsNone:
392 exposure.setWcs(None)
393 return Struct(multiplicative=medianError/stdevMeas, additive=bgMedian)
395 def detectFootprints(self, exposure, doSmooth=True, sigma=None, clearMask=True, expId=None,
396 background=None, backgroundToPhotometricRatio=None):
397 """Detect footprints with a dynamic threshold
399 This varies from the vanilla ``detectFootprints`` method because we
400 do detection three times: first with a high threshold to detect
401 "bright" (both positive and negative, the latter to identify very
402 over-subtracted regions) sources for which we grow the DETECTED and
403 DETECTED_NEGATIVE masks significantly to account for wings. Second,
404 with a low threshold to mask all non-empty regions of the image. These
405 two masks are combined and used to identify regions of sky
406 uncontaminated by objects. A final round of detection is then done
407 with the new calculated threshold.
409 Parameters
410 ----------
411 exposure : `lsst.afw.image.Exposure`
412 Exposure to process; DETECTED{,_NEGATIVE} mask plane will be
413 set in-place.
414 doSmooth : `bool`, optional
415 If True, smooth the image before detection using a Gaussian
416 of width ``sigma``.
417 sigma : `float`, optional
418 Gaussian Sigma of PSF (pixels); used for smoothing and to grow
419 detections; if `None` then measure the sigma of the PSF of the
420 ``exposure``.
421 clearMask : `bool`, optional
422 Clear both DETECTED and DETECTED_NEGATIVE planes before running
423 detection.
424 expId : `int`, optional
425 Exposure identifier, used as a seed for the random number
426 generator. If absent, the seed will be the sum of the image.
427 background : `lsst.afw.math.BackgroundList`, optional
428 Background that was already subtracted from the exposure; will be
429 modified in-place if ``reEstimateBackground=True``.
430 backgroundToPhotometricRatio : `lsst.afw.image.Image`, optional
431 Unused; if set will Raise.
433 Returns
434 -------
435 results : `lsst.pipe.base.Struct`
436 The results `~lsst.pipe.base.Struct` contains:
438 ``positive``
439 Positive polarity footprints.
440 (`lsst.afw.detection.FootprintSet` or `None`)
441 ``negative``
442 Negative polarity footprints.
443 (`lsst.afw.detection.FootprintSet` or `None`)
444 ``numPos``
445 Number of footprints in positive or 0 if detection polarity was
446 negative. (`int`)
447 ``numNeg``
448 Number of footprints in negative or 0 if detection polarity was
449 positive. (`int`)
450 ``background``
451 Re-estimated background. `None` or the input ``background``
452 if ``reEstimateBackground==False``.
453 (`lsst.afw.math.BackgroundList`)
454 ``factor``
455 Multiplication factor applied to the configured detection
456 threshold. (`float`)
457 ``prelim``
458 Results from preliminary detection pass.
459 (`lsst.pipe.base.Struct`)
460 """
461 if backgroundToPhotometricRatio is not None:
462 raise RuntimeError("DynamicDetectionTask does not support backgroundToPhotometricRatio.")
463 maskedImage = exposure.maskedImage
465 if clearMask:
466 self.clearMask(maskedImage.mask)
467 else:
468 oldDetected = maskedImage.mask.array & maskedImage.mask.getPlaneBitMask(["DETECTED",
469 "DETECTED_NEGATIVE"])
470 nPix = maskedImage.mask.array.size
471 badPixelMask = lsst.afw.image.Mask.getPlaneBitMask(["NO_DATA", "BAD"])
472 nGoodPix = np.sum(maskedImage.mask.array & badPixelMask == 0)
473 self.log.info("Number of good data pixels (i.e. not NO_DATA or BAD): {} ({:.2f}% of total)".
474 format(nGoodPix, 100*nGoodPix/nPix))
475 if nGoodPix/nPix < self.config.minGoodPixelFraction:
476 msg = (f"Image has a very low good pixel fraction ({nGoodPix} of {nPix}), so not worth further "
477 "consideration")
478 raise TooManyMaskedPixelsError(msg)
480 with self.tempWideBackgroundContext(exposure):
481 # Could potentially smooth with a wider kernel than the PSF in
482 # order to better pick up the wings of stars and galaxies, but for
483 # now sticking with the PSF as that's more simple.
484 psf = self.getPsf(exposure, sigma=sigma)
485 convolveResults = self.convolveImage(maskedImage, psf, doSmooth=doSmooth)
487 if self.config.doThresholdScaling:
488 if self.config.doBrightPrelimDetection:
489 brightDetectedMask, brightFactorNeg = self._computeBrightDetectionMask(
490 maskedImage, convolveResults)
491 # Scale the factor for negative polarity detections based
492 # on what was required in the bright detection pass to
493 # avoid too many pixels marked as DETECTED_NEGATIVE (but
494 # capping it at 20.0 as a guardrail).
495 factorNeg = min(20.0, brightFactorNeg/self.config.brightNegFactor)
496 else:
497 prelim = None
498 factor = 1.0
499 factorNeg = 1.0
501 # seed needs to fit in a C++ 'int' so pybind doesn't choke on it
502 seed = (expId if expId is not None else int(maskedImage.image.array.sum())) % (2**31 - 1)
504 middle = convolveResults.middle
505 sigma = convolveResults.sigma
506 if self.config.doThresholdScaling:
507 factorNeg *= self.config.prelimNegMultiplier*self.config.prelimThresholdFactor
508 prelim = self.applyThreshold(
509 middle, maskedImage.getBBox(), factor=self.config.prelimThresholdFactor,
510 factorNeg=factorNeg
511 )
512 self.finalizeFootprints(
513 maskedImage.mask, prelim, sigma, factor=self.config.prelimThresholdFactor,
514 factorNeg=factorNeg
515 )
516 if self.config.doBrightPrelimDetection:
517 # Combine prelim and bright detection masks for multiplier
518 # determination.
519 maskedImage.mask.array |= brightDetectedMask
521 # Calculate the proper threshold
522 threshResults = self.calculateThreshold(exposure, seed, sigma=sigma)
523 if (self.config.minThresholdScaleFactor
524 and threshResults.multiplicative < self.config.minThresholdScaleFactor):
525 self.log.warning("Measured threshold scaling factor (%.2f) is outside [min, max] "
526 "bounds [%.2f, %.2f]. Setting factor to lower limit: %.2f.",
527 threshResults.multiplicative, self.config.minThresholdScaleFactor,
528 self.config.maxThresholdScaleFactor, self.config.minThresholdScaleFactor)
529 factor = self.config.minThresholdScaleFactor
530 elif (self.config.maxThresholdScaleFactor
531 and threshResults.multiplicative > self.config.maxThresholdScaleFactor):
532 self.log.warning("Measured threshold scaling factor (%.2f) is outside [min, max] "
533 "bounds [%.2f, %.2f]. Setting factor to upper limit: %.2f.",
534 threshResults.multiplicative, self.config.minThresholdScaleFactor,
535 self.config.maxThresholdScaleFactor, self.config.maxThresholdScaleFactor)
536 factor = self.config.maxThresholdScaleFactor
537 else:
538 factor = threshResults.multiplicative
539 # Also scale the factor for negative polarity detections
540 factorNeg *= factor
541 self.log.info("Modifying configured detection threshold by factor %.2f to %.2f",
542 factor, factor*self.config.thresholdValue)
544 growOverride = None
545 inFinalize = True
546 while inFinalize:
547 inFinalize = False
548 # Blow away preliminary (low threshold) detection mask
549 self.clearMask(maskedImage.mask)
550 if not clearMask:
551 maskedImage.mask.array |= oldDetected
553 # Rinse and repeat thresholding with new calculated threshold
554 results = self.applyThreshold(
555 middle, maskedImage.getBBox(), factor=factor, factorNeg=factorNeg
556 )
557 results.prelim = prelim
558 results.background = background if background is not None else lsst.afw.math.BackgroundList()
559 if self.config.doTempLocalBackground:
560 self.applyTempLocalBackground(exposure, middle, results)
561 self.finalizeFootprints(maskedImage.mask, results, sigma, factor=factor, factorNeg=factorNeg,
562 growOverride=growOverride)
563 if results.numPos == 0:
564 msg = "No footprints were detected, so further processing would be moot"
565 raise ZeroFootprintError(msg)
566 else:
567 self.log.warning("nPeaks/nFootprint = %.2f (max is %.1f)",
568 results.numPosPeaks/results.numPos,
569 self.config.maxPeakToFootRatio)
570 if results.numPosPeaks/results.numPos > self.config.maxPeakToFootRatio:
571 if results.numPosPeaks/results.numPos > 3*self.config.maxPeakToFootRatio:
572 factor *= 1.4
573 else:
574 factor *= 1.2
575 if factor > 2.0:
576 if growOverride is None:
577 growOverride = 0.75*self.config.nSigmaToGrow
578 else:
579 growOverride *= 0.75
580 self.log.warning("Decreasing nSigmaToGrow to %.2f", growOverride)
581 if factor >= 5:
582 self.log.warning("New theshold value would be > 5 times the initially requested "
583 "one (%.2f > %.2f). Leaving inFinalize iteration without "
584 "getting the number of peaks per footprint below %.1f",
585 factor*self.config.thresholdValue, self.config.thresholdValue,
586 self.config.maxPeakToFootRatio)
587 inFinalize = False
588 else:
589 inFinalize = True
590 self.log.warning("numPosPeaks/numPos (%d) > maxPeakPerFootprint (%.1f). "
591 "Increasing threshold factor to %.2f and re-running,",
592 results.numPosPeaks/results.numPos,
593 self.config.maxPeakToFootRatio, factor)
595 self.clearUnwantedResults(maskedImage.mask, results)
597 if self.config.reEstimateBackground:
598 self.reEstimateBackground(maskedImage, results.background)
600 self.display(exposure, results, middle)
602 # Re-do the background tweak after any temporary backgrounds have
603 # been restored.
604 #
605 # But we want to keep any large-scale background (e.g., scattered
606 # light from bright stars) from being selected for sky objects in
607 # the calculation, so do another detection pass without either the
608 # local or wide temporary background subtraction; the DETECTED
609 # pixels will mark the area to ignore.
611 # The following if/else is to workaround the fact that it is
612 # currently not possible to persist an empty BackgroundList, so
613 # we instead set the value of the backround tweak to 0.0 if
614 # doBackgroundTweak is False and call the tweakBackground function
615 # regardless to get at least one background into the list (do we
616 # need a TODO here?).
617 if self.config.doBackgroundTweak:
618 originalMask = maskedImage.mask.array.copy()
619 try:
620 self.clearMask(exposure.mask)
621 convolveResults = self.convolveImage(maskedImage, psf, doSmooth=doSmooth)
622 # Don't use factorNeg if the image has had its background
623 # subtracted after source detection.
624 tweakFactorNeg = None if self.config.reEstimateBackground else factorNeg
625 tweakDetResults = self.applyThreshold(convolveResults.middle, maskedImage.getBBox(),
626 factor=factor, factorNeg=tweakFactorNeg)
627 self.finalizeFootprints(maskedImage.mask, tweakDetResults, sigma, factor=factor,
628 factorNeg=tweakFactorNeg)
629 bgLevel = self.calculateThreshold(exposure, seed, sigma=sigma, minFractionSourcesFactor=0.5,
630 isBgTweak=True).additive
631 if self.config.minBackgroundTweak:
632 minBackgroundTweak = self.config.minBackgroundTweak
633 # Increase the minimum bg tweak allowed if the factorNeg indicates image
634 # was very oversubtracted.
635 if factorNeg > 5*factor:
636 if factorNeg > 8*factor:
637 minBackgroundTweak *= 3.0
638 else:
639 minBackgroundTweak *= 2.0
640 self.log.warning("All evidence suggests the image is very oversubtracted. "
641 "Allowing for a larger corection (%.2f), than that set in "
642 "config.minBackgroundTweak (%.2f).", minBackgroundTweak,
643 self.config.minBackgroundTweak)
644 if bgLevel < minBackgroundTweak:
645 self.log.warning("Measured background tweak (%.2f) is outside [min, max] bounds "
646 "[%.2f, %.2f]. Setting tweak to lower limit: %.2f.", bgLevel,
647 minBackgroundTweak, self.config.maxBackgroundTweak,
648 minBackgroundTweak)
649 bgLevel = minBackgroundTweak
650 if self.config.maxBackgroundTweak and bgLevel > self.config.maxBackgroundTweak:
651 self.log.warning("Measured background tweak (%.2f) is outside [min, max] bounds "
652 "[%.2f, %.2f]. Setting tweak to upper limit: %.2f.", bgLevel,
653 self.config.minBackgroundTweak, self.config.maxBackgroundTweak,
654 self.config.maxBackgroundTweak)
655 bgLevel = self.config.maxBackgroundTweak
656 finally:
657 maskedImage.mask.array[:] = originalMask
658 else:
659 bgLevel = 0.0
660 self.tweakBackground(exposure, bgLevel, results.background)
662 return results
664 def tweakBackground(self, exposure, bgLevel, bgList=None):
665 """Modify the background by a constant value
667 Parameters
668 ----------
669 exposure : `lsst.afw.image.Exposure`
670 Exposure for which to tweak background.
671 bgLevel : `float`
672 Background level to remove
673 bgList : `lsst.afw.math.BackgroundList`, optional
674 List of backgrounds to append to.
676 Returns
677 -------
678 bg : `lsst.afw.math.BackgroundMI`
679 Constant background model.
680 """
681 if bgLevel != 0.0:
682 self.log.info("Tweaking background by %.3f to match sky photometry", bgLevel)
683 exposure.image -= bgLevel
684 bgStats = lsst.afw.image.MaskedImageF(1, 1)
685 bgStats.set(bgLevel, 0, bgLevel)
686 bg = lsst.afw.math.BackgroundMI(exposure.getBBox(), bgStats)
687 bgData = (bg, lsst.afw.math.Interpolate.LINEAR, lsst.afw.math.REDUCE_INTERP_ORDER,
688 lsst.afw.math.ApproximateControl.UNKNOWN, 0, 0, False)
689 if bgList is not None:
690 bgList.append(bgData)
691 return bg
693 def _computeBrightDetectionMask(self, maskedImage, convolveResults):
694 """Perform an initial bright source detection pass.
696 Perform an initial bright object detection pass using a high detection
697 threshold. The footprints in this pass are grown significantly more
698 than is typical to account for wings around bright sources. The
699 negative polarity detections in this pass help in masking severely
700 over-subtracted regions.
702 A maximum fraction of masked pixel from this pass is ensured via
703 the config ``brightMaskFractionMax``. If the masked pixel fraction is
704 above this value, the detection thresholds here are increased by
705 ``bisectFactor`` in a while loop until the detected masked fraction
706 falls below this value.
708 Parameters
709 ----------
710 maskedImage : `lsst.afw.image.MaskedImage`
711 Masked image on which to run the detection.
712 convolveResults : `lsst.pipe.base.Struct`
713 The results of the self.convolveImage function with attributes:
715 ``middle``
716 Convolved image, without the edges
717 (`lsst.afw.image.MaskedImage`).
718 ``sigma``
719 Gaussian sigma used for the convolution (`float`).
721 Returns
722 -------
723 brightDetectedMask : `numpy.ndarray`
724 Boolean array representing the union of the bright detection pass
725 DETECTED and DETECTED_NEGATIVE masks.
726 brightFactorNeg : `float`
727 Factor applied to the threshold for negative polarity detections.
728 This can get altered significantly in background over-subtracted
729 images and is useful to know for subsequent steps in the dynamic
730 detection.
731 """
732 # Initialize some parameters.
733 brightFactorPos = (
734 self.config.prelimThresholdFactor*self.config.brightMultiplier
735 )
736 brightFactorNeg = self.config.brightNegFactor
737 brightMaskFractionMax = self.config.brightMaskFractionMax
738 # Set a lower max value tolerated for negative detection mask fraction.
739 brightMaskNegFractionMax = max(0.3, 0.75*brightMaskFractionMax)
741 badPixelMask = lsst.afw.image.Mask.getPlaneBitMask(["NO_DATA", "BAD"])
742 nGoodPix = np.sum(maskedImage.mask.array & badPixelMask == 0)
744 # Loop until masked fraction is smaller than
745 # brightMaskFractionMax, increasing the thresholds by
746 # config.bisectFactor on each iteration (rarely necessary
747 # for current defaults).
748 for nIter in range(self.config.brightDetectionIterMax):
749 self.clearMask(maskedImage.mask)
750 prelimBright = self.applyThreshold(convolveResults.middle, maskedImage.getBBox(),
751 factor=brightFactorPos, factorNeg=brightFactorNeg)
752 self.finalizeFootprints(
753 maskedImage.mask, prelimBright, convolveResults.sigma*self.config.brightGrowFactor,
754 factor=brightFactorPos, factorNeg=brightFactorNeg
755 )
756 # Check that not too many pixels got masked.
757 detectedPosPixelMask = lsst.afw.image.Mask.getPlaneBitMask(["DETECTED"])
758 detectedNegPixelMask = lsst.afw.image.Mask.getPlaneBitMask(["DETECTED_NEGATIVE"])
759 badPixelMask = lsst.afw.image.Mask.getPlaneBitMask(["NO_DATA", "BAD"])
760 nPixDetPos = np.sum((maskedImage.mask.array & detectedPosPixelMask != 0)
761 & (maskedImage.mask.array & badPixelMask == 0))
762 nPixDetNeg = np.sum((maskedImage.mask.array & detectedNegPixelMask != 0)
763 & (maskedImage.mask.array & badPixelMask == 0))
764 self.log.info("Number (%) of bright DETECTED pix: {} ({:.1f}%)".
765 format(nPixDetPos, 100*nPixDetPos/nGoodPix))
766 self.log.info("Number (%) of bright DETECTED_NEGATIVE pix: {} ({:.1f}%)".
767 format(nPixDetNeg, 100*nPixDetNeg/nGoodPix))
769 if nPixDetPos/nGoodPix > brightMaskFractionMax or nPixDetNeg/nGoodPix > brightMaskNegFractionMax:
770 if nIter == self.config.brightDetectionIterMax - 1:
771 self.log.warning("Reached maximum number of iterations and still have too high "
772 "detected mask fractions in bright detection pass. Image is "
773 "likely mostly masked with BAD or NO_DATA or \"bad\" in some "
774 "other respect (so expected to likely fail further downstream).")
775 break
776 if nPixDetPos/nGoodPix > brightMaskFractionMax:
777 brightFactorPos *= self.config.bisectFactor
778 self.log.warning("Too high a fraction (%.2f > %.2f) of pixels were masked as "
779 "DETECTED with current \"bright\" detection round thresholds "
780 "(at nIter = %d). Increasing to a factor of %.2f and trying again.",
781 nPixDetPos/nGoodPix, brightMaskFractionMax, nIter, brightFactorPos)
782 if nPixDetNeg/nGoodPix > brightMaskNegFractionMax:
783 extraFactorNeg = 1.2
784 if nPixDetNeg/nGoodPix > min(0.98, 1.25*brightMaskNegFractionMax):
785 if nIter == 0:
786 extraFactorNeg = 3.0
787 elif nPixDetNeg/nGoodPix > 0.9999:
788 extraFactorNeg = 1.8
789 else:
790 extraFactorNeg = 1.5
791 brightFactorNeg *= self.config.bisectFactor*extraFactorNeg
792 self.log.warning("Too high a fraction (%.2f > %.2f) of pixels were masked as "
793 "DETECTED_NEGATIVE with current \"bright\" detection round thresholds "
794 "(at nIter = %d). Increasing to a factor of %.2f and trying again.",
795 nPixDetNeg/nGoodPix, brightMaskNegFractionMax, nIter, brightFactorNeg)
796 else:
797 break
799 # Save the mask planes from the "bright" detection round, then
800 # clear them before moving on to the "prelim" detection phase.
801 brightDetectedMask = (maskedImage.mask.array
802 & maskedImage.mask.getPlaneBitMask(["DETECTED", "DETECTED_NEGATIVE"]))
803 self.clearMask(maskedImage.mask)
804 return brightDetectedMask, brightFactorNeg
807def countMaskedPixels(maskedIm, maskPlane):
808 """Count the number of pixels in a given mask plane.
810 Parameters
811 ----------
812 maskedIm : `lsst.afw.image.MaskedImage`
813 Masked image to examine.
814 maskPlane : `str`
815 Name of the mask plane to examine.
817 Returns
818 -------
819 nPixMasked : `int`
820 Number of pixels with ``maskPlane`` bit set.
821 """
822 maskBit = maskedIm.mask.getPlaneBitMask(maskPlane)
823 nPixMasked = np.sum(np.bitwise_and(maskedIm.mask.array, maskBit))/maskBit
824 return nPixMasked