Coverage for python/lsst/meas/extensions/scarlet/scarletDeblendTask.py: 16%
521 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-10 03:23 -0700
1# This file is part of meas_extensions_scarlet.
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/>.
22from __future__ import annotations
24import logging
25from dataclasses import dataclass
26from functools import partial
27import time
28import traceback
29from typing import cast
31import lsst.afw.detection as afwDet
32import lsst.afw.geom.ellipses as afwEll
33import lsst.afw.image as afwImage
34import lsst.afw.table as afwTable
35import lsst.pex.config as pexConfig
36import lsst.pipe.base as pipeBase
37import lsst.scarlet.lite as scl
38import numpy as np
39from lsst.scarlet.lite.detect_pybind11 import get_footprints
40from scipy import ndimage
41from lsst.utils.logging import PeriodicLogger
42from lsst.utils.timer import timeMethod
44from . import io, utils
45from .footprint import scarletFootprintToAfw
46from .source import IsolatedSource
48__all__ = ["deblend", "ScarletDeblendContext", "ScarletDeblendConfig", "ScarletDeblendTask"]
50logger = logging.getLogger(__name__)
53class DeblenderError(Exception):
54 """Exception raised when the deblender fails.
56 This is used to catch errors in the deblender and set the appropriate flags
57 on the parent source.
58 """
60 def __init__(self, message: str, parentId: int, errorName: str):
61 super().__init__(message)
62 self.message = message
63 self.parentId = parentId
64 self.errorName = errorName
66 def __str__(self) -> str:
67 return f"DeblenderError: {self.args[0]} (parentId: {self.parentId})"
70class DeblenderSkippedError(Exception):
71 """Exception raised when a blend is skipped.
73 This is used to catch cases where the deblender does not process
74 a deconvolved parent because it is skipped for some reason.
75 """
76 def __init__(self, message: str, parentId: int, skipKey):
77 super().__init__(message)
78 self.message = message
79 self.parentId = parentId
80 self.skipKey = skipKey
82 def __str__(self) -> str:
83 return f"DeblenderSkippedError: {self.args[0]} (parentId: {self.parentId}, skipKey: {self.skipKey})"
86def _checkBlendConvergence(blend: scl.Blend, f_rel: float) -> bool:
87 """Check whether or not a blend has converged"""
88 if len(blend.loss) < 2:
89 # No fitting ran (or only one iteration): there is no delta
90 # to test against, so report vacuous convergence. The
91 # ``convergenceFailed`` flag is reserved for a fitted blend
92 # that did not converge.
93 return True
94 deltaLoss = np.abs(blend.loss[-2] - blend.loss[-1])
95 convergence = f_rel * np.abs(blend.loss[-1])
96 return deltaLoss < convergence
99def isPseudoSource(source: afwTable.source.SourceRecord, pseudoColumns: list[str]) -> bool:
100 """Check if a source is a pseudo source.
102 This is mostly for skipping sky objects,
103 but any other column can also be added to disable
104 deblending on a parent or individual source when
105 set to `True`.
107 Parameters
108 ----------
109 source :
110 The source to check for the pseudo bit.
111 pseudoColumns :
112 A list of columns to check for pseudo sources.
113 """
114 isPseudo = False
115 for col in pseudoColumns:
116 try:
117 isPseudo |= source[col]
118 except KeyError:
119 pass
120 return isPseudo
123def _getDeconvolvedFootprints(
124 mDeconvolved: afwImage.MultibandExposure,
125 sources: afwTable.SourceCatalog,
126 config: ScarletDeblendConfig,
127 sigma: np.ndarray,
128) -> tuple[list[scl.detect.Footprint], scl.Image]:
129 """Detect footprints in the deconvolved image
131 Parameters
132 ----------
133 mDeconvolved :
134 The deconvolved multiband exposure to detect footprints in.
135 sources :
136 The source catalog for the entire coadd.
137 config :
138 The configuration for the deblender.
139 sigma :
140 Per-band noise scale used to normalize the deconvolved image
141 before stacking. Must be sourced from the input coadd's
142 variance plane; the deconvolved exposure deliberately carries
143 an invalidated (zero-filled) variance.
145 Returns
146 -------
147 footprints :
148 The detected footprints in the deconvolved image.
149 footprintImage :
150 A footprint image as returned by `scarlet.detect.footprints_to_image`.
151 """
152 bbox = mDeconvolved.getBBox()
153 xmin, ymin = bbox.getMin()
154 detect = np.nansum(mDeconvolved.image.array/sigma[:, None, None], axis=0)
156 # We don't use the variance here because testing in DM-47738
157 # has shown that we get better results without it
158 # (we're detecting footprints without peaks in the deconvolved,
159 # noise reduced image).
160 footprints = scl.detect.detect_footprints(
161 images=np.array([detect]),
162 variance=np.ones((1, detect.shape[0], detect.shape[1]), dtype=detect.dtype),
163 scales=1,
164 min_area=config.minDeconvolvedArea,
165 footprint_thresh=config.footprintSNRThresh,
166 find_peaks=False,
167 origin=(ymin, xmin)
168 )
170 # Create an indexed image of the footprints so that the value of a pixel
171 # gives the index + 1 of the footprints that contain that pixel.
172 scarletBox = utils.bboxToScarletBox(bbox)
173 footprintImage = scl.detect.footprints_to_image(footprints, scarletBox)
175 # Ensure that there is a minimal footprint for all peaks
176 # in the detection catalog.
177 footprintArray = footprintImage.data > 0
178 for source in sources:
179 for peak in source.getFootprint().peaks:
180 x, y = peak.getI()
181 footprintArray[y - ymin, x - xmin] = True
183 # Grow the footprints
184 psfMinimalSize = config.growSize * 2 + 1
185 detectionArray = ndimage.binary_dilation(
186 footprintArray,
187 scl.utils.get_circle_mask(psfMinimalSize, bool)
188 )
190 # Ensure that all of the deconvolved footprints are contained within
191 # a single footprint from the input catalog.
192 # This should be unecessary, however creating entries in the output
193 # catalog will produce unexpected results if a deconvolved footprint
194 # is in more than one footprint from the source catalog or has
195 # flux outside of its parent footprint.
196 sourceImage = afwDet.footprintsToNumpy(sources, shape=detect.shape, xy0=(xmin, ymin))
197 detectionArray = detectionArray * sourceImage
199 footprints = get_footprints(
200 image=detectionArray.astype(np.float32),
201 min_separation=0,
202 min_area=1,
203 peak_thresh=0,
204 footprint_thresh=0,
205 find_peaks=False,
206 y0=ymin,
207 x0=xmin,
208 )
210 footprintImage = scl.detect.footprints_to_image(footprints, scarletBox)
212 return footprints, footprintImage
215@dataclass(kw_only=True)
216class ScarletDeblendContext:
217 """Context with parameters and config options for deblending
219 Attributes
220 ----------
221 monotonicity :
222 The monotonicity operator.
223 observation :
224 The observation for the entire coadd.
225 deconvolved :
226 The deconvolved image.
227 residual :
228 The residual image
229 (observation - deconvolved mode convolved with difference kernel).
230 footprints :
231 The footprints in the deconvolved image.
232 footprintImage :
233 An indexed image of the scarlet footprints so that the value
234 of a pixel gives the index + 1 of the footprints that
235 contain that pixel.
236 config :
237 The configuration for the deblender.
238 """
239 monotonicity: scl.operators.Monotonicity
240 observation: scl.Observation
241 deconvolved: scl.Image
242 footprints: list[scl.Footprint]
243 footprintImage: scl.Image
244 config: ScarletDeblendConfig
246 @staticmethod
247 def build(
248 mExposure: afwImage.MultibandExposure,
249 mDeconvolved: afwImage.MultibandExposure,
250 catalog: afwTable.SourceCatalog,
251 config: ScarletDeblendConfig,
252 ) -> ScarletDeblendContext:
253 """Build the context from a minimal set of inputs
255 Parameters
256 ----------
257 mExposure :
258 The multiband exposure for the entire coadd.
259 mDeconvolved :
260 The deconvolved multiband exposure for the entire coadd.
261 catalog :
262 The parent catalog for the entire coadd.
263 config :
264 The configuration for the deblender.
265 """
266 # The PSF of the model in the deconvolved space
267 modelPsf = scl.utils.integrated_circular_gaussian(
268 sigma=config.modelPsfSigma
269 )
270 # Initialize the monotonicity operator with a size of 101 x 101 pixels.
271 # Note: If a component is > 101x101 in either axis then the
272 # monotonicity operator will resize itself.
273 monotonicity = scl.operators.Monotonicity((101, 101))
274 # Build the observation for the entire coadd
275 observation = utils.buildObservation(
276 modelPsf=modelPsf,
277 psfCenter=mExposure.getBBox().getCenter(),
278 mExposure=mExposure,
279 badPixelMasks=config.badMask,
280 useWeights=config.useWeights,
281 convolutionType=config.convolutionType,
282 catalog=catalog,
283 )
285 # Create the deconvolved image
286 yx0 = observation.images.yx0
287 bands = observation.images.bands
288 deconvolved = scl.Image(mDeconvolved.image.array, bands=mDeconvolved.bands, yx0=yx0)
289 if len(bands) == 1:
290 deconvolved = deconvolved[bands[0]:]
291 else:
292 deconvolved = deconvolved[bands]
294 # Detect footprints in the deconvolved image. The deconvolved
295 # exposure carries an invalidated (zero-filled) variance plane,
296 # so the per-band noise scale must come from the input coadd.
297 sigma = np.nanmedian(np.sqrt(mExposure.variance.array), axis=(1, 2))
299 footprints, footprintImage = _getDeconvolvedFootprints(
300 mDeconvolved=mDeconvolved,
301 sources=catalog,
302 config=config,
303 sigma=sigma,
304 )
306 return ScarletDeblendContext(
307 monotonicity=monotonicity,
308 observation=observation,
309 deconvolved=deconvolved,
310 footprints=footprints,
311 footprintImage=footprintImage,
312 config=config,
313 )
316def deblend(
317 context: ScarletDeblendContext,
318 footprint: afwDet.Footprint,
319 config: ScarletDeblendConfig,
320 spectrumInit: bool = True,
321) -> scl.Blend:
322 """Deblend a parent footprint
324 Parameters
325 ----------
326 context :
327 Context with parameters and config options for deblending
328 footprint :
329 The parent footprint to deblend
330 config :
331 The configuration for the deblender
332 spectrumInit :
333 Whether or not to initialize the sources with their best-fit spectra
335 Returns
336 -------
337 blend : `scarlet.lite.Blend`
338 The blend this is to be deblended
339 skippedSources : `list[int]`
340 Indices of sources that were skipped due to no flux.
341 This usually means that a source was a spurrious detection in one
342 band that should not have been included in the merged catalog.
343 skippedBands : `list[str]`
344 Bands that were skipped because a PSF could not be generated for them.
345 """
346 # Define the bounding boxes in lsst.geom.Box2I and lsst.scarlet.lite.Box
347 footBox = footprint.getBBox()
348 bbox = utils.bboxToScarletBox(footBox)
350 # Extract the observation that covers the footprint and make
351 # a copy so that the changes don't affect the original observation.
352 observation = context.observation[:, bbox].copy()
353 footprintData = footprint.spans.asArray()
355 # Mask the pixels outside of the footprint
356 observation.weights.data[:] *= footprintData
358 # Drop pseudo peaks (e.g. sky objects) so that sources are
359 # initialized only for real detections. Keep the surviving peak
360 # records so the source-to-peak back-pointer below is indexed in
361 # the same filtered list that drives initialization.
362 non_pseudo_peaks = [
363 peak
364 for peak in footprint.peaks
365 if not isPseudoSource(peak, config.pseudoColumns)
366 ]
367 # Convert the peaks into an array
368 peaks = [
369 np.array([peak.getIy(), peak.getIx()], dtype=int)
370 for peak in non_pseudo_peaks
371 ]
373 detect_image = np.sum(context.deconvolved[:, bbox].data, axis=0)
375 # Initialize the sources
376 sources = scl.initialization.FactorizedInitialization(
377 observation=observation,
378 centers=peaks,
379 detect=detect_image,
380 min_snr=config.minSNR,
381 monotonicity=context.monotonicity,
382 bg_thresh=config.backgroundThresh,
383 initial_bg_thresh=config.initialBackgroundThresh,
384 ).sources
386 blend = scl.Blend(sources, observation)
388 # Initialize each source with its best fit spectrum
389 if spectrumInit:
390 try:
391 blend.fit_spectra()
392 except Exception as e:
393 # If the spectrum initialization fails, we will just skip it
394 # and use the default spectrum.
395 logger.warning(
396 "Spectrum initialization failed with error: %s", e, exc_info=True
397 )
399 # Set the optimizer
400 if config.optimizer == "adaprox":
401 blend.parameterize(
402 partial(
403 scl.component.default_adaprox_parameterization,
404 noise_rms=observation.noise_rms / 10,
405 )
406 )
407 elif config.optimizer == "fista":
408 blend.parameterize(scl.component.default_fista_parameterization)
409 else:
410 raise ValueError("Unrecognized optimizer. Must be either 'adaprox' or 'fista'.")
412 if config.maxIter > 0:
413 blend.fit(
414 max_iter=config.maxIter,
415 e_rel=config.relativeError,
416 min_iter=config.minIter,
417 resize=config.resizeFrequency,
418 )
420 # Attach the peak to all of the initialized sources
421 for k, center in enumerate(peaks):
422 # This is just to make sure that there isn't a coding bug
423 if len(sources[k].components) > 0 and np.any(sources[k].center != center):
424 raise ValueError(
425 f"Misaligned center, expected {center} but got {sources[k].center}"
426 )
427 # Store the record for the peak with the appropriate source
428 sources[k].detectedPeak = non_pseudo_peaks[k]
430 return blend
433class ScarletDeblendConfig(pexConfig.Config):
434 """MultibandDeblendConfig
436 Configuration for the multiband deblender.
437 The parameters are organized by the parameter types, which are
438 - Stopping Criteria: Used to determine if the fit has converged
439 - Position Fitting Criteria: Used to fit the positions of the peaks
440 - Constraints: Used to apply constraints to the peaks and their components
441 - Other: Parameters that don't fit into the above categories
442 """
444 # Stopping Criteria
445 minIter = pexConfig.Field[int](
446 default=5,
447 doc="Minimum number of iterations before the optimizer is allowed to stop.",
448 )
449 maxIter = pexConfig.Field[int](
450 default=300,
451 doc=("Maximum number of iterations to deblend a single parent"),
452 )
453 relativeError = pexConfig.Field[float](
454 default=1e-3,
455 doc=(
456 "Change in the loss function between iterations to exit fitter. "
457 "Typically this is `1e-3` if measurements will be made on the "
458 "flux re-distributed models and `1e-4` when making measurements "
459 "on the models themselves."
460 ),
461 )
462 resizeFrequency = pexConfig.Field[int](
463 default=3,
464 doc="Number of iterations between resizing sources.",
465 )
467 # Constraints
468 morphThresh = pexConfig.Field[float](
469 default=1,
470 doc="Fraction of background RMS a pixel must have"
471 "to be included in the initial morphology",
472 )
473 # Lite Parameters
474 # All of these parameters (except version) are only valid if version='lite'
475 optimizer = pexConfig.ChoiceField[str](
476 default="adaprox",
477 allowed={
478 "adaprox": "Proximal ADAM optimization",
479 "fista": "Accelerated proximal gradient method",
480 },
481 doc="The optimizer to use for fitting parameters.",
482 )
483 backgroundThresh = pexConfig.Field[float](
484 default=1.0,
485 doc="Fraction of background to use for a sparsity threshold. "
486 "This prevents sources from growing unrealistically outside "
487 "the parent footprint while still modeling flux correctly "
488 "for bright sources.",
489 )
490 initialBackgroundThresh = pexConfig.Field[float](
491 default=1.0,
492 doc="Same as `backgroundThresh` but used only for source initialization.",
493 )
494 maxProxIter = pexConfig.Field[int](
495 default=1,
496 doc="Maximum number of proximal operator iterations inside of each "
497 "iteration of the optimizer. "
498 "This config field is only used if version='lite' and optimizer='adaprox'.",
499 )
500 # Other scarlet paremeters
501 useWeights = pexConfig.Field[bool](
502 default=True,
503 doc=(
504 "Whether or not use use inverse variance weighting."
505 "If `useWeights` is `False` then flat weights are used"
506 ),
507 )
508 modelPsfSize = pexConfig.Field[int](
509 default=11, doc="Model PSF side length in pixels"
510 )
511 modelPsfSigma = pexConfig.Field[float](
512 default=0.8, doc="Define sigma for the model frame PSF"
513 )
514 minSNR = pexConfig.Field[float](
515 default=50,
516 doc="Minimum Signal to noise to accept the source."
517 "Sources with lower flux will be initialized with the PSF but updated "
518 "like an ordinary ExtendedSource (known in scarlet as a `CompactSource`).",
519 )
520 saveTemplates = pexConfig.Field[bool](
521 default=True, doc="Whether or not to save the SEDs and templates"
522 )
523 processSingles = pexConfig.Field[bool](
524 default=False,
525 doc="Whether or not to process isolated sources in the deblender",
526 )
527 persistIsolated = pexConfig.Field[bool](
528 default=True,
529 doc="Whether or not to persist isolated sources in the scarlet models",
530 )
531 convolutionType = pexConfig.Field[str](
532 default="fft",
533 doc="Type of convolution to render the model to the observations.\n"
534 "- 'fft': perform convolutions in Fourier space\n"
535 "- 'real': peform convolutions in real space.",
536 )
537 setSpectra = pexConfig.Field[bool](
538 default=True,
539 doc="Whether or not to solve for the best-fit spectra during initialization. "
540 "This makes initialization slightly longer, as it requires a convolution "
541 "to set the optimal spectra, but results in a much better initial log-likelihood "
542 "and reduced total runtime, with convergence in fewer iterations."
543 "This option is only used when "
544 "peaks*area < `maxSpectrumCutoff` will use the improved initialization.",
545 )
546 footprintSNRThresh = pexConfig.Field[float](
547 default=5.0,
548 doc="Minimum SNR for a pixel to be detected in a footprint.",
549 )
550 growSize = pexConfig.Field[int](
551 default=2,
552 doc="Number of pixels to grow the deconvolved footprints before final detection.",
553 )
555 # Mask-plane restrictions
556 badMask = pexConfig.ListField[str](
557 default=utils.defaultBadPixelMasks,
558 doc="Whether or not to process isolated sources in the deblender",
559 )
560 statsMask = pexConfig.ListField[str](
561 default=["SAT", "INTRP", "NO_DATA"],
562 doc="Mask planes to ignore when performing statistics",
563 )
564 maskLimits = pexConfig.DictField(
565 keytype=str,
566 itemtype=float,
567 default={},
568 doc=(
569 "Mask planes with the corresponding limit on the fraction of masked pixels. "
570 "Sources violating this limit will not be deblended. "
571 "If the fraction is `0` then the limit is a single pixel."
572 ),
573 )
574 minDeconvolvedArea = pexConfig.Field[int](
575 default=9,
576 doc="Minimum area for a single footprint in the deconvolved image. "
577 "Detected footprints smaller than this will not be created.",
578 )
580 # Size restrictions
581 maxNumberOfPeaks = pexConfig.Field[int](
582 default=600,
583 doc=(
584 "Only deblend the brightest maxNumberOfPeaks peaks in the parent"
585 " (<= 0: unlimited)"
586 ),
587 )
588 maxFootprintArea = pexConfig.Field[int](
589 default=2_000_000,
590 doc=(
591 "Maximum area for footprints before they are ignored as large; "
592 "non-positive means no threshold applied"
593 ),
594 )
595 maxAreaTimesPeaks = pexConfig.Field[int](
596 default=1_000_000_000,
597 doc=(
598 "Maximum rectangular footprint area * nPeaks in the footprint. "
599 "This was introduced in DM-33690 to prevent fields that are crowded or have a "
600 "LSB galaxy that causes memory intensive initialization in scarlet from dominating "
601 "the overall runtime and/or causing the task to run out of memory. "
602 "(<= 0: unlimited)"
603 ),
604 )
605 maxFootprintSize = pexConfig.Field[int](
606 default=0,
607 doc=(
608 "Maximum linear dimension for footprints before they are ignored "
609 "as large; non-positive means no threshold applied"
610 ),
611 )
612 minFootprintAxisRatio = pexConfig.Field[float](
613 default=0.0,
614 doc=(
615 "Minimum axis ratio for footprints before they are ignored "
616 "as large; non-positive means no threshold applied"
617 ),
618 )
619 maxSpectrumCutoff = pexConfig.Field[int](
620 default=1_000_000,
621 doc=(
622 "Maximum number of pixels * number of sources in a blend. "
623 "This is different than `maxFootprintArea` because this isn't "
624 "the footprint area but the area of the bounding box that "
625 "contains the footprint, and is also multiplied by the number of"
626 "sources in the footprint. This prevents large skinny blends with "
627 "a high density of sources from running out of memory. "
628 "If `maxSpectrumCutoff == -1` then there is no cutoff."
629 ),
630 )
631 # Failure modes
632 fallback = pexConfig.Field[bool](
633 default=True,
634 doc="Whether or not to fallback to a smaller number of components if a source does not initialize",
635 )
636 notDeblendedMask = pexConfig.Field[str](
637 default="NOT_DEBLENDED",
638 optional=True,
639 doc="Mask name for footprints not deblended, or None",
640 )
641 catchFailures = pexConfig.Field[bool](
642 default=True,
643 doc=(
644 "If True, catch exceptions thrown by the deblender, log them, "
645 "and set a flag on the parent, instead of letting them propagate up"
646 ),
647 )
649 # Other options
650 columnInheritance = pexConfig.DictField(
651 keytype=str,
652 itemtype=str,
653 default={
654 "deblend_nChild": "deblend_parentNChild",
655 "deblend_nPeaks": "deblend_parentNPeaks",
656 },
657 doc="Columns to pass from the parent to the child. "
658 "The child catalog records the number of peaks and children "
659 "in the parent footprint so downstream measurement consumers "
660 "can recover the parent context without a separate join.",
661 )
662 pseudoColumns = pexConfig.ListField[str](
663 default=["merge_peak_sky", "sky_source"],
664 doc="Names of flags which should never be deblended.",
665 )
666 measureParents = pexConfig.Field[bool](
667 default=False,
668 doc="Whether to add parents to the object catalog for measurement in downstream tasks.",
669 )
671 # Testing options
672 # Some obs packages and ci packages run the full pipeline on a small
673 # subset of data to test that the pipeline is functioning properly.
674 # This is not meant as scientific validation, so it can be useful
675 # to only run on a small subset of the data that is large enough to
676 # test the desired pipeline features but not so long that the deblender
677 # is the tall pole in terms of execution times.
678 useCiLimits = pexConfig.Field[bool](
679 default=False,
680 doc="Limit the number of sources deblended for CI to prevent long build times",
681 )
682 ciDeblendChildRange = pexConfig.ListField[int](
683 default=[5, 10],
684 doc="Only deblend parent Footprints with a number of peaks in the (inclusive) range indicated."
685 "If `useCiLimits==False` then this parameter is ignored.",
686 )
687 ciNumParentsToDeblend = pexConfig.Field[int](
688 default=10,
689 doc="Only use the first `ciNumParentsToDeblend` parent footprints with a total peak count "
690 "within `ciDebledChildRange`. "
691 "If `useCiLimits==False` then this parameter is ignored.",
692 )
695class ScarletDeblendTask(pipeBase.Task):
696 """ScarletDeblendTask
698 Split blended sources into individual sources.
700 This task has no return value; it only modifies the SourceCatalog in-place.
701 """
703 ConfigClass = ScarletDeblendConfig
704 _DefaultName = "scarletDeblend"
706 def __init__(
707 self,
708 schema: afwTable.Schema,
709 peakSchema: afwTable.Schema = None,
710 **kwargs
711 ):
712 """Create the task, adding necessary fields to the given schema.
714 Parameters
715 ----------
716 schema :
717 Schema object for measurement fields; will be modified in-place.
718 peakSchema :
719 Schema of Footprint Peaks that will be passed to the deblender.
720 Any fields beyond the PeakTable minimal schema will be transferred
721 to the main source Schema. If None, no fields will be transferred
722 from the Peaks.
723 **kwargs
724 Passed to Task.__init__.
725 """
726 pipeBase.Task.__init__(self, **kwargs)
728 peakMinimalSchema = afwDet.PeakTable.makeMinimalSchema()
729 if peakSchema is None:
730 # In this case, the peakSchemaMapper will transfer nothing, but
731 # we'll still have one
732 # to simplify downstream code
733 self.peakSchemaMapper = afwTable.SchemaMapper(peakMinimalSchema, schema)
734 else:
735 self.peakSchemaMapper = afwTable.SchemaMapper(peakSchema, schema)
736 for item in peakSchema:
737 if item.key not in peakMinimalSchema:
738 self.peakSchemaMapper.addMapping(item.key, item.field)
739 # Because SchemaMapper makes a copy of the output schema
740 # you give its ctor, it isn't updating this Schema in
741 # place. That's probably a design flaw, but in the
742 # meantime, we'll keep that schema in sync with the
743 # peakSchemaMapper.getOutputSchema() manually, by adding
744 # the same fields to both.
745 schema.addField(item.field)
746 assert (
747 schema == self.peakSchemaMapper.getOutputSchema()
748 ), "Logic bug mapping schemas"
750 # Add the parent keys to the parent catalog schema
751 self.parentSchemaMapper = afwTable.SchemaMapper(schema)
752 self.parentSchemaMapper.addMinimalSchema(schema, doMap=True)
753 parentOutSchema = self.parentSchemaMapper.editOutputSchema()
754 self._addParentSchemaKeys(parentOutSchema)
755 self.parentSchema = parentOutSchema
756 # Mirror peakSchemaMapper's extra-field handling so any
757 # merge_peak_* (or other non-minimal) peak fields are also
758 # carried onto parent records.
759 if peakSchema is None:
760 self.parentPeakSchemaMapper = afwTable.SchemaMapper(peakMinimalSchema, self.parentSchema)
761 else:
762 self.parentPeakSchemaMapper = afwTable.SchemaMapper(peakSchema, self.parentSchema)
763 for item in peakSchema:
764 if item.key not in peakMinimalSchema:
765 self.parentPeakSchemaMapper.addMapping(item.key, item.field)
767 # Add keys for isolated sources and deblended children to the schema.
768 self._addChildSchemaKeys(schema)
769 self.objectSchema = schema
770 self.objectPeakSchemaMapper = afwTable.SchemaMapper(peakMinimalSchema, self.objectSchema)
771 self.toCopyFromParent = [
772 name
773 for item in self.objectSchema
774 if ((name := item.field.getName()).startswith("merge_footprint")
775 or (name := item.field.getName()).startswith("merge_peak"))
776 ]
778 # Any pseudoColumn must be present on the deblender schema —
779 # the science pipeline always adds these upstream (e.g.
780 # ``SkyObjectsTask`` for ``sky_source``, the coadd merge step
781 # for ``merge_peak_*``). A missing column means a typo in
782 # config or a missing upstream task; warn so the user knows
783 # the corresponding filter is silently inactive.
784 schemaNames = set(self.objectSchema.getNames())
785 unknownPseudoColumns = [
786 col for col in self.config.pseudoColumns
787 if col not in schemaNames
788 ]
789 if unknownPseudoColumns:
790 self.log.warning(
791 "pseudoColumns %s not found on the deblender schema. "
792 "Pseudo-source filtering for these columns will be "
793 "silently skipped at runtime — check for config typos "
794 "or missing upstream tasks (e.g. SkyObjectsTask "
795 "provides 'sky_source', the coadd merge provides "
796 "'merge_peak_sky').",
797 unknownPseudoColumns,
798 )
800 def _addParentSchemaKeys(self, schema: afwTable.Schema):
801 """Add parent specific keys to the schema"""
802 # Parent (blend) fields
803 schema.addField(
804 "deblend_runtime", type=np.float32, doc="runtime in ms"
805 )
806 schema.addField(
807 "deblend_iterations", type=np.int32, doc="iterations to converge"
808 )
809 schema.addField(
810 "deblend_nChild",
811 type=np.int32,
812 doc="Number of children this object has (defaults to 0)",
813 )
814 schema.addField(
815 "deblend_nPeaks",
816 type=np.int32,
817 doc="Number of initial peaks in the blend. "
818 "This includes peaks that may have been culled "
819 "during deblending or failed to deblend",
820 )
821 schema.addField(
822 "deblend_spectrumInitFlag",
823 type="Flag",
824 doc="True when scarlet initializes sources "
825 "in the blend with a more accurate spectrum. "
826 "The algorithm uses a lot of memory, "
827 "so large dense blends will use "
828 "a less accurate initialization.",
829 )
830 # Skipped flags
831 schema.addField(
832 "deblend_skipped",
833 type="Flag",
834 doc=(
835 "The deblender skipped this source. On a deconvolved "
836 "sub-blend the flag means that sub-blend itself was "
837 "skipped. On a top-level parent the flag is the union "
838 "over its deconvolved sub-blends: it fires whenever at "
839 "least one sub-blend was skipped, even when the others "
840 "succeeded. The per-reason ``deblend_skipped_*`` "
841 "sub-flags follow the same union convention on a "
842 "top-level parent, so a parent record can carry "
843 "multiple sub-flags when different sub-blends were "
844 "skipped for different reasons."
845 ),
846 )
847 schema.addField(
848 "deblend_skipped_isolatedParent",
849 type="Flag",
850 doc="The source has only a single peak " "and was not deblended",
851 )
852 schema.addField(
853 "deblend_skipped_isPseudo",
854 type="Flag",
855 doc='The source is identified as a "pseudo" source and '
856 "was not deblended",
857 )
858 schema.addField(
859 "deblend_skipped_tooManyPeaks",
860 type="Flag",
861 doc="Source had too many peaks; " "only the brightest were included",
862 )
863 schema.addField(
864 "deblend_skipped_parentTooBig",
865 type="Flag",
866 doc="Parent footprint covered too many pixels",
867 )
868 schema.addField(
869 "deblend_skipped_masked",
870 type="Flag",
871 doc="Parent footprint had too many masked pixels",
872 )
873 # Convergence flags
874 schema.addField(
875 "deblend_blendConvergenceFailedFlag",
876 type="Flag",
877 doc="at least one source in the blend" "failed to converge",
878 )
879 # Error flags
880 schema.addField(
881 "deblend_failed", type="Flag", doc="Deblending failed on source"
882 )
883 schema.addField(
884 "deblend_childFailed",
885 type="Flag",
886 doc="Deblending failed on at least one child blend. "
887 " This is set in a parent when at least one of its children "
888 "is a blend that failed to deblend.",
889 )
890 schema.addField(
891 "deblend_error",
892 type="String",
893 size=25,
894 doc="Name of error if the blend failed",
895 )
896 schema.addField(
897 "deblend_incompleteData",
898 type="Flag",
899 doc="True when a blend has at least one band "
900 "that could not generate a PSF and was "
901 "not included in the model.",
902 )
903 schema.addField(
904 "deblend_nComponents",
905 type=np.int32,
906 doc="Number of components in a ScarletLiteSource.",
907 )
908 schema.addField(
909 "deblend_chi2",
910 type=np.float32,
911 doc="Final reduced chi2 (per pixel), used to identify goodness of fit.",
912 )
914 def _addChildSchemaKeys(self, schema: afwTable.Schema):
915 """Add deblender specific keys to the schema"""
916 afwTable.Point2IKey.addFields(
917 schema,
918 name="deblend_peak_center",
919 doc="Center used to apply constraints in scarlet",
920 unit="pixel",
921 )
922 schema.addField(
923 "deblend_nChild",
924 type=np.int32,
925 doc="Number of children this object has (defaults to 0)",
926 )
927 schema.addField(
928 "deblend_peakId",
929 type=np.int32,
930 doc="ID of the peak in the parent footprint. "
931 "This is not unique, but the combination of 'parent'"
932 "and 'peakId' should be for all child sources. "
933 "Top level blends with no parents have 'peakId=0'",
934 )
935 schema.addField(
936 "deblend_peak_instFlux",
937 type=float,
938 units="count",
939 doc="The instFlux at the peak position of deblended mode",
940 )
941 schema.addField(
942 "deblend_scarletFlux", type=np.float32, doc="Flux measurement from scarlet"
943 )
944 schema.addField(
945 "deblend_edgePixels",
946 type="Flag",
947 doc="Source had flux on the edge of the parent footprint",
948 )
949 schema.addField(
950 "deblend_dataCoverage",
951 type=np.float32,
952 doc="Fraction of pixels with data. "
953 "In other words, 1 - fraction of pixels with NO_DATA set.",
954 )
955 schema.addField(
956 "deblend_zeroFlux", type="Flag", doc="Source has zero flux."
957 )
958 # Blendedness/classification metrics
959 schema.addField(
960 "deblend_maxOverlap",
961 type=np.float32,
962 doc="Maximum overlap with all of the other neighbors flux "
963 "combined."
964 "This is useful as a metric for determining how blended a "
965 "source is because if it only overlaps with other sources "
966 "at or below the noise level, it is likely to be a mostly "
967 "isolated source in the deconvolved model frame.",
968 )
969 schema.addField(
970 "deblend_fluxOverlap",
971 type=np.float32,
972 doc="This is the total flux from neighboring objects that "
973 "overlaps with this source.",
974 )
975 schema.addField(
976 "deblend_fluxOverlapFraction",
977 type=np.float32,
978 doc="This is the fraction of "
979 "`flux from neighbors/source flux` "
980 "for a given source within the source's"
981 "footprint.",
982 )
983 schema.addField(
984 "deblend_blendedness",
985 type=np.float32,
986 doc="The Bosch et al. 2018 metric for 'blendedness.' ",
987 )
988 schema.addField(
989 "deblend_blendId",
990 type=np.int64,
991 doc="Parents in the catalog may be subdivided by deblending "
992 "into multiple deconvolved blends that are each "
993 "deblended separately. This is the ID of the "
994 "deconvolved blend in the catalog."
995 )
996 schema.addField(
997 "deblend_blendNChild",
998 type=np.int32,
999 doc="The number of children in the deconvolved blend."
1000 )
1001 schema.addField(
1002 "deblend_parentNPeaks",
1003 type=np.int32,
1004 doc="deblend_nPeaks from this records parent.",
1005 )
1006 schema.addField(
1007 "deblend_parentNChild",
1008 type=np.int32,
1009 doc="deblend_nChild from this records parent.",
1010 )
1011 schema.addField(
1012 "deblend_nComponents",
1013 type=np.int32,
1014 doc="Number of components in a ScarletLiteSource.",
1015 )
1016 schema.addField(
1017 "deblend_chi2",
1018 type=np.float32,
1019 doc="Final reduced chi2 (per pixel), used to identify goodness of fit.",
1020 )
1022 @timeMethod
1023 def run(
1024 self,
1025 mExposure: afwImage.MultibandExposure,
1026 mDeconvolved: afwImage.MultibandExposure,
1027 mergedSources: afwTable.SourceCatalog,
1028 ) -> pipeBase.Struct:
1029 """Get the psf from each exposure and then run deblend().
1031 Parameters
1032 ----------
1033 mExposure :
1034 The exposures should be co-added images of the same
1035 shape and region of the sky.
1036 mDeconvolved :
1037 The deconvolved images of the same shape and region of the sky.
1038 mergedSources :
1039 The merged `SourceCatalog` that contains parent footprints
1040 to (potentially) deblend.
1042 Returns
1043 -------
1044 templateCatalogs: dict
1045 Keys are the names of the bands and the values are
1046 `lsst.afw.table.source.source.SourceCatalog`'s.
1047 These are catalogs with heavy footprints that are the templates
1048 created by the multiband templates.
1049 """
1050 # Create a table to hold object records
1051 table = afwTable.SourceTable.make(self.objectSchema)
1052 objectCatalog = afwTable.SourceCatalog(table)
1054 # Create a table to hold parent records
1055 table = afwTable.SourceTable.make(self.parentSchema, mergedSources.getIdFactory())
1056 parentCatalog = afwTable.SourceCatalog(table)
1057 parentCatalog.extend(mergedSources, self.parentSchemaMapper)
1059 return self.deblend(mExposure, mDeconvolved, objectCatalog, parentCatalog)
1061 @timeMethod
1062 def deblend(
1063 self,
1064 mExposure: afwImage.MultibandExposure,
1065 mDeconvolved: afwImage.MultibandExposure,
1066 objectCatalog: afwTable.SourceCatalog,
1067 parentCatalog: afwTable.SourceCatalog,
1068 ) -> pipeBase.Struct:
1069 """Deblend a data cube of multiband images
1071 Deblending iterates over sources from the input catalog,
1072 which are blends of peaks with overlapping PSFs (depth 0 parents).
1073 In many cases those footprints can be subdived into multiple
1074 deconvolved footprints, which have an intermediate
1075 parent record added to the catalog and are be deblended separately.
1076 All deblended peaks have a source record added to the catalog,
1077 each of which has a depth one greater than the parent.
1079 Parameters
1080 ----------
1081 mExposure :
1082 The exposures should be co-added images of the same
1083 shape and region of the sky.
1084 mDeconvolved :
1085 The deconvolved images of the same shape and region of the sky.
1086 objectCatalog :
1087 An empty `SourceCatalog` with the schema to hold isolated sources
1088 and deblended children. This catalog is filled in place.
1089 parentCatalog :
1090 The merged `SourceCatalog` that contains parent footprints
1091 to (potentially) deblend. If a parent is subdivided into
1092 multiple deconvolved parents, the deconvolved parents are
1093 added to this catalog in place.
1095 Returns
1096 -------
1097 deblendedCatalog :
1098 The ``deblendedCatalog`` isolated and deblended child sources.
1099 scarletModelData :
1100 The persistable data model for the deblender.
1101 objectParents :
1102 The parent catalog with deconvolved parents added.
1103 """
1105 # Cull footprints if required by ci
1106 if self.config.useCiLimits:
1107 self.log.info(
1108 "Using CI catalog limits, the original number of sources to deblend was %d.",
1109 len(objectCatalog),
1110 )
1111 # Select parents with a number of children in the range
1112 # config.ciDeblendChildRange
1113 minChildren, maxChildren = self.config.ciDeblendChildRange
1114 nPeaks = np.array([len(src.getFootprint().peaks) for src in parentCatalog])
1115 childrenInRange = np.where(
1116 (nPeaks >= minChildren) & (nPeaks <= maxChildren)
1117 )[0]
1118 if len(childrenInRange) < self.config.ciNumParentsToDeblend:
1119 raise ValueError(
1120 "Fewer than ciNumParentsToDeblend children were contained in the range "
1121 "indicated by ciDeblendChildRange. Adjust this range to include more "
1122 "parents."
1123 )
1124 # Keep all of the isolated parents and the first
1125 # `ciNumParentsToDeblend` children
1126 parents = nPeaks == 1
1127 children = np.zeros((len(parentCatalog),), dtype=bool)
1128 children[childrenInRange[: self.config.ciNumParentsToDeblend]] = True
1129 parentCatalog = parentCatalog[parents | children]
1130 # We need to update the IdFactory, otherwise the the source ids
1131 # will not be sequential
1132 idFactory = parentCatalog.getIdFactory()
1133 maxId = np.max(parentCatalog["id"])
1134 idFactory.notify(maxId)
1135 del children
1137 self.log.info(
1138 "Deblending %d sources in %d exposure bands", len(parentCatalog), len(mExposure),
1139 )
1140 periodicLog = PeriodicLogger(self.log)
1142 # Add the NOT_DEBLENDED mask to the mask plane in each band
1143 if self.config.notDeblendedMask:
1144 for mask in mExposure.mask:
1145 mask.addMaskPlane(self.config.notDeblendedMask)
1147 # Create the context for the entire coadd
1148 context = ScarletDeblendContext.build(
1149 mExposure=mExposure,
1150 mDeconvolved=mDeconvolved,
1151 config=self.config,
1152 catalog=parentCatalog,
1153 )
1154 nBands = len(context.observation.bands)
1156 # Initialize the persistable ScarletModelData object
1157 modelData = io.LsstScarletModelData(metadata={
1158 "model_psf": context.observation.model_psf[0],
1159 "psf": context.observation.psfs,
1160 "bands": context.observation.bands,
1161 })
1163 if self.config.persistIsolated:
1164 # Add isolated sources to the model data
1165 for source in parentCatalog:
1166 if (
1167 len(source.getFootprint().peaks) == 1
1168 and not isPseudoSource(source, self.config.pseudoColumns)
1169 ):
1170 isolated = IsolatedSource.from_footprint(
1171 footprint=source.getFootprint(),
1172 mCoadd=mExposure,
1173 dtype=np.float32,
1174 )
1175 modelData.isolated[source.getId()] = isolated.to_data()
1177 # Attach full image objects to the task to simplify the API
1178 # and use for debugging.
1179 self.objectCatalog = objectCatalog
1180 self.parentCatalog = parentCatalog
1181 self.context = context
1182 self.modelData = modelData
1183 self.mExposure = mExposure
1185 # Subdivide the psf blended parents into deconvolved parents
1186 # using the deconvolved footprints stored in the context.
1187 nParents = len(parentCatalog)
1188 self._initializeCatalogs(
1189 objectCatalog,
1190 parentCatalog,
1191 context,
1192 )
1193 nBlends = len(parentCatalog)
1195 self.log.info(
1196 "Subdivided %d top level parents to create %d deconvolved parents.",
1197 nParents,
1198 nBlends-nParents,
1199 )
1201 # Deblend sources
1202 for parentIndex in range(nParents):
1203 # Log a message if it has been a while since the last log.
1204 periodicLog.log(
1205 "Deblended %d out of %d parents",
1206 parentIndex,
1207 nParents,
1208 )
1210 parentRecord = parentCatalog[parentIndex]
1211 blendRecords = parentCatalog[parentCatalog["parent"] == parentRecord.getId()]
1213 if len(blendRecords) == 0:
1214 # The source is isolated.
1215 if self.config.processSingles:
1216 blendRecords = [parentRecord]
1217 else:
1218 continue
1220 self.log.trace(
1221 "Split parent %d into %d deconvolved parents",
1222 parentRecord.getId(),
1223 len(blendRecords),
1224 )
1225 # Create an image to keep track of the cumulative model
1226 # for all sub blends in the parent footprint.
1227 parentFootprint = parentRecord.getFootprint()
1228 bbox = parentFootprint.getBBox()
1229 width, height = bbox.getDimensions()
1230 x0, y0 = bbox.getMin()
1231 emptyModel = np.zeros(
1232 (nBands, height, width),
1233 dtype=mExposure.image.array.dtype,
1234 )
1235 parentModel = scl.Image(
1236 emptyModel,
1237 bands=context.observation.images.bands,
1238 yx0=(y0, x0),
1239 )
1241 sourceRecords = []
1242 parentBlends = {}
1243 for subBlendIndex, blendRecord in enumerate(blendRecords):
1244 # Log a message if it has been a while since the last log.
1245 periodicLog.log(
1246 "Deblending sub-blend %d/%d of parent %d",
1247 subBlendIndex + 1,
1248 len(blendRecords),
1249 parentRecord.getId(),
1250 )
1251 try:
1252 blend, blendModel, chi2 = self._deblendParent(blendRecord)
1253 except DeblenderSkippedError as e:
1254 self._skipBlend(blendRecord, e.skipKey, e.message)
1255 parentRecord.set("deblend_skipped", True)
1256 parentRecord.set(e.skipKey, True)
1257 continue
1258 except DeblenderError as e:
1259 blendRecord.set("deblend_error", e.errorName)
1260 blendRecord.set("deblend_failed", True)
1261 self._skipBlend(blendRecord, "deblend_failed", e.message)
1262 parentRecord.set("deblend_childFailed", True)
1263 continue
1265 # Update the parent model
1266 parentModel.insert(blendModel)
1268 # Add each deblended source to the catalog
1269 for scarletSource in blend.sources:
1270 # Add all fields except the HeavyFootprint to the
1271 # source record
1272 scarletSource.metadata = {}
1273 scarletSource.metadata["peak_id"] = scarletSource.detectedPeak.getId()
1274 sourceRecord = self._addDeblendedSource(
1275 parentId=parentRecord.getId(),
1276 blendRecord=blendRecord,
1277 peak=scarletSource.detectedPeak,
1278 objectCatalog=objectCatalog,
1279 scarletSource=scarletSource,
1280 chi2=chi2,
1281 )
1282 scarletSource.metadata["id"] = sourceRecord.getId()
1283 sourceRecords.append(sourceRecord)
1285 # Store the blend information so that it can be persisted
1286 blendData = blend.to_data()
1287 parentBlends[blendRecord.getId()] = blendData
1289 if len(parentBlends) == 0:
1290 # All of the deconvolved blends failed to deblend
1291 self._updateParentRecord(
1292 parentRecord=parentRecord,
1293 nPeaks=len(parentRecord.getFootprint().peaks),
1294 nChild=0,
1295 nComponents=0,
1296 runtime=np.nan,
1297 iterations=0,
1298 logL=np.nan,
1299 chi2=np.nan,
1300 spectrumInit=False,
1301 # No blend was fit, so convergence does not apply.
1302 convergenceFailed=False,
1303 )
1304 continue
1306 # Calculate the reduced chi2 for the PSF parent
1307 parentFootprintImage = parentModel.data > 0
1308 chi2 = utils.calcChi2(parentModel, context.observation, parentFootprintImage)
1309 # Defensive: a detected parent's aggregate model has
1310 # positive flux somewhere, so the area should never be 0.
1311 parentArea = np.sum(parentFootprintImage)
1312 parentReducedChi2 = (
1313 np.sum(chi2.data) / parentArea if parentArea > 0 else np.nan
1314 )
1316 # Update the parent record with the deblending results
1317 self._updateParentRecord(
1318 parentRecord=parentRecord,
1319 nPeaks=len(parentFootprint.peaks),
1320 nChild=np.sum([child["deblend_nChild"] for child in blendRecords]),
1321 nComponents=np.sum([child["deblend_nComponents"] for child in blendRecords]),
1322 runtime=np.sum([child["deblend_runtime"] for child in blendRecords]),
1323 iterations=np.sum([child["deblend_iterations"] for child in blendRecords]),
1324 logL=np.nan,
1325 chi2=parentReducedChi2,
1326 spectrumInit=np.all([
1327 child["deblend_spectrumInitFlag"]
1328 for child in blendRecords
1329 ]), # type: ignore
1330 convergenceFailed=np.any([
1331 child["deblend_blendConvergenceFailedFlag"]
1332 for child in blendRecords
1333 ]), # type: ignore
1334 )
1335 # Persist parent columns to the children
1336 for child in sourceRecords:
1337 for key in self.toCopyFromParent:
1338 child.set(key, parentRecord.get(key))
1339 for parentColumn, childColumn in self.config.columnInheritance.items():
1340 child.set(childColumn, parentRecord.get(parentColumn))
1342 # Persist the blend data
1343 modelData.blends[parentRecord.getId()] = scl.io.HierarchicalBlendData(
1344 children=parentBlends,
1345 metadata={
1346 "spans": parentFootprint.getSpans().asArray(),
1347 "origin": (y0, x0),
1348 }
1349 )
1351 nDeblendedSources = np.sum(objectCatalog["parent"] != 0)
1352 self.log.info(
1353 "Deblender results: %d parent sources were "
1354 "split into %d deconvolved parents,"
1355 "resulting in %d deblended sources, "
1356 "for a total catalog size of %d sources",
1357 nParents,
1358 nBlends,
1359 nDeblendedSources,
1360 len(objectCatalog),
1361 )
1363 table = afwTable.SourceTable.make(self.objectSchema)
1364 sortedCatalog = afwTable.SourceCatalog(table)
1365 sortedCatalog.extend(objectCatalog, deep=True)
1366 table = afwTable.SourceTable.make(self.parentSchema)
1367 sortedBlendCatalog = afwTable.SourceCatalog(table)
1368 sortedBlendCatalog.extend(parentCatalog, deep=True)
1370 return pipeBase.Struct(
1371 deblendedCatalog=sortedCatalog,
1372 scarletModelData=modelData,
1373 objectParents=sortedBlendCatalog,
1374 )
1376 def _deblendParent(
1377 self,
1378 blendRecord: afwTable.SourceRecord,
1379 ) -> tuple[scl.Blend, scl.Image, scl.Image]:
1380 """Deblend a parent source record
1382 Parameters
1383 ----------
1384 blendRecord :
1385 The parent source record to deblend.
1387 Returns
1388 -------
1389 blend :
1390 The `scl.Blend` object that contains the deblended sources.
1391 blendModel :
1392 The `scl.Image` model of the blend.
1393 chi2 :
1394 The reduced chi2 of the blend model.
1395 """
1396 footprint = blendRecord.getFootprint()
1397 bbox = footprint.getBBox()
1398 peaks = footprint.getPeaks()
1400 # Skip the source if it meets the skipping criteria
1401 isSkipped = self._checkSkipped(blendRecord, self.mExposure)
1402 if isSkipped is not None:
1403 skipKey, skipMessage = isSkipped
1404 raise DeblenderSkippedError(
1405 skipMessage,
1406 blendRecord.getId(),
1407 skipKey,
1408 )
1410 self.log.trace(
1411 "Blend %d: deblending {%d} peaks",
1412 blendRecord.getId(),
1413 len(peaks),
1414 )
1415 # Choose whether or not to use improved spectral initialization.
1416 # This significantly cuts down on the number of iterations
1417 # that the optimizer needs and usually results in a better
1418 # fit.
1419 # But using least squares on a very large blend causes memory
1420 # issues, so it is not done for large blends
1421 if self.config.setSpectra:
1422 if self.config.maxSpectrumCutoff <= 0:
1423 spectrumInit = True
1424 else:
1425 spectrumInit = (
1426 len(footprint.peaks) * bbox.getArea() < self.config.maxSpectrumCutoff
1427 )
1428 else:
1429 spectrumInit = False
1431 try:
1432 t0 = time.monotonic()
1433 # Build the parameter lists with the same ordering
1434 blend = deblend(self.context, footprint, self.config, spectrumInit)
1435 tf = time.monotonic()
1436 runtime = (tf - t0) * 1000
1437 convergenceFailed = not _checkBlendConvergence(
1438 blend, self.config.relativeError
1439 )
1440 # Store the number of components in the blend
1441 nComponents = len(blend.components)
1442 nChild = len(blend.sources)
1443 # Catch all errors and filter out the ones that we know about
1444 except Exception as e:
1445 blendError = type(e).__name__
1446 if self.config.catchFailures:
1447 # Make it easy to find UnknownErrors in the log file
1448 self.log.warning("UnknownError")
1449 traceback.print_exc()
1450 else:
1451 raise
1453 raise DeblenderError(
1454 f"Unable to deblend parent {blendRecord.getId()}: {blendError}",
1455 blendRecord.getId(),
1456 blendError,
1457 )
1459 # Calculate the reduced chi2
1460 blendModel = blend.get_model(convolve=False)
1461 blendFootprintImage = blendModel.data > 0
1462 chi2 = utils.calcChi2(blendModel, self.context.observation, blendFootprintImage)
1463 # Defensive: same unreachable-in-practice guard as the
1464 # aggregate-parent branch above.
1465 blendArea = np.sum(blendFootprintImage)
1466 blendReducedChi2 = (
1467 np.sum(chi2.data) / blendArea if blendArea > 0 else np.nan
1468 )
1470 # Update the blend record with the deblending results
1471 self._updateParentRecord(
1472 parentRecord=blendRecord,
1473 nPeaks=len(peaks),
1474 nChild=nChild,
1475 nComponents=nComponents,
1476 runtime=runtime,
1477 iterations=len(blend.loss),
1478 logL=blend.log_likelihood,
1479 chi2=blendReducedChi2,
1480 spectrumInit=spectrumInit,
1481 convergenceFailed=convergenceFailed,
1482 )
1484 return blend, blendModel, chi2
1486 def _isLargeFootprint(self, footprint: afwDet.Footprint) -> bool:
1487 """Returns whether a Footprint is large
1489 'Large' is defined by thresholds on the area, size and axis ratio,
1490 and total area of the bounding box multiplied by
1491 the number of children.
1492 These may be disabled independently by configuring them to be
1493 non-positive.
1494 """
1495 if (
1496 self.config.maxFootprintArea > 0
1497 and footprint.getBBox().getArea() > self.config.maxFootprintArea
1498 ):
1499 return True
1500 if self.config.maxFootprintSize > 0:
1501 bbox = footprint.getBBox()
1502 if max(bbox.getWidth(), bbox.getHeight()) > self.config.maxFootprintSize:
1503 return True
1504 if self.config.minFootprintAxisRatio > 0:
1505 axes = afwEll.Axes(footprint.getShape())
1506 if axes.getB() < self.config.minFootprintAxisRatio * axes.getA():
1507 return True
1508 if self.config.maxAreaTimesPeaks > 0:
1509 if (
1510 footprint.getBBox().getArea() * len(footprint.peaks)
1511 > self.config.maxAreaTimesPeaks
1512 ):
1513 return True
1514 return False
1516 def _isMasked(self, footprint: afwDet.Footprint, mExposure: afwImage.MultibandExposure) -> bool:
1517 """Returns whether the footprint violates the mask limits
1519 Parameters
1520 ----------
1521 footprint :
1522 The footprint to check for masked pixels
1523 mExposure :
1524 The multiband exposure to check for masked pixels.
1526 Returns
1527 -------
1528 isMasked : `bool`
1529 `True` if `self.config.maskPlaneLimits` is less than the
1530 fraction of pixels for a given mask in
1531 `self.config.maskLimits`.
1532 """
1533 bbox = footprint.getBBox()
1534 # AND across bands: a pixel counts as masked only when the bit
1535 # is set in every band. Matches ``buildObservation``'s per-band
1536 # weight zeroing — a pixel is truly unconstrained only when
1537 # masked in all bands; otherwise the unmasked bands still
1538 # contribute.
1539 mask = np.bitwise_and.reduce(mExposure.mask[:, bbox].array, axis=0)
1540 size = float(footprint.getArea())
1541 for maskName, limit in self.config.maskLimits.items():
1542 maskVal = mExposure.mask.getPlaneBitMask(maskName)
1543 _mask = afwImage.MaskX(mask & maskVal, xy0=bbox.getMin())
1544 # spanset of masked pixels
1545 maskedSpan = footprint.spans.intersect(_mask, maskVal)
1546 if (maskedSpan.getArea()) / size > limit:
1547 return True
1548 return False
1550 def _skipBlend(
1551 self,
1552 blendRecord: afwTable.SourceRecord,
1553 skipKey: str,
1554 logMessage: str | None,
1555 ):
1556 """Update a parent record that is not being deblended.
1558 This is a fairly trivial function but is implemented to ensure
1559 that a skipped blend updates the appropriate columns
1560 consistently, and always has a flag to mark the reason that
1561 it is being skipped.
1563 Parameters
1564 ----------
1565 blendRecord :
1566 The blend record to flag as skipped.
1567 skipKey :
1568 The name of the flag to mark the reason for skipping.
1569 logMessage :
1570 The message to display in a log.trace when a source
1571 is skipped.
1572 """
1573 if logMessage:
1574 self.log.trace(logMessage)
1575 footprint = blendRecord.getFootprint()
1576 self._updateParentRecord(
1577 parentRecord=blendRecord,
1578 nPeaks=len(footprint.peaks),
1579 nChild=0,
1580 nComponents=0,
1581 runtime=np.nan,
1582 iterations=0,
1583 logL=np.nan,
1584 chi2=np.nan,
1585 spectrumInit=False,
1586 # A skipped blend was never fit, so convergence does not
1587 # apply; leave the convergence-failure flag unset.
1588 convergenceFailed=False,
1589 )
1591 # Mark the source as skipped by the deblender and
1592 # flag the reason why.
1593 blendRecord.set("deblend_skipped", True)
1594 blendRecord.set(skipKey, True)
1596 # Add the NOT_DEBLENDED mask to the mask plane in each band
1597 if self.config.notDeblendedMask:
1598 for mask in self.mExposure.mask:
1599 footprint.spans.setMask(
1600 mask, mask.getPlaneBitMask(self.config.notDeblendedMask)
1601 )
1603 def _checkSkipped(
1604 self,
1605 parent: afwTable.SourceRecord,
1606 mExposure: afwImage.MultibandExposure
1607 ) -> tuple[afwTable.Key, str] | None:
1608 """Update a parent record that is not being deblended.
1610 This is a fairly trivial function but is implemented to ensure
1611 that a skipped parent updates the appropriate columns
1612 consistently, and always has a flag to mark the reason that
1613 it is being skipped.
1615 Parameters
1616 ----------
1617 parent :
1618 The parent record to flag as skipped.
1619 mExposure :
1620 The exposures should be co-added images of the same
1621 shape and region of the sky.
1623 Returns
1624 -------
1625 skip :
1626 `True` if the deblender will skip the parent
1627 """
1628 skipKey = None
1629 skipMessage = ""
1630 footprint = parent.getFootprint()
1631 if isPseudoSource(parent, self.config.pseudoColumns):
1632 # We also skip pseudo sources, like sky objects, which
1633 # are intended to be skipped.
1634 skipKey = "deblend_skipped_isPseudo"
1635 elif self._isLargeFootprint(footprint):
1636 # The footprint is above the maximum footprint size limit
1637 skipKey = "deblend_skipped_parentTooBig"
1638 skipMessage = f"Parent {parent.getId()}: skipping large footprint"
1639 elif self._isMasked(footprint, mExposure):
1640 # The footprint exceeds the maximum number of masked pixels
1641 skipKey = "deblend_skipped_masked"
1642 skipMessage = f"Parent {parent.getId()}: skipping masked footprint"
1643 elif (
1644 self.config.maxNumberOfPeaks > 0
1645 and len(footprint.peaks) > self.config.maxNumberOfPeaks
1646 ):
1647 # Unlike meas_deblender, in scarlet we skip the entire blend
1648 # if the number of peaks exceeds max peaks, since neglecting
1649 # to model any peaks often results in catastrophic failure
1650 # of scarlet to generate models for the brighter sources.
1651 skipKey = "deblend_skipped_tooManyPeaks"
1652 skipMessage = f"Parent {parent.getId()}: skipping blend with too many peaks"
1653 if skipKey is not None:
1654 return (cast(afwTable.Key, skipKey), skipMessage)
1655 return None
1657 def _updateParentRecord(
1658 self,
1659 parentRecord: afwTable.SourceRecord,
1660 nPeaks: int,
1661 nChild: int,
1662 nComponents: int,
1663 runtime: float,
1664 iterations: int,
1665 logL: float,
1666 chi2: float,
1667 spectrumInit: bool,
1668 convergenceFailed: bool,
1669 ):
1670 """Update a parent record in all of the single band catalogs.
1672 Ensure that all locations that update a blend record,
1673 whether it is skipped or updated after deblending,
1674 update all of the appropriate columns.
1676 Parameters
1677 ----------
1678 blendRecord :
1679 The catalog record to update.
1680 nPeaks :
1681 Number of peaks in the parent footprint.
1682 nChild :
1683 Number of children deblended from the parent.
1684 This may differ from `nPeaks` if some of the peaks
1685 were culled and have no deblended model.
1686 nComponents :
1687 Total number of components in the parent.
1688 This is usually different than the number of children,
1689 since it is common for a single source to have multiple
1690 components.
1691 runtime :
1692 Total runtime for deblending.
1693 iterations :
1694 Total number of iterations in scarlet before convergence.
1695 logL :
1696 Final log likelihood of the blend.
1697 chi2 :
1698 Final reduced chi2 of the blend.
1699 spectrumInit :
1700 True when scarlet used `set_spectra` to initialize all
1701 sources with better initial intensities.
1702 convergenceFailed :
1703 True when the blend was fit but the optimizer reached the
1704 maximum number of iterations without converging. False both
1705 when the blend converged and when no fit was attempted
1706 (skipped or isolated parents), so the flag never fires for
1707 blends where convergence does not apply.
1708 """
1709 parentRecord.set("deblend_nPeaks", nPeaks)
1710 parentRecord.set("deblend_nChild", nChild)
1711 parentRecord.set("deblend_nComponents", nComponents)
1712 parentRecord.set("deblend_runtime", runtime)
1713 parentRecord.set("deblend_iterations", iterations)
1714 parentRecord.set("deblend_spectrumInitFlag", spectrumInit)
1715 parentRecord.set(
1716 "deblend_blendConvergenceFailedFlag", convergenceFailed
1717 )
1718 parentRecord.set("deblend_chi2", chi2)
1720 def _initializeCatalogs(
1721 self,
1722 objectCatalog: afwTable.SourceCatalog,
1723 parentCatalog: afwTable.SourceCatalog,
1724 context: ScarletDeblendContext,
1725 ) -> None:
1726 """Create footprints for deconvolved parents
1728 Each parent may be subdivided into multiple blends that are
1729 isolated in deconvolved space but still blended in the image.
1730 This method finds all of the deconvolved footprints that overlap
1731 with a single parent footprint from the input catalog and
1732 returns a dictionary to map the parent ids to a list of
1733 deconvolved footprints.
1735 Parameters
1736 ----------
1737 objectCatalog :
1738 The catalog of isolated and deblended sources.
1739 parentCatalog :
1740 The catalog of parents and deconvolved parents used for deblending.
1741 context :
1742 The context for the entire coadd.
1743 """
1744 nParents = len(parentCatalog)
1746 isolated = []
1747 for n in range(nParents):
1748 parent = parentCatalog[n]
1749 parentFoot = parent.getFootprint()
1750 # Since we use the first peak for the parent object, we should
1751 # propagate its flags to the parent source.
1752 # For example, this propagates `merge_peak_sky` to the parent
1753 parent.assign(parent.getFootprint().peaks[0], self.parentPeakSchemaMapper)
1755 if isPseudoSource(parent, self.config.pseudoColumns):
1756 # Skip pseudo sources
1757 # Note: this does not flag isolated sources as skipped or
1758 # set the NOT_DEBLENDED mask in the exposure,
1759 # since these aren't really any skipped blends.
1760 isolated.append(parent)
1761 continue
1763 if len(parentFoot.peaks) == 1:
1764 # Isolated source and we are not processing singles
1765 isolated.append(parent)
1766 parent.set("deblend_skipped_isolatedParent", True)
1767 continue
1769 # Find deconvolved footprints that intersect with the parent
1770 # and add them to the blend catalog.
1771 parentId = parent.getId()
1772 self._buildIntersectingFootprints(
1773 parentId,
1774 parentFoot,
1775 parentCatalog,
1776 context.footprints,
1777 context.footprintImage
1778 )
1780 # Update the IdFactory to account for all of the parents
1781 # and deconvolved parents added to the catalog.
1782 idFactory = objectCatalog.getIdFactory()
1783 maxId = np.max(parentCatalog["id"])
1784 idFactory.notify(maxId)
1786 # Add the isolated sources to the object catalog.
1787 # This has to be done before adding the other deblended children
1788 # because a SourceCatalog must be ordered by parentId and all isolated
1789 # sources have parentId == 0.
1790 for parent in isolated:
1791 self._addIsolatedSource(parent, objectCatalog)
1793 def _buildIntersectingFootprints(
1794 self,
1795 parentId: int,
1796 afwFootprint: afwDet.Footprint,
1797 parentCatalog: afwTable.SourceCatalog,
1798 sclFootprints: list[scl.detect.Footprint],
1799 footprintImage: scl.Image,
1800 ) -> None:
1801 """Get the intersection of two footprints
1803 Parameters
1804 ----------
1805 parentId :
1806 The parent id containing the footprints.
1807 afwFootprint :
1808 The afw footprint
1809 parentCatalog :
1810 The catalog of parents and deconvolved parents.
1811 sclFootprints :
1812 List of scarlet lite Footprints.
1813 footprintImage :
1814 An indexed image of the scarlet footprints so that the value
1815 of a pixel gives the index + 1 of the footprints that
1816 contain that pixel.
1818 Returns
1819 -------
1820 intersection :
1821 The intersection of the two footprints
1822 """
1823 footprintIndices = set()
1824 ymin, xmin = footprintImage.bbox.origin
1826 # Get the index of the deconvolved footprint at the peak location
1827 height, width = footprintImage.data.shape
1828 for peak in afwFootprint.peaks:
1829 x = peak["i_x"] - xmin
1830 y = peak["i_y"] - ymin
1831 # NumPy wraps negative indices, so a try/except on
1832 # IndexError would silently accept peaks west/south of
1833 # the bbox origin. Bounds-check both directions
1834 # explicitly.
1835 if x < 0 or y < 0 or x >= width or y >= height:
1836 raise RuntimeError(f"no footprint at ({y}, {x})")
1837 footprintIndex = footprintImage.data[y, x] - 1
1838 if footprintIndex >= 0:
1839 footprintIndices.add(footprintIndex)
1841 # Get the intersection of each deconvolved footprint with
1842 # the parent footprint.
1843 for index in footprintIndices:
1844 _sclFootprint = scarletFootprintToAfw(sclFootprints[index])
1845 intersection = afwFootprint.intersect(_sclFootprint, copyPeaks=True)
1846 if len(intersection.peaks) > 0:
1847 self._addBlendRecord(
1848 parentId=parentId,
1849 parentCatalog=parentCatalog,
1850 footprint=intersection,
1851 )
1853 def _addBlendRecord(
1854 self,
1855 parentId: int,
1856 parentCatalog: afwTable.SourceCatalog,
1857 footprint: afwDet.Footprint,
1858 ) -> None:
1859 """Add deconvolved parents to the parent catalog
1861 Each parent may be subdivided into multiple blends that are
1862 isolated in deconvolved space but still blended in the image.
1863 This function adds the sub-parents to the catalog.
1865 Parameters
1866 ----------
1867 parentId :
1868 The ID of the parent of the sub-parents.
1869 parentCatalog :
1870 The catalog of parents and deconvolved parents.
1871 footprint :
1872 The footprint of the deconvolved parent.
1873 """
1874 blendRecord = parentCatalog.addNew()
1875 blendRecord.setParent(parentId)
1876 blendRecord.setFootprint(footprint)
1877 # Propagate the first peak's schema fields onto the
1878 # sub-blend record, mirroring what ``_initializeCatalogs``
1879 # does for top-level parents.
1880 blendRecord.assign(footprint.peaks[0], self.parentPeakSchemaMapper)
1882 def _addDeblendedSource(
1883 self,
1884 parentId: int,
1885 blendRecord: afwTable.SourceRecord,
1886 peak: afwDet.PeakRecord,
1887 objectCatalog: afwTable.SourceCatalog,
1888 scarletSource: scl.Source,
1889 chi2: scl.Image,
1890 ):
1891 """Add a deblended source to a catalog.
1893 This creates a new child in the source catalog,
1894 assigning it a parent id, and adding all columns
1895 that are independent across all filter bands and not
1896 updated after deblending.
1898 Parameters
1899 ----------
1900 parent :
1901 The parent of the new child record.
1902 blendRecord :
1903 The deconvolved parent of the new child record.
1904 peak :
1905 The peak record for the peak from the parent peak catalog.
1906 catalog :
1907 The source catalog that the child is added to.
1908 scarletSource :
1909 The scarlet model for the new source record.
1910 chi2 :
1911 The chi2 for each pixel.
1913 Returns
1914 -------
1915 src :
1916 The new child source record.
1917 """
1918 src = objectCatalog.addNew()
1919 # The peak catalog is the same for all bands,
1920 # so we just use the first peak catalog
1921 src.assign(peak, self.peakSchemaMapper)
1922 src.setParent(parentId)
1923 src.set("deblend_blendId", blendRecord.getId())
1924 src.set("deblend_blendNChild", len(blendRecord.getFootprint().peaks))
1925 src.set("deblend_nChild", 0)
1927 # Set the position of the peak from the parent footprint
1928 # This will make it easier to match the same source across
1929 # deblenders and across observations, where the peak
1930 # position is unlikely to change unless enough time passes
1931 # for a source to move on the sky.
1932 src.set("deblend_peak_center_x", peak["i_x"])
1933 src.set("deblend_peak_center_y", peak["i_y"])
1934 src.set("deblend_peakId", peak["id"])
1936 # Store the number of components for the source
1937 src.set("deblend_nComponents", len(scarletSource.components))
1939 # Calculate the reduced chi2 for the source. The ``chi2``
1940 # image is the blend's, so within this source's bbox it
1941 # carries contributions wherever the *combined* blend model
1942 # was positive — including pixels where neighboring sources'
1943 # models extend into this bbox. Mask by this source's own
1944 # positive-model footprint so the chi2 sum and the area
1945 # normalization are over the same pixel set.
1946 # Defensive: a detected peak's source model has positive flux,
1947 # so ``area`` should never be 0.
1948 sourceFootprint = scarletSource.get_model().data > 0
1949 area = np.sum(sourceFootprint)
1950 sourceReducedChi2 = (
1951 np.sum(chi2[:, scarletSource.bbox].data * sourceFootprint) / area
1952 if area > 0 else np.nan
1953 )
1954 src.set("deblend_chi2", sourceReducedChi2)
1956 return src
1958 def _addIsolatedSource(
1959 self,
1960 parentRecord: afwTable.SourceRecord,
1961 objectCatalog: afwTable.SourceCatalog,
1962 ) -> afwTable.SourceRecord:
1963 """Add an isolated source to the catalog.
1965 This creates a new record in the object catalog,
1966 and assigns values for all fields not dependent on deblending.
1968 Parameters
1969 ----------
1970 parentRecord :
1971 The parent record to copy values from.
1972 objectCatalog :
1973 The source catalog that the child is added to.
1975 Returns
1976 -------
1977 src :
1978 The new isolated source record.
1979 """
1980 # The isolated source has the same peak as the parent
1981 footprint = parentRecord.getFootprint()
1982 peak = footprint.peaks[0]
1983 src = objectCatalog.addNew()
1984 src.assign(peak, self.peakSchemaMapper)
1985 src.setParent(0)
1986 # The id of the isolated source is the same as the parent
1987 src.setId(parentRecord.getId())
1988 src.setFootprint(footprint)
1989 src.set("deblend_blendId", 0)
1990 src.set("deblend_blendNChild", 0)
1991 src.set("deblend_peak_center_x", peak["i_x"])
1992 src.set("deblend_peak_center_y", peak["i_y"])
1993 src.set("deblend_peakId", peak["id"])
1995 # Isolated sources are not modeled, so they don't have components
1996 src.set("deblend_nComponents", 0)
1998 # There are no children so we need to update the parent
1999 # record and add the isolated source to the object catalog
2000 self._updateParentRecord(
2001 parentRecord=parentRecord,
2002 nPeaks=len(parentRecord.getFootprint().peaks),
2003 nChild=0,
2004 nComponents=0,
2005 runtime=0.0,
2006 iterations=0,
2007 logL=np.nan,
2008 chi2=np.nan,
2009 spectrumInit=False,
2010 # An isolated source is not fit, so convergence does not apply.
2011 convergenceFailed=False,
2012 )
2014 # Persist parent columns to the isolated source
2015 # This is usually detection flags, like merge_peak_sky
2016 # or merge_footprint_r, etc.
2017 for key in self.toCopyFromParent:
2018 src.set(key, parentRecord.get(key))
2020 return src