Coverage for python/lsst/meas/extensions/scarlet/scarletDeblendTask.py: 86%
521 statements
« prev ^ index » next coverage.py v7.14.3, created at 2026-07-09 02:24 -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: 289 ↛ 290line 289 didn't jump to line 290 because the condition on line 289 was never true
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": 400 ↛ 407line 400 didn't jump to line 407 because the condition on line 400 was always true
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): 423 ↛ 424line 423 didn't jump to line 424 because the condition on line 423 was never true
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: 789 ↛ exitline 789 didn't return from function '__init__' because the condition on line 789 was always true
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: 1106 ↛ 1107line 1106 didn't jump to line 1107 because the condition on line 1106 was never true
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: 1143 ↛ 1148line 1143 didn't jump to line 1148 because the condition on line 1143 was always true
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(
1158 bands=context.observation.bands,
1159 model_psf=context.observation.model_psf[0],
1160 psf=context.observation.psfs,
1161 )
1163 if self.config.persistIsolated: 1163 ↛ 1179line 1163 didn't jump to line 1179 because the condition on line 1163 was always true
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: 1215 ↛ 1216line 1215 didn't jump to line 1216 because the condition on line 1215 was never true
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: 1337 ↛ 1338line 1337 didn't jump to line 1338 because the loop on line 1337 never started
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()] = io.LsstHierarchicalBlendData(
1344 children=parentBlends,
1345 span_array=parentFootprint.getSpans().asArray(),
1346 origin=(y0, x0),
1347 legacy_spans=False,
1348 )
1350 nDeblendedSources = np.sum(objectCatalog["parent"] != 0)
1351 self.log.info(
1352 "Deblender results: %d parent sources were "
1353 "split into %d deconvolved parents,"
1354 "resulting in %d deblended sources, "
1355 "for a total catalog size of %d sources",
1356 nParents,
1357 nBlends,
1358 nDeblendedSources,
1359 len(objectCatalog),
1360 )
1362 table = afwTable.SourceTable.make(self.objectSchema)
1363 sortedCatalog = afwTable.SourceCatalog(table)
1364 sortedCatalog.extend(objectCatalog, deep=True)
1365 table = afwTable.SourceTable.make(self.parentSchema)
1366 sortedBlendCatalog = afwTable.SourceCatalog(table)
1367 sortedBlendCatalog.extend(parentCatalog, deep=True)
1369 return pipeBase.Struct(
1370 deblendedCatalog=sortedCatalog,
1371 scarletModelData=modelData,
1372 objectParents=sortedBlendCatalog,
1373 )
1375 def _deblendParent(
1376 self,
1377 blendRecord: afwTable.SourceRecord,
1378 ) -> tuple[scl.Blend, scl.Image, scl.Image]:
1379 """Deblend a parent source record
1381 Parameters
1382 ----------
1383 blendRecord :
1384 The parent source record to deblend.
1386 Returns
1387 -------
1388 blend :
1389 The `scl.Blend` object that contains the deblended sources.
1390 blendModel :
1391 The `scl.Image` model of the blend.
1392 chi2 :
1393 The reduced chi2 of the blend model.
1394 """
1395 footprint = blendRecord.getFootprint()
1396 bbox = footprint.getBBox()
1397 peaks = footprint.getPeaks()
1399 # Skip the source if it meets the skipping criteria
1400 isSkipped = self._checkSkipped(blendRecord, self.mExposure)
1401 if isSkipped is not None:
1402 skipKey, skipMessage = isSkipped
1403 raise DeblenderSkippedError(
1404 skipMessage,
1405 blendRecord.getId(),
1406 skipKey,
1407 )
1409 self.log.trace(
1410 "Blend %d: deblending {%d} peaks",
1411 blendRecord.getId(),
1412 len(peaks),
1413 )
1414 # Choose whether or not to use improved spectral initialization.
1415 # This significantly cuts down on the number of iterations
1416 # that the optimizer needs and usually results in a better
1417 # fit.
1418 # But using least squares on a very large blend causes memory
1419 # issues, so it is not done for large blends
1420 if self.config.setSpectra: 1420 ↛ 1428line 1420 didn't jump to line 1428 because the condition on line 1420 was always true
1421 if self.config.maxSpectrumCutoff <= 0: 1421 ↛ 1422line 1421 didn't jump to line 1422 because the condition on line 1421 was never true
1422 spectrumInit = True
1423 else:
1424 spectrumInit = (
1425 len(footprint.peaks) * bbox.getArea() < self.config.maxSpectrumCutoff
1426 )
1427 else:
1428 spectrumInit = False
1430 try:
1431 t0 = time.monotonic()
1432 # Build the parameter lists with the same ordering
1433 blend = deblend(self.context, footprint, self.config, spectrumInit)
1434 tf = time.monotonic()
1435 runtime = (tf - t0) * 1000
1436 convergenceFailed = not _checkBlendConvergence(
1437 blend, self.config.relativeError
1438 )
1439 # Store the number of components in the blend
1440 nComponents = len(blend.components)
1441 nChild = len(blend.sources)
1442 # Catch all errors and filter out the ones that we know about
1443 except Exception as e:
1444 blendError = type(e).__name__
1445 if self.config.catchFailures:
1446 # Make it easy to find UnknownErrors in the log file
1447 self.log.warning("UnknownError")
1448 traceback.print_exc()
1449 else:
1450 raise
1452 raise DeblenderError(
1453 f"Unable to deblend parent {blendRecord.getId()}: {blendError}",
1454 blendRecord.getId(),
1455 blendError,
1456 )
1458 # Calculate the reduced chi2
1459 blendModel = blend.get_model(convolve=False)
1460 blendFootprintImage = blendModel.data > 0
1461 chi2 = utils.calcChi2(blendModel, self.context.observation, blendFootprintImage)
1462 # Defensive: same unreachable-in-practice guard as the
1463 # aggregate-parent branch above.
1464 blendArea = np.sum(blendFootprintImage)
1465 blendReducedChi2 = (
1466 np.sum(chi2.data) / blendArea if blendArea > 0 else np.nan
1467 )
1469 # Update the blend record with the deblending results
1470 self._updateParentRecord(
1471 parentRecord=blendRecord,
1472 nPeaks=len(peaks),
1473 nChild=nChild,
1474 nComponents=nComponents,
1475 runtime=runtime,
1476 iterations=len(blend.loss),
1477 logL=blend.log_likelihood,
1478 chi2=blendReducedChi2,
1479 spectrumInit=spectrumInit,
1480 convergenceFailed=convergenceFailed,
1481 )
1483 return blend, blendModel, chi2
1485 def _isLargeFootprint(self, footprint: afwDet.Footprint) -> bool:
1486 """Returns whether a Footprint is large
1488 'Large' is defined by thresholds on the area, size and axis ratio,
1489 and total area of the bounding box multiplied by
1490 the number of children.
1491 These may be disabled independently by configuring them to be
1492 non-positive.
1493 """
1494 if (
1495 self.config.maxFootprintArea > 0
1496 and footprint.getBBox().getArea() > self.config.maxFootprintArea
1497 ):
1498 return True
1499 if self.config.maxFootprintSize > 0: 1499 ↛ 1500line 1499 didn't jump to line 1500 because the condition on line 1499 was never true
1500 bbox = footprint.getBBox()
1501 if max(bbox.getWidth(), bbox.getHeight()) > self.config.maxFootprintSize:
1502 return True
1503 if self.config.minFootprintAxisRatio > 0: 1503 ↛ 1504line 1503 didn't jump to line 1504 because the condition on line 1503 was never true
1504 axes = afwEll.Axes(footprint.getShape())
1505 if axes.getB() < self.config.minFootprintAxisRatio * axes.getA():
1506 return True
1507 if self.config.maxAreaTimesPeaks > 0: 1507 ↛ 1513line 1507 didn't jump to line 1513 because the condition on line 1507 was always true
1508 if ( 1508 ↛ 1512line 1508 didn't jump to line 1512 because the condition on line 1508 was never true
1509 footprint.getBBox().getArea() * len(footprint.peaks)
1510 > self.config.maxAreaTimesPeaks
1511 ):
1512 return True
1513 return False
1515 def _isMasked(self, footprint: afwDet.Footprint, mExposure: afwImage.MultibandExposure) -> bool:
1516 """Returns whether the footprint violates the mask limits
1518 Parameters
1519 ----------
1520 footprint :
1521 The footprint to check for masked pixels
1522 mExposure :
1523 The multiband exposure to check for masked pixels.
1525 Returns
1526 -------
1527 isMasked : `bool`
1528 `True` if `self.config.maskPlaneLimits` is less than the
1529 fraction of pixels for a given mask in
1530 `self.config.maskLimits`.
1531 """
1532 bbox = footprint.getBBox()
1533 # AND across bands: a pixel counts as masked only when the bit
1534 # is set in every band. Matches ``buildObservation``'s per-band
1535 # weight zeroing — a pixel is truly unconstrained only when
1536 # masked in all bands; otherwise the unmasked bands still
1537 # contribute.
1538 mask = np.bitwise_and.reduce(mExposure.mask[:, bbox].array, axis=0)
1539 size = float(footprint.getArea())
1540 for maskName, limit in self.config.maskLimits.items():
1541 maskVal = mExposure.mask.getPlaneBitMask(maskName)
1542 _mask = afwImage.MaskX(mask & maskVal, xy0=bbox.getMin())
1543 # spanset of masked pixels
1544 maskedSpan = footprint.spans.intersect(_mask, maskVal)
1545 if (maskedSpan.getArea()) / size > limit:
1546 return True
1547 return False
1549 def _skipBlend(
1550 self,
1551 blendRecord: afwTable.SourceRecord,
1552 skipKey: str,
1553 logMessage: str | None,
1554 ):
1555 """Update a parent record that is not being deblended.
1557 This is a fairly trivial function but is implemented to ensure
1558 that a skipped blend updates the appropriate columns
1559 consistently, and always has a flag to mark the reason that
1560 it is being skipped.
1562 Parameters
1563 ----------
1564 blendRecord :
1565 The blend record to flag as skipped.
1566 skipKey :
1567 The name of the flag to mark the reason for skipping.
1568 logMessage :
1569 The message to display in a log.trace when a source
1570 is skipped.
1571 """
1572 if logMessage: 1572 ↛ 1574line 1572 didn't jump to line 1574 because the condition on line 1572 was always true
1573 self.log.trace(logMessage)
1574 footprint = blendRecord.getFootprint()
1575 self._updateParentRecord(
1576 parentRecord=blendRecord,
1577 nPeaks=len(footprint.peaks),
1578 nChild=0,
1579 nComponents=0,
1580 runtime=np.nan,
1581 iterations=0,
1582 logL=np.nan,
1583 chi2=np.nan,
1584 spectrumInit=False,
1585 # A skipped blend was never fit, so convergence does not
1586 # apply; leave the convergence-failure flag unset.
1587 convergenceFailed=False,
1588 )
1590 # Mark the source as skipped by the deblender and
1591 # flag the reason why.
1592 blendRecord.set("deblend_skipped", True)
1593 blendRecord.set(skipKey, True)
1595 # Add the NOT_DEBLENDED mask to the mask plane in each band
1596 if self.config.notDeblendedMask: 1596 ↛ exitline 1596 didn't return from function '_skipBlend' because the condition on line 1596 was always true
1597 for mask in self.mExposure.mask:
1598 footprint.spans.setMask(
1599 mask, mask.getPlaneBitMask(self.config.notDeblendedMask)
1600 )
1602 def _checkSkipped(
1603 self,
1604 parent: afwTable.SourceRecord,
1605 mExposure: afwImage.MultibandExposure
1606 ) -> tuple[afwTable.Key, str] | None:
1607 """Update a parent record that is not being deblended.
1609 This is a fairly trivial function but is implemented to ensure
1610 that a skipped parent updates the appropriate columns
1611 consistently, and always has a flag to mark the reason that
1612 it is being skipped.
1614 Parameters
1615 ----------
1616 parent :
1617 The parent record to flag as skipped.
1618 mExposure :
1619 The exposures should be co-added images of the same
1620 shape and region of the sky.
1622 Returns
1623 -------
1624 skip :
1625 `True` if the deblender will skip the parent
1626 """
1627 skipKey = None
1628 skipMessage = ""
1629 footprint = parent.getFootprint()
1630 if isPseudoSource(parent, self.config.pseudoColumns): 1630 ↛ 1633line 1630 didn't jump to line 1633 because the condition on line 1630 was never true
1631 # We also skip pseudo sources, like sky objects, which
1632 # are intended to be skipped.
1633 skipKey = "deblend_skipped_isPseudo"
1634 elif self._isLargeFootprint(footprint):
1635 # The footprint is above the maximum footprint size limit
1636 skipKey = "deblend_skipped_parentTooBig"
1637 skipMessage = f"Parent {parent.getId()}: skipping large footprint"
1638 elif self._isMasked(footprint, mExposure): 1638 ↛ 1640line 1638 didn't jump to line 1640 because the condition on line 1638 was never true
1639 # The footprint exceeds the maximum number of masked pixels
1640 skipKey = "deblend_skipped_masked"
1641 skipMessage = f"Parent {parent.getId()}: skipping masked footprint"
1642 elif (
1643 self.config.maxNumberOfPeaks > 0
1644 and len(footprint.peaks) > self.config.maxNumberOfPeaks
1645 ):
1646 # Unlike meas_deblender, in scarlet we skip the entire blend
1647 # if the number of peaks exceeds max peaks, since neglecting
1648 # to model any peaks often results in catastrophic failure
1649 # of scarlet to generate models for the brighter sources.
1650 skipKey = "deblend_skipped_tooManyPeaks"
1651 skipMessage = f"Parent {parent.getId()}: skipping blend with too many peaks"
1652 if skipKey is not None:
1653 return (cast(afwTable.Key, skipKey), skipMessage)
1654 return None
1656 def _updateParentRecord(
1657 self,
1658 parentRecord: afwTable.SourceRecord,
1659 nPeaks: int,
1660 nChild: int,
1661 nComponents: int,
1662 runtime: float,
1663 iterations: int,
1664 logL: float,
1665 chi2: float,
1666 spectrumInit: bool,
1667 convergenceFailed: bool,
1668 ):
1669 """Update a parent record in all of the single band catalogs.
1671 Ensure that all locations that update a blend record,
1672 whether it is skipped or updated after deblending,
1673 update all of the appropriate columns.
1675 Parameters
1676 ----------
1677 blendRecord :
1678 The catalog record to update.
1679 nPeaks :
1680 Number of peaks in the parent footprint.
1681 nChild :
1682 Number of children deblended from the parent.
1683 This may differ from `nPeaks` if some of the peaks
1684 were culled and have no deblended model.
1685 nComponents :
1686 Total number of components in the parent.
1687 This is usually different than the number of children,
1688 since it is common for a single source to have multiple
1689 components.
1690 runtime :
1691 Total runtime for deblending.
1692 iterations :
1693 Total number of iterations in scarlet before convergence.
1694 logL :
1695 Final log likelihood of the blend.
1696 chi2 :
1697 Final reduced chi2 of the blend.
1698 spectrumInit :
1699 True when scarlet used `set_spectra` to initialize all
1700 sources with better initial intensities.
1701 convergenceFailed :
1702 True when the blend was fit but the optimizer reached the
1703 maximum number of iterations without converging. False both
1704 when the blend converged and when no fit was attempted
1705 (skipped or isolated parents), so the flag never fires for
1706 blends where convergence does not apply.
1707 """
1708 parentRecord.set("deblend_nPeaks", nPeaks)
1709 parentRecord.set("deblend_nChild", nChild)
1710 parentRecord.set("deblend_nComponents", nComponents)
1711 parentRecord.set("deblend_runtime", runtime)
1712 parentRecord.set("deblend_iterations", iterations)
1713 parentRecord.set("deblend_spectrumInitFlag", spectrumInit)
1714 parentRecord.set(
1715 "deblend_blendConvergenceFailedFlag", convergenceFailed
1716 )
1717 parentRecord.set("deblend_chi2", chi2)
1719 def _initializeCatalogs(
1720 self,
1721 objectCatalog: afwTable.SourceCatalog,
1722 parentCatalog: afwTable.SourceCatalog,
1723 context: ScarletDeblendContext,
1724 ) -> None:
1725 """Create footprints for deconvolved parents
1727 Each parent may be subdivided into multiple blends that are
1728 isolated in deconvolved space but still blended in the image.
1729 This method finds all of the deconvolved footprints that overlap
1730 with a single parent footprint from the input catalog and
1731 returns a dictionary to map the parent ids to a list of
1732 deconvolved footprints.
1734 Parameters
1735 ----------
1736 objectCatalog :
1737 The catalog of isolated and deblended sources.
1738 parentCatalog :
1739 The catalog of parents and deconvolved parents used for deblending.
1740 context :
1741 The context for the entire coadd.
1742 """
1743 nParents = len(parentCatalog)
1745 isolated = []
1746 for n in range(nParents):
1747 parent = parentCatalog[n]
1748 parentFoot = parent.getFootprint()
1749 # Since we use the first peak for the parent object, we should
1750 # propagate its flags to the parent source.
1751 # For example, this propagates `merge_peak_sky` to the parent
1752 parent.assign(parent.getFootprint().peaks[0], self.parentPeakSchemaMapper)
1754 if isPseudoSource(parent, self.config.pseudoColumns): 1754 ↛ 1759line 1754 didn't jump to line 1759 because the condition on line 1754 was never true
1755 # Skip pseudo sources
1756 # Note: this does not flag isolated sources as skipped or
1757 # set the NOT_DEBLENDED mask in the exposure,
1758 # since these aren't really any skipped blends.
1759 isolated.append(parent)
1760 continue
1762 if len(parentFoot.peaks) == 1:
1763 # Isolated source and we are not processing singles
1764 isolated.append(parent)
1765 parent.set("deblend_skipped_isolatedParent", True)
1766 continue
1768 # Find deconvolved footprints that intersect with the parent
1769 # and add them to the blend catalog.
1770 parentId = parent.getId()
1771 self._buildIntersectingFootprints(
1772 parentId,
1773 parentFoot,
1774 parentCatalog,
1775 context.footprints,
1776 context.footprintImage
1777 )
1779 # Update the IdFactory to account for all of the parents
1780 # and deconvolved parents added to the catalog.
1781 idFactory = objectCatalog.getIdFactory()
1782 maxId = np.max(parentCatalog["id"])
1783 idFactory.notify(maxId)
1785 # Add the isolated sources to the object catalog.
1786 # This has to be done before adding the other deblended children
1787 # because a SourceCatalog must be ordered by parentId and all isolated
1788 # sources have parentId == 0.
1789 for parent in isolated:
1790 self._addIsolatedSource(parent, objectCatalog)
1792 def _buildIntersectingFootprints(
1793 self,
1794 parentId: int,
1795 afwFootprint: afwDet.Footprint,
1796 parentCatalog: afwTable.SourceCatalog,
1797 sclFootprints: list[scl.detect.Footprint],
1798 footprintImage: scl.Image,
1799 ) -> None:
1800 """Get the intersection of two footprints
1802 Parameters
1803 ----------
1804 parentId :
1805 The parent id containing the footprints.
1806 afwFootprint :
1807 The afw footprint
1808 parentCatalog :
1809 The catalog of parents and deconvolved parents.
1810 sclFootprints :
1811 List of scarlet lite Footprints.
1812 footprintImage :
1813 An indexed image of the scarlet footprints so that the value
1814 of a pixel gives the index + 1 of the footprints that
1815 contain that pixel.
1817 Returns
1818 -------
1819 intersection :
1820 The intersection of the two footprints
1821 """
1822 footprintIndices = set()
1823 ymin, xmin = footprintImage.bbox.origin
1825 # Get the index of the deconvolved footprint at the peak location
1826 height, width = footprintImage.data.shape
1827 for peak in afwFootprint.peaks:
1828 x = peak["i_x"] - xmin
1829 y = peak["i_y"] - ymin
1830 # NumPy wraps negative indices, so a try/except on
1831 # IndexError would silently accept peaks west/south of
1832 # the bbox origin. Bounds-check both directions
1833 # explicitly.
1834 if x < 0 or y < 0 or x >= width or y >= height:
1835 raise RuntimeError(f"no footprint at ({y}, {x})")
1836 footprintIndex = footprintImage.data[y, x] - 1
1837 if footprintIndex >= 0: 1837 ↛ 1827line 1837 didn't jump to line 1827 because the condition on line 1837 was always true
1838 footprintIndices.add(footprintIndex)
1840 # Get the intersection of each deconvolved footprint with
1841 # the parent footprint.
1842 for index in footprintIndices:
1843 _sclFootprint = scarletFootprintToAfw(sclFootprints[index])
1844 intersection = afwFootprint.intersect(_sclFootprint, copyPeaks=True)
1845 if len(intersection.peaks) > 0: 1845 ↛ 1842line 1845 didn't jump to line 1842 because the condition on line 1845 was always true
1846 self._addBlendRecord(
1847 parentId=parentId,
1848 parentCatalog=parentCatalog,
1849 footprint=intersection,
1850 )
1852 def _addBlendRecord(
1853 self,
1854 parentId: int,
1855 parentCatalog: afwTable.SourceCatalog,
1856 footprint: afwDet.Footprint,
1857 ) -> None:
1858 """Add deconvolved parents to the parent catalog
1860 Each parent may be subdivided into multiple blends that are
1861 isolated in deconvolved space but still blended in the image.
1862 This function adds the sub-parents to the catalog.
1864 Parameters
1865 ----------
1866 parentId :
1867 The ID of the parent of the sub-parents.
1868 parentCatalog :
1869 The catalog of parents and deconvolved parents.
1870 footprint :
1871 The footprint of the deconvolved parent.
1872 """
1873 blendRecord = parentCatalog.addNew()
1874 blendRecord.setParent(parentId)
1875 blendRecord.setFootprint(footprint)
1876 # Propagate the first peak's schema fields onto the
1877 # sub-blend record, mirroring what ``_initializeCatalogs``
1878 # does for top-level parents.
1879 blendRecord.assign(footprint.peaks[0], self.parentPeakSchemaMapper)
1881 def _addDeblendedSource(
1882 self,
1883 parentId: int,
1884 blendRecord: afwTable.SourceRecord,
1885 peak: afwDet.PeakRecord,
1886 objectCatalog: afwTable.SourceCatalog,
1887 scarletSource: scl.Source,
1888 chi2: scl.Image,
1889 ):
1890 """Add a deblended source to a catalog.
1892 This creates a new child in the source catalog,
1893 assigning it a parent id, and adding all columns
1894 that are independent across all filter bands and not
1895 updated after deblending.
1897 Parameters
1898 ----------
1899 parent :
1900 The parent of the new child record.
1901 blendRecord :
1902 The deconvolved parent of the new child record.
1903 peak :
1904 The peak record for the peak from the parent peak catalog.
1905 catalog :
1906 The source catalog that the child is added to.
1907 scarletSource :
1908 The scarlet model for the new source record.
1909 chi2 :
1910 The chi2 for each pixel.
1912 Returns
1913 -------
1914 src :
1915 The new child source record.
1916 """
1917 src = objectCatalog.addNew()
1918 # The peak catalog is the same for all bands,
1919 # so we just use the first peak catalog
1920 src.assign(peak, self.peakSchemaMapper)
1921 src.setParent(parentId)
1922 src.set("deblend_blendId", blendRecord.getId())
1923 src.set("deblend_blendNChild", len(blendRecord.getFootprint().peaks))
1924 src.set("deblend_nChild", 0)
1926 # Set the position of the peak from the parent footprint
1927 # This will make it easier to match the same source across
1928 # deblenders and across observations, where the peak
1929 # position is unlikely to change unless enough time passes
1930 # for a source to move on the sky.
1931 src.set("deblend_peak_center_x", peak["i_x"])
1932 src.set("deblend_peak_center_y", peak["i_y"])
1933 src.set("deblend_peakId", peak["id"])
1935 # Store the number of components for the source
1936 src.set("deblend_nComponents", len(scarletSource.components))
1938 # Calculate the reduced chi2 for the source. The ``chi2``
1939 # image is the blend's, so within this source's bbox it
1940 # carries contributions wherever the *combined* blend model
1941 # was positive — including pixels where neighboring sources'
1942 # models extend into this bbox. Mask by this source's own
1943 # positive-model footprint so the chi2 sum and the area
1944 # normalization are over the same pixel set.
1945 # Defensive: a detected peak's source model has positive flux,
1946 # so ``area`` should never be 0.
1947 sourceFootprint = scarletSource.get_model().data > 0
1948 area = np.sum(sourceFootprint)
1949 sourceReducedChi2 = (
1950 np.sum(chi2[:, scarletSource.bbox].data * sourceFootprint) / area
1951 if area > 0 else np.nan
1952 )
1953 src.set("deblend_chi2", sourceReducedChi2)
1955 return src
1957 def _addIsolatedSource(
1958 self,
1959 parentRecord: afwTable.SourceRecord,
1960 objectCatalog: afwTable.SourceCatalog,
1961 ) -> afwTable.SourceRecord:
1962 """Add an isolated source to the catalog.
1964 This creates a new record in the object catalog,
1965 and assigns values for all fields not dependent on deblending.
1967 Parameters
1968 ----------
1969 parentRecord :
1970 The parent record to copy values from.
1971 objectCatalog :
1972 The source catalog that the child is added to.
1974 Returns
1975 -------
1976 src :
1977 The new isolated source record.
1978 """
1979 # The isolated source has the same peak as the parent
1980 footprint = parentRecord.getFootprint()
1981 peak = footprint.peaks[0]
1982 src = objectCatalog.addNew()
1983 src.assign(peak, self.peakSchemaMapper)
1984 src.setParent(0)
1985 # The id of the isolated source is the same as the parent
1986 src.setId(parentRecord.getId())
1987 src.setFootprint(footprint)
1988 src.set("deblend_blendId", 0)
1989 src.set("deblend_blendNChild", 0)
1990 src.set("deblend_peak_center_x", peak["i_x"])
1991 src.set("deblend_peak_center_y", peak["i_y"])
1992 src.set("deblend_peakId", peak["id"])
1994 # Isolated sources are not modeled, so they don't have components
1995 src.set("deblend_nComponents", 0)
1997 # There are no children so we need to update the parent
1998 # record and add the isolated source to the object catalog
1999 self._updateParentRecord(
2000 parentRecord=parentRecord,
2001 nPeaks=len(parentRecord.getFootprint().peaks),
2002 nChild=0,
2003 nComponents=0,
2004 runtime=0.0,
2005 iterations=0,
2006 logL=np.nan,
2007 chi2=np.nan,
2008 spectrumInit=False,
2009 # An isolated source is not fit, so convergence does not apply.
2010 convergenceFailed=False,
2011 )
2013 # Persist parent columns to the isolated source
2014 # This is usually detection flags, like merge_peak_sky
2015 # or merge_footprint_r, etc.
2016 for key in self.toCopyFromParent: 2016 ↛ 2017line 2016 didn't jump to line 2017 because the loop on line 2016 never started
2017 src.set(key, parentRecord.get(key))
2019 return src