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