118 isBgTweak=False, nPixMaskErode=None, maxMaskErodeIter=10):
119 """Calculate new threshold
121 This is the main functional addition to the vanilla
122 `SourceDetectionTask`.
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.
131 exposure : `lsst.afw.image.Exposure`
132 Exposure on which we're detecting sources.
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
146 Set to ``True`` for the background tweak pass (for more helpful
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.
156 result : `lsst.pipe.base.Struct`
157 Result struct with components:
160 Multiplicative factor to be applied to the
161 configured detection threshold (`float`).
163 Additive factor to be applied to the background
169 Raised if the number of good sky sources found is less than the
171 (``self.config.minFractionSources``*``minFractionSourcesFactor``)
172 of the number requested (``self.skyObjects.config.nSources``).
174 wcsIsNone = exposure.getWcs()
is None
176 self.log.info(
"WCS for exposure is None. Setting a dummy WCS for dynamic detection.")
179 cdMatrix=makeCdMatrix(scale=1e-5*geom.degrees)))
180 minNumSources = int(self.config.minFractionSources*self.skyObjects.config.nSources)
183 if minFractionSourcesFactor != 1.0:
184 minNumSources = max(3, int(minNumSources*minFractionSourcesFactor))
185 fp = self.skyObjects.run(exposure.maskedImage.mask, seed)
187 if self.config.allowMaskErode:
188 detectedMaskPlanes = [
"DETECTED",
"DETECTED_NEGATIVE"]
189 mask = exposure.maskedImage.mask
190 for nIter
in range(maxMaskErodeIter):
192 fp = self.skyObjects.run(mask, seed)
193 if len(fp) < int(2*minNumSources):
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:
202 elif len(fp) < int(0.75*minNumSources):
206 for maskName
in detectedMaskPlanes:
208 detectedMaskBit = mask.getPlaneBitMask(maskName)
209 detectedMaskSpanSet = SpanSet.fromMask(mask, detectedMaskBit)
210 detectedMaskSpanSet = detectedMaskSpanSet.eroded(nPixMaskErode)
212 detectedMask = mask.getMaskPlane(maskName)
213 mask.clearMaskPlane(detectedMask)
215 detectedMaskSpanSet.setMask(mask, detectedMaskBit)
220 skyFootprints.setFootprints(fp)
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())
231 source.updateCoord(exposure.getWcs(), include_covariance=
False)
234 self.
skyMeasurement.run(catalog, exposure, catalog, exposure.getWcs())
237 fluxes = catalog[
"base_PsfFlux_instFlux"]
238 area = catalog[
"base_PsfFlux_area"]
239 bg = catalog[
"base_LocalBackground_instFlux"]
241 good = (~catalog[
"base_PsfFlux_flag"] & ~catalog[
"base_LocalBackground_flag"]
242 & np.isfinite(fluxes) & np.isfinite(area) & np.isfinite(bg))
244 if good.sum() < minNumSources:
246 msg = (f
"Insufficient good sky source flux measurements ({good.sum()} < "
247 f
"{minNumSources}) for dynamic threshold calculation.")
249 msg = (f
"Insufficient good sky source flux measurements ({good.sum()} < "
250 f
"{minNumSources}) for background tweak calculation.")
252 nPix = exposure.mask.array.size
254 nGoodPix = np.sum(exposure.mask.array & badPixelMask == 0)
255 if nGoodPix/nPix > 0.2:
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)
266 self.log.info(
"Number of good sky sources used for dynamic detection: %d (of %d requested).",
267 good.sum(), self.skyObjects.config.nSources)
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])
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])
277 exposure.setWcs(
None)
278 return Struct(multiplicative=medianError/stdevMeas, additive=bgMedian)
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
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.
296 exposure : `lsst.afw.image.Exposure`
297 Exposure to process; DETECTED{,_NEGATIVE} mask plane will be
299 doSmooth : `bool`, optional
300 If True, smooth the image before detection using a Gaussian
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
306 clearMask : `bool`, optional
307 Clear both DETECTED and DETECTED_NEGATIVE planes before running
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.
320 results : `lsst.pipe.base.Struct`
321 The results `~lsst.pipe.base.Struct` contains:
324 Positive polarity footprints.
325 (`lsst.afw.detection.FootprintSet` or `None`)
327 Negative polarity footprints.
328 (`lsst.afw.detection.FootprintSet` or `None`)
330 Number of footprints in positive or 0 if detection polarity was
333 Number of footprints in negative or 0 if detection polarity was
336 Re-estimated background. `None` or the input ``background``
337 if ``reEstimateBackground==False``.
338 (`lsst.afw.math.BackgroundList`)
340 Multiplication factor applied to the configured detection
343 Results from preliminary detection pass.
344 (`lsst.pipe.base.Struct`)
346 if backgroundToPhotometricRatio
is not None:
347 raise RuntimeError(
"DynamicDetectionTask does not support backgroundToPhotometricRatio.")
348 maskedImage = exposure.maskedImage
353 oldDetected = maskedImage.mask.array & maskedImage.mask.getPlaneBitMask([
"DETECTED",
354 "DETECTED_NEGATIVE"])
355 nPix = maskedImage.mask.array.size
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))
365 psf = self.
getPsf(exposure, sigma=sigma)
366 convolveResults = self.
convolveImage(maskedImage, psf, doSmooth=doSmooth)
368 if self.config.doBrightPrelimDetection:
371 middle = convolveResults.middle
372 sigma = convolveResults.sigma
374 middle, maskedImage.getBBox(), factor=self.config.prelimThresholdFactor,
375 factorNeg=self.config.prelimNegMultiplier*self.config.prelimThresholdFactor
378 maskedImage.mask, prelim, sigma, factor=self.config.prelimThresholdFactor,
379 factorNeg=self.config.prelimNegMultiplier*self.config.prelimThresholdFactor
381 if self.config.doBrightPrelimDetection:
384 maskedImage.mask.array |= brightDetectedMask
388 seed = (expId
if expId
is not None else int(maskedImage.image.array.sum())) % (2**31 - 1)
390 factor = threshResults.multiplicative
391 self.log.info(
"Modifying configured detection threshold by factor %f to %f",
392 factor, factor*self.config.thresholdValue)
397 maskedImage.mask.array |= oldDetected
400 results = self.
applyThreshold(middle, maskedImage.getBBox(), factor)
401 results.prelim = prelim
403 if self.config.doTempLocalBackground:
409 if self.config.reEstimateBackground:
412 self.
display(exposure, results, middle)
414 if self.config.doBackgroundTweak:
423 originalMask = maskedImage.mask.array.copy()
426 convolveResults = self.
convolveImage(maskedImage, psf, doSmooth=doSmooth)
427 tweakDetResults = self.
applyThreshold(convolveResults.middle, maskedImage.getBBox(), factor)
429 bgLevel = self.
calculateThreshold(exposure, seed, sigma=sigma, minFractionSourcesFactor=0.5,
430 isBgTweak=
True).additive
432 maskedImage.mask.array[:] = originalMask
438 """Modify the background by a constant value
442 exposure : `lsst.afw.image.Exposure`
443 Exposure for which to tweak background.
445 Background level to remove
446 bgList : `lsst.afw.math.BackgroundList`, optional
447 List of backgrounds to append to.
451 bg : `lsst.afw.math.BackgroundMI`
452 Constant background model.
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)
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)
466 """Perform an initial bright source detection pass.
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.
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.
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:
488 Convolved image, without the edges
489 (`lsst.afw.image.MaskedImage`).
491 Gaussian sigma used for the convolution (`float`).
495 brightDetectedMask : `numpy.ndarray`
496 Boolean array representing the union of the bright detection pass
497 DETECTED and DETECTED_NEGATIVE masks.
501 self.config.prelimThresholdFactor*self.config.brightMultiplier/self.config.bisectFactor
503 brightNegFactor = self.config.brightNegFactor/self.config.bisectFactor
507 brightMaskFractionMax = self.config.brightMaskFractionMax
513 while nPixDetNeg/nPix > brightMaskFractionMax
or nPixDet/nPix > brightMaskFractionMax:
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
524 nPix = maskedImage.mask.array.size
526 self.log.info(
"Number (%) of bright DETECTED pix: {} ({:.1f}%)".
527 format(nPixDet, 100*nPixDet/nPix))
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)
539 brightDetectedMask = (maskedImage.mask.array
540 & maskedImage.mask.getPlaneBitMask([
"DETECTED",
"DETECTED_NEGATIVE"]))
542 return brightDetectedMask