LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
detection.py
Go to the documentation of this file.
2# LSST Data Management System
3#
4# Copyright 2008-2017 AURA/LSST.
5#
6# This product includes software developed by the
7# LSST Project (http://www.lsst.org/).
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the LSST License Statement and
20# the GNU General Public License along with this program. If not,
21# see <https://www.lsstcorp.org/LegalNotices/>.
22#
23
24__all__ = ("SourceDetectionConfig", "SourceDetectionTask", "addExposures")
25
26from contextlib import contextmanager
27
28import numpy as np
29
30import lsst.geom
31import lsst.afw.display as afwDisplay
32import lsst.afw.detection as afwDet
33import lsst.afw.geom as afwGeom
34import lsst.afw.image as afwImage
35import lsst.afw.math as afwMath
36import lsst.afw.table as afwTable
37import lsst.pex.config as pexConfig
38import lsst.pipe.base as pipeBase
39from lsst.utils.timer import timeMethod
40from .subtractBackground import SubtractBackgroundTask
41
42
43class SourceDetectionConfig(pexConfig.Config):
44 """Configuration parameters for the SourceDetectionTask
45 """
46 minPixels = pexConfig.RangeField(
47 doc="detected sources with fewer than the specified number of pixels will be ignored",
48 dtype=int, optional=False, default=1, min=0,
49 )
50 isotropicGrow = pexConfig.Field(
51 doc="Grow pixels as isotropically as possible? If False, use a Manhattan metric instead.",
52 dtype=bool, default=True,
53 )
54 combinedGrow = pexConfig.Field(
55 doc="Grow all footprints at the same time? This allows disconnected footprints to merge.",
56 dtype=bool, default=True,
57 )
58 nSigmaToGrow = pexConfig.Field(
59 doc="Grow detections by nSigmaToGrow * [PSF RMS width]; if 0 then do not grow",
60 dtype=float, default=2.4, # 2.4 pixels/sigma is roughly one pixel/FWHM
61 )
62 returnOriginalFootprints = pexConfig.Field(
63 doc="Grow detections to set the image mask bits, but return the original (not-grown) footprints",
64 dtype=bool, optional=False, default=False,
65 )
66 thresholdValue = pexConfig.RangeField(
67 doc="Threshold for detecting footprints; exact meaning and units depend on thresholdType.",
68 dtype=float, optional=False, default=5.0, min=0.0,
69 )
70 includeThresholdMultiplier = pexConfig.RangeField(
71 doc="Multiplier on thresholdValue for whether a source is included in the output catalog."
72 " For example, thresholdValue=5, includeThresholdMultiplier=10, thresholdType='pixel_stdev' "
73 "results in a catalog of sources at >50 sigma with the detection mask and footprints "
74 "including pixels >5 sigma.",
75 dtype=float, default=1.0, min=0.0,
76 )
77 thresholdType = pexConfig.ChoiceField(
78 doc="Specifies the meaning of thresholdValue.",
79 dtype=str, optional=False, default="pixel_stdev",
80 allowed={
81 "variance": "threshold applied to image variance",
82 "stdev": "threshold applied to image std deviation",
83 "value": "threshold applied to image value",
84 "pixel_stdev": "threshold applied to per-pixel std deviation",
85 },
86 )
87 thresholdPolarity = pexConfig.ChoiceField(
88 doc="Specifies whether to detect positive, or negative sources, or both.",
89 dtype=str, optional=False, default="positive",
90 allowed={
91 "positive": "detect only positive sources",
92 "negative": "detect only negative sources",
93 "both": "detect both positive and negative sources",
94 },
95 )
96 adjustBackground = pexConfig.Field(
97 dtype=float,
98 doc="Fiddle factor to add to the background; debugging only",
99 default=0.0,
100 )
101 reEstimateBackground = pexConfig.Field(
102 dtype=bool,
103 doc="Estimate the background again after final source detection?",
104 default=True, optional=False,
105 )
106 background = pexConfig.ConfigurableField(
107 doc="Background re-estimation; ignored if reEstimateBackground false",
108 target=SubtractBackgroundTask,
109 )
110 tempLocalBackground = pexConfig.ConfigurableField(
111 doc=("A local (small-scale), temporary background estimation step run between "
112 "detecting above-threshold regions and detecting the peaks within "
113 "them; used to avoid detecting spuerious peaks in the wings."),
114 target=SubtractBackgroundTask,
115 )
116 doTempLocalBackground = pexConfig.Field(
117 dtype=bool,
118 doc="Enable temporary local background subtraction? (see tempLocalBackground)",
119 default=True,
120 )
121 tempWideBackground = pexConfig.ConfigurableField(
122 doc=("A wide (large-scale) background estimation and removal before footprint and peak detection. "
123 "It is added back into the image after detection. The purpose is to suppress very large "
124 "footprints (e.g., from large artifacts) that the deblender may choke on."),
125 target=SubtractBackgroundTask,
126 )
127 doTempWideBackground = pexConfig.Field(
128 dtype=bool,
129 doc="Do temporary wide (large-scale) background subtraction before footprint detection?",
130 default=False,
131 )
132 nPeaksMaxSimple = pexConfig.Field(
133 dtype=int,
134 doc=("The maximum number of peaks in a Footprint before trying to "
135 "replace its peaks using the temporary local background"),
136 default=1,
137 )
138 nSigmaForKernel = pexConfig.Field(
139 dtype=float,
140 doc=("Multiple of PSF RMS size to use for convolution kernel bounding box size; "
141 "note that this is not a half-size. The size will be rounded up to the nearest odd integer"),
142 default=7.0,
143 )
144 statsMask = pexConfig.ListField(
145 dtype=str,
146 doc="Mask planes to ignore when calculating statistics of image (for thresholdType=stdev)",
147 default=['BAD', 'SAT', 'EDGE', 'NO_DATA'],
148 )
149 excludeMaskPlanes = lsst.pex.config.ListField(
150 dtype=str,
151 default=[],
152 doc="Mask planes to exclude when detecting sources."
153 )
154
155 def setDefaults(self):
156 self.tempLocalBackground.binSize = 64
157 self.tempLocalBackground.algorithm = "AKIMA_SPLINE"
158 self.tempLocalBackground.useApprox = False
159 # Background subtraction to remove a large-scale background (e.g., scattered light); restored later.
160 # Want to keep it from exceeding the deblender size limit of 1 Mpix, so half that is reasonable.
161 self.tempWideBackground.binSize = 512
162 self.tempWideBackground.algorithm = "AKIMA_SPLINE"
163 self.tempWideBackground.useApprox = False
164 # Ensure we can remove even bright scattered light that is DETECTED
165 for maskPlane in ("DETECTED", "DETECTED_NEGATIVE"):
166 if maskPlane in self.tempWideBackground.ignoredPixelMask:
167 self.tempWideBackground.ignoredPixelMask.remove(maskPlane)
168
169
170class SourceDetectionTask(pipeBase.Task):
171 """Detect peaks and footprints of sources in an image.
172
173 This task expects the image to have been background subtracted first.
174 Running detection on images with a non-zero-centered background may result
175 in a single source detected on the entire image containing thousands of
176 peaks, or other pathological outputs.
177
178 Parameters
179 ----------
180 schema : `lsst.afw.table.Schema`
181 Schema object used to create the output `lsst.afw.table.SourceCatalog`
182 **kwds
183 Keyword arguments passed to `lsst.pipe.base.Task.__init__`
184
185 If schema is not None and configured for 'both' detections,
186 a 'flags.negative' field will be added to label detections made with a
187 negative threshold.
188
189 Notes
190 -----
191 This task convolves the image with a Gaussian approximation to the PSF,
192 matched to the sigma of the input exposure, because this is separable and
193 fast. The PSF would have to be very non-Gaussian or non-circular for this
194 approximation to have a significant impact on the signal-to-noise of the
195 detected sources.
196 """
197 ConfigClass = SourceDetectionConfig
198 _DefaultName = "sourceDetection"
199
200 def __init__(self, schema=None, **kwds):
201 pipeBase.Task.__init__(self, **kwds)
202 if schema is not None and self.config.thresholdPolarity == "both":
203 self.negativeFlagKey = schema.addField(
204 "flags_negative", type="Flag",
205 doc="set if source was detected as significantly negative"
206 )
207 else:
208 if self.config.thresholdPolarity == "both":
209 self.log.warning("Detection polarity set to 'both', but no flag will be "
210 "set to distinguish between positive and negative detections")
211 self.negativeFlagKey = None
212 if self.config.reEstimateBackground:
213 self.makeSubtask("background")
214 if self.config.doTempLocalBackground:
215 self.makeSubtask("tempLocalBackground")
216 if self.config.doTempWideBackground:
217 self.makeSubtask("tempWideBackground")
218
219 @timeMethod
220 def run(self, table, exposure, doSmooth=True, sigma=None, clearMask=True, expId=None,
221 background=None):
222 r"""Detect sources and return catalog(s) of detections.
223
224 Parameters
225 ----------
226 table : `lsst.afw.table.SourceTable`
227 Table object that will be used to create the SourceCatalog.
228 exposure : `lsst.afw.image.Exposure`
229 Exposure to process; DETECTED mask plane will be set in-place.
230 doSmooth : `bool`, optional
231 If True, smooth the image before detection using a Gaussian of width
232 ``sigma``, or the measured PSF width. Set to False when running on
233 e.g. a pre-convolved image, or a mask plane.
234 sigma : `float`, optional
235 Sigma of PSF (pixels); used for smoothing and to grow detections;
236 if None then measure the sigma of the PSF of the exposure
237 clearMask : `bool`, optional
238 Clear DETECTED{,_NEGATIVE} planes before running detection.
239 expId : `int`, optional
240 Exposure identifier; unused by this implementation, but used for
241 RNG seed by subclasses.
242 background : `lsst.afw.math.BackgroundList`, optional
243 Background that was already subtracted from the exposure; will be
244 modified in-place if ``reEstimateBackground=True``.
245
246 Returns
247 -------
248 result : `lsst.pipe.base.Struct`
249 The `~lsst.pipe.base.Struct` contains:
250
251 ``sources``
252 Detected sources on the exposure.
253 (`lsst.afw.table.SourceCatalog`)
254 ``positive``
255 Positive polarity footprints.
256 (`lsst.afw.detection.FootprintSet` or `None`)
257 ``negative``
258 Negative polarity footprints.
259 (`lsst.afw.detection.FootprintSet` or `None`)
260 ``numPos``
261 Number of footprints in positive or 0 if detection polarity was
262 negative. (`int`)
263 ``numNeg``
264 Number of footprints in negative or 0 if detection polarity was
265 positive. (`int`)
266 ``background``
267 Re-estimated background. `None` if
268 ``reEstimateBackground==False``.
269 (`lsst.afw.math.BackgroundList`)
270 ``factor``
271 Multiplication factor applied to the configured detection
272 threshold. (`float`)
273
274 Raises
275 ------
276 ValueError
277 Raised if flags.negative is needed, but isn't in table's schema.
278 lsst.pipe.base.TaskError
279 Raised if sigma=None, doSmooth=True and the exposure has no PSF.
280
281 Notes
282 -----
283 If you want to avoid dealing with Sources and Tables, you can use
284 `detectFootprints()` to just get the
285 `~lsst.afw.detection.FootprintSet`\s.
286 """
287 if self.negativeFlagKey is not None and self.negativeFlagKey not in table.getSchema():
288 raise ValueError("Table has incorrect Schema")
289 results = self.detectFootprints(exposure=exposure, doSmooth=doSmooth, sigma=sigma,
290 clearMask=clearMask, expId=expId, background=background)
291 sources = afwTable.SourceCatalog(table)
292 sources.reserve(results.numPos + results.numNeg)
293 if results.negative:
294 results.negative.makeSources(sources)
295 if self.negativeFlagKey:
296 for record in sources:
297 record.set(self.negativeFlagKey, True)
298 if results.positive:
299 results.positive.makeSources(sources)
300 results.sources = sources
301 return results
302
303 def display(self, exposure, results, convolvedImage=None):
304 """Display detections if so configured
305
306 Displays the ``exposure`` in frame 0, overlays the detection peaks.
307
308 Requires that ``lsstDebug`` has been set up correctly, so that
309 ``lsstDebug.Info("lsst.meas.algorithms.detection")`` evaluates `True`.
310
311 If the ``convolvedImage`` is non-`None` and
312 ``lsstDebug.Info("lsst.meas.algorithms.detection") > 1``, the
313 ``convolvedImage`` will be displayed in frame 1.
314
315 Parameters
316 ----------
317 exposure : `lsst.afw.image.Exposure`
318 Exposure to display, on which will be plotted the detections.
319 results : `lsst.pipe.base.Struct`
320 Results of the 'detectFootprints' method, containing positive and
321 negative footprints (which contain the peak positions that we will
322 plot). This is a `Struct` with ``positive`` and ``negative``
323 elements that are of type `lsst.afw.detection.FootprintSet`.
324 convolvedImage : `lsst.afw.image.Image`, optional
325 Convolved image used for thresholding.
326 """
327 try:
328 import lsstDebug
329 display = lsstDebug.Info(__name__).display
330 except ImportError:
331 try:
332 display
333 except NameError:
334 display = False
335 if not display:
336 return
337
338 afwDisplay.setDefaultMaskTransparency(75)
339
340 disp0 = afwDisplay.Display(frame=0)
341 disp0.mtv(exposure, title="detection")
342
343 def plotPeaks(fps, ctype):
344 if fps is None:
345 return
346 with disp0.Buffering():
347 for fp in fps.getFootprints():
348 for pp in fp.getPeaks():
349 disp0.dot("+", pp.getFx(), pp.getFy(), ctype=ctype)
350 plotPeaks(results.positive, "yellow")
351 plotPeaks(results.negative, "red")
352
353 if convolvedImage and display > 1:
354 disp1 = afwDisplay.Display(frame=1)
355 disp1.mtv(convolvedImage, title="PSF smoothed")
356
357 disp2 = afwDisplay.Display(frame=2)
358 disp2.mtv(afwImage.ImageF(np.sqrt(exposure.variance.array)), title="stddev")
359
360 def applyTempLocalBackground(self, exposure, middle, results):
361 """Apply a temporary local background subtraction
362
363 This temporary local background serves to suppress noise fluctuations
364 in the wings of bright objects.
365
366 Peaks in the footprints will be updated.
367
368 Parameters
369 ----------
370 exposure : `lsst.afw.image.Exposure`
371 Exposure for which to fit local background.
372 middle : `lsst.afw.image.MaskedImage`
373 Convolved image on which detection will be performed
374 (typically smaller than ``exposure`` because the
375 half-kernel has been removed around the edges).
376 results : `lsst.pipe.base.Struct`
377 Results of the 'detectFootprints' method, containing positive and
378 negative footprints (which contain the peak positions that we will
379 plot). This is a `Struct` with ``positive`` and ``negative``
380 elements that are of type `lsst.afw.detection.FootprintSet`.
381 """
382 # Subtract the local background from the smoothed image. Since we
383 # never use the smoothed again we don't need to worry about adding
384 # it back in.
385 bg = self.tempLocalBackground.fitBackground(exposure.getMaskedImage())
386 bgImage = bg.getImageF(self.tempLocalBackground.config.algorithm,
387 self.tempLocalBackground.config.undersampleStyle)
388 middle -= bgImage.Factory(bgImage, middle.getBBox())
389 if self.config.thresholdPolarity != "negative":
390 results.positiveThreshold = self.makeThreshold(middle, "positive")
391 self.updatePeaks(results.positive, middle, results.positiveThreshold)
392 if self.config.thresholdPolarity != "positive":
393 results.negativeThreshold = self.makeThreshold(middle, "negative")
394 self.updatePeaks(results.negative, middle, results.negativeThreshold)
395
396 def clearMask(self, mask):
397 """Clear the DETECTED and DETECTED_NEGATIVE mask planes.
398
399 Removes any previous detection mask in preparation for a new
400 detection pass.
401
402 Parameters
403 ----------
404 mask : `lsst.afw.image.Mask`
405 Mask to be cleared.
406 """
407 mask &= ~(mask.getPlaneBitMask("DETECTED") | mask.getPlaneBitMask("DETECTED_NEGATIVE"))
408
409 def calculateKernelSize(self, sigma):
410 """Calculate the size of the smoothing kernel.
411
412 Uses the ``nSigmaForKernel`` configuration parameter. Note
413 that that is the full width of the kernel bounding box
414 (so a value of 7 means 3.5 sigma on either side of center).
415 The value will be rounded up to the nearest odd integer.
416
417 Parameters
418 ----------
419 sigma : `float`
420 Gaussian sigma of smoothing kernel.
421
422 Returns
423 -------
424 size : `int`
425 Size of the smoothing kernel.
426 """
427 return (int(sigma * self.config.nSigmaForKernel + 0.5)//2)*2 + 1 # make sure it is odd
428
429 def getPsf(self, exposure, sigma=None):
430 """Create a single Gaussian PSF for an exposure.
431
432 If ``sigma`` is provided, we make a `~lsst.afw.detection.GaussianPsf`
433 with that, otherwise use the sigma from the psf of the ``exposure`` to
434 make the `~lsst.afw.detection.GaussianPsf`.
435
436 Parameters
437 ----------
438 exposure : `lsst.afw.image.Exposure`
439 Exposure from which to retrieve the PSF.
440 sigma : `float`, optional
441 Gaussian sigma to use if provided.
442
443 Returns
444 -------
445 psf : `lsst.afw.detection.GaussianPsf`
446 PSF to use for detection.
447
448 Raises
449 ------
450 RuntimeError
451 Raised if ``sigma`` is not provided and ``exposure`` does not
452 contain a ``Psf`` object.
453 """
454 if sigma is None:
455 psf = exposure.getPsf()
456 if psf is None:
457 raise RuntimeError("Unable to determine PSF to use for detection: no sigma provided")
458 sigma = psf.computeShape(psf.getAveragePosition()).getDeterminantRadius()
459 size = self.calculateKernelSize(sigma)
460 psf = afwDet.GaussianPsf(size, size, sigma)
461 return psf
462
463 def convolveImage(self, maskedImage, psf, doSmooth=True):
464 """Convolve the image with the PSF.
465
466 We convolve the image with a Gaussian approximation to the PSF,
467 because this is separable and therefore fast. It's technically a
468 correlation rather than a convolution, but since we use a symmetric
469 Gaussian there's no difference.
470
471 The convolution can be disabled with ``doSmooth=False``. If we do
472 convolve, we mask the edges as ``EDGE`` and return the convolved image
473 with the edges removed. This is because we can't convolve the edges
474 because the kernel would extend off the image.
475
476 Parameters
477 ----------
478 maskedImage : `lsst.afw.image.MaskedImage`
479 Image to convolve.
480 psf : `lsst.afw.detection.Psf`
481 PSF to convolve with (actually with a Gaussian approximation
482 to it).
483 doSmooth : `bool`
484 Actually do the convolution? Set to False when running on
485 e.g. a pre-convolved image, or a mask plane.
486
487 Returns
488 -------
489 results : `lsst.pipe.base.Struct`
490 The `~lsst.pipe.base.Struct` contains:
491
492 ``middle``
493 Convolved image, without the edges. (`lsst.afw.image.MaskedImage`)
494 ``sigma``
495 Gaussian sigma used for the convolution. (`float`)
496 """
497 self.metadata["doSmooth"] = doSmooth
498 sigma = psf.computeShape(psf.getAveragePosition()).getDeterminantRadius()
499 self.metadata["sigma"] = sigma
500
501 if not doSmooth:
502 middle = maskedImage.Factory(maskedImage, deep=True)
503 return pipeBase.Struct(middle=middle, sigma=sigma)
504
505 # Smooth using a Gaussian (which is separable, hence fast) of width sigma
506 # Make a SingleGaussian (separable) kernel with the 'sigma'
507 kWidth = self.calculateKernelSize(sigma)
508 self.metadata["smoothingKernelWidth"] = kWidth
509 gaussFunc = afwMath.GaussianFunction1D(sigma)
510 gaussKernel = afwMath.SeparableKernel(kWidth, kWidth, gaussFunc, gaussFunc)
511
512 convolvedImage = maskedImage.Factory(maskedImage.getBBox())
513
514 afwMath.convolve(convolvedImage, maskedImage, gaussKernel, afwMath.ConvolutionControl())
515
516 # Only search psf-smoothed part of frame
517 goodBBox = gaussKernel.shrinkBBox(convolvedImage.getBBox())
518 middle = convolvedImage.Factory(convolvedImage, goodBBox, afwImage.PARENT, False)
519
520 # Mark the parts of the image outside goodBBox as EDGE
521 self.setEdgeBits(maskedImage, goodBBox, maskedImage.getMask().getPlaneBitMask("EDGE"))
522
523 return pipeBase.Struct(middle=middle, sigma=sigma)
524
525 def applyThreshold(self, middle, bbox, factor=1.0, factorNeg=None):
526 r"""Apply thresholds to the convolved image
527
528 Identifies `~lsst.afw.detection.Footprint`\s, both positive and negative.
529 The threshold can be modified by the provided multiplication
530 ``factor``.
531
532 Parameters
533 ----------
534 middle : `lsst.afw.image.MaskedImage`
535 Convolved image to threshold.
536 bbox : `lsst.geom.Box2I`
537 Bounding box of unconvolved image.
538 factor : `float`
539 Multiplier for the configured threshold.
540 factorNeg : `float` or `None`
541 Multiplier for the configured threshold for negative detection polarity.
542 If `None`, will be set equal to ``factor`` (i.e. equal to the factor used
543 for positive detection polarity).
544
545 Returns
546 -------
547 results : `lsst.pipe.base.Struct`
548 The `~lsst.pipe.base.Struct` contains:
549
550 ``positive``
551 Positive detection footprints, if configured.
552 (`lsst.afw.detection.FootprintSet` or `None`)
553 ``negative``
554 Negative detection footprints, if configured.
555 (`lsst.afw.detection.FootprintSet` or `None`)
556 ``factor``
557 Multiplier for the configured threshold.
558 (`float`)
559 ``factorNeg``
560 Multiplier for the configured threshold for negative detection polarity.
561 (`float`)
562 """
563 if factorNeg is None:
564 factorNeg = factor
565 self.log.info("Setting factor for negative detections equal to that for positive "
566 "detections: %f", factor)
567 results = pipeBase.Struct(positive=None, negative=None, factor=factor, factorNeg=factorNeg,
568 positiveThreshold=None, negativeThreshold=None)
569 # Detect the Footprints (peaks may be replaced if doTempLocalBackground)
570 if self.config.reEstimateBackground or self.config.thresholdPolarity != "negative":
571 results.positiveThreshold = self.makeThreshold(middle, "positive", factor=factor)
572 results.positive = afwDet.FootprintSet(
573 middle,
574 results.positiveThreshold,
575 "DETECTED",
576 self.config.minPixels
577 )
578 results.positive.setRegion(bbox)
579 if self.config.reEstimateBackground or self.config.thresholdPolarity != "positive":
580 results.negativeThreshold = self.makeThreshold(middle, "negative", factor=factorNeg)
581 results.negative = afwDet.FootprintSet(
582 middle,
583 results.negativeThreshold,
584 "DETECTED_NEGATIVE",
585 self.config.minPixels
586 )
587 results.negative.setRegion(bbox)
588
589 return results
590
591 def finalizeFootprints(self, mask, results, sigma, factor=1.0, factorNeg=None):
592 """Finalize the detected footprints.
593
594 Grow the footprints, set the ``DETECTED`` and ``DETECTED_NEGATIVE``
595 mask planes, and log the results.
596
597 ``numPos`` (number of positive footprints), ``numPosPeaks`` (number
598 of positive peaks), ``numNeg`` (number of negative footprints),
599 ``numNegPeaks`` (number of negative peaks) entries are added to the
600 ``results`` struct.
601
602 Parameters
603 ----------
604 mask : `lsst.afw.image.Mask`
605 Mask image on which to flag detected pixels.
606 results : `lsst.pipe.base.Struct`
607 Struct of detection results, including ``positive`` and
608 ``negative`` entries; modified.
609 sigma : `float`
610 Gaussian sigma of PSF.
611 factor : `float`
612 Multiplier for the configured threshold. Note that this is only
613 used here for logging purposes.
614 factorNeg : `float` or `None`
615 Multiplier used for the negative detection polarity threshold.
616 If `None`, a factor equal to ``factor`` (i.e. equal to the one used
617 for positive detection polarity) is assumed. Note that this is only
618 used here for logging purposes.
619 """
620 factorNeg = factor if factorNeg is None else factorNeg
621 for polarity, maskName in (("positive", "DETECTED"), ("negative", "DETECTED_NEGATIVE")):
622 fpSet = getattr(results, polarity)
623 if fpSet is None:
624 continue
625 if self.config.nSigmaToGrow > 0:
626 nGrow = int((self.config.nSigmaToGrow * sigma) + 0.5)
627 self.metadata["nGrow"] = nGrow
628 if self.config.combinedGrow:
629 fpSet = afwDet.FootprintSet(fpSet, nGrow, self.config.isotropicGrow)
630 else:
631 stencil = (afwGeom.Stencil.CIRCLE if self.config.isotropicGrow else
632 afwGeom.Stencil.MANHATTAN)
633 for fp in fpSet:
634 fp.dilate(nGrow, stencil)
635 fpSet.setMask(mask, maskName)
636 if not self.config.returnOriginalFootprints:
637 setattr(results, polarity, fpSet)
638
639 results.numPos = 0
640 results.numPosPeaks = 0
641 results.numNeg = 0
642 results.numNegPeaks = 0
643 positive = ""
644 negative = ""
645
646 if results.positive is not None:
647 results.numPos = len(results.positive.getFootprints())
648 results.numPosPeaks = sum(len(fp.getPeaks()) for fp in results.positive.getFootprints())
649 positive = " %d positive peaks in %d footprints" % (results.numPosPeaks, results.numPos)
650 if results.negative is not None:
651 results.numNeg = len(results.negative.getFootprints())
652 results.numNegPeaks = sum(len(fp.getPeaks()) for fp in results.negative.getFootprints())
653 negative = " %d negative peaks in %d footprints" % (results.numNegPeaks, results.numNeg)
654
655 self.log.info("Detected%s%s%s to %g +ve and %g -ve %s",
656 positive, " and" if positive and negative else "", negative,
657 self.config.thresholdValue*self.config.includeThresholdMultiplier*factor,
658 self.config.thresholdValue*self.config.includeThresholdMultiplier*factorNeg,
659 "DN" if self.config.thresholdType == "value" else "sigma")
660
661 def reEstimateBackground(self, maskedImage, backgrounds):
662 """Estimate the background after detection
663
664 Parameters
665 ----------
666 maskedImage : `lsst.afw.image.MaskedImage`
667 Image on which to estimate the background.
668 backgrounds : `lsst.afw.math.BackgroundList`
669 List of backgrounds; modified.
670
671 Returns
672 -------
673 bg : `lsst.afw.math.backgroundMI`
674 Empirical background model.
675 """
676 bg = self.background.fitBackground(maskedImage)
677 if self.config.adjustBackground:
678 self.log.warning("Fiddling the background by %g", self.config.adjustBackground)
679 bg += self.config.adjustBackground
680 self.log.info("Resubtracting the background after object detection")
681 maskedImage -= bg.getImageF(self.background.config.algorithm,
682 self.background.config.undersampleStyle)
683
684 actrl = bg.getBackgroundControl().getApproximateControl()
685 backgrounds.append((bg, getattr(afwMath.Interpolate, self.background.config.algorithm),
686 bg.getAsUsedUndersampleStyle(), actrl.getStyle(), actrl.getOrderX(),
687 actrl.getOrderY(), actrl.getWeighting()))
688 return bg
689
690 def clearUnwantedResults(self, mask, results):
691 """Clear unwanted results from the Struct of results
692
693 If we specifically want only positive or only negative detections,
694 drop the ones we don't want, and its associated mask plane.
695
696 Parameters
697 ----------
698 mask : `lsst.afw.image.Mask`
699 Mask image.
700 results : `lsst.pipe.base.Struct`
701 Detection results, with ``positive`` and ``negative`` elements;
702 modified.
703 """
704 if self.config.thresholdPolarity == "positive":
705 if self.config.reEstimateBackground:
706 mask &= ~mask.getPlaneBitMask("DETECTED_NEGATIVE")
707 results.negative = None
708 elif self.config.thresholdPolarity == "negative":
709 if self.config.reEstimateBackground:
710 mask &= ~mask.getPlaneBitMask("DETECTED")
711 results.positive = None
712
713 @timeMethod
714 def detectFootprints(self, exposure, doSmooth=True, sigma=None, clearMask=True, expId=None,
715 background=None):
716 """Detect footprints on an exposure.
717
718 Parameters
719 ----------
720 exposure : `lsst.afw.image.Exposure`
721 Exposure to process; DETECTED{,_NEGATIVE} mask plane will be
722 set in-place.
723 doSmooth : `bool`, optional
724 If True, smooth the image before detection using a Gaussian
725 of width ``sigma``, or the measured PSF width of ``exposure``.
726 Set to False when running on e.g. a pre-convolved image, or a mask
727 plane.
728 sigma : `float`, optional
729 Gaussian Sigma of PSF (pixels); used for smoothing and to grow
730 detections; if `None` then measure the sigma of the PSF of the
731 ``exposure``.
732 clearMask : `bool`, optional
733 Clear both DETECTED and DETECTED_NEGATIVE planes before running
734 detection.
735 expId : `dict`, optional
736 Exposure identifier; unused by this implementation, but used for
737 RNG seed by subclasses.
738 background : `lsst.afw.math.BackgroundList`, optional
739 Background that was already subtracted from the exposure; will be
740 modified in-place if ``reEstimateBackground=True``.
741
742 Returns
743 -------
744 results : `lsst.pipe.base.Struct`
745 A `~lsst.pipe.base.Struct` containing:
746
747 ``positive``
748 Positive polarity footprints.
749 (`lsst.afw.detection.FootprintSet` or `None`)
750 ``negative``
751 Negative polarity footprints.
752 (`lsst.afw.detection.FootprintSet` or `None`)
753 ``numPos``
754 Number of footprints in positive or 0 if detection polarity was
755 negative. (`int`)
756 ``numNeg``
757 Number of footprints in negative or 0 if detection polarity was
758 positive. (`int`)
759 ``background``
760 Re-estimated background. `None` or the input ``background``
761 if ``reEstimateBackground==False``.
762 (`lsst.afw.math.BackgroundList`)
763 ``factor``
764 Multiplication factor applied to the configured detection
765 threshold. (`float`)
766 """
767 maskedImage = exposure.maskedImage
768
769 if clearMask:
770 self.clearMask(maskedImage.getMask())
771
772 psf = self.getPsf(exposure, sigma=sigma)
773 with self.tempWideBackgroundContext(exposure):
774 convolveResults = self.convolveImage(maskedImage, psf, doSmooth=doSmooth)
775 middle = convolveResults.middle
776 sigma = convolveResults.sigma
777 self.removeBadPixels(middle)
778
779 results = self.applyThreshold(middle, maskedImage.getBBox())
780 results.background = background if background is not None else afwMath.BackgroundList()
781
782 if self.config.doTempLocalBackground:
783 self.applyTempLocalBackground(exposure, middle, results)
784 self.finalizeFootprints(maskedImage.mask, results, sigma)
785
786 # Compute the significance of peaks after the peaks have been
787 # finalized and after local background correction/updatePeaks, so
788 # that the significance represents the "final" detection S/N.
789 results.positive = self.setPeakSignificance(middle, results.positive, results.positiveThreshold)
790 results.negative = self.setPeakSignificance(middle, results.negative, results.negativeThreshold,
791 negative=True)
792
793 if self.config.reEstimateBackground:
794 self.reEstimateBackground(maskedImage, results.background)
795
796 self.clearUnwantedResults(maskedImage.getMask(), results)
797
798 self.display(exposure, results, middle)
799
800 return results
801
802 def removeBadPixels(self, middle):
803 """Set the significance of flagged pixels to zero.
804
805 Parameters
806 ----------
807 middle : `lsst.afw.image.ExposureF`
808 Score or maximum likelihood difference image.
809 The image plane will be modified in place.
810 """
811 badPixelMask = lsst.afw.image.Mask.getPlaneBitMask(self.config.excludeMaskPlanes)
812 badPixels = middle.mask.array & badPixelMask > 0
813 middle.image.array[badPixels] = 0
814
815 def setPeakSignificance(self, exposure, footprints, threshold, negative=False):
816 """Set the significance of each detected peak to the pixel value divided
817 by the appropriate standard-deviation for ``config.thresholdType``.
818
819 Only sets significance for "stdev" and "pixel_stdev" thresholdTypes;
820 we leave it undefined for "value" and "variance" as it does not have a
821 well-defined meaning in those cases.
822
823 Parameters
824 ----------
825 exposure : `lsst.afw.image.Exposure`
826 Exposure that footprints were detected on, likely the convolved,
827 local background-subtracted image.
828 footprints : `lsst.afw.detection.FootprintSet`
829 Footprints detected on the image.
830 threshold : `lsst.afw.detection.Threshold`
831 Threshold used to find footprints.
832 negative : `bool`, optional
833 Are we calculating for negative sources?
834 """
835 if footprints is None or footprints.getFootprints() == []:
836 return footprints
837 polarity = -1 if negative else 1
838
839 # All incoming footprints have the same schema.
840 mapper = afwTable.SchemaMapper(footprints.getFootprints()[0].peaks.schema)
841 mapper.addMinimalSchema(footprints.getFootprints()[0].peaks.schema)
842 mapper.addOutputField("significance", type=float,
843 doc="Ratio of peak value to configured standard deviation.")
844
845 # Copy the old peaks to the new ones with a significance field.
846 # Do this independent of the threshold type, so we always have a
847 # significance field.
848 newFootprints = afwDet.FootprintSet(footprints)
849 for old, new in zip(footprints.getFootprints(), newFootprints.getFootprints()):
850 newPeaks = afwDet.PeakCatalog(mapper.getOutputSchema())
851 newPeaks.extend(old.peaks, mapper=mapper)
852 new.getPeaks().clear()
853 new.setPeakCatalog(newPeaks)
854
855 # Compute the significance values.
856 if self.config.thresholdType == "pixel_stdev":
857 for footprint in newFootprints.getFootprints():
858 footprint.updatePeakSignificance(exposure.variance, polarity)
859 elif self.config.thresholdType == "stdev":
860 sigma = threshold.getValue() / self.config.thresholdValue
861 for footprint in newFootprints.getFootprints():
862 footprint.updatePeakSignificance(polarity*sigma)
863 else:
864 for footprint in newFootprints.getFootprints():
865 for peak in footprint.peaks:
866 peak["significance"] = 0
867
868 return newFootprints
869
870 def makeThreshold(self, image, thresholdParity, factor=1.0):
871 """Make an afw.detection.Threshold object corresponding to the task's
872 configuration and the statistics of the given image.
873
874 Parameters
875 ----------
876 image : `afw.image.MaskedImage`
877 Image to measure noise statistics from if needed.
878 thresholdParity: `str`
879 One of "positive" or "negative", to set the kind of fluctuations
880 the Threshold will detect.
881 factor : `float`
882 Factor by which to multiply the configured detection threshold.
883 This is useful for tweaking the detection threshold slightly.
884
885 Returns
886 -------
887 threshold : `lsst.afw.detection.Threshold`
888 Detection threshold.
889 """
890 parity = False if thresholdParity == "negative" else True
891 thresholdValue = self.config.thresholdValue
892 thresholdType = self.config.thresholdType
893 if self.config.thresholdType == 'stdev':
894 bad = image.getMask().getPlaneBitMask(self.config.statsMask)
896 sctrl.setAndMask(bad)
897 stats = afwMath.makeStatistics(image, afwMath.STDEVCLIP, sctrl)
898 thresholdValue *= stats.getValue(afwMath.STDEVCLIP)
899 thresholdType = 'value'
900
901 threshold = afwDet.createThreshold(thresholdValue*factor, thresholdType, parity)
902 threshold.setIncludeMultiplier(self.config.includeThresholdMultiplier)
903 self.log.debug("Detection threshold: %s", threshold)
904 return threshold
905
906 def updatePeaks(self, fpSet, image, threshold):
907 """Update the Peaks in a FootprintSet by detecting new Footprints and
908 Peaks in an image and using the new Peaks instead of the old ones.
909
910 Parameters
911 ----------
912 fpSet : `afw.detection.FootprintSet`
913 Set of Footprints whose Peaks should be updated.
914 image : `afw.image.MaskedImage`
915 Image to detect new Footprints and Peak in.
916 threshold : `afw.detection.Threshold`
917 Threshold object for detection.
918
919 Input Footprints with fewer Peaks than self.config.nPeaksMaxSimple
920 are not modified, and if no new Peaks are detected in an input
921 Footprint, the brightest original Peak in that Footprint is kept.
922 """
923 for footprint in fpSet.getFootprints():
924 oldPeaks = footprint.getPeaks()
925 if len(oldPeaks) <= self.config.nPeaksMaxSimple:
926 continue
927 # We detect a new FootprintSet within each non-simple Footprint's
928 # bbox to avoid a big O(N^2) comparison between the two sets of
929 # Footprints.
930 sub = image.Factory(image, footprint.getBBox())
931 fpSetForPeaks = afwDet.FootprintSet(
932 sub,
933 threshold,
934 "", # don't set a mask plane
935 self.config.minPixels
936 )
937 newPeaks = afwDet.PeakCatalog(oldPeaks.getTable())
938 for fpForPeaks in fpSetForPeaks.getFootprints():
939 for peak in fpForPeaks.getPeaks():
940 if footprint.contains(peak.getI()):
941 newPeaks.append(peak)
942 if len(newPeaks) > 0:
943 del oldPeaks[:]
944 oldPeaks.extend(newPeaks)
945 else:
946 del oldPeaks[1:]
947
948 @staticmethod
949 def setEdgeBits(maskedImage, goodBBox, edgeBitmask):
950 """Set the edgeBitmask bits for all of maskedImage outside goodBBox
951
952 Parameters
953 ----------
954 maskedImage : `lsst.afw.image.MaskedImage`
955 Image on which to set edge bits in the mask.
956 goodBBox : `lsst.geom.Box2I`
957 Bounding box of good pixels, in ``LOCAL`` coordinates.
958 edgeBitmask : `lsst.afw.image.MaskPixel`
959 Bit mask to OR with the existing mask bits in the region
960 outside ``goodBBox``.
961 """
962 msk = maskedImage.getMask()
963
964 mx0, my0 = maskedImage.getXY0()
965 for x0, y0, w, h in ([0, 0,
966 msk.getWidth(), goodBBox.getBeginY() - my0],
967 [0, goodBBox.getEndY() - my0, msk.getWidth(),
968 maskedImage.getHeight() - (goodBBox.getEndY() - my0)],
969 [0, 0,
970 goodBBox.getBeginX() - mx0, msk.getHeight()],
971 [goodBBox.getEndX() - mx0, 0,
972 maskedImage.getWidth() - (goodBBox.getEndX() - mx0), msk.getHeight()],
973 ):
974 edgeMask = msk.Factory(msk, lsst.geom.BoxI(lsst.geom.PointI(x0, y0),
975 lsst.geom.ExtentI(w, h)), afwImage.LOCAL)
976 edgeMask |= edgeBitmask
977
978 @contextmanager
979 def tempWideBackgroundContext(self, exposure):
980 """Context manager for removing wide (large-scale) background
981
982 Removing a wide (large-scale) background helps to suppress the
983 detection of large footprints that may overwhelm the deblender.
984 It does, however, set a limit on the maximum scale of objects.
985
986 The background that we remove will be restored upon exit from
987 the context manager.
988
989 Parameters
990 ----------
991 exposure : `lsst.afw.image.Exposure`
992 Exposure on which to remove large-scale background.
993
994 Returns
995 -------
996 context : context manager
997 Context manager that will ensure the temporary wide background
998 is restored.
999 """
1000 doTempWideBackground = self.config.doTempWideBackground
1001 if doTempWideBackground:
1002 self.log.info("Applying temporary wide background subtraction")
1003 original = exposure.maskedImage.image.array[:].copy()
1004 self.tempWideBackground.run(exposure).background
1005 # Remove NO_DATA regions (e.g., edge of the field-of-view); these can cause detections after
1006 # subtraction because of extrapolation of the background model into areas with no constraints.
1007 image = exposure.maskedImage.image
1008 mask = exposure.maskedImage.mask
1009 noData = mask.array & mask.getPlaneBitMask("NO_DATA") > 0
1010 isGood = mask.array & mask.getPlaneBitMask(self.config.statsMask) == 0
1011 image.array[noData] = np.median(image.array[~noData & isGood])
1012 try:
1013 yield
1014 finally:
1015 if doTempWideBackground:
1016 exposure.maskedImage.image.array[:] = original
1017
1018
1019def addExposures(exposureList):
1020 """Add a set of exposures together.
1021
1022 Parameters
1023 ----------
1024 exposureList : `list` of `lsst.afw.image.Exposure`
1025 Sequence of exposures to add.
1026
1027 Returns
1028 -------
1029 addedExposure : `lsst.afw.image.Exposure`
1030 An exposure of the same size as each exposure in ``exposureList``,
1031 with the metadata from ``exposureList[0]`` and a masked image equal
1032 to the sum of all the exposure's masked images.
1033 """
1034 exposure0 = exposureList[0]
1035 image0 = exposure0.getMaskedImage()
1036
1037 addedImage = image0.Factory(image0, True)
1038 addedImage.setXY0(image0.getXY0())
1039
1040 for exposure in exposureList[1:]:
1041 image = exposure.getMaskedImage()
1042 addedImage += image
1043
1044 addedExposure = exposure0.Factory(addedImage, exposure0.getWcs())
1045 return addedExposure
A set of Footprints, associated with a MaskedImage.
A circularly symmetric Gaussian Psf class with no spatial variation, intended mostly for testing purp...
Definition GaussianPsf.h:42
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
Parameters to control convolution.
A kernel described by a pair of functions: func(x, y) = colFunc(x) * rowFunc(y)
Definition Kernel.h:860
Pass parameters to a Statistics object.
Definition Statistics.h:83
A mapping between the keys of two Schemas, used to copy data between them.
An integer coordinate rectangle.
Definition Box.h:55
makeThreshold(self, image, thresholdParity, factor=1.0)
Definition detection.py:870
applyTempLocalBackground(self, exposure, middle, results)
Definition detection.py:360
run(self, table, exposure, doSmooth=True, sigma=None, clearMask=True, expId=None, background=None)
Definition detection.py:221
setEdgeBits(maskedImage, goodBBox, edgeBitmask)
Definition detection.py:949
setPeakSignificance(self, exposure, footprints, threshold, negative=False)
Definition detection.py:815
detectFootprints(self, exposure, doSmooth=True, sigma=None, clearMask=True, expId=None, background=None)
Definition detection.py:715
applyThreshold(self, middle, bbox, factor=1.0, factorNeg=None)
Definition detection.py:525
convolveImage(self, maskedImage, psf, doSmooth=True)
Definition detection.py:463
display(self, exposure, results, convolvedImage=None)
Definition detection.py:303
finalizeFootprints(self, mask, results, sigma, factor=1.0, factorNeg=None)
Definition detection.py:591
reEstimateBackground(self, maskedImage, backgrounds)
Definition detection.py:661
Threshold createThreshold(const double value, const std::string type="value", const bool polarity=true)
Factory method for creating Threshold objects.
Definition Threshold.cc:109
Statistics makeStatistics(lsst::afw::image::Image< Pixel > const &img, lsst::afw::image::Mask< image::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Handle a watered-down front-end to the constructor (no variance)
Definition Statistics.h:361
void convolve(OutImageT &convolvedImage, InImageT const &inImage, KernelT const &kernel, ConvolutionControl const &convolutionControl=ConvolutionControl())
Convolve an Image or MaskedImage with a Kernel, setting pixels of an existing output image.