Coverage for python / lsst / meas / deblender / baseline.py: 15%
260 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-28 08:57 +0000
1# This file is part of meas_deblender.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
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 GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
22__all__ = ["DEFAULT_PLUGINS", "DeblenderResult", "DeblendedParent", "MultiColorPeak",
23 "DeblendedPeak", "deblend", "newDeblend", "CachingPsf"]
25from collections import OrderedDict
26import numpy as np
28import lsst.pex.exceptions
29import lsst.afw.image as afwImage
30import lsst.afw.detection as afwDet
31import lsst.geom as geom
32import lsst.afw.math as afwMath
33import lsst.utils.logging
35from . import plugins
37DEFAULT_PLUGINS = [
38 plugins.DeblenderPlugin(plugins.fitPsfs),
39 plugins.DeblenderPlugin(plugins.buildSymmetricTemplates),
40 plugins.DeblenderPlugin(plugins.medianSmoothTemplates),
41 plugins.DeblenderPlugin(plugins.makeTemplatesMonotonic),
42 plugins.DeblenderPlugin(plugins.clipFootprintsToNonzero),
43 plugins.DeblenderPlugin(plugins.reconstructTemplates, onReset=5),
44 plugins.DeblenderPlugin(plugins.apportionFlux),
45]
48class DeblenderResult:
49 r"""Collection of objects in multiple bands for a single parent footprint.
51 Parameters
52 ----------
53 footprint: `afw.detection.Footprint`
54 Parent footprint to deblend. While `maskedImages`, `psfs`,
55 and `psffwhms` are lists of objects, one for each band,
56 `footprint` is a single parent footprint (from a `mergedDet`)
57 this is used for all bands.
58 mMaskedImage: `MaskedImage`\ s or `MultibandMaskedImage`
59 Masked image containing the ``footprint`` in each band.
60 psfs: list of `afw.detection.Psf`s
61 Psf of the ``maskedImage`` for each band.
62 psffwhm: list of `float`s
63 FWHM of the ``maskedImage``'s ``psf`` in each band.
64 maxNumberOfPeaks: `int`, optional
65 If positive, the maximum number of peaks to deblend.
66 If the total number of peaks is greater than ``maxNumberOfPeaks``,
67 then only the first ``maxNumberOfPeaks`` sources are deblended.
68 The default is 0, which deblends all of the peaks.
69 avgNoise: `float`or list of `float`s, optional
70 Average noise level in each ``maskedImage``.
71 The default is ``None``, which estimates the noise from the median value of the
72 variance plane of ``maskedImage`` for each filter.
73 """
75 def __init__(self, footprint, mMaskedImage, psfs, psffwhms, log,
76 maxNumberOfPeaks=0, avgNoise=None):
77 # Check if this is collection of footprints in multiple bands or a single footprint
78 if not isinstance(mMaskedImage, afwImage.MultibandMaskedImage):
79 mMaskedImage = [mMaskedImage]
80 self.filters = [0]
81 else:
82 self.filters = mMaskedImage.filters
83 try:
84 len(psfs)
85 except TypeError:
86 psfs = [psfs]
87 try:
88 len(psffwhms)
89 except TypeError:
90 psffwhms = [psffwhms]
91 try:
92 len(avgNoise)
93 except TypeError:
94 if avgNoise is None:
95 avgNoise = [None]*len(psfs)
96 else:
97 avgNoise = [avgNoise]
98 # Now check that all of the parameters have the same number of entries
99 if any([len(self.filters) != len(p) for p in [psfs, psffwhms, avgNoise]]):
100 raise ValueError("To use the multi-color deblender, "
101 "'maskedImage', 'psf', 'psffwhm', 'avgNoise'"
102 "must have the same length, but instead have lengths: "
103 "{0}".format([len(p) for p in [mMaskedImage,
104 psfs,
105 psffwhms,
106 avgNoise]]))
108 self.log = log
109 self.mMaskedImage = mMaskedImage
110 self.footprint = footprint
111 self.psfs = psfs
113 self.peakCount = len(footprint.getPeaks())
114 if maxNumberOfPeaks > 0 and maxNumberOfPeaks < self.peakCount:
115 self.peakCount = maxNumberOfPeaks
117 # Create a DeblendedParent for the Footprint in every filter
118 self.deblendedParents = OrderedDict([])
119 for n, f in enumerate(self.filters):
120 f = self.filters[n]
121 dp = DeblendedParent(f, footprint, mMaskedImage[f], psfs[n],
122 psffwhms[n], avgNoise[n], maxNumberOfPeaks, self)
123 self.deblendedParents[self.filters[n]] = dp
125 # Group the peaks in each color
126 self.peaks = []
127 for idx in range(self.peakCount):
128 peakDict = OrderedDict([(f, dp.peaks[idx]) for f, dp in self.deblendedParents.items()])
129 multiPeak = MultiColorPeak(peakDict, idx, self)
130 self.peaks.append(multiPeak)
132 # Result from multiband debender (if used)
133 self.blend = None
134 self.failed = False
136 def getParentProperty(self, propertyName):
137 """Get the footprint in each filter
138 """
139 return [getattr(dp, propertyName) for dp in self.deblendedParents]
141 def setTemplateSums(self, templateSums, fidx=None):
142 if fidx is not None:
143 self.templateSums[fidx] = templateSums
144 else:
145 for f, templateSum in templateSums.items():
146 self.deblendedParents[f].templateSum = templateSum
149class DeblendedParent:
150 """Deblender result of a single parent footprint, in a single band
152 Convenience class to store useful objects used by the deblender for a single band,
153 such as the maskedImage, psf, etc as well as the results from the deblender.
155 Parameters
156 ----------
157 filterName: `str`
158 Name of the filter used for `maskedImage`
159 footprint: list of `afw.detection.Footprint`s
160 Parent footprint to deblend in each band.
161 maskedImages: list of `afw.image.MaskedImageF`s
162 Masked image containing the ``footprint`` in each band.
163 psf: list of `afw.detection.Psf`s
164 Psf of the ``maskedImage`` for each band.
165 psffwhm: list of `float`s
166 FWHM of the ``maskedImage``'s ``psf`` in each band.
167 avgNoise: `float`or list of `float`s, optional
168 Average noise level in each ``maskedImage``.
169 The default is ``None``, which estimates the noise from the median value of the
170 variance plane of ``maskedImage`` for each filter.
171 maxNumberOfPeaks: `int`, optional
172 If positive, the maximum number of peaks to deblend.
173 If the total number of peaks is greater than ``maxNumberOfPeaks``,
174 then only the first ``maxNumberOfPeaks`` sources are deblended.
175 The default is 0, which deblends all of the peaks.
176 debResult: `DeblenderResult`
177 The ``DeblenderResult`` that contains this ``DeblendedParent``.
178 """
179 def __init__(self, filterName, footprint, maskedImage, psf, psffwhm, avgNoise,
180 maxNumberOfPeaks, debResult):
181 self.filter = filterName
182 self.fp = footprint
183 self.maskedImage = maskedImage
184 self.psf = psf
185 self.psffwhm = psffwhm
186 self.img = maskedImage.getImage()
187 self.imbb = self.img.getBBox()
188 self.varimg = maskedImage.getVariance()
189 self.mask = maskedImage.getMask()
190 self.avgNoise = avgNoise
191 self.updateFootprintBbox()
192 self.debResult = debResult
193 self.peakCount = debResult.peakCount
194 self.templateSum = None
196 # avgNoise is an estiamte of the average noise level for the image in this filter
197 if avgNoise is None:
198 stats = afwMath.makeStatistics(self.varimg, self.mask, afwMath.MEDIAN)
199 avgNoise = np.sqrt(stats.getValue(afwMath.MEDIAN))
200 debResult.log.trace('Estimated avgNoise for filter %s = %f', self.filter, avgNoise)
201 self.avgNoise = avgNoise
203 # Store all of the peak information in a single object
204 self.peaks = []
205 peaks = self.fp.getPeaks()
206 for idx in range(self.peakCount):
207 deblendedPeak = DeblendedPeak(peaks[idx], idx, self)
208 self.peaks.append(deblendedPeak)
210 def updateFootprintBbox(self):
211 """Update the bounding box of the parent footprint
213 If the parent Footprint is resized it will be necessary to update the bounding box of the footprint.
214 """
215 # Pull out the image bounds of the parent Footprint
216 self.bb = self.fp.getBBox()
217 if not self.imbb.contains(self.bb):
218 raise ValueError(('Footprint bounding-box %s extends outside image bounding-box %s') %
219 (str(self.bb), str(self.imbb)))
220 self.W, self.H = self.bb.getWidth(), self.bb.getHeight()
221 self.x0, self.y0 = self.bb.getMinX(), self.bb.getMinY()
222 self.x1, self.y1 = self.bb.getMaxX(), self.bb.getMaxY()
225class MultiColorPeak:
226 """Collection of single peak deblender results in multiple bands.
228 There is one of these objects for each Peak in the footprint.
230 Parameters
231 ----------
232 peaks: `dict` of `afw.detection.PeakRecord`
233 Each entry in ``peaks`` is a peak record for the same peak in a different band
234 pki: int
235 Index of the peak in `parent.peaks`
236 parent: `DeblendedParent`
237 DeblendedParent object that contains the peak as a child.
238 """
239 def __init__(self, peaks, pki, parent):
240 self.filters = list(peaks.keys())
241 self.deblendedPeaks = peaks
242 self.parent = parent
243 for filter, peak in self.deblendedPeaks.items():
244 peak.multiColorPeak = self
246 # Fields common to the peak in all bands that will be set by the deblender
247 # In the future this is likely to contain information about the probability of the peak
248 # being a point source, color-information about templateFootprints, etc.
249 self.pki = pki
250 self.skip = False
251 self.deblendedAsPsf = False
252 self.x = self.deblendedPeaks[self.filters[0]].peak.getFx()
253 self.y = self.deblendedPeaks[self.filters[0]].peak.getFy()
256class DeblendedPeak:
257 """Result of deblending a single Peak within a parent Footprint.
259 There is one of these objects for each Peak in the Footprint.
261 Parameters
262 ----------
263 peak: `afw.detection.PeakRecord`
264 Peak object in a single band from a peak record
265 pki: `int`
266 Index of the peak in `multiColorPeak.parent.peaks`
267 parent: `DeblendedParent`
268 Parent in the same filter that contains the peak
269 multiColorPeak: `MultiColorPeak`
270 Object containing the same peak in multiple bands
271 """
272 def __init__(self, peak, pki, parent, multiColorPeak=None):
273 # Peak object
274 self.peak = peak
275 # int, peak index number
276 self.pki = pki
277 self.parent = parent
278 self.multiColorPeak = multiColorPeak
279 # union of all the ways of failing...
280 self.skip = False
282 self.outOfBounds = False
283 self.tinyFootprint = False
284 self.noValidPixels = False
285 self.deblendedAsPsf = False
286 self.degenerate = False
288 # Field set during _fitPsf:
289 self.psfFitFailed = False
290 self.psfFitBadDof = False
291 # (chisq, dof) for PSF fit without decenter
292 self.psfFit1 = None
293 # (chisq, dof) for PSF fit with decenter
294 self.psfFit2 = None
295 # (chisq, dof) for PSF fit after applying decenter
296 self.psfFit3 = None
297 # decentered PSF fit wanted to move the center too much
298 self.psfFitBigDecenter = False
299 # was the fit with decenter better?
300 self.psfFitWithDecenter = False
301 #
302 self.psfFitR0 = None
303 self.psfFitR1 = None
304 self.psfFitStampExtent = None
305 self.psfFitCenter = None
306 self.psfFitBest = None
307 self.psfFitParams = None
308 self.psfFitFlux = None
309 self.psfFitNOthers = None
311 # Things only set in _fitPsf when debugging is turned on:
312 self.psfFitDebugPsf0Img = None
313 self.psfFitDebugPsfImg = None
314 self.psfFitDebugPsfDerivImg = None
315 self.psfFitDebugPsfModel = None
317 self.failedSymmetricTemplate = False
319 # The actual template Image and Footprint
320 self.templateImage = None
321 self.templateFootprint = None
323 # The flux assigned to this template -- a MaskedImage
324 self.fluxPortion = None
326 # The stray flux assigned to this template (may be None), a HeavyFootprint
327 self.strayFlux = None
329 self.hasRampedTemplate = False
331 self.patched = False
333 # debug -- a copy of the original symmetric template
334 self.origTemplate = None
335 self.origFootprint = None
336 # MaskedImage
337 self.rampedTemplate = None
338 # MaskedImage
339 self.medianFilteredTemplate = None
341 # when least-squares fitting templates, the template weight.
342 self.templateWeight = 1.0
344 def __str__(self):
345 return (('deblend result: outOfBounds: %s, deblendedAsPsf: %s') %
346 (self.outOfBounds, self.deblendedAsPsf))
348 @property
349 def psfFitChisq(self):
350 chisq, dof = self.psfFitBest
351 return chisq
353 @property
354 def psfFitDof(self):
355 chisq, dof = self.psfFitBest
356 return dof
358 def getFluxPortion(self, strayFlux=True):
359 """
360 Return a HeavyFootprint containing the flux apportioned to this peak.
362 @param[in] strayFlux include stray flux also?
363 """
364 if self.templateFootprint is None or self.fluxPortion is None:
365 return None
366 heavy = afwDet.makeHeavyFootprint(self.templateFootprint, self.fluxPortion)
367 if strayFlux:
368 if self.strayFlux is not None:
369 heavy = afwDet.mergeHeavyFootprints(heavy, self.strayFlux)
371 return heavy
373 def setStrayFlux(self, stray):
374 self.strayFlux = stray
376 def setFluxPortion(self, mimg):
377 self.fluxPortion = mimg
379 def setTemplateWeight(self, w):
380 self.templateWeight = w
382 def setPatched(self):
383 self.patched = True
385 # DEBUG
386 def setOrigTemplate(self, t, tfoot):
387 self.origTemplate = t.Factory(t, True)
388 self.origFootprint = tfoot
390 def setRampedTemplate(self, t, tfoot):
391 self.hasRampedTemplate = True
392 self.rampedTemplate = t.Factory(t, True)
394 def setMedianFilteredTemplate(self, t, tfoot):
395 self.medianFilteredTemplate = t.Factory(t, True)
397 def setPsfTemplate(self, tim, tfoot):
398 self.psfFootprint = afwDet.Footprint(tfoot)
399 self.psfTemplate = tim.Factory(tim, True)
401 def setOutOfBounds(self):
402 self.outOfBounds = True
403 self.skip = True
405 def setTinyFootprint(self):
406 self.tinyFootprint = True
407 self.skip = True
409 def setNoValidPixels(self):
410 self.noValidPixels = True
411 self.skip = True
413 def setPsfFitFailed(self):
414 self.psfFitFailed = True
416 def setBadPsfDof(self):
417 self.psfFitBadDof = True
419 def setDeblendedAsPsf(self):
420 self.deblendedAsPsf = True
422 def setFailedSymmetricTemplate(self):
423 self.failedSymmetricTemplate = True
424 self.skip = True
426 def setTemplate(self, image, footprint):
427 self.templateImage = image
428 self.templateFootprint = footprint
431def deblend(footprint, maskedImage, psf, psffwhm,
432 psfChisqCut1=1.5, psfChisqCut2=1.5, psfChisqCut2b=1.5, fitPsfs=True,
433 medianSmoothTemplate=True, medianFilterHalfsize=2,
434 monotonicTemplate=True, weightTemplates=False,
435 log=None, verbose=False, sigma1=None, maxNumberOfPeaks=0,
436 assignStrayFlux=True, strayFluxToPointSources='necessary', strayFluxAssignment='r-to-peak',
437 rampFluxAtEdge=False, patchEdges=False, tinyFootprintSize=2,
438 getTemplateSum=False, clipStrayFluxFraction=0.001, clipFootprintToNonzero=True,
439 removeDegenerateTemplates=False, maxTempDotProd=0.5
440 ):
441 r"""Deblend a parent ``Footprint`` in a ``MaskedImageF``.
443 Deblending assumes that ``footprint`` has multiple peaks, as it will still create a
444 `PerFootprint` object with a list of peaks even if there is only one peak in the list.
445 It is recommended to first check that ``footprint`` has more than one peak, similar to the
446 execution of `lsst.meas.deblender.deblend.SourceDeblendTask`.
448 .. note::
449 This is the API for the old deblender, however the function builds the plugins necessary
450 to use the new deblender to perform identically to the old deblender.
451 To test out newer functionality use ``newDeblend`` instead.
453 Deblending involves several mandatory and optional steps:
455 # Optional: If ``fitPsfs`` is True, find all peaks that are well-fit by a PSF + background model
457 * Peaks that pass the cuts have their footprints modified to the PSF + background model
458 and their ``deblendedAsPsf`` property set to ``True``.
459 * Relevant parameters: ``psfChisqCut1``, ``psfChisqCut2``, ``psfChisqCut2b``,
460 ``tinyFootprintSize``.
461 * See the parameter descriptions for more.
463 # Build a symmetric template for each peak not well-fit by the PSF model
465 * Given ``maskedImageF``, ``footprint``, and a ``DeblendedPeak``, creates a symmetric
466 template (``templateImage`` and ``templateFootprint``) around the peak
467 for all peaks not flagged as ``skip`` or ``deblendedAsPsf``.
468 * If ``patchEdges=True`` and if ``footprint`` touches pixels with the
469 ``EDGE`` bit set, then ``footprint`` is grown to include spans whose
470 symmetric mirror is outside of the image.
471 * Relevant parameters: ``sigma1`` and ``patchEdges``.
473 # Optional: If ``rampFluxAtEdge`` is True, adjust flux on the edges of the template footprints
475 * Using the PSF, a peak ``Footprint`` with pixels on the edge of of ``footprint``
476 is grown by the psffwhm*1.5 and filled in with zeros.
477 * The result is a new symmetric footprint template for the peaks near the edge.
478 * Relevant parameter: ``patchEdges``.
480 # Optionally (``medianSmoothTemplate=True``) filter the template images
482 * Apply a median smoothing filter to all of the template images.
483 * Relevant parameters: ``medianFilterHalfSize``
485 # Optional: If ``monotonicTemplate`` is True, make the templates monotonic.
487 * The pixels in the templates are modified such that pixels
488 further from the peak will have values smaller than those closer to the peak.
490 # Optional: If ``clipFootprintToNonzero`` is True, clip non-zero spans in the template footprints
492 * Peak ``Footprint``\ s are clipped to the region in the image containing non-zero values
493 by dropping spans that are completely zero and moving endpoints to non-zero pixels
494 (but does not split spans that have internal zeros).
496 # Optional: If ``weightTemplates`` is True, weight the templates to best fit the observed image
498 * Re-weight the templates so that their linear combination
499 best represents the observed ``maskedImage``
501 # Optional: If ``removeDegenerateTempaltes`` is True, reconstruct shredded galaxies
503 * If galaxies have substructure, such as face-on spirals, the process of identifying peaks can
504 "shred" the galaxy into many pieces. The templates of shredded galaxies are typically quite
505 similar because they represent the same galaxy, so we try to identify these "degenerate" peaks
506 by looking at the inner product (in pixel space) of pairs of templates.
507 * If they are nearly parallel, we only keep one of the peaks and reject the other.
508 * If only one of the peaks is a PSF template, the other template is used,
509 otherwise the one with the maximum template value is kept.
510 * Relevant parameters: ``maxTempDotProduct``
512 # Apportion flux to all of the peak templates
514 * Divide the ``maskedImage`` flux amongst all of the templates based on the fraction of
515 flux assigned to each ``tempalteFootprint``.
516 * Leftover "stray flux" is assigned to peaks based on the other parameters.
517 * Relevant parameters: ``clipStrayFluxFraction``, ``strayFluxAssignment``,
518 ``strayFluxToPointSources``, ``assignStrayFlux``
520 Parameters
521 ----------
522 footprint: `afw.detection.Footprint`
523 Parent footprint to deblend
524 maskedImage: `afw.image.MaskedImageF`
525 Masked image containing the ``footprint``
526 psf: `afw.detection.Psf`
527 Psf of the ``maskedImage``
528 psffwhm: `float`
529 FWHM of the ``maskedImage``\'s ``psf``
530 psfChisqCut*: `float`, optional
531 If ``fitPsfs==True``, all of the peaks are fit to the image PSF.
532 ``psfChisqCut1`` is the maximum chi-squared-per-degree-of-freedom allowed for a peak to
533 be considered a PSF match without recentering.
534 A fit is also made that includes terms to recenter the PSF.
535 ``psfChisqCut2`` is the same as ``psfChisqCut1`` except it determines the restriction on the
536 fit that includes recentering terms.
537 If the peak is a match for a re-centered PSF, the PSF is repositioned at the new center and
538 the peak footprint is fit again, this time to the new PSF.
539 If the resulting chi-squared-per-degree-of-freedom is less than ``psfChisqCut2b`` then it
540 passes the re-centering algorithm.
541 If the peak passes both the re-centered and fixed position cuts, the better of the two is accepted,
542 but parameters for all three psf fits are stored in the ``DeblendedPeak``.
543 The default for ``psfChisqCut1``, ``psfChisqCut2``, and ``psfChisqCut2b`` is ``1.5``.
544 fitPsfs: `bool`, optional
545 If True then all of the peaks will be compared to the image PSF to
546 distinguish stars from galaxies.
547 medianSmoothTemplate: ``bool``, optional
548 If ``medianSmoothTemplate==True`` it a median smoothing filter is applied to the ``maskedImage``.
549 The default is ``True``.
550 medianFilterHalfSize: `int`, optional
551 Half the box size of the median filter, ie a ``medianFilterHalfSize`` of 50 means that
552 each output pixel will be the median of the pixels in a 101 x 101-pixel box in the input image.
553 This parameter is only used when ``medianSmoothTemplate==True``, otherwise it is ignored.
554 The default value is 2.
555 monotonicTempalte: `bool`, optional
556 If True then make the template monotonic.
557 The default is True.
558 weightTemplates: `bool`, optional
559 If True, re-weight the templates so that their linear combination best represents
560 the observed ``maskedImage``.
561 The default is False.
562 log:`lsst.log.Logger` or `lsst.utils.logging.LsstLogAdapter`, optional
563 LSST logger for logging purposes.
564 If `None`, a default logger will be created named after this module.
565 verbose: `bool`, optional
566 Whether or not to show a more verbose output. This option only affects
567 the logger creeated internally and will not change the reporting level
568 of an externally-supplied logger.
569 The default is ``False``.
570 sigma1: `float`, optional
571 Average noise level in ``maskedImage``.
572 The default is ``None``, which estimates the noise from the median value of ``maskedImage``.
573 maxNumberOfPeaks: `int`, optional
574 If nonzero, the maximum number of peaks to deblend.
575 If the total number of peaks is greater than ``maxNumberOfPeaks``,
576 then only the first ``maxNumberOfPeaks`` sources are deblended.
577 The default is 0, which deblends all of the peaks.
578 assignStrayFlux: `bool`, optional
579 If True then flux in the parent footprint that is not covered by any of the
580 template footprints is assigned to templates based on their 1/(1+r^2) distance.
581 How the flux is apportioned is determined by ``strayFluxAssignment``.
582 The default is True.
583 strayFluxToPointSources: `string`
584 Determines how stray flux is apportioned to point sources
586 * ``never``: never apportion stray flux to point sources
587 * ``necessary`` (default): point sources are included only if there are no extended sources nearby
588 * ``always``: point sources are always included in the 1/(1+r^2) splitting
590 strayFluxAssignment: `string`, optional
591 Determines how stray flux is apportioned.
593 * ``trim``: Trim stray flux and do not include in any footprints
594 * ``r-to-peak`` (default): Stray flux is assigned based on (1/(1+r^2) from the peaks
595 * ``r-to-footprint``: Stray flux is distributed to the footprints based on 1/(1+r^2) of the
596 minimum distance from the stray flux to footprint
597 * ``nearest-footprint``: Stray flux is assigned to the footprint with lowest L-1 (Manhattan)
598 distance to the stray flux
600 rampFluxAtEdge: `bool`, optional
601 If True then extend footprints with excessive flux on the edges as described above.
602 The default is False.
603 patchEdges: `bool`, optional
604 If True and if the footprint touches pixels with the ``EDGE`` bit set,
605 then grow the footprint to include all symmetric templates.
606 The default is ``False``.
607 tinyFootprintSize: `float`, optional
608 The PSF model is shrunk to the size that contains the original footprint.
609 If the bbox of the clipped PSF model for a peak is smaller than ``max(tinyFootprintSize,2)``
610 then ``tinyFootprint`` for the peak is set to ``True`` and the peak is not fit.
611 The default is 2.
612 getTemplateSum: `bool`, optional
613 As part of the flux calculation, the sum of the templates is calculated.
614 If ``getTemplateSum==True`` then the sum of the templates is stored in the result (a `PerFootprint`).
615 The default is False.
616 clipStrayFluxFraction: `float`, optional
617 Minimum stray-flux portion.
618 Any stray-flux portion less than ``clipStrayFluxFraction`` is clipped to zero.
619 The default is 0.001.
620 clipFootprintToNonzero: `bool`, optional
621 If True then clip non-zero spans in the template footprints. See above for more.
622 The default is True.
623 removeDegenerateTemplates: `bool`, optional
624 If True then we try to identify "degenerate" peaks by looking at the inner product
625 (in pixel space) of pairs of templates.
626 The default is False.
627 maxTempDotProduct: `float`, optional
628 All dot products between templates greater than ``maxTempDotProduct`` will result in one
629 of the templates removed. This parameter is only used when ``removeDegenerateTempaltes==True``.
630 The default is 0.5.
632 Returns
633 -------
634 res: `PerFootprint`
635 Deblender result that contains a list of ``DeblendedPeak``\ s for each peak and (optionally)
636 the template sum.
637 """
638 avgNoise = sigma1
640 debPlugins = []
642 # Add activated deblender plugins
643 if fitPsfs:
644 debPlugins.append(plugins.DeblenderPlugin(plugins.fitPsfs,
645 psfChisqCut1=psfChisqCut1,
646 psfChisqCut2=psfChisqCut2,
647 psfChisqCut2b=psfChisqCut2b,
648 tinyFootprintSize=tinyFootprintSize))
649 debPlugins.append(plugins.DeblenderPlugin(plugins.buildSymmetricTemplates, patchEdges=patchEdges))
650 if rampFluxAtEdge:
651 debPlugins.append(plugins.DeblenderPlugin(plugins.rampFluxAtEdge, patchEdges=patchEdges))
652 if medianSmoothTemplate:
653 debPlugins.append(plugins.DeblenderPlugin(plugins.medianSmoothTemplates,
654 medianFilterHalfsize=medianFilterHalfsize))
655 if monotonicTemplate:
656 debPlugins.append(plugins.DeblenderPlugin(plugins.makeTemplatesMonotonic))
657 if clipFootprintToNonzero:
658 debPlugins.append(plugins.DeblenderPlugin(plugins.clipFootprintsToNonzero))
659 if weightTemplates:
660 debPlugins.append(plugins.DeblenderPlugin(plugins.weightTemplates))
661 if removeDegenerateTemplates:
662 if weightTemplates:
663 onReset = len(debPlugins)-1
664 else:
665 onReset = len(debPlugins)
666 debPlugins.append(plugins.DeblenderPlugin(plugins.reconstructTemplates,
667 onReset=onReset,
668 maxTempDotProd=maxTempDotProd))
669 debPlugins.append(plugins.DeblenderPlugin(plugins.apportionFlux,
670 clipStrayFluxFraction=clipStrayFluxFraction,
671 assignStrayFlux=assignStrayFlux,
672 strayFluxAssignment=strayFluxAssignment,
673 strayFluxToPointSources=strayFluxToPointSources,
674 getTemplateSum=getTemplateSum))
676 debResult = newDeblend(debPlugins, footprint, maskedImage, psf, psffwhm, log, verbose, avgNoise)
678 return debResult
681def newDeblend(debPlugins, footprint, mMaskedImage, psfs, psfFwhms,
682 log=None, verbose=False, avgNoise=None, maxNumberOfPeaks=0):
683 r"""Deblend a parent ``Footprint`` in a ``MaskedImageF``.
685 Deblending assumes that ``footprint`` has multiple peaks, as it will still create a
686 `PerFootprint` object with a list of peaks even if there is only one peak in the list.
687 It is recommended to first check that ``footprint`` has more than one peak, similar to the
688 execution of `lsst.meas.deblender.deblend.SourceDeblendTask`.
690 This version of the deblender accepts a list of plugins to execute, with the option to re-run parts
691 of the deblender if templates are changed during any of the steps.
693 Parameters
694 ----------
695 debPlugins: list of `meas.deblender.plugins.DeblenderPlugins`
696 Plugins to execute (in order of execution) for the deblender.
697 footprint: `afw.detection.Footprint` or list of Footprints
698 Parent footprint to deblend.
699 mMaskedImage: `MultibandMaskedImage` or `MaskedImage`
700 Masked image in each band.
701 psfs: `afw.detection.Psf` or list of Psfs
702 Psf of the ``maskedImage``.
703 psfFwhms: `float` or list of floats
704 FWHM of the ``maskedImage``\'s ``psf``.
705 log: `lsst.log.Logger` or `lsst.utils.logging.LsstLogAdapter`, optional
706 Logger for logging purposes. Must support a ``trace`` method.
707 If `None`, a default logger will be created named after this module.
708 verbose: `bool`, optional
709 Whether or not to show a more verbose output. This option only affects
710 the logger creeated internally and will not change the reporting level
711 of an externally-supplied logger.
712 The default is ``False``.
713 avgNoise: `float`or list of `float`\ s, optional
714 Average noise level in each ``maskedImage``.
715 The default is ``None``, which estimates the noise from the median value of the
716 variance plane of ``maskedImage`` for each filter.
717 maxNumberOfPeaks: `int`, optional
718 If nonzero, the maximum number of peaks to deblend.
719 If the total number of peaks is greater than ``maxNumberOfPeaks``,
720 then only the first ``maxNumberOfPeaks`` sources are deblended.
722 Returns
723 -------
724 debResult: `DeblendedParent`
725 Deblender result that contains a list of ``MultiColorPeak``\ s for each peak and
726 information that describes the footprint in all filters.
727 """
728 # Import C++ routines
730 if log is None:
731 log = lsst.utils.logging.getLogger(__name__)
733 if verbose:
734 log.setLevel(log.TRACE)
735 else:
736 log.setLevel(log.INFO)
738 # get object that will hold our results
739 debResult = DeblenderResult(footprint, mMaskedImage, psfs, psfFwhms, log,
740 maxNumberOfPeaks=maxNumberOfPeaks, avgNoise=avgNoise)
742 step = 0
743 while step < len(debPlugins):
744 # If a failure occurs at any step,
745 # the result is flagged as `failed`
746 # and the remaining steps are skipped
747 if not debResult.failed:
748 reset = debPlugins[step].run(debResult, log)
749 else:
750 log.warning("Skipping steps %s", debPlugins[step:])
751 return debResult
752 if reset:
753 step = debPlugins[step].onReset
754 else:
755 step += 1
757 return debResult
760class CachingPsf:
761 """Cache the PSF models
763 In the PSF fitting code, we request PSF models for all peaks near
764 the one being fit. This was turning out to be quite expensive in
765 some cases. Here, we cache the PSF models to bring the cost down
766 closer to O(N) rather than O(N^2).
767 """
769 def __init__(self, psf):
770 self.cache = {}
771 self.psf = psf
773 def computeImage(self, cx, cy):
774 im = self.cache.get((cx, cy), None)
775 if im is not None:
776 return im
777 try:
778 im = self.psf.computeImage(geom.Point2D(cx, cy))
779 except lsst.pex.exceptions.Exception:
780 im = self.psf.computeImage(self.psf.getAveragePosition())
781 self.cache[(cx, cy)] = im
782 return im