Coverage for python/lsst/meas/extensions/scarlet/scarletDeblendTask.py: 86%
521 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-23 09:08 +0000
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 ]
777 # Any pseudoColumn must be present on the deblender schema —
778 # the science pipeline always adds these upstream (e.g.
779 # ``SkyObjectsTask`` for ``sky_source``, the coadd merge step
780 # for ``merge_peak_*``). A missing column means a typo in
781 # config or a missing upstream task; warn so the user knows
782 # the corresponding filter is silently inactive.
783 schemaNames = set(self.objectSchema.getNames())
784 unknownPseudoColumns = [
785 col for col in self.config.pseudoColumns
786 if col not in schemaNames
787 ]
788 if unknownPseudoColumns: 788 ↛ exitline 788 didn't return from function '__init__' because the condition on line 788 was always true
789 self.log.warning(
790 "pseudoColumns %s not found on the deblender schema. "
791 "Pseudo-source filtering for these columns will be "
792 "silently skipped at runtime — check for config typos "
793 "or missing upstream tasks (e.g. SkyObjectsTask "
794 "provides 'sky_source', the coadd merge provides "
795 "'merge_peak_sky').",
796 unknownPseudoColumns,
797 )
799 def _addParentSchemaKeys(self, schema: afwTable.Schema):
800 """Add parent specific keys to the schema"""
801 # Parent (blend) fields
802 schema.addField(
803 "deblend_runtime", type=np.float32, doc="runtime in ms"
804 )
805 schema.addField(
806 "deblend_iterations", type=np.int32, doc="iterations to converge"
807 )
808 schema.addField(
809 "deblend_nChild",
810 type=np.int32,
811 doc="Number of children this object has (defaults to 0)",
812 )
813 schema.addField(
814 "deblend_nPeaks",
815 type=np.int32,
816 doc="Number of initial peaks in the blend. "
817 "This includes peaks that may have been culled "
818 "during deblending or failed to deblend",
819 )
820 schema.addField(
821 "deblend_spectrumInitFlag",
822 type="Flag",
823 doc="True when scarlet initializes sources "
824 "in the blend with a more accurate spectrum. "
825 "The algorithm uses a lot of memory, "
826 "so large dense blends will use "
827 "a less accurate initialization.",
828 )
829 # Skipped flags
830 schema.addField(
831 "deblend_skipped",
832 type="Flag",
833 doc=(
834 "The deblender skipped this source. On a deconvolved "
835 "sub-blend the flag means that sub-blend itself was "
836 "skipped. On a top-level parent the flag is the union "
837 "over its deconvolved sub-blends: it fires whenever at "
838 "least one sub-blend was skipped, even when the others "
839 "succeeded. The per-reason ``deblend_skipped_*`` "
840 "sub-flags follow the same union convention on a "
841 "top-level parent, so a parent record can carry "
842 "multiple sub-flags when different sub-blends were "
843 "skipped for different reasons."
844 ),
845 )
846 schema.addField(
847 "deblend_skipped_isolatedParent",
848 type="Flag",
849 doc="The source has only a single peak " "and was not deblended",
850 )
851 schema.addField(
852 "deblend_skipped_isPseudo",
853 type="Flag",
854 doc='The source is identified as a "pseudo" source and '
855 "was not deblended",
856 )
857 schema.addField(
858 "deblend_skipped_tooManyPeaks",
859 type="Flag",
860 doc="Source had too many peaks; " "only the brightest were included",
861 )
862 schema.addField(
863 "deblend_skipped_parentTooBig",
864 type="Flag",
865 doc="Parent footprint covered too many pixels",
866 )
867 schema.addField(
868 "deblend_skipped_masked",
869 type="Flag",
870 doc="Parent footprint had too many masked pixels",
871 )
872 # Convergence flags
873 schema.addField(
874 "deblend_blendConvergenceFailedFlag",
875 type="Flag",
876 doc="at least one source in the blend" "failed to converge",
877 )
878 # Error flags
879 schema.addField(
880 "deblend_failed", type="Flag", doc="Deblending failed on source"
881 )
882 schema.addField(
883 "deblend_childFailed",
884 type="Flag",
885 doc="Deblending failed on at least one child blend. "
886 " This is set in a parent when at least one of its children "
887 "is a blend that failed to deblend.",
888 )
889 schema.addField(
890 "deblend_error",
891 type="String",
892 size=25,
893 doc="Name of error if the blend failed",
894 )
895 schema.addField(
896 "deblend_incompleteData",
897 type="Flag",
898 doc="True when a blend has at least one band "
899 "that could not generate a PSF and was "
900 "not included in the model.",
901 )
902 schema.addField(
903 "deblend_nComponents",
904 type=np.int32,
905 doc="Number of components in a ScarletLiteSource.",
906 )
907 schema.addField(
908 "deblend_chi2",
909 type=np.float32,
910 doc="Final reduced chi2 (per pixel), used to identify goodness of fit.",
911 )
913 def _addChildSchemaKeys(self, schema: afwTable.Schema):
914 """Add deblender specific keys to the schema"""
915 afwTable.Point2IKey.addFields(
916 schema,
917 name="deblend_peak_center",
918 doc="Center used to apply constraints in scarlet",
919 unit="pixel",
920 )
921 schema.addField(
922 "deblend_nChild",
923 type=np.int32,
924 doc="Number of children this object has (defaults to 0)",
925 )
926 schema.addField(
927 "deblend_peakId",
928 type=np.int32,
929 doc="ID of the peak in the parent footprint. "
930 "This is not unique, but the combination of 'parent'"
931 "and 'peakId' should be for all child sources. "
932 "Top level blends with no parents have 'peakId=0'",
933 )
934 schema.addField(
935 "deblend_peak_instFlux",
936 type=float,
937 units="count",
938 doc="The instFlux at the peak position of deblended mode",
939 )
940 schema.addField(
941 "deblend_scarletFlux", type=np.float32, doc="Flux measurement from scarlet"
942 )
943 schema.addField(
944 "deblend_edgePixels",
945 type="Flag",
946 doc="Source had flux on the edge of the parent footprint",
947 )
948 schema.addField(
949 "deblend_dataCoverage",
950 type=np.float32,
951 doc="Fraction of pixels with data. "
952 "In other words, 1 - fraction of pixels with NO_DATA set.",
953 )
954 schema.addField(
955 "deblend_zeroFlux", type="Flag", doc="Source has zero flux."
956 )
957 # Blendedness/classification metrics
958 schema.addField(
959 "deblend_maxOverlap",
960 type=np.float32,
961 doc="Maximum overlap with all of the other neighbors flux "
962 "combined."
963 "This is useful as a metric for determining how blended a "
964 "source is because if it only overlaps with other sources "
965 "at or below the noise level, it is likely to be a mostly "
966 "isolated source in the deconvolved model frame.",
967 )
968 schema.addField(
969 "deblend_fluxOverlap",
970 type=np.float32,
971 doc="This is the total flux from neighboring objects that "
972 "overlaps with this source.",
973 )
974 schema.addField(
975 "deblend_fluxOverlapFraction",
976 type=np.float32,
977 doc="This is the fraction of "
978 "`flux from neighbors/source flux` "
979 "for a given source within the source's"
980 "footprint.",
981 )
982 schema.addField(
983 "deblend_blendedness",
984 type=np.float32,
985 doc="The Bosch et al. 2018 metric for 'blendedness.' ",
986 )
987 schema.addField(
988 "deblend_blendId",
989 type=np.int64,
990 doc="Parents in the catalog may be subdivided by deblending "
991 "into multiple deconvolved blends that are each "
992 "deblended separately. This is the ID of the "
993 "deconvolved blend in the catalog."
994 )
995 schema.addField(
996 "deblend_blendNChild",
997 type=np.int32,
998 doc="The number of children in the deconvolved blend."
999 )
1000 schema.addField(
1001 "deblend_parentNPeaks",
1002 type=np.int32,
1003 doc="deblend_nPeaks from this records parent.",
1004 )
1005 schema.addField(
1006 "deblend_parentNChild",
1007 type=np.int32,
1008 doc="deblend_nChild from this records parent.",
1009 )
1010 schema.addField(
1011 "deblend_nComponents",
1012 type=np.int32,
1013 doc="Number of components in a ScarletLiteSource.",
1014 )
1015 schema.addField(
1016 "deblend_chi2",
1017 type=np.float32,
1018 doc="Final reduced chi2 (per pixel), used to identify goodness of fit.",
1019 )
1021 @timeMethod
1022 def run(
1023 self,
1024 mExposure: afwImage.MultibandExposure,
1025 mDeconvolved: afwImage.MultibandExposure,
1026 mergedSources: afwTable.SourceCatalog,
1027 ) -> pipeBase.Struct:
1028 """Get the psf from each exposure and then run deblend().
1030 Parameters
1031 ----------
1032 mExposure :
1033 The exposures should be co-added images of the same
1034 shape and region of the sky.
1035 mDeconvolved :
1036 The deconvolved images of the same shape and region of the sky.
1037 mergedSources :
1038 The merged `SourceCatalog` that contains parent footprints
1039 to (potentially) deblend.
1041 Returns
1042 -------
1043 result : `lsst.pipe.base.Struct`
1044 See `deblend`; its output is returned unmodified.
1045 """
1046 # Create a table to hold object records
1047 table = afwTable.SourceTable.make(self.objectSchema)
1048 objectCatalog = afwTable.SourceCatalog(table)
1050 # Create a table to hold parent records
1051 table = afwTable.SourceTable.make(self.parentSchema, mergedSources.getIdFactory())
1052 parentCatalog = afwTable.SourceCatalog(table)
1053 parentCatalog.extend(mergedSources, self.parentSchemaMapper)
1055 return self.deblend(mExposure, mDeconvolved, objectCatalog, parentCatalog)
1057 @timeMethod
1058 def deblend(
1059 self,
1060 mExposure: afwImage.MultibandExposure,
1061 mDeconvolved: afwImage.MultibandExposure,
1062 objectCatalog: afwTable.SourceCatalog,
1063 parentCatalog: afwTable.SourceCatalog,
1064 ) -> pipeBase.Struct:
1065 """Deblend a data cube of multiband images
1067 Deblending iterates over sources from the input catalog,
1068 which are blends of peaks with overlapping PSFs (depth 0 parents).
1069 In many cases those footprints can be subdived into multiple
1070 deconvolved footprints, which have an intermediate
1071 parent record added to the catalog and are be deblended separately.
1072 All deblended peaks have a source record added to the catalog,
1073 each of which has a depth one greater than the parent.
1075 Parameters
1076 ----------
1077 mExposure :
1078 The exposures should be co-added images of the same
1079 shape and region of the sky.
1080 mDeconvolved :
1081 The deconvolved images of the same shape and region of the sky.
1082 objectCatalog :
1083 An empty `SourceCatalog` with the schema to hold isolated sources
1084 and deblended children. This catalog is filled in place.
1085 parentCatalog :
1086 The merged `SourceCatalog` that contains parent footprints
1087 to (potentially) deblend. If a parent is subdivided into
1088 multiple deconvolved parents, the deconvolved parents are
1089 added to this catalog in place.
1091 Returns
1092 -------
1093 result : `lsst.pipe.base.Struct`
1094 A struct with the following attributes (all
1095 `lsst.afw.table.SourceCatalog` instances):
1097 - ``deblendedCatalog``: isolated and deblended child sources.
1098 - ``scarletModelData``: persistable deblended model data for the
1099 deblender.
1100 - ``objectParents``: the parent catalog with deconvolved parents
1101 added.
1102 """
1104 # Cull footprints if required by ci
1105 if self.config.useCiLimits: 1105 ↛ 1106line 1105 didn't jump to line 1106 because the condition on line 1105 was never true
1106 self.log.info(
1107 "Using CI catalog limits, the original number of sources to deblend was %d.",
1108 len(objectCatalog),
1109 )
1110 # Select parents with a number of children in the range
1111 # config.ciDeblendChildRange
1112 minChildren, maxChildren = self.config.ciDeblendChildRange
1113 nPeaks = np.array([len(src.getFootprint().peaks) for src in parentCatalog])
1114 childrenInRange = np.where(
1115 (nPeaks >= minChildren) & (nPeaks <= maxChildren)
1116 )[0]
1117 if len(childrenInRange) < self.config.ciNumParentsToDeblend:
1118 raise ValueError(
1119 "Fewer than ciNumParentsToDeblend children were contained in the range "
1120 "indicated by ciDeblendChildRange. Adjust this range to include more "
1121 "parents."
1122 )
1123 # Keep all of the isolated parents and the first
1124 # `ciNumParentsToDeblend` children
1125 parents = nPeaks == 1
1126 children = np.zeros((len(parentCatalog),), dtype=bool)
1127 children[childrenInRange[: self.config.ciNumParentsToDeblend]] = True
1128 parentCatalog = parentCatalog[parents | children]
1129 # We need to update the IdFactory, otherwise the the source ids
1130 # will not be sequential
1131 idFactory = parentCatalog.getIdFactory()
1132 maxId = np.max(parentCatalog["id"])
1133 idFactory.notify(maxId)
1134 del children
1136 self.log.info(
1137 "Deblending %d sources in %d exposure bands", len(parentCatalog), len(mExposure),
1138 )
1139 periodicLog = PeriodicLogger(self.log)
1141 # Add the NOT_DEBLENDED mask to the mask plane in each band
1142 if self.config.notDeblendedMask: 1142 ↛ 1147line 1142 didn't jump to line 1147 because the condition on line 1142 was always true
1143 for mask in mExposure.mask:
1144 mask.addMaskPlane(self.config.notDeblendedMask)
1146 # Create the context for the entire coadd
1147 context = ScarletDeblendContext.build(
1148 mExposure=mExposure,
1149 mDeconvolved=mDeconvolved,
1150 config=self.config,
1151 catalog=parentCatalog,
1152 )
1153 nBands = len(context.observation.bands)
1155 # Initialize the persistable ScarletModelData object
1156 modelData = io.LsstScarletModelData(
1157 bands=context.observation.bands,
1158 model_psf=context.observation.model_psf[0],
1159 psf=context.observation.psfs,
1160 )
1162 if self.config.persistIsolated: 1162 ↛ 1178line 1162 didn't jump to line 1178 because the condition on line 1162 was always true
1163 # Add isolated sources to the model data
1164 for source in parentCatalog:
1165 if (
1166 len(source.getFootprint().peaks) == 1
1167 and not isPseudoSource(source, self.config.pseudoColumns)
1168 ):
1169 isolated = IsolatedSource.from_footprint(
1170 footprint=source.getFootprint(),
1171 mCoadd=mExposure,
1172 dtype=np.float32,
1173 )
1174 modelData.isolated[source.getId()] = isolated.to_data()
1176 # Attach full image objects to the task to simplify the API
1177 # and use for debugging.
1178 self.objectCatalog = objectCatalog
1179 self.parentCatalog = parentCatalog
1180 self.context = context
1181 self.modelData = modelData
1182 self.mExposure = mExposure
1184 # Subdivide the psf blended parents into deconvolved parents
1185 # using the deconvolved footprints stored in the context.
1186 nParents = len(parentCatalog)
1187 self._initializeCatalogs(
1188 objectCatalog,
1189 parentCatalog,
1190 context,
1191 )
1192 nBlends = len(parentCatalog)
1194 self.log.info(
1195 "Subdivided %d top level parents to create %d deconvolved parents.",
1196 nParents,
1197 nBlends-nParents,
1198 )
1200 # Deblend sources
1201 for parentIndex in range(nParents):
1202 # Log a message if it has been a while since the last log.
1203 periodicLog.log(
1204 "Deblended %d out of %d parents",
1205 parentIndex,
1206 nParents,
1207 )
1209 parentRecord = parentCatalog[parentIndex]
1210 blendRecords = parentCatalog[parentCatalog["parent"] == parentRecord.getId()]
1212 if len(blendRecords) == 0:
1213 # The source is isolated.
1214 if self.config.processSingles: 1214 ↛ 1215line 1214 didn't jump to line 1215 because the condition on line 1214 was never true
1215 blendRecords = [parentRecord]
1216 else:
1217 continue
1219 self.log.trace(
1220 "Split parent %d into %d deconvolved parents",
1221 parentRecord.getId(),
1222 len(blendRecords),
1223 )
1224 # Create an image to keep track of the cumulative model
1225 # for all sub blends in the parent footprint.
1226 parentFootprint = parentRecord.getFootprint()
1227 bbox = parentFootprint.getBBox()
1228 width, height = bbox.getDimensions()
1229 x0, y0 = bbox.getMin()
1230 emptyModel = np.zeros(
1231 (nBands, height, width),
1232 dtype=mExposure.image.array.dtype,
1233 )
1234 parentModel = scl.Image(
1235 emptyModel,
1236 bands=context.observation.images.bands,
1237 yx0=(y0, x0),
1238 )
1240 sourceRecords = []
1241 parentBlends = {}
1242 for subBlendIndex, blendRecord in enumerate(blendRecords):
1243 # Log a message if it has been a while since the last log.
1244 periodicLog.log(
1245 "Deblending sub-blend %d/%d of parent %d",
1246 subBlendIndex + 1,
1247 len(blendRecords),
1248 parentRecord.getId(),
1249 )
1250 try:
1251 blend, blendModel, chi2 = self._deblendParent(blendRecord)
1252 except DeblenderSkippedError as e:
1253 self._skipBlend(blendRecord, e.skipKey, e.message)
1254 parentRecord.set("deblend_skipped", True)
1255 parentRecord.set(e.skipKey, True)
1256 continue
1257 except DeblenderError as e:
1258 blendRecord.set("deblend_error", e.errorName)
1259 blendRecord.set("deblend_failed", True)
1260 self._skipBlend(blendRecord, "deblend_failed", e.message)
1261 parentRecord.set("deblend_childFailed", True)
1262 continue
1264 # Update the parent model
1265 parentModel.insert(blendModel)
1267 # Add each deblended source to the catalog
1268 for scarletSource in blend.sources:
1269 # Add all fields except the HeavyFootprint to the
1270 # source record
1271 scarletSource.metadata = {}
1272 scarletSource.metadata["peak_id"] = scarletSource.detectedPeak.getId()
1273 sourceRecord = self._addDeblendedSource(
1274 parentId=parentRecord.getId(),
1275 blendRecord=blendRecord,
1276 peak=scarletSource.detectedPeak,
1277 objectCatalog=objectCatalog,
1278 scarletSource=scarletSource,
1279 chi2=chi2,
1280 )
1281 scarletSource.metadata["id"] = sourceRecord.getId()
1282 sourceRecords.append(sourceRecord)
1284 # Store the blend information so that it can be persisted
1285 blendData = blend.to_data()
1286 parentBlends[blendRecord.getId()] = blendData
1288 if len(parentBlends) == 0:
1289 # All of the deconvolved blends failed to deblend
1290 self._updateParentRecord(
1291 parentRecord=parentRecord,
1292 nPeaks=len(parentRecord.getFootprint().peaks),
1293 nChild=0,
1294 nComponents=0,
1295 runtime=np.nan,
1296 iterations=0,
1297 logL=np.nan,
1298 chi2=np.nan,
1299 spectrumInit=False,
1300 # No blend was fit, so convergence does not apply.
1301 convergenceFailed=False,
1302 )
1303 continue
1305 # Calculate the reduced chi2 for the PSF parent
1306 parentFootprintImage = parentModel.data > 0
1307 chi2 = utils.calcChi2(parentModel, context.observation, parentFootprintImage)
1308 # Defensive: a detected parent's aggregate model has
1309 # positive flux somewhere, so the area should never be 0.
1310 parentArea = np.sum(parentFootprintImage)
1311 parentReducedChi2 = (
1312 np.sum(chi2.data) / parentArea if parentArea > 0 else np.nan
1313 )
1315 # Update the parent record with the deblending results
1316 self._updateParentRecord(
1317 parentRecord=parentRecord,
1318 nPeaks=len(parentFootprint.peaks),
1319 nChild=np.sum([child["deblend_nChild"] for child in blendRecords]),
1320 nComponents=np.sum([child["deblend_nComponents"] for child in blendRecords]),
1321 runtime=np.sum([child["deblend_runtime"] for child in blendRecords]),
1322 iterations=np.sum([child["deblend_iterations"] for child in blendRecords]),
1323 logL=np.nan,
1324 chi2=parentReducedChi2,
1325 spectrumInit=np.all([
1326 child["deblend_spectrumInitFlag"]
1327 for child in blendRecords
1328 ]), # type: ignore
1329 convergenceFailed=np.any([
1330 child["deblend_blendConvergenceFailedFlag"]
1331 for child in blendRecords
1332 ]), # type: ignore
1333 )
1334 # Persist parent columns to the children
1335 for child in sourceRecords:
1336 for key in self.toCopyFromParent: 1336 ↛ 1337line 1336 didn't jump to line 1337 because the loop on line 1336 never started
1337 child.set(key, parentRecord.get(key))
1338 for parentColumn, childColumn in self.config.columnInheritance.items():
1339 child.set(childColumn, parentRecord.get(parentColumn))
1341 # Persist the blend data
1342 modelData.blends[parentRecord.getId()] = io.LsstHierarchicalBlendData(
1343 children=parentBlends,
1344 span_array=parentFootprint.getSpans().asArray(),
1345 origin=(y0, x0),
1346 legacy_spans=False,
1347 )
1349 nDeblendedSources = np.sum(objectCatalog["parent"] != 0)
1350 self.log.info(
1351 "Deblender results: %d parent sources were "
1352 "split into %d deconvolved parents,"
1353 "resulting in %d deblended sources, "
1354 "for a total catalog size of %d sources",
1355 nParents,
1356 nBlends,
1357 nDeblendedSources,
1358 len(objectCatalog),
1359 )
1361 table = afwTable.SourceTable.make(self.objectSchema)
1362 sortedCatalog = afwTable.SourceCatalog(table)
1363 sortedCatalog.extend(objectCatalog, deep=True)
1364 table = afwTable.SourceTable.make(self.parentSchema)
1365 sortedBlendCatalog = afwTable.SourceCatalog(table)
1366 sortedBlendCatalog.extend(parentCatalog, deep=True)
1368 return pipeBase.Struct(
1369 deblendedCatalog=sortedCatalog,
1370 scarletModelData=modelData,
1371 objectParents=sortedBlendCatalog,
1372 )
1374 def _deblendParent(
1375 self,
1376 blendRecord: afwTable.SourceRecord,
1377 ) -> tuple[scl.Blend, scl.Image, scl.Image]:
1378 """Deblend a parent source record
1380 Parameters
1381 ----------
1382 blendRecord :
1383 The parent source record to deblend.
1385 Returns
1386 -------
1387 blend :
1388 The `scl.Blend` object that contains the deblended sources.
1389 blendModel :
1390 The `scl.Image` model of the blend.
1391 chi2 :
1392 The reduced chi2 of the blend model.
1393 """
1394 footprint = blendRecord.getFootprint()
1395 bbox = footprint.getBBox()
1396 peaks = footprint.getPeaks()
1398 # Skip the source if it meets the skipping criteria
1399 isSkipped = self._checkSkipped(blendRecord, self.mExposure)
1400 if isSkipped is not None:
1401 skipKey, skipMessage = isSkipped
1402 raise DeblenderSkippedError(
1403 skipMessage,
1404 blendRecord.getId(),
1405 skipKey,
1406 )
1408 self.log.trace(
1409 "Blend %d: deblending {%d} peaks",
1410 blendRecord.getId(),
1411 len(peaks),
1412 )
1413 # Choose whether or not to use improved spectral initialization.
1414 # This significantly cuts down on the number of iterations
1415 # that the optimizer needs and usually results in a better
1416 # fit.
1417 # But using least squares on a very large blend causes memory
1418 # issues, so it is not done for large blends
1419 if self.config.setSpectra: 1419 ↛ 1427line 1419 didn't jump to line 1427 because the condition on line 1419 was always true
1420 if self.config.maxSpectrumCutoff <= 0: 1420 ↛ 1421line 1420 didn't jump to line 1421 because the condition on line 1420 was never true
1421 spectrumInit = True
1422 else:
1423 spectrumInit = (
1424 len(footprint.peaks) * bbox.getArea() < self.config.maxSpectrumCutoff
1425 )
1426 else:
1427 spectrumInit = False
1429 try:
1430 t0 = time.monotonic()
1431 # Build the parameter lists with the same ordering
1432 blend = deblend(self.context, footprint, self.config, spectrumInit)
1433 tf = time.monotonic()
1434 runtime = (tf - t0) * 1000
1435 convergenceFailed = not _checkBlendConvergence(
1436 blend, self.config.relativeError
1437 )
1438 # Store the number of components in the blend
1439 nComponents = len(blend.components)
1440 nChild = len(blend.sources)
1441 # Catch all errors and filter out the ones that we know about
1442 except Exception as e:
1443 blendError = type(e).__name__
1444 if self.config.catchFailures:
1445 # Make it easy to find UnknownErrors in the log file
1446 self.log.warning("UnknownError")
1447 traceback.print_exc()
1448 else:
1449 raise
1451 raise DeblenderError(
1452 f"Unable to deblend parent {blendRecord.getId()}: {blendError}",
1453 blendRecord.getId(),
1454 blendError,
1455 )
1457 # Calculate the reduced chi2
1458 blendModel = blend.get_model(convolve=False)
1459 blendFootprintImage = blendModel.data > 0
1460 chi2 = utils.calcChi2(blendModel, self.context.observation, blendFootprintImage)
1461 # Defensive: same unreachable-in-practice guard as the
1462 # aggregate-parent branch above.
1463 blendArea = np.sum(blendFootprintImage)
1464 blendReducedChi2 = (
1465 np.sum(chi2.data) / blendArea if blendArea > 0 else np.nan
1466 )
1468 # Update the blend record with the deblending results
1469 self._updateParentRecord(
1470 parentRecord=blendRecord,
1471 nPeaks=len(peaks),
1472 nChild=nChild,
1473 nComponents=nComponents,
1474 runtime=runtime,
1475 iterations=len(blend.loss),
1476 logL=blend.log_likelihood,
1477 chi2=blendReducedChi2,
1478 spectrumInit=spectrumInit,
1479 convergenceFailed=convergenceFailed,
1480 )
1482 return blend, blendModel, chi2
1484 def _isLargeFootprint(self, footprint: afwDet.Footprint) -> bool:
1485 """Returns whether a Footprint is large
1487 'Large' is defined by thresholds on the area, size and axis ratio,
1488 and total area of the bounding box multiplied by
1489 the number of children.
1490 These may be disabled independently by configuring them to be
1491 non-positive.
1492 """
1493 if (
1494 self.config.maxFootprintArea > 0
1495 and footprint.getBBox().getArea() > self.config.maxFootprintArea
1496 ):
1497 return True
1498 if self.config.maxFootprintSize > 0: 1498 ↛ 1499line 1498 didn't jump to line 1499 because the condition on line 1498 was never true
1499 bbox = footprint.getBBox()
1500 if max(bbox.getWidth(), bbox.getHeight()) > self.config.maxFootprintSize:
1501 return True
1502 if self.config.minFootprintAxisRatio > 0: 1502 ↛ 1503line 1502 didn't jump to line 1503 because the condition on line 1502 was never true
1503 axes = afwEll.Axes(footprint.getShape())
1504 if axes.getB() < self.config.minFootprintAxisRatio * axes.getA():
1505 return True
1506 if self.config.maxAreaTimesPeaks > 0: 1506 ↛ 1512line 1506 didn't jump to line 1512 because the condition on line 1506 was always true
1507 if ( 1507 ↛ 1511line 1507 didn't jump to line 1511 because the condition on line 1507 was never true
1508 footprint.getBBox().getArea() * len(footprint.peaks)
1509 > self.config.maxAreaTimesPeaks
1510 ):
1511 return True
1512 return False
1514 def _isMasked(self, footprint: afwDet.Footprint, mExposure: afwImage.MultibandExposure) -> bool:
1515 """Returns whether the footprint violates the mask limits
1517 Parameters
1518 ----------
1519 footprint :
1520 The footprint to check for masked pixels
1521 mExposure :
1522 The multiband exposure to check for masked pixels.
1524 Returns
1525 -------
1526 isMasked : `bool`
1527 `True` if `self.config.maskPlaneLimits` is less than the
1528 fraction of pixels for a given mask in
1529 `self.config.maskLimits`.
1530 """
1531 bbox = footprint.getBBox()
1532 # AND across bands: a pixel counts as masked only when the bit
1533 # is set in every band. Matches ``buildObservation``'s per-band
1534 # weight zeroing — a pixel is truly unconstrained only when
1535 # masked in all bands; otherwise the unmasked bands still
1536 # contribute.
1537 mask = np.bitwise_and.reduce(mExposure.mask[:, bbox].array, axis=0)
1538 size = float(footprint.getArea())
1539 for maskName, limit in self.config.maskLimits.items():
1540 maskVal = mExposure.mask.getPlaneBitMask(maskName)
1541 _mask = afwImage.MaskX(mask & maskVal, xy0=bbox.getMin())
1542 # spanset of masked pixels
1543 maskedSpan = footprint.spans.intersect(_mask, maskVal)
1544 if (maskedSpan.getArea()) / size > limit:
1545 return True
1546 return False
1548 def _skipBlend(
1549 self,
1550 blendRecord: afwTable.SourceRecord,
1551 skipKey: str,
1552 logMessage: str | None,
1553 ):
1554 """Update a parent record that is not being deblended.
1556 This is a fairly trivial function but is implemented to ensure
1557 that a skipped blend updates the appropriate columns
1558 consistently, and always has a flag to mark the reason that
1559 it is being skipped.
1561 Parameters
1562 ----------
1563 blendRecord :
1564 The blend record to flag as skipped.
1565 skipKey :
1566 The name of the flag to mark the reason for skipping.
1567 logMessage :
1568 The message to display in a log.trace when a source
1569 is skipped.
1570 """
1571 if logMessage: 1571 ↛ 1573line 1571 didn't jump to line 1573 because the condition on line 1571 was always true
1572 self.log.trace(logMessage)
1573 footprint = blendRecord.getFootprint()
1574 self._updateParentRecord(
1575 parentRecord=blendRecord,
1576 nPeaks=len(footprint.peaks),
1577 nChild=0,
1578 nComponents=0,
1579 runtime=np.nan,
1580 iterations=0,
1581 logL=np.nan,
1582 chi2=np.nan,
1583 spectrumInit=False,
1584 # A skipped blend was never fit, so convergence does not
1585 # apply; leave the convergence-failure flag unset.
1586 convergenceFailed=False,
1587 )
1589 # Mark the source as skipped by the deblender and
1590 # flag the reason why.
1591 blendRecord.set("deblend_skipped", True)
1592 blendRecord.set(skipKey, True)
1594 # Add the NOT_DEBLENDED mask to the mask plane in each band
1595 if self.config.notDeblendedMask: 1595 ↛ exitline 1595 didn't return from function '_skipBlend' because the condition on line 1595 was always true
1596 for mask in self.mExposure.mask:
1597 footprint.spans.setMask(
1598 mask, mask.getPlaneBitMask(self.config.notDeblendedMask)
1599 )
1601 def _checkSkipped(
1602 self,
1603 parent: afwTable.SourceRecord,
1604 mExposure: afwImage.MultibandExposure
1605 ) -> tuple[afwTable.Key, str] | None:
1606 """Update a parent record that is not being deblended.
1608 This is a fairly trivial function but is implemented to ensure
1609 that a skipped parent updates the appropriate columns
1610 consistently, and always has a flag to mark the reason that
1611 it is being skipped.
1613 Parameters
1614 ----------
1615 parent :
1616 The parent record to flag as skipped.
1617 mExposure :
1618 The exposures should be co-added images of the same
1619 shape and region of the sky.
1621 Returns
1622 -------
1623 skip :
1624 `True` if the deblender will skip the parent
1625 """
1626 skipKey = None
1627 skipMessage = ""
1628 footprint = parent.getFootprint()
1629 if isPseudoSource(parent, self.config.pseudoColumns): 1629 ↛ 1632line 1629 didn't jump to line 1632 because the condition on line 1629 was never true
1630 # We also skip pseudo sources, like sky objects, which
1631 # are intended to be skipped.
1632 skipKey = "deblend_skipped_isPseudo"
1633 elif self._isLargeFootprint(footprint):
1634 # The footprint is above the maximum footprint size limit
1635 skipKey = "deblend_skipped_parentTooBig"
1636 skipMessage = f"Parent {parent.getId()}: skipping large footprint"
1637 elif self._isMasked(footprint, mExposure): 1637 ↛ 1639line 1637 didn't jump to line 1639 because the condition on line 1637 was never true
1638 # The footprint exceeds the maximum number of masked pixels
1639 skipKey = "deblend_skipped_masked"
1640 skipMessage = f"Parent {parent.getId()}: skipping masked footprint"
1641 elif (
1642 self.config.maxNumberOfPeaks > 0
1643 and len(footprint.peaks) > self.config.maxNumberOfPeaks
1644 ):
1645 # Unlike meas_deblender, in scarlet we skip the entire blend
1646 # if the number of peaks exceeds max peaks, since neglecting
1647 # to model any peaks often results in catastrophic failure
1648 # of scarlet to generate models for the brighter sources.
1649 skipKey = "deblend_skipped_tooManyPeaks"
1650 skipMessage = f"Parent {parent.getId()}: skipping blend with too many peaks"
1651 if skipKey is not None:
1652 return (cast(afwTable.Key, skipKey), skipMessage)
1653 return None
1655 def _updateParentRecord(
1656 self,
1657 parentRecord: afwTable.SourceRecord,
1658 nPeaks: int,
1659 nChild: int,
1660 nComponents: int,
1661 runtime: float,
1662 iterations: int,
1663 logL: float,
1664 chi2: float,
1665 spectrumInit: bool,
1666 convergenceFailed: bool,
1667 ):
1668 """Update a parent record in all of the single band catalogs.
1670 Ensure that all locations that update a blend record,
1671 whether it is skipped or updated after deblending,
1672 update all of the appropriate columns.
1674 Parameters
1675 ----------
1676 blendRecord :
1677 The catalog record to update.
1678 nPeaks :
1679 Number of peaks in the parent footprint.
1680 nChild :
1681 Number of children deblended from the parent.
1682 This may differ from `nPeaks` if some of the peaks
1683 were culled and have no deblended model.
1684 nComponents :
1685 Total number of components in the parent.
1686 This is usually different than the number of children,
1687 since it is common for a single source to have multiple
1688 components.
1689 runtime :
1690 Total runtime for deblending.
1691 iterations :
1692 Total number of iterations in scarlet before convergence.
1693 logL :
1694 Final log likelihood of the blend.
1695 chi2 :
1696 Final reduced chi2 of the blend.
1697 spectrumInit :
1698 True when scarlet used `set_spectra` to initialize all
1699 sources with better initial intensities.
1700 convergenceFailed :
1701 True when the blend was fit but the optimizer reached the
1702 maximum number of iterations without converging. False both
1703 when the blend converged and when no fit was attempted
1704 (skipped or isolated parents), so the flag never fires for
1705 blends where convergence does not apply.
1706 """
1707 parentRecord.set("deblend_nPeaks", nPeaks)
1708 parentRecord.set("deblend_nChild", nChild)
1709 parentRecord.set("deblend_nComponents", nComponents)
1710 parentRecord.set("deblend_runtime", runtime)
1711 parentRecord.set("deblend_iterations", iterations)
1712 parentRecord.set("deblend_spectrumInitFlag", spectrumInit)
1713 parentRecord.set(
1714 "deblend_blendConvergenceFailedFlag", convergenceFailed
1715 )
1716 parentRecord.set("deblend_chi2", chi2)
1718 def _initializeCatalogs(
1719 self,
1720 objectCatalog: afwTable.SourceCatalog,
1721 parentCatalog: afwTable.SourceCatalog,
1722 context: ScarletDeblendContext,
1723 ) -> None:
1724 """Create footprints for deconvolved parents
1726 Each parent may be subdivided into multiple blends that are
1727 isolated in deconvolved space but still blended in the image.
1728 This method finds all of the deconvolved footprints that overlap
1729 with a single parent footprint from the input catalog and
1730 returns a dictionary to map the parent ids to a list of
1731 deconvolved footprints.
1733 Parameters
1734 ----------
1735 objectCatalog :
1736 The catalog of isolated and deblended sources.
1737 parentCatalog :
1738 The catalog of parents and deconvolved parents used for deblending.
1739 context :
1740 The context for the entire coadd.
1741 """
1742 nParents = len(parentCatalog)
1744 isolated = []
1745 for n in range(nParents):
1746 parent = parentCatalog[n]
1747 parentFoot = parent.getFootprint()
1748 # Since we use the first peak for the parent object, we should
1749 # propagate its flags to the parent source.
1750 # For example, this propagates `merge_peak_sky` to the parent
1751 parent.assign(parent.getFootprint().peaks[0], self.parentPeakSchemaMapper)
1753 if isPseudoSource(parent, self.config.pseudoColumns): 1753 ↛ 1758line 1753 didn't jump to line 1758 because the condition on line 1753 was never true
1754 # Skip pseudo sources
1755 # Note: this does not flag isolated sources as skipped or
1756 # set the NOT_DEBLENDED mask in the exposure,
1757 # since these aren't really any skipped blends.
1758 isolated.append(parent)
1759 continue
1761 if len(parentFoot.peaks) == 1:
1762 # Isolated source and we are not processing singles
1763 isolated.append(parent)
1764 parent.set("deblend_skipped_isolatedParent", True)
1765 continue
1767 # Find deconvolved footprints that intersect with the parent
1768 # and add them to the blend catalog.
1769 parentId = parent.getId()
1770 self._buildIntersectingFootprints(
1771 parentId,
1772 parentFoot,
1773 parentCatalog,
1774 context.footprints,
1775 context.footprintImage
1776 )
1778 # Update the IdFactory to account for all of the parents
1779 # and deconvolved parents added to the catalog.
1780 idFactory = objectCatalog.getIdFactory()
1781 maxId = np.max(parentCatalog["id"])
1782 idFactory.notify(maxId)
1784 # Add the isolated sources to the object catalog.
1785 # This has to be done before adding the other deblended children
1786 # because a SourceCatalog must be ordered by parentId and all isolated
1787 # sources have parentId == 0.
1788 for parent in isolated:
1789 self._addIsolatedSource(parent, objectCatalog)
1791 def _buildIntersectingFootprints(
1792 self,
1793 parentId: int,
1794 afwFootprint: afwDet.Footprint,
1795 parentCatalog: afwTable.SourceCatalog,
1796 sclFootprints: list[scl.detect.Footprint],
1797 footprintImage: scl.Image,
1798 ) -> None:
1799 """Get the intersection of two footprints
1801 Parameters
1802 ----------
1803 parentId :
1804 The parent id containing the footprints.
1805 afwFootprint :
1806 The afw footprint
1807 parentCatalog :
1808 The catalog of parents and deconvolved parents.
1809 sclFootprints :
1810 List of scarlet lite Footprints.
1811 footprintImage :
1812 An indexed image of the scarlet footprints so that the value
1813 of a pixel gives the index + 1 of the footprints that
1814 contain that pixel.
1816 Returns
1817 -------
1818 intersection :
1819 The intersection of the two footprints
1820 """
1821 footprintIndices = set()
1822 ymin, xmin = footprintImage.bbox.origin
1824 # Get the index of the deconvolved footprint at the peak location
1825 height, width = footprintImage.data.shape
1826 for peak in afwFootprint.peaks:
1827 x = peak["i_x"] - xmin
1828 y = peak["i_y"] - ymin
1829 # NumPy wraps negative indices, so a try/except on
1830 # IndexError would silently accept peaks west/south of
1831 # the bbox origin. Bounds-check both directions
1832 # explicitly.
1833 if x < 0 or y < 0 or x >= width or y >= height:
1834 raise RuntimeError(f"no footprint at ({y}, {x})")
1835 footprintIndex = footprintImage.data[y, x] - 1
1836 if footprintIndex >= 0: 1836 ↛ 1826line 1836 didn't jump to line 1826 because the condition on line 1836 was always true
1837 footprintIndices.add(footprintIndex)
1839 # Get the intersection of each deconvolved footprint with
1840 # the parent footprint.
1841 for index in footprintIndices:
1842 _sclFootprint = scarletFootprintToAfw(sclFootprints[index])
1843 intersection = afwFootprint.intersect(_sclFootprint, copyPeaks=True)
1844 if len(intersection.peaks) > 0: 1844 ↛ 1841line 1844 didn't jump to line 1841 because the condition on line 1844 was always true
1845 self._addBlendRecord(
1846 parentId=parentId,
1847 parentCatalog=parentCatalog,
1848 footprint=intersection,
1849 )
1851 def _addBlendRecord(
1852 self,
1853 parentId: int,
1854 parentCatalog: afwTable.SourceCatalog,
1855 footprint: afwDet.Footprint,
1856 ) -> None:
1857 """Add deconvolved parents to the parent catalog
1859 Each parent may be subdivided into multiple blends that are
1860 isolated in deconvolved space but still blended in the image.
1861 This function adds the sub-parents to the catalog.
1863 Parameters
1864 ----------
1865 parentId :
1866 The ID of the parent of the sub-parents.
1867 parentCatalog :
1868 The catalog of parents and deconvolved parents.
1869 footprint :
1870 The footprint of the deconvolved parent.
1871 """
1872 blendRecord = parentCatalog.addNew()
1873 blendRecord.setParent(parentId)
1874 blendRecord.setFootprint(footprint)
1875 # Propagate the first peak's schema fields onto the
1876 # sub-blend record, mirroring what ``_initializeCatalogs``
1877 # does for top-level parents.
1878 blendRecord.assign(footprint.peaks[0], self.parentPeakSchemaMapper)
1880 def _addDeblendedSource(
1881 self,
1882 parentId: int,
1883 blendRecord: afwTable.SourceRecord,
1884 peak: afwDet.PeakRecord,
1885 objectCatalog: afwTable.SourceCatalog,
1886 scarletSource: scl.Source,
1887 chi2: scl.Image,
1888 ):
1889 """Add a deblended source to a catalog.
1891 This creates a new child in the source catalog,
1892 assigning it a parent id, and adding all columns
1893 that are independent across all filter bands and not
1894 updated after deblending.
1896 Parameters
1897 ----------
1898 parent :
1899 The parent of the new child record.
1900 blendRecord :
1901 The deconvolved parent of the new child record.
1902 peak :
1903 The peak record for the peak from the parent peak catalog.
1904 catalog :
1905 The source catalog that the child is added to.
1906 scarletSource :
1907 The scarlet model for the new source record.
1908 chi2 :
1909 The chi2 for each pixel.
1911 Returns
1912 -------
1913 src :
1914 The new child source record.
1915 """
1916 src = objectCatalog.addNew()
1917 # The peak catalog is the same for all bands,
1918 # so we just use the first peak catalog
1919 src.assign(peak, self.peakSchemaMapper)
1920 src.setParent(parentId)
1921 src.set("deblend_blendId", blendRecord.getId())
1922 src.set("deblend_blendNChild", len(blendRecord.getFootprint().peaks))
1923 src.set("deblend_nChild", 0)
1925 # Set the position of the peak from the parent footprint
1926 # This will make it easier to match the same source across
1927 # deblenders and across observations, where the peak
1928 # position is unlikely to change unless enough time passes
1929 # for a source to move on the sky.
1930 src.set("deblend_peak_center_x", peak["i_x"])
1931 src.set("deblend_peak_center_y", peak["i_y"])
1932 src.set("deblend_peakId", peak["id"])
1934 # Store the number of components for the source
1935 src.set("deblend_nComponents", len(scarletSource.components))
1937 # Calculate the reduced chi2 for the source. The ``chi2``
1938 # image is the blend's, so within this source's bbox it
1939 # carries contributions wherever the *combined* blend model
1940 # was positive — including pixels where neighboring sources'
1941 # models extend into this bbox. Mask by this source's own
1942 # positive-model footprint so the chi2 sum and the area
1943 # normalization are over the same pixel set.
1944 # Defensive: a detected peak's source model has positive flux,
1945 # so ``area`` should never be 0.
1946 sourceFootprint = scarletSource.get_model().data > 0
1947 area = np.sum(sourceFootprint)
1948 sourceReducedChi2 = (
1949 np.sum(chi2[:, scarletSource.bbox].data * sourceFootprint) / area
1950 if area > 0 else np.nan
1951 )
1952 src.set("deblend_chi2", sourceReducedChi2)
1954 return src
1956 def _addIsolatedSource(
1957 self,
1958 parentRecord: afwTable.SourceRecord,
1959 objectCatalog: afwTable.SourceCatalog,
1960 ) -> afwTable.SourceRecord:
1961 """Add an isolated source to the catalog.
1963 This creates a new record in the object catalog,
1964 and assigns values for all fields not dependent on deblending.
1966 Parameters
1967 ----------
1968 parentRecord :
1969 The parent record to copy values from.
1970 objectCatalog :
1971 The source catalog that the child is added to.
1973 Returns
1974 -------
1975 src :
1976 The new isolated source record.
1977 """
1978 # The isolated source has the same peak as the parent
1979 footprint = parentRecord.getFootprint()
1980 peak = footprint.peaks[0]
1981 src = objectCatalog.addNew()
1982 src.assign(peak, self.peakSchemaMapper)
1983 src.setParent(0)
1984 # The id of the isolated source is the same as the parent
1985 src.setId(parentRecord.getId())
1986 src.setFootprint(footprint)
1987 src.set("deblend_blendId", 0)
1988 src.set("deblend_blendNChild", 0)
1989 src.set("deblend_peak_center_x", peak["i_x"])
1990 src.set("deblend_peak_center_y", peak["i_y"])
1991 src.set("deblend_peakId", peak["id"])
1993 # Isolated sources are not modeled, so they don't have components
1994 src.set("deblend_nComponents", 0)
1996 # There are no children so we need to update the parent
1997 # record and add the isolated source to the object catalog
1998 self._updateParentRecord(
1999 parentRecord=parentRecord,
2000 nPeaks=len(parentRecord.getFootprint().peaks),
2001 nChild=0,
2002 nComponents=0,
2003 runtime=0.0,
2004 iterations=0,
2005 logL=np.nan,
2006 chi2=np.nan,
2007 spectrumInit=False,
2008 # An isolated source is not fit, so convergence does not apply.
2009 convergenceFailed=False,
2010 )
2012 # Persist parent columns to the isolated source
2013 # This is usually detection flags, like merge_peak_sky
2014 # or merge_footprint_r, etc.
2015 for key in self.toCopyFromParent: 2015 ↛ 2016line 2015 didn't jump to line 2016 because the loop on line 2015 never started
2016 src.set(key, parentRecord.get(key))
2018 return src