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