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