Coverage for python/lsst/meas/algorithms/dynamicDetection.py: 15%
222 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-04 02:01 -0800
2__all__ = ["DynamicDetectionConfig", "DynamicDetectionTask", "InsufficientSourcesError"]
4import numpy as np
6from lsst.pex.config import Field, ConfigurableField
7from lsst.pipe.base import Struct
9from .detection import SourceDetectionConfig, SourceDetectionTask
10from .skyObjects import SkyObjectsTask
12from lsst.afw.detection import FootprintSet
13from lsst.afw.geom import makeCdMatrix, makeSkyWcs, SpanSet
14from lsst.afw.table import SourceCatalog, SourceTable
15from lsst.meas.base import ForcedMeasurementTask
17import lsst.afw.image
18import lsst.afw.math
19import lsst.geom as geom
22class InsufficientSourcesError(Exception):
23 """Raised if an insufficient number of sky sources are found for
24 dynamic detection.
26 Parameters
27 ----------
28 msg : `str`
29 Error message.
30 nGoodPix : `int`
31 Number of good pixels (i.e. not NO_DATA or BAD).
32 nPix : `int`
33 Total number of pixels.
34 **kwargs : `dict`, optional
35 Additional keyword arguments to initialize the Exception base class.
36 """
37 def __init__(self, msg, nGoodPix, nPix, **kwargs):
38 self.msg = msg
39 self._metadata = kwargs
40 super().__init__(msg, **kwargs)
41 self._metadata["nGoodPix"] = int(nGoodPix)
42 self._metadata["nPix"] = int(nPix)
44 def __str__(self):
45 # Exception doesn't handle **kwargs, so we need a custom str.
46 return f"{self.msg}: {self.metadata}"
48 @property
49 def metadata(self):
50 for key, value in self._metadata.items():
51 if not isinstance(value, (int, float, str)):
52 raise TypeError(f"{key} is of type {type(value)}, but only (int, float, str) are allowed.")
53 return self._metadata
56class DynamicDetectionConfig(SourceDetectionConfig):
57 """Configuration for DynamicDetectionTask
58 """
59 prelimThresholdFactor = Field(dtype=float, default=0.5,
60 doc="Factor by which to multiply the main detection threshold "
61 "(thresholdValue) to use for first pass (to find sky objects).")
62 prelimNegMultiplier = Field(dtype=float, default=2.5,
63 doc="Multiplier for the negative (relative to positive) polarity "
64 "detections threshold to use for first pass (to find sky objects).")
65 skyObjects = ConfigurableField(target=SkyObjectsTask, doc="Generate sky objects.")
66 doBackgroundTweak = Field(dtype=bool, default=True,
67 doc="Tweak background level so median PSF flux of sky objects is zero?")
68 minFractionSources = Field(dtype=float, default=0.02,
69 doc="Minimum fraction of the requested number of sky sources for dynamic "
70 "detection to be considered a success. If the number of good sky sources "
71 "identified falls below this threshold, a NoWorkFound error is raised so "
72 "that this dataId is no longer considered in downstream processing.")
73 doBrightPrelimDetection = Field(dtype=bool, default=True,
74 doc="Do initial bright detection pass where footprints are grown "
75 "by brightGrowFactor?")
76 brightMultiplier = Field(dtype=float, default=2000.0,
77 doc="Multiplier to apply to the prelimThresholdFactor for the "
78 "\"bright\" detections stage (want this to be large to only "
79 "detect the brightest sources).")
80 brightNegFactor = Field(dtype=float, default=2.2,
81 doc="Factor by which to multiply the threshold for the negative polatiry "
82 "detections for the \"bright\" detections stage (this needs to be fairly "
83 "low given the nature of the negative polarity detections in the very "
84 "large positive polarity threshold).")
85 brightGrowFactor = Field(dtype=int, default=40,
86 doc="Factor by which to grow the footprints of sources detected in the "
87 "\"bright\" detections stage (want this to be large to mask wings of "
88 "bright sources).")
89 brightMaskFractionMax = Field(dtype=float, default=0.95,
90 doc="Maximum allowed fraction of masked pixes from the \"bright\" "
91 "detection stage (to mask regions unsuitable for sky sourcess). "
92 "If this fraction is exeeded, the detection threshold for this stage "
93 "will be increased by bisectFactor until the fraction of masked "
94 "pixels drops below this threshold.")
95 bisectFactor = Field(dtype=float, default=1.2,
96 doc="Factor by which to increase thresholds in brightMaskFractionMax loop.")
97 allowMaskErode = Field(dtype=bool, default=True,
98 doc="Crowded/large fill-factor fields make it difficult to find "
99 "suitable locations to lay down sky objects. To allow for best effort "
100 "sky source placement, if True, this allows for a slight erosion of "
101 "the detection masks.")
103 def setDefaults(self):
104 SourceDetectionConfig.setDefaults(self)
105 self.skyObjects.nSources = 1000 # For good statistics
106 for maskStr in ["SAT"]:
107 if maskStr not in self.skyObjects.avoidMask:
108 self.skyObjects.avoidMask.append(maskStr)
110 def validate(self):
111 super().validate()
113 if self.doApplyFlatBackgroundRatio:
114 raise ValueError("DynamicDetectionTask does not support doApplyFlatBackgroundRatio.")
117class DynamicDetectionTask(SourceDetectionTask):
118 """Detection of sources on an image with a dynamic threshold
120 We first detect sources using a lower threshold than normal (see config
121 parameter ``prelimThresholdFactor``) in order to identify good sky regions
122 (configurable ``skyObjects``). Then we perform forced PSF photometry on
123 those sky regions. Using those PSF flux measurements and estimated errors,
124 we set the threshold so that the stdev of the measurements matches the
125 median estimated error.
127 Besides the usual initialisation of configurables, we also set up
128 the forced measurement which is deliberately not represented in
129 this Task's configuration parameters because we're using it as
130 part of the algorithm and we don't want to allow it to be modified.
131 """
132 ConfigClass = DynamicDetectionConfig
133 _DefaultName = "dynamicDetection"
135 def __init__(self, *args, **kwargs):
137 SourceDetectionTask.__init__(self, *args, **kwargs)
138 self.makeSubtask("skyObjects")
140 # Set up forced measurement.
141 config = ForcedMeasurementTask.ConfigClass()
142 config.plugins.names = ['base_TransformedCentroid', 'base_PsfFlux', 'base_LocalBackground']
143 # We'll need the "centroid" and "psfFlux" slots
144 for slot in ("shape", "psfShape", "apFlux", "modelFlux", "gaussianFlux", "calibFlux"):
145 setattr(config.slots, slot, None)
146 config.copyColumns = {}
147 self.skySchema = SourceTable.makeMinimalSchema()
148 self.skyMeasurement = ForcedMeasurementTask(config=config, name="skyMeasurement", parentTask=self,
149 refSchema=self.skySchema)
151 def calculateThreshold(self, exposure, seed, sigma=None, minFractionSourcesFactor=1.0,
152 isBgTweak=False, nPixMaskErode=None, maxMaskErodeIter=10):
153 """Calculate new threshold
155 This is the main functional addition to the vanilla
156 `SourceDetectionTask`.
158 We identify sky objects and perform forced PSF photometry on
159 them. Using those PSF flux measurements and estimated errors,
160 we set the threshold so that the stdev of the measurements
161 matches the median estimated error.
163 Parameters
164 ----------
165 exposure : `lsst.afw.image.Exposure`
166 Exposure on which we're detecting sources.
167 seed : `int`
168 RNG seed to use for finding sky objects.
169 sigma : `float`, optional
170 Gaussian sigma of smoothing kernel; if not provided,
171 will be deduced from the exposure's PSF.
172 minFractionSourcesFactor : `float`
173 Change the fraction of required sky sources from that set in
174 ``self.config.minFractionSources`` by this factor. NOTE: this
175 is intended for use in the background tweak pass (the detection
176 threshold is much lower there, so many more pixels end up marked
177 as DETECTED or DETECTED_NEGATIVE, leaving less room for sky
178 object placement).
179 isBgTweak : `bool`
180 Set to ``True`` for the background tweak pass (for more helpful
181 log messages).
182 nPixMaskErode : `int`, optional
183 Number of pixels by which to erode the detection masks on each
184 iteration of best-effort sky object placement.
185 maxMaskErodeIter : `int`, optional
186 Maximum number of iterations for the detection mask erosion.
188 Returns
189 -------
190 result : `lsst.pipe.base.Struct`
191 Result struct with components:
193 ``multiplicative``
194 Multiplicative factor to be applied to the
195 configured detection threshold (`float`).
196 ``additive``
197 Additive factor to be applied to the background
198 level (`float`).
200 Raises
201 ------
202 NoWorkFound
203 Raised if the number of good sky sources found is less than the
204 minimum fraction
205 (``self.config.minFractionSources``*``minFractionSourcesFactor``)
206 of the number requested (``self.skyObjects.config.nSources``).
207 """
208 wcsIsNone = exposure.getWcs() is None
209 if wcsIsNone: # create a dummy WCS as needed by ForcedMeasurementTask
210 self.log.info("WCS for exposure is None. Setting a dummy WCS for dynamic detection.")
211 exposure.setWcs(makeSkyWcs(crpix=geom.Point2D(0, 0),
212 crval=geom.SpherePoint(0, 0, geom.degrees),
213 cdMatrix=makeCdMatrix(scale=1e-5*geom.degrees)))
214 minNumSources = int(self.config.minFractionSources*self.skyObjects.config.nSources)
215 # Reduce the number of sky sources required if requested, but ensure
216 # a minumum of 3.
217 if minFractionSourcesFactor != 1.0:
218 minNumSources = max(3, int(minNumSources*minFractionSourcesFactor))
219 fp = self.skyObjects.run(exposure.maskedImage.mask, seed)
221 if self.config.allowMaskErode:
222 detectedMaskPlanes = ["DETECTED", "DETECTED_NEGATIVE"]
223 mask = exposure.maskedImage.mask
224 for nIter in range(maxMaskErodeIter):
225 if nIter > 0:
226 fp = self.skyObjects.run(mask, seed)
227 if len(fp) < int(2*minNumSources): # Allow for measurement failures
228 self.log.info("Current number of sky sources is below 2*minimum required "
229 "(%d < %d, allowing for some subsequent measurement failures). "
230 "Allowing erosion of detected mask planes for sky placement "
231 "nIter: %d [of %d max]",
232 len(fp), 2*minNumSources, nIter, maxMaskErodeIter)
233 if nPixMaskErode is None:
234 if len(fp) == 0:
235 nPixMaskErode = 4
236 elif len(fp) < int(0.75*minNumSources):
237 nPixMaskErode = 2
238 else:
239 nPixMaskErode = 1
240 for maskName in detectedMaskPlanes:
241 # Compute the eroded detection mask plane using SpanSet
242 detectedMaskBit = mask.getPlaneBitMask(maskName)
243 detectedMaskSpanSet = SpanSet.fromMask(mask, detectedMaskBit)
244 detectedMaskSpanSet = detectedMaskSpanSet.eroded(nPixMaskErode)
245 # Clear the detected mask plane
246 detectedMask = mask.getMaskPlane(maskName)
247 mask.clearMaskPlane(detectedMask)
248 # Set the mask plane to the eroded one
249 detectedMaskSpanSet.setMask(mask, detectedMaskBit)
250 else:
251 break
253 skyFootprints = FootprintSet(exposure.getBBox())
254 skyFootprints.setFootprints(fp)
255 table = SourceTable.make(self.skyMeasurement.schema)
256 catalog = SourceCatalog(table)
257 catalog.reserve(len(skyFootprints.getFootprints()))
258 skyFootprints.makeSources(catalog)
259 key = catalog.getCentroidSlot().getMeasKey()
260 for source in catalog:
261 peaks = source.getFootprint().getPeaks()
262 assert len(peaks) == 1
263 source.set(key, peaks[0].getF())
264 # Coordinate covariance is not used, so don't bother calulating it.
265 source.updateCoord(exposure.getWcs(), include_covariance=False)
267 # Forced photometry on sky objects
268 self.skyMeasurement.run(catalog, exposure, catalog, exposure.getWcs())
270 # Calculate new threshold
271 fluxes = catalog["base_PsfFlux_instFlux"]
272 area = catalog["base_PsfFlux_area"]
273 bg = catalog["base_LocalBackground_instFlux"]
275 good = (~catalog["base_PsfFlux_flag"] & ~catalog["base_LocalBackground_flag"]
276 & np.isfinite(fluxes) & np.isfinite(area) & np.isfinite(bg))
278 if good.sum() < minNumSources:
279 if not isBgTweak:
280 msg = (f"Insufficient good sky source flux measurements ({good.sum()} < "
281 f"{minNumSources}) for dynamic threshold calculation.")
282 else:
283 msg = (f"Insufficient good sky source flux measurements ({good.sum()} < "
284 f"{minNumSources}) for background tweak calculation.")
286 nPix = exposure.mask.array.size
287 badPixelMask = lsst.afw.image.Mask.getPlaneBitMask(["NO_DATA", "BAD"])
288 nGoodPix = np.sum(exposure.mask.array & badPixelMask == 0)
289 if nGoodPix/nPix > 0.2:
290 detectedPixelMask = lsst.afw.image.Mask.getPlaneBitMask(["DETECTED", "DETECTED_NEGATIVE"])
291 nDetectedPix = np.sum(exposure.mask.array & detectedPixelMask != 0)
292 msg += (f" However, {nGoodPix}/{nPix} pixels are not marked NO_DATA or BAD, "
293 "so there should be sufficient area to locate suitable sky sources. "
294 f"Note that {nDetectedPix} of {nGoodPix} \"good\" pixels were marked "
295 "as DETECTED or DETECTED_NEGATIVE.")
296 raise InsufficientSourcesError(msg, nGoodPix, nPix)
297 raise InsufficientSourcesError(msg, nGoodPix, nPix)
299 if not isBgTweak:
300 self.log.info("Number of good sky sources used for dynamic detection: %d (of %d requested).",
301 good.sum(), self.skyObjects.config.nSources)
302 else:
303 self.log.info("Number of good sky sources used for dynamic detection background tweak:"
304 " %d (of %d requested).", good.sum(), self.skyObjects.config.nSources)
305 bgMedian = np.median((fluxes/area)[good])
307 lq, uq = np.percentile((fluxes - bg*area)[good], [25.0, 75.0])
308 stdevMeas = 0.741*(uq - lq)
309 medianError = np.median(catalog["base_PsfFlux_instFluxErr"][good])
310 if wcsIsNone:
311 exposure.setWcs(None)
312 return Struct(multiplicative=medianError/stdevMeas, additive=bgMedian)
314 def detectFootprints(self, exposure, doSmooth=True, sigma=None, clearMask=True, expId=None,
315 background=None, backgroundToPhotometricRatio=None):
316 """Detect footprints with a dynamic threshold
318 This varies from the vanilla ``detectFootprints`` method because we
319 do detection three times: first with a high threshold to detect
320 "bright" (both positive and negative, the latter to identify very
321 over-subtracted regions) sources for which we grow the DETECTED and
322 DETECTED_NEGATIVE masks significantly to account for wings. Second,
323 with a low threshold to mask all non-empty regions of the image. These
324 two masks are combined and used to identify regions of sky
325 uncontaminated by objects. A final round of detection is then done
326 with the new calculated threshold.
328 Parameters
329 ----------
330 exposure : `lsst.afw.image.Exposure`
331 Exposure to process; DETECTED{,_NEGATIVE} mask plane will be
332 set in-place.
333 doSmooth : `bool`, optional
334 If True, smooth the image before detection using a Gaussian
335 of width ``sigma``.
336 sigma : `float`, optional
337 Gaussian Sigma of PSF (pixels); used for smoothing and to grow
338 detections; if `None` then measure the sigma of the PSF of the
339 ``exposure``.
340 clearMask : `bool`, optional
341 Clear both DETECTED and DETECTED_NEGATIVE planes before running
342 detection.
343 expId : `int`, optional
344 Exposure identifier, used as a seed for the random number
345 generator. If absent, the seed will be the sum of the image.
346 background : `lsst.afw.math.BackgroundList`, optional
347 Background that was already subtracted from the exposure; will be
348 modified in-place if ``reEstimateBackground=True``.
349 backgroundToPhotometricRatio : `lsst.afw.image.Image`, optional
350 Unused; if set will Raise.
352 Returns
353 -------
354 results : `lsst.pipe.base.Struct`
355 The results `~lsst.pipe.base.Struct` contains:
357 ``positive``
358 Positive polarity footprints.
359 (`lsst.afw.detection.FootprintSet` or `None`)
360 ``negative``
361 Negative polarity footprints.
362 (`lsst.afw.detection.FootprintSet` or `None`)
363 ``numPos``
364 Number of footprints in positive or 0 if detection polarity was
365 negative. (`int`)
366 ``numNeg``
367 Number of footprints in negative or 0 if detection polarity was
368 positive. (`int`)
369 ``background``
370 Re-estimated background. `None` or the input ``background``
371 if ``reEstimateBackground==False``.
372 (`lsst.afw.math.BackgroundList`)
373 ``factor``
374 Multiplication factor applied to the configured detection
375 threshold. (`float`)
376 ``prelim``
377 Results from preliminary detection pass.
378 (`lsst.pipe.base.Struct`)
379 """
380 if backgroundToPhotometricRatio is not None:
381 raise RuntimeError("DynamicDetectionTask does not support backgroundToPhotometricRatio.")
382 maskedImage = exposure.maskedImage
384 if clearMask:
385 self.clearMask(maskedImage.mask)
386 else:
387 oldDetected = maskedImage.mask.array & maskedImage.mask.getPlaneBitMask(["DETECTED",
388 "DETECTED_NEGATIVE"])
389 nPix = maskedImage.mask.array.size
390 badPixelMask = lsst.afw.image.Mask.getPlaneBitMask(["NO_DATA", "BAD"])
391 nGoodPix = np.sum(maskedImage.mask.array & badPixelMask == 0)
392 self.log.info("Number of good data pixels (i.e. not NO_DATA or BAD): {} ({:.1f}% of total)".
393 format(nGoodPix, 100*nGoodPix/nPix))
395 with self.tempWideBackgroundContext(exposure):
396 # Could potentially smooth with a wider kernel than the PSF in
397 # order to better pick up the wings of stars and galaxies, but for
398 # now sticking with the PSF as that's more simple.
399 psf = self.getPsf(exposure, sigma=sigma)
400 convolveResults = self.convolveImage(maskedImage, psf, doSmooth=doSmooth)
402 if self.config.doBrightPrelimDetection:
403 brightDetectedMask = self._computeBrightDetectionMask(maskedImage, convolveResults)
405 middle = convolveResults.middle
406 sigma = convolveResults.sigma
407 prelim = self.applyThreshold(
408 middle, maskedImage.getBBox(), factor=self.config.prelimThresholdFactor,
409 factorNeg=self.config.prelimNegMultiplier*self.config.prelimThresholdFactor
410 )
411 self.finalizeFootprints(
412 maskedImage.mask, prelim, sigma, factor=self.config.prelimThresholdFactor,
413 factorNeg=self.config.prelimNegMultiplier*self.config.prelimThresholdFactor
414 )
415 if self.config.doBrightPrelimDetection:
416 # Combine prelim and bright detection masks for multiplier
417 # determination.
418 maskedImage.mask.array |= brightDetectedMask
420 # Calculate the proper threshold
421 # seed needs to fit in a C++ 'int' so pybind doesn't choke on it
422 seed = (expId if expId is not None else int(maskedImage.image.array.sum())) % (2**31 - 1)
423 threshResults = self.calculateThreshold(exposure, seed, sigma=sigma)
424 factor = threshResults.multiplicative
425 self.log.info("Modifying configured detection threshold by factor %f to %f",
426 factor, factor*self.config.thresholdValue)
428 # Blow away preliminary (low threshold) detection mask
429 self.clearMask(maskedImage.mask)
430 if not clearMask:
431 maskedImage.mask.array |= oldDetected
433 # Rinse and repeat thresholding with new calculated threshold
434 results = self.applyThreshold(middle, maskedImage.getBBox(), factor)
435 results.prelim = prelim
436 results.background = background if background is not None else lsst.afw.math.BackgroundList()
437 if self.config.doTempLocalBackground:
438 self.applyTempLocalBackground(exposure, middle, results)
439 self.finalizeFootprints(maskedImage.mask, results, sigma, factor=factor)
441 self.clearUnwantedResults(maskedImage.mask, results)
443 if self.config.reEstimateBackground:
444 self.reEstimateBackground(maskedImage, results.background)
446 self.display(exposure, results, middle)
448 if self.config.doBackgroundTweak:
449 # Re-do the background tweak after any temporary backgrounds have
450 # been restored.
451 #
452 # But we want to keep any large-scale background (e.g., scattered
453 # light from bright stars) from being selected for sky objects in
454 # the calculation, so do another detection pass without either the
455 # local or wide temporary background subtraction; the DETECTED
456 # pixels will mark the area to ignore.
457 originalMask = maskedImage.mask.array.copy()
458 try:
459 self.clearMask(exposure.mask)
460 convolveResults = self.convolveImage(maskedImage, psf, doSmooth=doSmooth)
461 tweakDetResults = self.applyThreshold(convolveResults.middle, maskedImage.getBBox(), factor)
462 self.finalizeFootprints(maskedImage.mask, tweakDetResults, sigma, factor=factor)
463 bgLevel = self.calculateThreshold(exposure, seed, sigma=sigma, minFractionSourcesFactor=0.5,
464 isBgTweak=True).additive
465 finally:
466 maskedImage.mask.array[:] = originalMask
467 self.tweakBackground(exposure, bgLevel, results.background)
469 return results
471 def tweakBackground(self, exposure, bgLevel, bgList=None):
472 """Modify the background by a constant value
474 Parameters
475 ----------
476 exposure : `lsst.afw.image.Exposure`
477 Exposure for which to tweak background.
478 bgLevel : `float`
479 Background level to remove
480 bgList : `lsst.afw.math.BackgroundList`, optional
481 List of backgrounds to append to.
483 Returns
484 -------
485 bg : `lsst.afw.math.BackgroundMI`
486 Constant background model.
487 """
488 self.log.info("Tweaking background by %f to match sky photometry", bgLevel)
489 exposure.image -= bgLevel
490 bgStats = lsst.afw.image.MaskedImageF(1, 1)
491 bgStats.set(bgLevel, 0, bgLevel)
492 bg = lsst.afw.math.BackgroundMI(exposure.getBBox(), bgStats)
493 bgData = (bg, lsst.afw.math.Interpolate.LINEAR, lsst.afw.math.REDUCE_INTERP_ORDER,
494 lsst.afw.math.ApproximateControl.UNKNOWN, 0, 0, False)
495 if bgList is not None:
496 bgList.append(bgData)
497 return bg
499 def _computeBrightDetectionMask(self, maskedImage, convolveResults):
500 """Perform an initial bright source detection pass.
502 Perform an initial bright object detection pass using a high detection
503 threshold. The footprints in this pass are grown significantly more
504 than is typical to account for wings around bright sources. The
505 negative polarity detections in this pass help in masking severely
506 over-subtracted regions.
508 A maximum fraction of masked pixel from this pass is ensured via
509 the config ``brightMaskFractionMax``. If the masked pixel fraction is
510 above this value, the detection thresholds here are increased by
511 ``bisectFactor`` in a while loop until the detected masked fraction
512 falls below this value.
514 Parameters
515 ----------
516 maskedImage : `lsst.afw.image.MaskedImage`
517 Masked image on which to run the detection.
518 convolveResults : `lsst.pipe.base.Struct`
519 The results of the self.convolveImage function with attributes:
521 ``middle``
522 Convolved image, without the edges
523 (`lsst.afw.image.MaskedImage`).
524 ``sigma``
525 Gaussian sigma used for the convolution (`float`).
527 Returns
528 -------
529 brightDetectedMask : `numpy.ndarray`
530 Boolean array representing the union of the bright detection pass
531 DETECTED and DETECTED_NEGATIVE masks.
532 """
533 # Initialize some parameters.
534 brightPosFactor = (
535 self.config.prelimThresholdFactor*self.config.brightMultiplier/self.config.bisectFactor
536 )
537 brightNegFactor = self.config.brightNegFactor/self.config.bisectFactor
538 nPix = 1
539 nPixDet = 1
540 nPixDetNeg = 1
541 brightMaskFractionMax = self.config.brightMaskFractionMax
543 # Loop until masked fraction is smaller than
544 # brightMaskFractionMax, increasing the thresholds by
545 # config.bisectFactor on each iteration (rarely necessary
546 # for current defaults).
547 while nPixDetNeg/nPix > brightMaskFractionMax or nPixDet/nPix > brightMaskFractionMax:
548 self.clearMask(maskedImage.mask)
549 brightPosFactor *= self.config.bisectFactor
550 brightNegFactor *= self.config.bisectFactor
551 prelimBright = self.applyThreshold(convolveResults.middle, maskedImage.getBBox(),
552 factor=brightPosFactor, factorNeg=brightNegFactor)
553 self.finalizeFootprints(
554 maskedImage.mask, prelimBright, convolveResults.sigma*self.config.brightGrowFactor,
555 factor=brightPosFactor, factorNeg=brightNegFactor
556 )
557 # Check that not too many pixels got masked.
558 nPix = maskedImage.mask.array.size
559 nPixDet = countMaskedPixels(maskedImage, "DETECTED")
560 self.log.info("Number (%) of bright DETECTED pix: {} ({:.1f}%)".
561 format(nPixDet, 100*nPixDet/nPix))
562 nPixDetNeg = countMaskedPixels(maskedImage, "DETECTED_NEGATIVE")
563 self.log.info("Number (%) of bright DETECTED_NEGATIVE pix: {} ({:.1f}%)".
564 format(nPixDetNeg, 100*nPixDetNeg/nPix))
565 if nPixDetNeg/nPix > brightMaskFractionMax or nPixDet/nPix > brightMaskFractionMax:
566 self.log.warn("Too high a fraction (%.1f > %.1f) of pixels were masked with current "
567 "\"bright\" detection round thresholds. Increasing by a factor of %f "
568 "and trying again.", max(nPixDetNeg, nPixDet)/nPix,
569 brightMaskFractionMax, self.config.bisectFactor)
571 # Save the mask planes from the "bright" detection round, then
572 # clear them before moving on to the "prelim" detection phase.
573 brightDetectedMask = (maskedImage.mask.array
574 & maskedImage.mask.getPlaneBitMask(["DETECTED", "DETECTED_NEGATIVE"]))
575 self.clearMask(maskedImage.mask)
576 return brightDetectedMask
579def countMaskedPixels(maskedIm, maskPlane):
580 """Count the number of pixels in a given mask plane.
582 Parameters
583 ----------
584 maskedIm : `lsst.afw.image.MaskedImage`
585 Masked image to examine.
586 maskPlane : `str`
587 Name of the mask plane to examine.
589 Returns
590 -------
591 nPixMasked : `int`
592 Number of pixels with ``maskPlane`` bit set.
593 """
594 maskBit = maskedIm.mask.getPlaneBitMask(maskPlane)
595 nPixMasked = np.sum(np.bitwise_and(maskedIm.mask.array, maskBit))/maskBit
596 return nPixMasked