22__all__ = [
'SourceDeblendConfig',
'SourceDeblendTask']
35from lsst.utils.timer
import timeMethod
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',
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.',
50 strayFluxToPointSources = pexConfig.ChoiceField(
51 doc=
'When the deblender should attribute stray flux to point sources',
52 dtype=str, default=
'necessary',
54 'necessary':
'When there is not an extended object in the footprint',
56 'never': (
'Never; stray flux will not be attributed to any deblended child '
57 'if the deblender thinks all peaks look like point sources'),
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.')
65 strayFluxRule = pexConfig.ChoiceField(
66 doc=
'How to split flux among peaks',
67 dtype=str, default=
'trim',
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')
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")
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 '
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(
122 doc=(
"Mask planes with the corresponding limit on the fraction of masked pixels. "
123 "Sources violating this limit will not be deblended."),
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 "
137 medianSmoothTemplate = pexConfig.Field(dtype=bool, default=
True,
138 doc=
"Apply a smoothing filter to all of the template images")
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.")
162 """Split blended sources into individual sources.
164 This task has no return value; it only modifies the SourceCatalog
in-place.
166 ConfigClass = SourceDeblendConfig
167 _DefaultName = "sourceDeblend"
169 def __init__(self, schema, peakSchema=None, **kwargs):
170 """Create the task, adding necessary fields to the given 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
182 Additional keyword arguments passed to ~lsst.pipe.base.task
184 pipeBase.Task.__init__(self, **kwargs)
187 if item.field.getName().startswith(
"merge_footprint")]
188 peakMinimalSchema = afwDet.PeakTable.makeMinimalSchema()
189 if peakSchema
is None:
195 for item
in peakSchema:
196 if item.key
not in peakMinimalSchema:
202 schema.addField(item.field)
203 assert schema == self.
peakSchemaMapperpeakSchemaMapper.getOutputSchema(),
"Logic bug mapping schemas"
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')
223 if self.config.catchFailures:
225 doc=
"Deblending failed on source")
228 doc=
"Deblender skipped this source")
231 'deblend_rampedTemplate', type=
'Flag',
232 doc=(
'This source was near an image edge and the deblender used '
233 '"ramp" edge-handling.'))
236 'deblend_patchedTemplate', type=
'Flag',
237 doc=(
'This source was near an image edge and the deblender used '
238 '"patched" edge-handling.'))
241 'deblend_hasStrayFlux', type=
'Flag',
242 doc=(
'This source was assigned some stray flux'))
244 self.log.
trace(
'Added keys to schema: %s',
", ".join(
str(x)
for x
in (
247 self.
peakCenterpeakCenter = afwTable.Point2IKey.addFields(schema, name=
"deblend_peak_center",
248 doc=
"Center used to apply constraints in scarlet",
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")
264 def run(self, exposure, sources):
265 """Get the PSF from the provided exposure and then run deblend.
270 Exposure to be processed
272 SourceCatalog containing sources detected on this exposure.
274 psf = exposure.getPsf()
275 assert sources.getSchema() == self.
schemaschema
276 self.
deblenddeblend(exposure, sources, psf)
278 def _getPsfFwhm(self, psf, position):
279 return psf.computeShape(position).getDeterminantRadius() * 2.35
288 Exposure to be processed
290 SourceCatalog containing sources detected on this exposure
292 Point source function
299 if self.config.useCiLimits:
300 self.log.
info(f
"Using CI catalog limits, "
301 f
"the original number of sources to deblend was {len(srcs)}.")
304 minChildren, maxChildren = self.config.ciDeblendChildRange
305 nPeaks = np.array([len(src.getFootprint().peaks)
for src
in srcs])
306 childrenInRange = np.where((nPeaks >= minChildren) & (nPeaks <= maxChildren))[0]
307 if len(childrenInRange) < self.config.ciNumParentsToDeblend:
308 raise ValueError(
"Fewer than ciNumParentsToDeblend children were contained in the range "
309 "indicated by ciDeblendChildRange. Adjust this range to include more "
313 parents = nPeaks == 1
314 children = np.zeros((len(srcs),), dtype=bool)
315 children[childrenInRange[:self.config.ciNumParentsToDeblend]] =
True
316 srcs = srcs[parents | children]
319 idFactory = srcs.getIdFactory()
320 maxId = np.max(srcs[
"id"])
321 idFactory.notify(maxId)
323 self.log.
info(
"Deblending %d sources", len(srcs))
328 mi = exposure.getMaskedImage()
330 statsCtrl.setAndMask(mi.getMask().getPlaneBitMask(self.config.maskPlanes))
332 sigma1 = math.sqrt(stats.getValue(afwMath.MEDIAN))
333 self.log.
trace(
'sigma1: %g', sigma1)
337 for i, src
in enumerate(srcs):
340 fp = src.getFootprint()
353 self.log.
warning(
'Parent %i: skipping large footprint (area: %i)',
354 int(src.getId()),
int(fp.getArea()))
356 if self.
isMaskedisMasked(fp, exposure.getMaskedImage().getMask()):
359 self.log.
warning(
'Parent %i: skipping masked footprint (area: %i)',
360 int(src.getId()),
int(fp.getArea()))
364 center = fp.getCentroid()
365 psf_fwhm = self.
_getPsfFwhm_getPsfFwhm(psf, center)
367 self.log.
trace(
'Parent %i: deblending %i peaks',
int(src.getId()), len(pks))
373 src.set(self.
tooManyPeaksKeytooManyPeaksKey, len(fp.getPeaks()) > self.config.maxNumberOfPeaks)
377 fp, mi, psf, psf_fwhm, sigma1=sigma1,
378 psfChisqCut1=self.config.psfChisq1,
379 psfChisqCut2=self.config.psfChisq2,
380 psfChisqCut2b=self.config.psfChisq2b,
381 maxNumberOfPeaks=self.config.maxNumberOfPeaks,
382 strayFluxToPointSources=self.config.strayFluxToPointSources,
383 assignStrayFlux=self.config.assignStrayFlux,
384 strayFluxAssignment=self.config.strayFluxRule,
385 rampFluxAtEdge=(self.config.edgeHandling ==
'ramp'),
386 patchEdges=(self.config.edgeHandling ==
'noclip'),
387 tinyFootprintSize=self.config.tinyFootprintSize,
388 clipStrayFluxFraction=self.config.clipStrayFluxFraction,
389 weightTemplates=self.config.weightTemplates,
390 removeDegenerateTemplates=self.config.removeDegenerateTemplates,
391 maxTempDotProd=self.config.maxTempDotProd,
392 medianSmoothTemplate=self.config.medianSmoothTemplate
394 if self.config.catchFailures:
396 except Exception
as e:
397 if self.config.catchFailures:
398 self.log.
warning(
"Unable to deblend source %d: %s", src.getId(), e)
401 traceback.print_exc()
408 for j, peak
in enumerate(res.deblendedParents[0].peaks):
409 heavy = peak.getFluxPortion()
410 if heavy
is None or peak.skip:
412 if not self.config.propagateAllPeaks:
417 self.log.
trace(
"Peak at (%i,%i) failed. Using minimal default info for child.",
418 pks[j].getIx(), pks[j].getIy())
422 peakList = foot.getPeaks()
424 peakList.append(peak.peak)
425 zeroMimg = afwImage.MaskedImageF(foot.getBBox())
427 if peak.deblendedAsPsf:
428 if peak.psfFitFlux
is None:
429 peak.psfFitFlux = 0.0
430 if peak.psfFitCenter
is None:
431 peak.psfFitCenter = (peak.peak.getIx(), peak.peak.getIy())
433 assert(len(heavy.getPeaks()) == 1)
436 child = srcs.addNew()
439 child.set(key, src.get(key))
441 child.setParent(src.getId())
442 child.setFootprint(heavy)
443 child.set(self.
psfKeypsfKey, peak.deblendedAsPsf)
444 child.set(self.
hasStrayFluxKeyhasStrayFluxKey, peak.strayFlux
is not None)
445 if peak.deblendedAsPsf:
446 (cx, cy) = peak.psfFitCenter
448 child.set(self.
psfFluxKeypsfFluxKey, peak.psfFitFlux)
458 child.set(self.
peakIdKeypeakIdKey, pks[j].getId())
472 spans = src.getFootprint().spans
474 spans = spans.union(child.getFootprint().spans)
475 src.getFootprint().setSpans(spans)
479 self.
postSingleDeblendHookpostSingleDeblendHook(exposure, srcs, i, npre, kids, fp, psf, psf_fwhm, sigma1, res)
483 self.log.
info(
'Deblended: of %i sources, %i were deblended, creating %i children, total %i sources',
484 n0, nparents, n1-n0, n1)
489 def postSingleDeblendHook(self, exposure, srcs, i, npre, kids, fp, psf, psf_fwhm, sigma1, res):
493 """Returns whether a Footprint is large
495 'Large' is defined by thresholds on the area, size
and axis ratio.
496 These may be disabled independently by configuring them to be non-positive.
498 This
is principally intended to get rid of satellite streaks, which the
499 deblender
or other downstream processing can have trouble dealing
with
500 (e.g., multiple large HeavyFootprints can chew up memory).
502 if self.config.maxFootprintArea > 0
and footprint.getArea() > self.config.maxFootprintArea:
504 if self.config.maxFootprintSize > 0:
505 bbox = footprint.getBBox()
506 if max(bbox.getWidth(), bbox.getHeight()) > self.config.maxFootprintSize:
508 if self.config.minFootprintAxisRatio > 0:
509 axes = afwEll.Axes(footprint.getShape())
510 if axes.getB() < self.config.minFootprintAxisRatio*axes.getA():
515 """Returns whether the footprint violates the mask limits
517 size = float(footprint.getArea())
518 for maskName, limit
in self.config.maskLimits.items():
519 maskVal = mask.getPlaneBitMask(maskName)
520 unmaskedSpan = footprint.spans.intersectNot(mask, maskVal)
521 if (size - unmaskedSpan.getArea())/size > limit:
526 """Indicate that the parent source is not being deblended
528 We set the appropriate flags and mask.
533 The source to flag
as skipped
537 fp = source.getFootprint()
539 if self.config.notDeblendedMask:
540 mask.addMaskPlane(self.config.notDeblendedMask)
541 fp.spans.setMask(mask, mask.getPlaneBitMask(self.config.notDeblendedMask))
545 centerX =
int(bbox.getMinX()+bbox.getWidth()/2)
546 centerY =
int(bbox.getMinY()+bbox.getHeight()/2)
552 source.set(self.
nPeaksKeynPeaksKey, len(fp.peaks))
A polymorphic base class for representing an image's Point Spread Function.
A class to contain the data, WCS, and other information needed to describe an image of the sky.
Represent a 2-dimensional array of bitmask pixels.
Pass parameters to a Statistics object.
Defines the fields and offsets for a table.
A mapping between the keys of two Schemas, used to copy data between them.
Record class that contains measurements made on a single exposure.
def isLargeFootprint(self, footprint)
def deblend(self, exposure, srcs, psf)
def _getPsfFwhm(self, psf, position)
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)
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)