22 __all__ = [
'SourceDeblendConfig',
'SourceDeblendTask']
35 from 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.
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
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.
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.
274 psf = exposure.getPsf()
275 assert sources.getSchema() == self.
schemaschema
276 self.
deblenddeblend(exposure, sources, psf)
278 def _getPsfFwhm(self, psf, bbox):
281 return psf.computeShape().getDeterminantRadius() * 2.35
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
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)}.")
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 "
315 parents = nPeaks == 1
316 children = np.zeros((len(srcs),), dtype=bool)
317 children[childrenInRange[:self.config.ciNumParentsToDeblend]] =
True
318 srcs = srcs[parents | children]
321 idFactory = srcs.getIdFactory()
322 maxId = np.max(srcs[
"id"])
323 idFactory.notify(maxId)
325 self.log.
info(
"Deblending %d sources", len(srcs))
330 mi = exposure.getMaskedImage()
332 statsCtrl.setAndMask(mi.getMask().getPlaneBitMask(self.config.maskPlanes))
334 sigma1 = math.sqrt(stats.getValue(afwMath.MEDIAN))
335 self.log.
trace(
'sigma1: %g', sigma1)
339 for i, src
in enumerate(srcs):
342 fp = src.getFootprint()
355 self.log.
warning(
'Parent %i: skipping large footprint (area: %i)',
356 int(src.getId()), int(fp.getArea()))
358 if self.
isMaskedisMasked(fp, exposure.getMaskedImage().getMask()):
361 self.log.
warning(
'Parent %i: skipping masked footprint (area: %i)',
362 int(src.getId()), int(fp.getArea()))
369 self.log.
trace(
'Parent %i: deblending %i peaks', int(src.getId()), len(pks))
375 src.set(self.
tooManyPeaksKeytooManyPeaksKey, len(fp.getPeaks()) > self.config.maxNumberOfPeaks)
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
396 if self.config.catchFailures:
398 except Exception
as e:
399 if self.config.catchFailures:
400 self.log.
warning(
"Unable to deblend source %d: %s", src.getId(), e)
403 traceback.print_exc()
410 for j, peak
in enumerate(res.deblendedParents[0].peaks):
411 heavy = peak.getFluxPortion()
412 if heavy
is None or peak.skip:
414 if not self.config.propagateAllPeaks:
419 self.log.
trace(
"Peak at (%i,%i) failed. Using minimal default info for child.",
420 pks[j].getIx(), pks[j].getIy())
424 peakList = foot.getPeaks()
426 peakList.append(peak.peak)
427 zeroMimg = afwImage.MaskedImageF(foot.getBBox())
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())
435 assert(len(heavy.getPeaks()) == 1)
438 child = srcs.addNew()
441 child.set(key, src.get(key))
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
450 child.set(self.
psfFluxKeypsfFluxKey, peak.psfFitFlux)
460 child.set(self.
peakIdKeypeakIdKey, pks[j].getId())
474 spans = src.getFootprint().spans
476 spans = spans.union(child.getFootprint().spans)
477 src.getFootprint().setSpans(spans)
481 self.
postSingleDeblendHookpostSingleDeblendHook(exposure, srcs, i, npre, kids, fp, psf, psf_fwhm, sigma1, res)
485 self.log.
info(
'Deblended: of %i sources, %i were deblended, creating %i children, total %i sources',
486 n0, nparents, n1-n0, n1)
491 def postSingleDeblendHook(self, exposure, srcs, i, npre, kids, fp, psf, psf_fwhm, sigma1, res):
495 """Returns whether a Footprint is large
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.
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).
504 if self.config.maxFootprintArea > 0
and footprint.getArea() > self.config.maxFootprintArea:
506 if self.config.maxFootprintSize > 0:
507 bbox = footprint.getBBox()
508 if max(bbox.getWidth(), bbox.getHeight()) > self.config.maxFootprintSize:
510 if self.config.minFootprintAxisRatio > 0:
511 axes = afwEll.Axes(footprint.getShape())
512 if axes.getB() < self.config.minFootprintAxisRatio*axes.getA():
517 """Returns whether the footprint violates the mask limits
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)
523 if (size - unmaskedSpan.getArea())/size > limit:
528 """Indicate that the parent source is not being deblended
530 We set the appropriate flags and mask.
534 source : `lsst.afw.table.SourceRecord`
535 The source to flag as skipped
536 mask : `lsst.afw.image.Mask`
539 fp = source.getFootprint()
541 if self.config.notDeblendedMask:
542 mask.addMaskPlane(self.config.notDeblendedMask)
543 fp.spans.setMask(mask, mask.getPlaneBitMask(self.config.notDeblendedMask))
547 centerX = int(bbox.getMinX()+bbox.getWidth()/2)
548 centerY = int(bbox.getMinY()+bbox.getHeight()/2)
554 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.
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)