LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
sourceDeblendTask.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 __all__ = ['SourceDeblendConfig', 'SourceDeblendTask']
23 
24 import math
25 import numpy as np
26 
27 import lsst.pex.config as pexConfig
28 import lsst.pipe.base as pipeBase
29 import lsst.afw.math as afwMath
30 import lsst.geom as geom
31 import lsst.afw.geom.ellipses as afwEll
32 import lsst.afw.image as afwImage
33 import lsst.afw.detection as afwDet
34 import lsst.afw.table as afwTable
35 from lsst.utils.timer import timeMethod
36 
37 
38 class SourceDeblendConfig(pexConfig.Config):
39 
40  edgeHandling = pexConfig.ChoiceField(
41  doc='What to do when a peak to be deblended is close to the edge of the image',
42  dtype=str, default='ramp',
43  allowed={
44  'clip': 'Clip the template at the edge AND the mirror of the edge.',
45  'ramp': 'Ramp down flux at the image edge by the PSF',
46  'noclip': 'Ignore the edge when building the symmetric template.',
47  }
48  )
49 
50  strayFluxToPointSources = pexConfig.ChoiceField(
51  doc='When the deblender should attribute stray flux to point sources',
52  dtype=str, default='necessary',
53  allowed={
54  'necessary': 'When there is not an extended object in the footprint',
55  'always': 'Always',
56  'never': ('Never; stray flux will not be attributed to any deblended child '
57  'if the deblender thinks all peaks look like point sources'),
58  }
59  )
60 
61  assignStrayFlux = pexConfig.Field(dtype=bool, default=True,
62  doc='Assign stray flux (not claimed by any child in the deblender) '
63  'to deblend children.')
64 
65  strayFluxRule = pexConfig.ChoiceField(
66  doc='How to split flux among peaks',
67  dtype=str, default='trim',
68  allowed={
69  'r-to-peak': '~ 1/(1+R^2) to the peak',
70  'r-to-footprint': ('~ 1/(1+R^2) to the closest pixel in the footprint. '
71  'CAUTION: this can be computationally expensive on large footprints!'),
72  'nearest-footprint': ('Assign 100% to the nearest footprint (using L-1 norm aka '
73  'Manhattan distance)'),
74  'trim': ('Shrink the parent footprint to pixels that are not assigned to children')
75  }
76  )
77 
78  clipStrayFluxFraction = pexConfig.Field(dtype=float, default=0.001,
79  doc=('When splitting stray flux, clip fractions below '
80  'this value to zero.'))
81  psfChisq1 = pexConfig.Field(dtype=float, default=1.5, optional=False,
82  doc=('Chi-squared per DOF cut for deciding a source is '
83  'a PSF during deblending (un-shifted PSF model)'))
84  psfChisq2 = pexConfig.Field(dtype=float, default=1.5, optional=False,
85  doc=('Chi-squared per DOF cut for deciding a source is '
86  'PSF during deblending (shifted PSF model)'))
87  psfChisq2b = pexConfig.Field(dtype=float, default=1.5, optional=False,
88  doc=('Chi-squared per DOF cut for deciding a source is '
89  'a PSF during deblending (shifted PSF model #2)'))
90  maxNumberOfPeaks = pexConfig.Field(dtype=int, default=0,
91  doc=("Only deblend the brightest maxNumberOfPeaks peaks in the parent"
92  " (<= 0: unlimited)"))
93  maxFootprintArea = pexConfig.Field(dtype=int, default=1000000,
94  doc=("Maximum area for footprints before they are ignored as large; "
95  "non-positive means no threshold applied"))
96  maxFootprintSize = pexConfig.Field(dtype=int, default=0,
97  doc=("Maximum linear dimension for footprints before they are ignored "
98  "as large; non-positive means no threshold applied"))
99  minFootprintAxisRatio = pexConfig.Field(dtype=float, default=0.0,
100  doc=("Minimum axis ratio for footprints before they are ignored "
101  "as large; non-positive means no threshold applied"))
102  notDeblendedMask = pexConfig.Field(dtype=str, default="NOT_DEBLENDED", optional=True,
103  doc="Mask name for footprints not deblended, or None")
104 
105  tinyFootprintSize = pexConfig.RangeField(dtype=int, default=2, min=2, inclusiveMin=True,
106  doc=('Footprints smaller in width or height than this value '
107  'will be ignored; minimum of 2 due to PSF gradient '
108  'calculation.'))
109 
110  propagateAllPeaks = pexConfig.Field(dtype=bool, default=False,
111  doc=('Guarantee that all peaks produce a child source.'))
112  catchFailures = pexConfig.Field(
113  dtype=bool, default=False,
114  doc=("If True, catch exceptions thrown by the deblender, log them, "
115  "and set a flag on the parent, instead of letting them propagate up"))
116  maskPlanes = pexConfig.ListField(dtype=str, default=["SAT", "INTRP", "NO_DATA"],
117  doc="Mask planes to ignore when performing statistics")
118  maskLimits = pexConfig.DictField(
119  keytype=str,
120  itemtype=float,
121  default={},
122  doc=("Mask planes with the corresponding limit on the fraction of masked pixels. "
123  "Sources violating this limit will not be deblended."),
124  )
125  weightTemplates = pexConfig.Field(
126  dtype=bool, default=False,
127  doc=("If true, a least-squares fit of the templates will be done to the "
128  "full image. The templates will be re-weighted based on this fit."))
129  removeDegenerateTemplates = pexConfig.Field(dtype=bool, default=False,
130  doc=("Try to remove similar templates?"))
131  maxTempDotProd = pexConfig.Field(
132  dtype=float, default=0.5,
133  doc=("If the dot product between two templates is larger than this value, we consider them to be "
134  "describing the same object (i.e. they are degenerate). If one of the objects has been "
135  "labeled as a PSF it will be removed, otherwise the template with the lowest value will "
136  "be removed."))
137  medianSmoothTemplate = pexConfig.Field(dtype=bool, default=True,
138  doc="Apply a smoothing filter to all of the template images")
139 
140  # Testing options
141  # Some obs packages and ci packages run the full pipeline on a small
142  # subset of data to test that the pipeline is functioning properly.
143  # This is not meant as scientific validation, so it can be useful
144  # to only run on a small subset of the data that is large enough to
145  # test the desired pipeline features but not so long that the deblender
146  # is the tall pole in terms of execution times.
147  useCiLimits = pexConfig.Field(
148  dtype=bool, default=False,
149  doc="Limit the number of sources deblended for CI to prevent long build times")
150  ciDeblendChildRange = pexConfig.ListField(
151  dtype=int, default=[2, 10],
152  doc="Only deblend parent Footprints with a number of peaks in the (inclusive) range indicated."
153  "If `useCiLimits==False` then this parameter is ignored.")
154  ciNumParentsToDeblend = pexConfig.Field(
155  dtype=int, default=10,
156  doc="Only use the first `ciNumParentsToDeblend` parent footprints with a total peak count "
157  "within `ciDebledChildRange`. "
158  "If `useCiLimits==False` then this parameter is ignored.")
159 
160 
161 class SourceDeblendTask(pipeBase.Task):
162  """Split blended sources into individual sources.
163 
164  This task has no return value; it only modifies the SourceCatalog in-place.
165  """
166  ConfigClass = SourceDeblendConfig
167  _DefaultName = "sourceDeblend"
168 
169  def __init__(self, schema, peakSchema=None, **kwargs):
170  """Create the task, adding necessary fields to the given schema.
171 
172  Parameters
173  ----------
174  schema : `lsst.afw.table.Schema`
175  Schema object for measurement fields; will be modified in-place.
176  peakSchema : `lsst.afw.table.peakSchema`
177  Schema of Footprint Peaks that will be passed to the deblender.
178  Any fields beyond the PeakTable minimal schema will be transferred
179  to the main source Schema. If None, no fields will be transferred
180  from the Peaks
181  **kwargs
182  Additional keyword arguments passed to ~lsst.pipe.base.task
183  """
184  pipeBase.Task.__init__(self, **kwargs)
185  self.schemaschema = schema
186  self.toCopyFromParenttoCopyFromParent = [item.key for item in self.schemaschema
187  if item.field.getName().startswith("merge_footprint")]
188  peakMinimalSchema = afwDet.PeakTable.makeMinimalSchema()
189  if peakSchema is None:
190  # In this case, the peakSchemaMapper will transfer nothing, but we'll still have one
191  # to simplify downstream code
192  self.peakSchemaMapperpeakSchemaMapper = afwTable.SchemaMapper(peakMinimalSchema, schema)
193  else:
194  self.peakSchemaMapperpeakSchemaMapper = afwTable.SchemaMapper(peakSchema, schema)
195  for item in peakSchema:
196  if item.key not in peakMinimalSchema:
197  self.peakSchemaMapperpeakSchemaMapper.addMapping(item.key, item.field)
198  # Because SchemaMapper makes a copy of the output schema you give its ctor, it isn't
199  # updating this Schema in place. That's probably a design flaw, but in the meantime,
200  # we'll keep that schema in sync with the peakSchemaMapper.getOutputSchema() manually,
201  # by adding the same fields to both.
202  schema.addField(item.field)
203  assert schema == self.peakSchemaMapperpeakSchemaMapper.getOutputSchema(), "Logic bug mapping schemas"
204  self.addSchemaKeysaddSchemaKeys(schema)
205 
206  def addSchemaKeys(self, schema):
207  self.nChildKeynChildKey = schema.addField('deblend_nChild', type=np.int32,
208  doc='Number of children this object has (defaults to 0)')
209  self.psfKeypsfKey = schema.addField('deblend_deblendedAsPsf', type='Flag',
210  doc='Deblender thought this source looked like a PSF')
211  self.psfCenterKeypsfCenterKey = afwTable.Point2DKey.addFields(schema, 'deblend_psfCenter',
212  'If deblended-as-psf, the PSF centroid', "pixel")
213  self.psfFluxKeypsfFluxKey = schema.addField('deblend_psf_instFlux', type='D',
214  doc='If deblended-as-psf, the instrumental PSF flux', units='count')
215  self.tooManyPeaksKeytooManyPeaksKey = schema.addField('deblend_tooManyPeaks', type='Flag',
216  doc='Source had too many peaks; '
217  'only the brightest were included')
218  self.tooBigKeytooBigKey = schema.addField('deblend_parentTooBig', type='Flag',
219  doc='Parent footprint covered too many pixels')
220  self.maskedKeymaskedKey = schema.addField('deblend_masked', type='Flag',
221  doc='Parent footprint was predominantly masked')
222 
223  if self.config.catchFailures:
224  self.deblendFailedKeydeblendFailedKey = schema.addField('deblend_failed', type='Flag',
225  doc="Deblending failed on source")
226 
227  self.deblendSkippedKeydeblendSkippedKey = schema.addField('deblend_skipped', type='Flag',
228  doc="Deblender skipped this source")
229 
230  self.deblendRampedTemplateKeydeblendRampedTemplateKey = schema.addField(
231  'deblend_rampedTemplate', type='Flag',
232  doc=('This source was near an image edge and the deblender used '
233  '"ramp" edge-handling.'))
234 
235  self.deblendPatchedTemplateKeydeblendPatchedTemplateKey = schema.addField(
236  'deblend_patchedTemplate', type='Flag',
237  doc=('This source was near an image edge and the deblender used '
238  '"patched" edge-handling.'))
239 
240  self.hasStrayFluxKeyhasStrayFluxKey = schema.addField(
241  'deblend_hasStrayFlux', type='Flag',
242  doc=('This source was assigned some stray flux'))
243 
244  self.log.trace('Added keys to schema: %s', ", ".join(str(x) for x in (
245  self.nChildKeynChildKey, self.psfKeypsfKey, self.psfCenterKeypsfCenterKey, self.psfFluxKeypsfFluxKey,
246  self.tooManyPeaksKeytooManyPeaksKey, self.tooBigKeytooBigKey)))
247  self.peakCenterpeakCenter = afwTable.Point2IKey.addFields(schema, name="deblend_peak_center",
248  doc="Center used to apply constraints in scarlet",
249  unit="pixel")
250  self.peakIdKeypeakIdKey = schema.addField("deblend_peakId", type=np.int32,
251  doc="ID of the peak in the parent footprint. "
252  "This is not unique, but the combination of 'parent'"
253  "and 'peakId' should be for all child sources. "
254  "Top level blends with no parents have 'peakId=0'")
255  self.nPeaksKeynPeaksKey = schema.addField("deblend_nPeaks", type=np.int32,
256  doc="Number of initial peaks in the blend. "
257  "This includes peaks that may have been culled "
258  "during deblending or failed to deblend")
259  self.parentNPeaksKeyparentNPeaksKey = schema.addField("deblend_parentNPeaks", type=np.int32,
260  doc="Same as deblend_n_peaks, but the number of peaks "
261  "in the parent footprint")
262 
263  @timeMethod
264  def run(self, exposure, sources):
265  """Get the PSF from the provided exposure and then run deblend.
266 
267  Parameters
268  ----------
269  exposure : `lsst.afw.image.Exposure`
270  Exposure to be processed
271  sources : `lsst.afw.table.SourceCatalog`
272  SourceCatalog containing sources detected on this exposure.
273  """
274  psf = exposure.getPsf()
275  assert sources.getSchema() == self.schemaschema
276  self.deblenddeblend(exposure, sources, psf)
277 
278  def _getPsfFwhm(self, psf, bbox):
279  # It should be easier to get a PSF's fwhm;
280  # https://dev.lsstcorp.org/trac/ticket/3030
281  return psf.computeShape().getDeterminantRadius() * 2.35
282 
283  @timeMethod
284  def deblend(self, exposure, srcs, psf):
285  """Deblend.
286 
287  Parameters
288  ----------
289  exposure : `lsst.afw.image.Exposure`
290  Exposure to be processed
291  srcs : `lsst.afw.table.SourceCatalog`
292  SourceCatalog containing sources detected on this exposure
293  psf : `lsst.afw.detection.Psf`
294  Point source function
295 
296  Returns
297  -------
298  None
299  """
300  # Cull footprints if required by ci
301  if self.config.useCiLimits:
302  self.log.info(f"Using CI catalog limits, "
303  f"the original number of sources to deblend was {len(srcs)}.")
304  # Select parents with a number of children in the range
305  # config.ciDeblendChildRange
306  minChildren, maxChildren = self.config.ciDeblendChildRange
307  nPeaks = np.array([len(src.getFootprint().peaks) for src in srcs])
308  childrenInRange = np.where((nPeaks >= minChildren) & (nPeaks <= maxChildren))[0]
309  if len(childrenInRange) < self.config.ciNumParentsToDeblend:
310  raise ValueError("Fewer than ciNumParentsToDeblend children were contained in the range "
311  "indicated by ciDeblendChildRange. Adjust this range to include more "
312  "parents.")
313  # Keep all of the isolated parents and the first
314  # `ciNumParentsToDeblend` children
315  parents = nPeaks == 1
316  children = np.zeros((len(srcs),), dtype=bool)
317  children[childrenInRange[:self.config.ciNumParentsToDeblend]] = True
318  srcs = srcs[parents | children]
319  # We need to update the IdFactory, otherwise the the source ids
320  # will not be sequential
321  idFactory = srcs.getIdFactory()
322  maxId = np.max(srcs["id"])
323  idFactory.notify(maxId)
324 
325  self.log.info("Deblending %d sources", len(srcs))
326 
327  from lsst.meas.deblender.baseline import deblend
328 
329  # find the median stdev in the image...
330  mi = exposure.getMaskedImage()
331  statsCtrl = afwMath.StatisticsControl()
332  statsCtrl.setAndMask(mi.getMask().getPlaneBitMask(self.config.maskPlanes))
333  stats = afwMath.makeStatistics(mi.getVariance(), mi.getMask(), afwMath.MEDIAN, statsCtrl)
334  sigma1 = math.sqrt(stats.getValue(afwMath.MEDIAN))
335  self.log.trace('sigma1: %g', sigma1)
336 
337  n0 = len(srcs)
338  nparents = 0
339  for i, src in enumerate(srcs):
340  # t0 = time.clock()
341 
342  fp = src.getFootprint()
343  pks = fp.getPeaks()
344 
345  # Since we use the first peak for the parent object, we should propagate its flags
346  # to the parent source.
347  src.assign(pks[0], self.peakSchemaMapperpeakSchemaMapper)
348 
349  if len(pks) < 2:
350  continue
351 
352  if self.isLargeFootprintisLargeFootprint(fp):
353  src.set(self.tooBigKeytooBigKey, True)
354  self.skipParentskipParent(src, mi.getMask())
355  self.log.warning('Parent %i: skipping large footprint (area: %i)',
356  int(src.getId()), int(fp.getArea()))
357  continue
358  if self.isMaskedisMasked(fp, exposure.getMaskedImage().getMask()):
359  src.set(self.maskedKeymaskedKey, True)
360  self.skipParentskipParent(src, mi.getMask())
361  self.log.warning('Parent %i: skipping masked footprint (area: %i)',
362  int(src.getId()), int(fp.getArea()))
363  continue
364 
365  nparents += 1
366  bb = fp.getBBox()
367  psf_fwhm = self._getPsfFwhm_getPsfFwhm(psf, bb)
368 
369  self.log.trace('Parent %i: deblending %i peaks', int(src.getId()), len(pks))
370 
371  self.preSingleDeblendHookpreSingleDeblendHook(exposure, srcs, i, fp, psf, psf_fwhm, sigma1)
372  npre = len(srcs)
373 
374  # This should really be set in deblend, but deblend doesn't have access to the src
375  src.set(self.tooManyPeaksKeytooManyPeaksKey, len(fp.getPeaks()) > self.config.maxNumberOfPeaks)
376 
377  try:
378  res = deblend(
379  fp, mi, psf, psf_fwhm, sigma1=sigma1,
380  psfChisqCut1=self.config.psfChisq1,
381  psfChisqCut2=self.config.psfChisq2,
382  psfChisqCut2b=self.config.psfChisq2b,
383  maxNumberOfPeaks=self.config.maxNumberOfPeaks,
384  strayFluxToPointSources=self.config.strayFluxToPointSources,
385  assignStrayFlux=self.config.assignStrayFlux,
386  strayFluxAssignment=self.config.strayFluxRule,
387  rampFluxAtEdge=(self.config.edgeHandling == 'ramp'),
388  patchEdges=(self.config.edgeHandling == 'noclip'),
389  tinyFootprintSize=self.config.tinyFootprintSize,
390  clipStrayFluxFraction=self.config.clipStrayFluxFraction,
391  weightTemplates=self.config.weightTemplates,
392  removeDegenerateTemplates=self.config.removeDegenerateTemplates,
393  maxTempDotProd=self.config.maxTempDotProd,
394  medianSmoothTemplate=self.config.medianSmoothTemplate
395  )
396  if self.config.catchFailures:
397  src.set(self.deblendFailedKeydeblendFailedKey, False)
398  except Exception as e:
399  if self.config.catchFailures:
400  self.log.warning("Unable to deblend source %d: %s", src.getId(), e)
401  src.set(self.deblendFailedKeydeblendFailedKey, True)
402  import traceback
403  traceback.print_exc()
404  continue
405  else:
406  raise
407 
408  kids = []
409  nchild = 0
410  for j, peak in enumerate(res.deblendedParents[0].peaks):
411  heavy = peak.getFluxPortion()
412  if heavy is None or peak.skip:
413  src.set(self.deblendSkippedKeydeblendSkippedKey, True)
414  if not self.config.propagateAllPeaks:
415  # Don't care
416  continue
417  # We need to preserve the peak: make sure we have enough info to create a minimal
418  # child src
419  self.log.trace("Peak at (%i,%i) failed. Using minimal default info for child.",
420  pks[j].getIx(), pks[j].getIy())
421  if heavy is None:
422  # copy the full footprint and strip out extra peaks
423  foot = afwDet.Footprint(src.getFootprint())
424  peakList = foot.getPeaks()
425  peakList.clear()
426  peakList.append(peak.peak)
427  zeroMimg = afwImage.MaskedImageF(foot.getBBox())
428  heavy = afwDet.makeHeavyFootprint(foot, zeroMimg)
429  if peak.deblendedAsPsf:
430  if peak.psfFitFlux is None:
431  peak.psfFitFlux = 0.0
432  if peak.psfFitCenter is None:
433  peak.psfFitCenter = (peak.peak.getIx(), peak.peak.getIy())
434 
435  assert(len(heavy.getPeaks()) == 1)
436 
437  src.set(self.deblendSkippedKeydeblendSkippedKey, False)
438  child = srcs.addNew()
439  nchild += 1
440  for key in self.toCopyFromParenttoCopyFromParent:
441  child.set(key, src.get(key))
442  child.assign(heavy.getPeaks()[0], self.peakSchemaMapperpeakSchemaMapper)
443  child.setParent(src.getId())
444  child.setFootprint(heavy)
445  child.set(self.psfKeypsfKey, peak.deblendedAsPsf)
446  child.set(self.hasStrayFluxKeyhasStrayFluxKey, peak.strayFlux is not None)
447  if peak.deblendedAsPsf:
448  (cx, cy) = peak.psfFitCenter
449  child.set(self.psfCenterKeypsfCenterKey, geom.Point2D(cx, cy))
450  child.set(self.psfFluxKeypsfFluxKey, peak.psfFitFlux)
451  child.set(self.deblendRampedTemplateKeydeblendRampedTemplateKey, peak.hasRampedTemplate)
452  child.set(self.deblendPatchedTemplateKeydeblendPatchedTemplateKey, peak.patched)
453 
454  # Set the position of the peak from the parent footprint
455  # This will make it easier to match the same source across
456  # deblenders and across observations, where the peak
457  # position is unlikely to change unless enough time passes
458  # for a source to move on the sky.
459  child.set(self.peakCenterpeakCenter, geom.Point2I(pks[j].getIx(), pks[j].getIy()))
460  child.set(self.peakIdKeypeakIdKey, pks[j].getId())
461 
462  # The children have a single peak
463  child.set(self.nPeaksKeynPeaksKey, 1)
464  # Set the number of peaks in the parent
465  child.set(self.parentNPeaksKeyparentNPeaksKey, len(pks))
466 
467  kids.append(child)
468 
469  # Child footprints may extend beyond the full extent of their parent's which
470  # results in a failure of the replace-by-noise code to reinstate these pixels
471  # to their original values. The following updates the parent footprint
472  # in-place to ensure it contains the full union of itself and all of its
473  # children's footprints.
474  spans = src.getFootprint().spans
475  for child in kids:
476  spans = spans.union(child.getFootprint().spans)
477  src.getFootprint().setSpans(spans)
478 
479  src.set(self.nChildKeynChildKey, nchild)
480 
481  self.postSingleDeblendHookpostSingleDeblendHook(exposure, srcs, i, npre, kids, fp, psf, psf_fwhm, sigma1, res)
482  # print('Deblending parent id', src.getId(), 'took', time.clock() - t0)
483 
484  n1 = len(srcs)
485  self.log.info('Deblended: of %i sources, %i were deblended, creating %i children, total %i sources',
486  n0, nparents, n1-n0, n1)
487 
488  def preSingleDeblendHook(self, exposure, srcs, i, fp, psf, psf_fwhm, sigma1):
489  pass
490 
491  def postSingleDeblendHook(self, exposure, srcs, i, npre, kids, fp, psf, psf_fwhm, sigma1, res):
492  pass
493 
494  def isLargeFootprint(self, footprint):
495  """Returns whether a Footprint is large
496 
497  'Large' is defined by thresholds on the area, size and axis ratio.
498  These may be disabled independently by configuring them to be non-positive.
499 
500  This is principally intended to get rid of satellite streaks, which the
501  deblender or other downstream processing can have trouble dealing with
502  (e.g., multiple large HeavyFootprints can chew up memory).
503  """
504  if self.config.maxFootprintArea > 0 and footprint.getArea() > self.config.maxFootprintArea:
505  return True
506  if self.config.maxFootprintSize > 0:
507  bbox = footprint.getBBox()
508  if max(bbox.getWidth(), bbox.getHeight()) > self.config.maxFootprintSize:
509  return True
510  if self.config.minFootprintAxisRatio > 0:
511  axes = afwEll.Axes(footprint.getShape())
512  if axes.getB() < self.config.minFootprintAxisRatio*axes.getA():
513  return True
514  return False
515 
516  def isMasked(self, footprint, mask):
517  """Returns whether the footprint violates the mask limits
518  """
519  size = float(footprint.getArea())
520  for maskName, limit in self.config.maskLimits.items():
521  maskVal = mask.getPlaneBitMask(maskName)
522  unmaskedSpan = footprint.spans.intersectNot(mask, maskVal) # spanset of unmasked pixels
523  if (size - unmaskedSpan.getArea())/size > limit:
524  return True
525  return False
526 
527  def skipParent(self, source, mask):
528  """Indicate that the parent source is not being deblended
529 
530  We set the appropriate flags and mask.
531 
532  Parameters
533  ----------
534  source : `lsst.afw.table.SourceRecord`
535  The source to flag as skipped
536  mask : `lsst.afw.image.Mask`
537  The mask to update
538  """
539  fp = source.getFootprint()
540  source.set(self.deblendSkippedKeydeblendSkippedKey, True)
541  if self.config.notDeblendedMask:
542  mask.addMaskPlane(self.config.notDeblendedMask)
543  fp.spans.setMask(mask, mask.getPlaneBitMask(self.config.notDeblendedMask))
544 
545  # Set the center of the parent
546  bbox = fp.getBBox()
547  centerX = int(bbox.getMinX()+bbox.getWidth()/2)
548  centerY = int(bbox.getMinY()+bbox.getHeight()/2)
549  source.set(self.peakCenterpeakCenter, geom.Point2I(centerX, centerY))
550  # There are no deblended children, so nChild = 0
551  source.set(self.nChildKeynChildKey, 0)
552  # But we also want to know how many peaks that we would have
553  # deblended if the parent wasn't skipped.
554  source.set(self.nPeaksKeynPeaksKey, len(fp.peaks))
555  # Top level parents are not a detected peak, so they have no peakId
556  source.set(self.peakIdKeypeakIdKey, 0)
557  # Top level parents also have no parentNPeaks
558  source.set(self.parentNPeaksKeyparentNPeaksKey, 0)
int max
Class to describe the properties of a detected object from an image.
Definition: Footprint.h:63
Pass parameters to a Statistics object.
Definition: Statistics.h:92
A mapping between the keys of two Schemas, used to copy data between them.
Definition: SchemaMapper.h:21
def postSingleDeblendHook(self, exposure, srcs, i, npre, kids, fp, psf, psf_fwhm, sigma1, res)
def preSingleDeblendHook(self, exposure, srcs, i, fp, psf, psf_fwhm, sigma1)
def __init__(self, schema, peakSchema=None, **kwargs)
HeavyFootprint< ImagePixelT, MaskPixelT, VariancePixelT > makeHeavyFootprint(Footprint const &foot, lsst::afw::image::MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > const &img, HeavyFootprintCtrl const *ctrl=nullptr)
Create a HeavyFootprint with footprint defined by the given Footprint and pixel values from the given...
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
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:359