22 __all__ = [
'SourceDeblendConfig',
'SourceDeblendTask']
42 edgeHandling = pexConfig.ChoiceField(
43 doc=
'What to do when a peak to be deblended is close to the edge of the image',
44 dtype=str, default=
'ramp',
46 'clip':
'Clip the template at the edge AND the mirror of the edge.',
47 'ramp':
'Ramp down flux at the image edge by the PSF',
48 'noclip':
'Ignore the edge when building the symmetric template.',
52 strayFluxToPointSources = pexConfig.ChoiceField(
53 doc=
'When the deblender should attribute stray flux to point sources',
54 dtype=str, default=
'necessary',
56 'necessary':
'When there is not an extended object in the footprint',
58 'never': (
'Never; stray flux will not be attributed to any deblended child '
59 'if the deblender thinks all peaks look like point sources'),
63 assignStrayFlux = pexConfig.Field(dtype=bool, default=
True,
64 doc=
'Assign stray flux (not claimed by any child in the deblender) '
65 'to deblend children.')
67 strayFluxRule = pexConfig.ChoiceField(
68 doc=
'How to split flux among peaks',
69 dtype=str, default=
'trim',
71 'r-to-peak':
'~ 1/(1+R^2) to the peak',
72 'r-to-footprint': (
'~ 1/(1+R^2) to the closest pixel in the footprint. '
73 'CAUTION: this can be computationally expensive on large footprints!'),
74 'nearest-footprint': (
'Assign 100% to the nearest footprint (using L-1 norm aka '
75 'Manhattan distance)'),
76 'trim': (
'Shrink the parent footprint to pixels that are not assigned to children')
80 clipStrayFluxFraction = pexConfig.Field(dtype=float, default=0.001,
81 doc=(
'When splitting stray flux, clip fractions below '
82 'this value to zero.'))
83 psfChisq1 = pexConfig.Field(dtype=float, default=1.5, optional=
False,
84 doc=(
'Chi-squared per DOF cut for deciding a source is '
85 'a PSF during deblending (un-shifted PSF model)'))
86 psfChisq2 = pexConfig.Field(dtype=float, default=1.5, optional=
False,
87 doc=(
'Chi-squared per DOF cut for deciding a source is '
88 'PSF during deblending (shifted PSF model)'))
89 psfChisq2b = pexConfig.Field(dtype=float, default=1.5, optional=
False,
90 doc=(
'Chi-squared per DOF cut for deciding a source is '
91 'a PSF during deblending (shifted PSF model #2)'))
92 maxNumberOfPeaks = pexConfig.Field(dtype=int, default=0,
93 doc=(
"Only deblend the brightest maxNumberOfPeaks peaks in the parent"
94 " (<= 0: unlimited)"))
95 maxFootprintArea = pexConfig.Field(dtype=int, default=1000000,
96 doc=(
"Maximum area for footprints before they are ignored as large; "
97 "non-positive means no threshold applied"))
98 maxFootprintSize = pexConfig.Field(dtype=int, default=0,
99 doc=(
"Maximum linear dimension for footprints before they are ignored "
100 "as large; non-positive means no threshold applied"))
101 minFootprintAxisRatio = pexConfig.Field(dtype=float, default=0.0,
102 doc=(
"Minimum axis ratio for footprints before they are ignored "
103 "as large; non-positive means no threshold applied"))
104 notDeblendedMask = pexConfig.Field(dtype=str, default=
"NOT_DEBLENDED", optional=
True,
105 doc=
"Mask name for footprints not deblended, or None")
107 tinyFootprintSize = pexConfig.RangeField(dtype=int, default=2, min=2, inclusiveMin=
True,
108 doc=(
'Footprints smaller in width or height than this value '
109 'will be ignored; minimum of 2 due to PSF gradient '
112 propagateAllPeaks = pexConfig.Field(dtype=bool, default=
False,
113 doc=(
'Guarantee that all peaks produce a child source.'))
114 catchFailures = pexConfig.Field(
115 dtype=bool, default=
False,
116 doc=(
"If True, catch exceptions thrown by the deblender, log them, "
117 "and set a flag on the parent, instead of letting them propagate up"))
118 maskPlanes = pexConfig.ListField(dtype=str, default=[
"SAT",
"INTRP",
"NO_DATA"],
119 doc=
"Mask planes to ignore when performing statistics")
120 maskLimits = pexConfig.DictField(
124 doc=(
"Mask planes with the corresponding limit on the fraction of masked pixels. "
125 "Sources violating this limit will not be deblended."),
127 weightTemplates = pexConfig.Field(
128 dtype=bool, default=
False,
129 doc=(
"If true, a least-squares fit of the templates will be done to the "
130 "full image. The templates will be re-weighted based on this fit."))
131 removeDegenerateTemplates = pexConfig.Field(dtype=bool, default=
False,
132 doc=(
"Try to remove similar templates?"))
133 maxTempDotProd = pexConfig.Field(
134 dtype=float, default=0.5,
135 doc=(
"If the dot product between two templates is larger than this value, we consider them to be "
136 "describing the same object (i.e. they are degenerate). If one of the objects has been "
137 "labeled as a PSF it will be removed, otherwise the template with the lowest value will "
139 medianSmoothTemplate = pexConfig.Field(dtype=bool, default=
True,
140 doc=
"Apply a smoothing filter to all of the template images")
149 useCiLimits = pexConfig.Field(
150 dtype=bool, default=
False,
151 doc=
"Limit the number of sources deblended for CI to prevent long build times")
152 ciDeblendChildRange = pexConfig.ListField(
153 dtype=int, default=[2, 10],
154 doc=
"Only deblend parent Footprints with a number of peaks in the (inclusive) range indicated."
155 "If `useCiLimits==False` then this parameter is ignored.")
156 ciNumParentsToDeblend = pexConfig.Field(
157 dtype=int, default=10,
158 doc=
"Only use the first `ciNumParentsToDeblend` parent footprints with a total peak count "
159 "within `ciDebledChildRange`. "
160 "If `useCiLimits==False` then this parameter is ignored.")
164 """Split blended sources into individual sources.
166 This task has no return value; it only modifies the SourceCatalog in-place.
168 ConfigClass = SourceDeblendConfig
169 _DefaultName =
"sourceDeblend"
171 def __init__(self, schema, peakSchema=None, **kwargs):
172 """Create the task, adding necessary fields to the given schema.
176 schema : `lsst.afw.table.Schema`
177 Schema object for measurement fields; will be modified in-place.
178 peakSchema : `lsst.afw.table.peakSchema`
179 Schema of Footprint Peaks that will be passed to the deblender.
180 Any fields beyond the PeakTable minimal schema will be transferred
181 to the main source Schema. If None, no fields will be transferred
184 Additional keyword arguments passed to ~lsst.pipe.base.task
186 pipeBase.Task.__init__(self, **kwargs)
189 if item.field.getName().startswith(
"merge_footprint")]
190 peakMinimalSchema = afwDet.PeakTable.makeMinimalSchema()
191 if peakSchema
is None:
197 for item
in peakSchema:
198 if item.key
not in peakMinimalSchema:
204 schema.addField(item.field)
205 assert schema == self.
peakSchemaMapperpeakSchemaMapper.getOutputSchema(),
"Logic bug mapping schemas"
209 self.
nChildKeynChildKey = schema.addField(
'deblend_nChild', type=np.int32,
210 doc=
'Number of children this object has (defaults to 0)')
211 self.
psfKeypsfKey = schema.addField(
'deblend_deblendedAsPsf', type=
'Flag',
212 doc=
'Deblender thought this source looked like a PSF')
213 self.
psfCenterKeypsfCenterKey = afwTable.Point2DKey.addFields(schema,
'deblend_psfCenter',
214 'If deblended-as-psf, the PSF centroid',
"pixel")
215 self.
psfFluxKeypsfFluxKey = schema.addField(
'deblend_psf_instFlux', type=
'D',
216 doc=
'If deblended-as-psf, the instrumental PSF flux', units=
'count')
217 self.
tooManyPeaksKeytooManyPeaksKey = schema.addField(
'deblend_tooManyPeaks', type=
'Flag',
218 doc=
'Source had too many peaks; '
219 'only the brightest were included')
220 self.
tooBigKeytooBigKey = schema.addField(
'deblend_parentTooBig', type=
'Flag',
221 doc=
'Parent footprint covered too many pixels')
222 self.
maskedKeymaskedKey = schema.addField(
'deblend_masked', type=
'Flag',
223 doc=
'Parent footprint was predominantly masked')
225 if self.config.catchFailures:
227 doc=
"Deblending failed on source")
230 doc=
"Deblender skipped this source")
233 'deblend_rampedTemplate', type=
'Flag',
234 doc=(
'This source was near an image edge and the deblender used '
235 '"ramp" edge-handling.'))
238 'deblend_patchedTemplate', type=
'Flag',
239 doc=(
'This source was near an image edge and the deblender used '
240 '"patched" edge-handling.'))
243 'deblend_hasStrayFlux', type=
'Flag',
244 doc=(
'This source was assigned some stray flux'))
246 self.log.
trace(
'Added keys to schema: %s',
", ".join(str(x)
for x
in (
249 self.
peakCenterpeakCenter = afwTable.Point2IKey.addFields(schema, name=
"deblend_peak_center",
250 doc=
"Center used to apply constraints in scarlet",
252 self.
peakIdKeypeakIdKey = schema.addField(
"deblend_peakId", type=np.int32,
253 doc=
"ID of the peak in the parent footprint. "
254 "This is not unique, but the combination of 'parent'"
255 "and 'peakId' should be for all child sources. "
256 "Top level blends with no parents have 'peakId=0'")
257 self.
nPeaksKeynPeaksKey = schema.addField(
"deblend_nPeaks", type=np.int32,
258 doc=
"Number of initial peaks in the blend. "
259 "This includes peaks that may have been culled "
260 "during deblending or failed to deblend")
261 self.
parentNPeaksKeyparentNPeaksKey = schema.addField(
"deblend_parentNPeaks", type=np.int32,
262 doc=
"Same as deblend_n_peaks, but the number of peaks "
263 "in the parent footprint")
266 def run(self, exposure, sources):
267 """Get the PSF from the provided exposure and then run deblend.
271 exposure : `lsst.afw.image.Exposure`
272 Exposure to be processed
273 sources : `lsst.afw.table.SourceCatalog`
274 SourceCatalog containing sources detected on this exposure.
276 psf = exposure.getPsf()
277 assert sources.getSchema() == self.
schemaschema
278 self.
deblenddeblend(exposure, sources, psf)
280 def _getPsfFwhm(self, psf, bbox):
283 return psf.computeShape().getDeterminantRadius() * 2.35
291 exposure : `lsst.afw.image.Exposure`
292 Exposure to be processed
293 srcs : `lsst.afw.table.SourceCatalog`
294 SourceCatalog containing sources detected on this exposure
295 psf : `lsst.afw.detection.Psf`
296 Point source function
303 if self.config.useCiLimits:
304 self.log.
info(f
"Using CI catalog limits, "
305 f
"the original number of sources to deblend was {len(srcs)}.")
308 minChildren, maxChildren = self.config.ciDeblendChildRange
309 nPeaks = np.array([len(src.getFootprint().peaks)
for src
in srcs])
310 childrenInRange = np.where((nPeaks >= minChildren) & (nPeaks <= maxChildren))[0]
311 if len(childrenInRange) < self.config.ciNumParentsToDeblend:
312 raise ValueError(
"Fewer than ciNumParentsToDeblend children were contained in the range "
313 "indicated by ciDeblendChildRange. Adjust this range to include more "
317 parents = nPeaks == 1
318 children = np.zeros((len(srcs),), dtype=bool)
319 children[childrenInRange[:self.config.ciNumParentsToDeblend]] =
True
320 srcs = srcs[parents | children]
323 idFactory = srcs.getIdFactory()
324 maxId = np.max(srcs[
"id"])
325 idFactory.notify(maxId)
327 self.log.
info(
"Deblending %d sources" % len(srcs))
332 mi = exposure.getMaskedImage()
334 statsCtrl.setAndMask(mi.getMask().getPlaneBitMask(self.config.maskPlanes))
336 sigma1 = math.sqrt(stats.getValue(afwMath.MEDIAN))
337 self.log.
trace(
'sigma1: %g', sigma1)
341 for i, src
in enumerate(srcs):
344 fp = src.getFootprint()
357 self.log.
warn(
'Parent %i: skipping large footprint (area: %i)',
358 int(src.getId()), int(fp.getArea()))
360 if self.
isMaskedisMasked(fp, exposure.getMaskedImage().getMask()):
363 self.log.
warn(
'Parent %i: skipping masked footprint (area: %i)',
364 int(src.getId()), int(fp.getArea()))
371 self.log.
trace(
'Parent %i: deblending %i peaks', int(src.getId()), len(pks))
377 src.set(self.
tooManyPeaksKeytooManyPeaksKey, len(fp.getPeaks()) > self.config.maxNumberOfPeaks)
381 fp, mi, psf, psf_fwhm, sigma1=sigma1,
382 psfChisqCut1=self.config.psfChisq1,
383 psfChisqCut2=self.config.psfChisq2,
384 psfChisqCut2b=self.config.psfChisq2b,
385 maxNumberOfPeaks=self.config.maxNumberOfPeaks,
386 strayFluxToPointSources=self.config.strayFluxToPointSources,
387 assignStrayFlux=self.config.assignStrayFlux,
388 strayFluxAssignment=self.config.strayFluxRule,
389 rampFluxAtEdge=(self.config.edgeHandling ==
'ramp'),
390 patchEdges=(self.config.edgeHandling ==
'noclip'),
391 tinyFootprintSize=self.config.tinyFootprintSize,
392 clipStrayFluxFraction=self.config.clipStrayFluxFraction,
393 weightTemplates=self.config.weightTemplates,
394 removeDegenerateTemplates=self.config.removeDegenerateTemplates,
395 maxTempDotProd=self.config.maxTempDotProd,
396 medianSmoothTemplate=self.config.medianSmoothTemplate
398 if self.config.catchFailures:
400 except Exception
as e:
401 if self.config.catchFailures:
402 self.log.
warn(
"Unable to deblend source %d: %s" % (src.getId(), e))
405 traceback.print_exc()
412 for j, peak
in enumerate(res.deblendedParents[0].peaks):
413 heavy = peak.getFluxPortion()
414 if heavy
is None or peak.skip:
416 if not self.config.propagateAllPeaks:
421 self.log.
trace(
"Peak at (%i,%i) failed. Using minimal default info for child.",
422 pks[j].getIx(), pks[j].getIy())
426 peakList = foot.getPeaks()
428 peakList.append(peak.peak)
429 zeroMimg = afwImage.MaskedImageF(foot.getBBox())
431 if peak.deblendedAsPsf:
432 if peak.psfFitFlux
is None:
433 peak.psfFitFlux = 0.0
434 if peak.psfFitCenter
is None:
435 peak.psfFitCenter = (peak.peak.getIx(), peak.peak.getIy())
437 assert(len(heavy.getPeaks()) == 1)
440 child = srcs.addNew()
443 child.set(key, src.get(key))
445 child.setParent(src.getId())
446 child.setFootprint(heavy)
447 child.set(self.
psfKeypsfKey, peak.deblendedAsPsf)
448 child.set(self.
hasStrayFluxKeyhasStrayFluxKey, peak.strayFlux
is not None)
449 if peak.deblendedAsPsf:
450 (cx, cy) = peak.psfFitCenter
452 child.set(self.
psfFluxKeypsfFluxKey, peak.psfFitFlux)
462 child.set(self.
peakIdKeypeakIdKey, pks[j].getId())
476 spans = src.getFootprint().spans
478 spans = spans.union(child.getFootprint().spans)
479 src.getFootprint().setSpans(spans)
483 self.
postSingleDeblendHookpostSingleDeblendHook(exposure, srcs, i, npre, kids, fp, psf, psf_fwhm, sigma1, res)
487 self.log.
info(
'Deblended: of %i sources, %i were deblended, creating %i children, total %i sources'
488 % (n0, nparents, n1-n0, n1))
493 def postSingleDeblendHook(self, exposure, srcs, i, npre, kids, fp, psf, psf_fwhm, sigma1, res):
497 """Returns whether a Footprint is large
499 'Large' is defined by thresholds on the area, size and axis ratio.
500 These may be disabled independently by configuring them to be non-positive.
502 This is principally intended to get rid of satellite streaks, which the
503 deblender or other downstream processing can have trouble dealing with
504 (e.g., multiple large HeavyFootprints can chew up memory).
506 if self.config.maxFootprintArea > 0
and footprint.getArea() > self.config.maxFootprintArea:
508 if self.config.maxFootprintSize > 0:
509 bbox = footprint.getBBox()
510 if max(bbox.getWidth(), bbox.getHeight()) > self.config.maxFootprintSize:
512 if self.config.minFootprintAxisRatio > 0:
513 axes = afwEll.Axes(footprint.getShape())
514 if axes.getB() < self.config.minFootprintAxisRatio*axes.getA():
519 """Returns whether the footprint violates the mask limits
521 size = float(footprint.getArea())
522 for maskName, limit
in self.config.maskLimits.items():
523 maskVal = mask.getPlaneBitMask(maskName)
524 unmaskedSpan = footprint.spans.intersectNot(mask, maskVal)
525 if (size - unmaskedSpan.getArea())/size > limit:
530 """Indicate that the parent source is not being deblended
532 We set the appropriate flags and mask.
536 source : `lsst.afw.table.SourceRecord`
537 The source to flag as skipped
538 mask : `lsst.afw.image.Mask`
541 fp = source.getFootprint()
543 if self.config.notDeblendedMask:
544 mask.addMaskPlane(self.config.notDeblendedMask)
545 fp.spans.setMask(mask, mask.getPlaneBitMask(self.config.notDeblendedMask))
549 centerX = int(bbox.getMinX()+bbox.getWidth()/2)
550 centerY = int(bbox.getMinY()+bbox.getHeight()/2)
556 source.set(self.
nPeaksKeynPeaksKey, len(fp.peaks))
Pass parameters to a Statistics object.
A mapping between the keys of two Schemas, used to copy data between them.
static Log getLogger(Log const &logger)
def isLargeFootprint(self, footprint)
def deblend(self, exposure, srcs, psf)
def addSchemaKeys(self, schema)
def postSingleDeblendHook(self, exposure, srcs, i, npre, kids, fp, psf, psf_fwhm, sigma1, res)
def isMasked(self, footprint, mask)
def skipParent(self, source, mask)
def preSingleDeblendHook(self, exposure, srcs, i, fp, psf, psf_fwhm, sigma1)
def __init__(self, schema, peakSchema=None, **kwargs)
def _getPsfFwhm(self, psf, bbox)
deblendPatchedTemplateKey
def run(self, exposure, sources)
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)