Coverage for python/lsst/meas/extensions/scarlet/scarletDeblendTask.py: 19%
472 statements
« prev ^ index » next coverage.py v7.10.6, created at 2025-10-11 01:54 -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
28from typing import cast
30import lsst.afw.detection as afwDet
31import lsst.afw.geom.ellipses as afwEll
32import lsst.afw.image as afwImage
33import lsst.afw.table as afwTable
34import lsst.pex.config as pexConfig
35import lsst.pipe.base as pipeBase
36import lsst.scarlet.lite as scl
37import numpy as np
38from lsst.scarlet.lite.detect_pybind11 import get_footprints
39from scipy import ndimage
40from lsst.utils.logging import PeriodicLogger
41from lsst.utils.timer import timeMethod
43from . import utils
44from .footprint import scarletFootprintToAfw
46# Scarlet and proxmin have a different definition of log levels than the stack,
47# so even "warnings" occur far more often than we would like.
48# So for now we only display scarlet and proxmin errors, as all other
49# scarlet outputs would be considered "TRACE" by our standards.
50scarletLogger = logging.getLogger("scarlet")
51scarletLogger.setLevel(logging.ERROR)
52proxminLogger = logging.getLogger("proxmin")
53proxminLogger.setLevel(logging.ERROR)
55__all__ = ["deblend", "ScarletDeblendContext", "ScarletDeblendConfig", "ScarletDeblendTask"]
57logger = logging.getLogger(__name__)
60class DeblenderError(Exception):
61 """Exception raised when the deblender fails.
63 This is used to catch errors in the deblender and set the appropriate flags
64 on the parent source.
65 """
67 def __init__(
68 self,
69 message: str,
70 parent: afwTable.source.SourceRecord,
71 errorName: str,
72 ):
73 super().__init__(message)
74 self.message = message
75 self.parent = parent
76 self.errorName = errorName
78 def __str__(self) -> str:
79 return f"DeblenderError: {self.args[0]} (parent: {self.parent})"
82class DeblenderSkippedError(Exception):
83 """Exception raised when a blend is skipped.
85 This is used to catch cases where the deblender does not process
86 a deconvolved parent because it is skipped for some reason.
87 """
88 def __init__(self, message: str, parent: afwTable.source.SourceRecord, skipKey):
89 super().__init__(message)
90 self.message = message
91 self.parent = parent
92 self.skipKey = skipKey
94 def __str__(self) -> str:
95 return f"DeblenderSkippedError: {self.args[0]} (parent: {self.parent}, skipKey: {self.skipKey})"
98def _checkBlendConvergence(blend: scl.Blend, f_rel: float) -> bool:
99 """Check whether or not a blend has converged"""
100 deltaLoss = np.abs(blend.loss[-2] - blend.loss[-1])
101 convergence = f_rel * np.abs(blend.loss[-1])
102 return deltaLoss < convergence
105def isPseudoSource(source: afwTable.source.SourceRecord, pseudoColumns: list[str]) -> bool:
106 """Check if a source is a pseudo source.
108 This is mostly for skipping sky objects,
109 but any other column can also be added to disable
110 deblending on a parent or individual source when
111 set to `True`.
113 Parameters
114 ----------
115 source :
116 The source to check for the pseudo bit.
117 pseudoColumns :
118 A list of columns to check for pseudo sources.
119 """
120 isPseudo = False
121 for col in pseudoColumns:
122 try:
123 isPseudo |= source[col]
124 except KeyError:
125 pass
126 return isPseudo
129def _getDeconvolvedFootprints(
130 mDeconvolved: afwImage.MultibandExposure,
131 sources: afwTable.SourceCatalog,
132 config: ScarletDeblendConfig,
133) -> tuple[list[scl.detect.Footprint], scl.Image]:
134 """Detect footprints in the deconvolved image
136 Parameters
137 ----------
138 mDeconvolved :
139 The deconvolved multiband exposure to detect footprints in.
140 sources :
141 The source catalog for the entire coadd.
142 config :
143 The configuration for the deblender.
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 sigma = np.median(np.sqrt(mDeconvolved.variance.array), axis=(1, 2))
155 detect = np.sum(mDeconvolved.image.array/sigma[:, None, None], axis=0)
157 # We don't use the variance here because testing in DM-47738
158 # has shown that we get better results without it
159 # (we're detecting footprints without peaks in the deconvolved,
160 # noise reduced image).
161 footprints = scl.detect.detect_footprints(
162 images=np.array([detect]),
163 variance=np.ones((1, detect.shape[0], detect.shape[1]), dtype=detect.dtype),
164 scales=1,
165 min_area=config.minDeconvolvedArea,
166 footprint_thresh=config.footprintSNRThresh,
167 find_peaks=False,
168 origin=(ymin, xmin)
169 )
171 # Create an indexed image of the footprints so that the value of a pixel
172 # gives the index + 1 of the footprints that contain that pixel.
173 scarletBox = utils.bboxToScarletBox(bbox)
174 footprintImage = scl.detect.footprints_to_image(footprints, scarletBox)
176 # Ensure that there is a minimal footprint for all peaks
177 # in the detection catalog.
178 footprintArray = footprintImage.data > 0
179 for source in sources:
180 for peak in source.getFootprint().peaks:
181 x, y = peak.getI()
182 footprintArray[y - ymin, x - xmin] = True
184 # Grow the footprints
185 psfMinimalSize = config.growSize * 2 + 1
186 detectionArray = ndimage.binary_dilation(
187 footprintArray,
188 scl.utils.get_circle_mask(psfMinimalSize, bool)
189 )
191 # Ensure that all of the deconvolved footprints are contained within
192 # a single footprint from the input catalog.
193 # This should be unecessary, however creating entries in the output
194 # catalog will produce unexpected results if a deconvolved footprint
195 # is in more than one footprint from the source catalog or has
196 # flux outside of its parent footprint.
197 sourceImage = afwDet.footprintsToNumpy(sources, shape=detect.shape, xy0=(xmin, ymin))
198 detectionArray = detectionArray * sourceImage
200 footprints = get_footprints(
201 image=detectionArray.astype(np.float32),
202 min_separation=0,
203 min_area=1,
204 peak_thresh=0,
205 footprint_thresh=0,
206 find_peaks=False,
207 y0=ymin,
208 x0=xmin,
209 )
211 footprintImage = scl.detect.footprints_to_image(footprints, scarletBox)
213 return footprints, footprintImage
216@dataclass(kw_only=True)
217class ScarletDeblendContext:
218 """Context with parameters and config options for deblending
220 Attributes
221 ----------
222 monotonicity :
223 The monotonicity operator.
224 observation :
225 The observation for the entire coadd.
226 deconvolved :
227 The deconvolved image.
228 residual :
229 The residual image
230 (observation - deconvolved mode convolved with difference kernel).
231 footprints :
232 The footprints in the deconvolved image.
233 footprintImage :
234 An indexed image of the scarlet footprints so that the value
235 of a pixel gives the index + 1 of the footprints that
236 contain that pixel.
237 config :
238 The configuration for the deblender.
239 """
240 monotonicity: scl.operators.Monotonicity
241 observation: scl.Observation
242 deconvolved: scl.Image
243 footprints: list[scl.Footprint]
244 footprintImage: scl.Image
245 config: ScarletDeblendConfig
247 @staticmethod
248 def build(
249 mExposure: afwImage.MultibandExposure,
250 mDeconvolved: afwImage.MultibandExposure,
251 catalog: afwTable.SourceCatalog,
252 config: ScarletDeblendConfig,
253 ) -> ScarletDeblendContext:
254 """Build the context from a minimal set of inputs
256 Parameters
257 ----------
258 mExposure :
259 The multiband exposure for the entire coadd.
260 mDeconvolved :
261 The deconvolved multiband exposure for the entire coadd.
262 catalog :
263 The source catalog for the entire coadd.
264 config :
265 The configuration for the deblender.
266 """
267 # The PSF of the model in the deconvolved space
268 modelPsf = scl.utils.integrated_circular_gaussian(
269 sigma=config.modelPsfSigma
270 )
271 # Initialize the monotonicity operator with a size of 101 x 101 pixels.
272 # Note: If a component is > 101x101 in either axis then the
273 # monotonicity operator will resize itself.
274 monotonicity = scl.operators.Monotonicity((101, 101))
275 # Build the observation for the entire coadd
276 observation = utils.buildObservation(
277 modelPsf=modelPsf,
278 psfCenter=mExposure.getBBox().getCenter(),
279 mExposure=mExposure,
280 badPixelMasks=config.badMask,
281 useWeights=config.useWeights,
282 convolutionType=config.convolutionType,
283 catalog=catalog,
284 )
286 # Create the deconvolved image
287 yx0 = observation.images.yx0
288 bands = observation.images.bands
289 deconvolved = scl.Image(mDeconvolved.image.array, bands=mDeconvolved.bands, yx0=yx0)
290 if len(bands) == 1:
291 deconvolved = deconvolved[bands[0]:]
292 else:
293 deconvolved = deconvolved[bands]
295 # Detect footprints in the deconvolved image
296 footprints, footprintImage = _getDeconvolvedFootprints(
297 mDeconvolved=mDeconvolved,
298 sources=catalog,
299 config=config,
300 )
302 return ScarletDeblendContext(
303 monotonicity=monotonicity,
304 observation=observation,
305 deconvolved=deconvolved,
306 footprints=footprints,
307 footprintImage=footprintImage,
308 config=config,
309 )
312def deblend(
313 context: ScarletDeblendContext,
314 footprint: afwDet.Footprint,
315 config: ScarletDeblendConfig,
316 spectrumInit: bool = True,
317) -> scl.Blend:
318 """Deblend a parent footprint
320 Parameters
321 ----------
322 context :
323 Context with parameters and config options for deblending
324 footprint :
325 The parent footprint to deblend
326 config :
327 The configuration for the deblender
328 spectrumInit :
329 Whether or not to initialize the sources with their best-fit spectra
331 Returns
332 -------
333 blend : `scarlet.lite.Blend`
334 The blend this is to be deblended
335 skippedSources : `list[int]`
336 Indices of sources that were skipped due to no flux.
337 This usually means that a source was a spurrious detection in one
338 band that should not have been included in the merged catalog.
339 skippedBands : `list[str]`
340 Bands that were skipped because a PSF could not be generated for them.
341 """
342 # Define the bounding boxes in lsst.geom.Box2I and lsst.scarlet.lite.Box
343 footBox = footprint.getBBox()
344 bbox = utils.bboxToScarletBox(footBox)
346 # Extract the observation that covers the footprint and make
347 # a copy so that the changes don't affect the original observation.
348 observation = context.observation[:, bbox].copy()
349 footprintData = footprint.spans.asArray()
351 # Mask the pixels outside of the footprint
352 observation.weights.data[:] *= footprintData
354 # Convert the peaks into an array
355 peaks = [
356 np.array([peak.getIy(), peak.getIx()], dtype=int)
357 for peak in footprint.peaks
358 if not isPseudoSource(peak, config.pseudoColumns)
359 ]
361 detect_image = np.sum(context.deconvolved[:, bbox].data, axis=0)
363 # Initialize the sources
364 sources = scl.initialization.FactorizedInitialization(
365 observation=observation,
366 centers=peaks,
367 detect=detect_image,
368 min_snr=config.minSNR,
369 monotonicity=context.monotonicity,
370 bg_thresh=config.backgroundThresh,
371 initial_bg_thresh=config.initialBackgroundThresh,
372 ).sources
374 blend = scl.Blend(sources, observation)
376 # Initialize each source with its best fit spectrum
377 if spectrumInit:
378 try:
379 blend.fit_spectra()
380 except Exception as e:
381 # If the spectrum initialization fails, we will just skip it
382 # and use the default spectrum.
383 logger.warning(
384 "Spectrum initialization failed with error: %s", e, exc_info=True
385 )
387 # Set the optimizer
388 if config.optimizer == "adaprox":
389 blend.parameterize(
390 partial(
391 scl.component.default_adaprox_parameterization,
392 noise_rms=observation.noise_rms / 10,
393 )
394 )
395 elif config.optimizer == "fista":
396 blend.parameterize(scl.component.default_fista_parameterization)
397 else:
398 raise ValueError("Unrecognized optimizer. Must be either 'adaprox' or 'fista'.")
400 if config.maxIter > 0:
401 blend.fit(
402 max_iter=config.maxIter,
403 e_rel=config.relativeError,
404 min_iter=config.minIter,
405 resize=config.resizeFrequency,
406 )
407 else:
408 loss = (blend.observation.images - blend.get_model(convolve=True)).data
409 blend.loss = [np.sum(loss), np.sum(loss)]
411 # Attach the peak to all of the initialized sources
412 for k, center in enumerate(peaks):
413 # This is just to make sure that there isn't a coding bug
414 if len(sources[k].components) > 0 and np.any(sources[k].center != center):
415 raise ValueError(
416 f"Misaligned center, expected {center} but got {sources[k].center}"
417 )
418 # Store the record for the peak with the appropriate source
419 sources[k].detectedPeak = footprint.peaks[k]
421 return blend
424class ScarletDeblendConfig(pexConfig.Config):
425 """MultibandDeblendConfig
427 Configuration for the multiband deblender.
428 The parameters are organized by the parameter types, which are
429 - Stopping Criteria: Used to determine if the fit has converged
430 - Position Fitting Criteria: Used to fit the positions of the peaks
431 - Constraints: Used to apply constraints to the peaks and their components
432 - Other: Parameters that don't fit into the above categories
433 """
435 # Stopping Criteria
436 minIter = pexConfig.Field[int](
437 default=5,
438 doc="Minimum number of iterations before the optimizer is allowed to stop.",
439 )
440 maxIter = pexConfig.Field[int](
441 default=300,
442 doc=("Maximum number of iterations to deblend a single parent"),
443 )
444 relativeError = pexConfig.Field[float](
445 default=1e-3,
446 doc=(
447 "Change in the loss function between iterations to exit fitter. "
448 "Typically this is `1e-3` if measurements will be made on the "
449 "flux re-distributed models and `1e-4` when making measurements "
450 "on the models themselves."
451 ),
452 )
453 resizeFrequency = pexConfig.Field[int](
454 default=3,
455 doc="Number of iterations between resizing sources.",
456 )
458 # Constraints
459 morphThresh = pexConfig.Field[float](
460 default=1,
461 doc="Fraction of background RMS a pixel must have"
462 "to be included in the initial morphology",
463 )
464 # Lite Parameters
465 # All of these parameters (except version) are only valid if version='lite'
466 version = pexConfig.ChoiceField[str](
467 default="lite",
468 allowed={
469 "lite": "LSST optimized version of scarlet for survey data from a single instrument",
470 },
471 doc="The version of scarlet to use.",
472 deprecated="This field is deprecated since the ony available `version` is `lite` "
473 "and will be removed after v29.0",
474 )
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 morphImage = pexConfig.ChoiceField[str](
484 default="chi2",
485 allowed={
486 "chi2": "Initialize sources on a chi^2 image made from all available bands",
487 },
488 doc="The type of image to use for initializing the morphology. "
489 "Must be either 'chi2' or 'wavelet'. ",
490 deprecated="This field is deprecated since testing has shown that only 'chi2' should be used "
491 "and 'wavelet' has been broken since v27.0. "
492 "This field will be removed in v29.0",
493 )
494 backgroundThresh = pexConfig.Field[float](
495 default=1.0,
496 doc="Fraction of background to use for a sparsity threshold. "
497 "This prevents sources from growing unrealistically outside "
498 "the parent footprint while still modeling flux correctly "
499 "for bright sources.",
500 )
501 initialBackgroundThresh = pexConfig.Field[float](
502 default=1.0,
503 doc="Same as `backgroundThresh` but used only for source initialization.",
504 )
505 maxProxIter = pexConfig.Field[int](
506 default=1,
507 doc="Maximum number of proximal operator iterations inside of each "
508 "iteration of the optimizer. "
509 "This config field is only used if version='lite' and optimizer='adaprox'.",
510 )
511 waveletScales = pexConfig.Field[int](
512 default=5,
513 doc="Number of wavelet scales to use for wavelet initialization. "
514 "This field is only used when `version`='lite' and `morphImage`='wavelet'.",
515 deprecated="This field is deprecated along with `morphImage` and will be removed in v29.0.",
516 )
518 # Other scarlet paremeters
519 useWeights = pexConfig.Field[bool](
520 default=True,
521 doc=(
522 "Whether or not use use inverse variance weighting."
523 "If `useWeights` is `False` then flat weights are used"
524 ),
525 )
526 modelPsfSize = pexConfig.Field[int](
527 default=11, doc="Model PSF side length in pixels"
528 )
529 modelPsfSigma = pexConfig.Field[float](
530 default=0.8, doc="Define sigma for the model frame PSF"
531 )
532 minSNR = pexConfig.Field[float](
533 default=50,
534 doc="Minimum Signal to noise to accept the source."
535 "Sources with lower flux will be initialized with the PSF but updated "
536 "like an ordinary ExtendedSource (known in scarlet as a `CompactSource`).",
537 )
538 saveTemplates = pexConfig.Field[bool](
539 default=True, doc="Whether or not to save the SEDs and templates"
540 )
541 processSingles = pexConfig.Field[bool](
542 default=True,
543 doc="Whether or not to process isolated sources in the deblender",
544 )
545 convolutionType = pexConfig.Field[str](
546 default="fft",
547 doc="Type of convolution to render the model to the observations.\n"
548 "- 'fft': perform convolutions in Fourier space\n"
549 "- 'real': peform convolutions in real space.",
550 )
551 sourceModel = pexConfig.Field[str](
552 default="double",
553 doc=(
554 "How to determine which model to use for sources, from\n"
555 "- 'single': use a single component for all sources\n"
556 "- 'double': use a bulge disk model for all sources\n"
557 "- 'compact': use a single component model, initialzed with a point source morphology, "
558 " for all sources\n"
559 "- 'point': use a point-source model for all sources\n"
560 "- 'fit: use a PSF fitting model to determine the number of components (not yet implemented)"
561 ),
562 deprecated="This field will be deprecated when the default for `version` is changed to `lite`.",
563 )
564 setSpectra = pexConfig.Field[bool](
565 default=True,
566 doc="Whether or not to solve for the best-fit spectra during initialization. "
567 "This makes initialization slightly longer, as it requires a convolution "
568 "to set the optimal spectra, but results in a much better initial log-likelihood "
569 "and reduced total runtime, with convergence in fewer iterations."
570 "This option is only used when "
571 "peaks*area < `maxSpectrumCutoff` will use the improved initialization.",
572 )
573 footprintSNRThresh = pexConfig.Field[float](
574 default=5.0,
575 doc="Minimum SNR for a pixel to be detected in a footprint.",
576 )
577 growSize = pexConfig.Field[int](
578 default=2,
579 doc="Number of pixels to grow the deconvolved footprints before final detection.",
580 )
582 # Mask-plane restrictions
583 badMask = pexConfig.ListField[str](
584 default=utils.defaultBadPixelMasks,
585 doc="Whether or not to process isolated sources in the deblender",
586 )
587 statsMask = pexConfig.ListField[str](
588 default=["SAT", "INTRP", "NO_DATA"],
589 doc="Mask planes to ignore when performing statistics",
590 )
591 maskLimits = pexConfig.DictField(
592 keytype=str,
593 itemtype=float,
594 default={},
595 doc=(
596 "Mask planes with the corresponding limit on the fraction of masked pixels. "
597 "Sources violating this limit will not be deblended. "
598 "If the fraction is `0` then the limit is a single pixel."
599 ),
600 )
601 minDeconvolvedArea = pexConfig.Field[int](
602 default=9,
603 doc="Minimum area for a single footprint in the deconvolved image. "
604 "Detected footprints smaller than this will not be created.",
605 )
607 # Size restrictions
608 maxNumberOfPeaks = pexConfig.Field[int](
609 default=600,
610 doc=(
611 "Only deblend the brightest maxNumberOfPeaks peaks in the parent"
612 " (<= 0: unlimited)"
613 ),
614 )
615 maxFootprintArea = pexConfig.Field[int](
616 default=2_000_000,
617 doc=(
618 "Maximum area for footprints before they are ignored as large; "
619 "non-positive means no threshold applied"
620 ),
621 )
622 maxAreaTimesPeaks = pexConfig.Field[int](
623 default=1_000_000_000,
624 doc=(
625 "Maximum rectangular footprint area * nPeaks in the footprint. "
626 "This was introduced in DM-33690 to prevent fields that are crowded or have a "
627 "LSB galaxy that causes memory intensive initialization in scarlet from dominating "
628 "the overall runtime and/or causing the task to run out of memory. "
629 "(<= 0: unlimited)"
630 ),
631 )
632 maxFootprintSize = pexConfig.Field[int](
633 default=0,
634 doc=(
635 "Maximum linear dimension for footprints before they are ignored "
636 "as large; non-positive means no threshold applied"
637 ),
638 )
639 minFootprintAxisRatio = pexConfig.Field[float](
640 default=0.0,
641 doc=(
642 "Minimum axis ratio for footprints before they are ignored "
643 "as large; non-positive means no threshold applied"
644 ),
645 )
646 maxSpectrumCutoff = pexConfig.Field[int](
647 default=1_000_000,
648 doc=(
649 "Maximum number of pixels * number of sources in a blend. "
650 "This is different than `maxFootprintArea` because this isn't "
651 "the footprint area but the area of the bounding box that "
652 "contains the footprint, and is also multiplied by the number of"
653 "sources in the footprint. This prevents large skinny blends with "
654 "a high density of sources from running out of memory. "
655 "If `maxSpectrumCutoff == -1` then there is no cutoff."
656 ),
657 )
658 # Failure modes
659 fallback = pexConfig.Field[bool](
660 default=True,
661 doc="Whether or not to fallback to a smaller number of components if a source does not initialize",
662 )
663 notDeblendedMask = pexConfig.Field[str](
664 default="NOT_DEBLENDED",
665 optional=True,
666 doc="Mask name for footprints not deblended, or None",
667 )
668 catchFailures = pexConfig.Field[bool](
669 default=True,
670 doc=(
671 "If True, catch exceptions thrown by the deblender, log them, "
672 "and set a flag on the parent, instead of letting them propagate up"
673 ),
674 )
676 # Other options
677 columnInheritance = pexConfig.DictField(
678 keytype=str,
679 itemtype=str,
680 default={
681 "deblend_nChild": "deblend_parentNChild",
682 "deblend_nPeaks": "deblend_parentNPeaks",
683 "deblend_spectrumInitFlag": "deblend_spectrumInitFlag",
684 "deblend_blendConvergenceFailedFlag": "deblend_blendConvergenceFailedFlag",
685 },
686 doc="Columns to pass from the parent to the child. "
687 "The key is the name of the column for the parent record, "
688 "the value is the name of the column to use for the child.",
689 )
690 pseudoColumns = pexConfig.ListField[str](
691 default=["merge_peak_sky", "sky_source"],
692 doc="Names of flags which should never be deblended.",
693 )
695 # Testing options
696 # Some obs packages and ci packages run the full pipeline on a small
697 # subset of data to test that the pipeline is functioning properly.
698 # This is not meant as scientific validation, so it can be useful
699 # to only run on a small subset of the data that is large enough to
700 # test the desired pipeline features but not so long that the deblender
701 # is the tall pole in terms of execution times.
702 useCiLimits = pexConfig.Field[bool](
703 default=False,
704 doc="Limit the number of sources deblended for CI to prevent long build times",
705 )
706 ciDeblendChildRange = pexConfig.ListField[int](
707 default=[5, 10],
708 doc="Only deblend parent Footprints with a number of peaks in the (inclusive) range indicated."
709 "If `useCiLimits==False` then this parameter is ignored.",
710 )
711 ciNumParentsToDeblend = pexConfig.Field[int](
712 default=10,
713 doc="Only use the first `ciNumParentsToDeblend` parent footprints with a total peak count "
714 "within `ciDebledChildRange`. "
715 "If `useCiLimits==False` then this parameter is ignored.",
716 )
719class ScarletDeblendTask(pipeBase.Task):
720 """ScarletDeblendTask
722 Split blended sources into individual sources.
724 This task has no return value; it only modifies the SourceCatalog in-place.
725 """
727 ConfigClass = ScarletDeblendConfig
728 _DefaultName = "scarletDeblend"
730 def __init__(
731 self,
732 schema: afwTable.Schema,
733 peakSchema: afwTable.Schema = None,
734 **kwargs
735 ):
736 """Create the task, adding necessary fields to the given schema.
738 Parameters
739 ----------
740 schema :
741 Schema object for measurement fields; will be modified in-place.
742 peakSchema :
743 Schema of Footprint Peaks that will be passed to the deblender.
744 Any fields beyond the PeakTable minimal schema will be transferred
745 to the main source Schema. If None, no fields will be transferred
746 from the Peaks.
747 **kwargs
748 Passed to Task.__init__.
749 """
750 pipeBase.Task.__init__(self, **kwargs)
752 peakMinimalSchema = afwDet.PeakTable.makeMinimalSchema()
753 if peakSchema is None:
754 # In this case, the peakSchemaMapper will transfer nothing, but
755 # we'll still have one
756 # to simplify downstream code
757 self.peakSchemaMapper = afwTable.SchemaMapper(peakMinimalSchema, schema)
758 else:
759 self.peakSchemaMapper = afwTable.SchemaMapper(peakSchema, schema)
760 for item in peakSchema:
761 if item.key not in peakMinimalSchema:
762 self.peakSchemaMapper.addMapping(item.key, item.field)
763 # Because SchemaMapper makes a copy of the output schema
764 # you give its ctor, it isn't updating this Schema in
765 # place. That's probably a design flaw, but in the
766 # meantime, we'll keep that schema in sync with the
767 # peakSchemaMapper.getOutputSchema() manually, by adding
768 # the same fields to both.
769 schema.addField(item.field)
770 assert (
771 schema == self.peakSchemaMapper.getOutputSchema()
772 ), "Logic bug mapping schemas"
774 # Add the parent keys to the parent catalog schema
775 blendSchema = afwTable.SourceTable.makeMinimalSchema()
776 self._addParentSchemaKeys(blendSchema)
777 self._addSharedSchemaKeys(blendSchema)
778 self.blendSchema = blendSchema
779 self.blendPeakSchemaMapper = afwTable.SchemaMapper(peakMinimalSchema, blendSchema)
781 # For now source records have parent schema keys and deblended
782 # source keys.
783 # DM-51670 will add the top level parents to the parent catalog
784 # and remove the parent keys from the source catalog.
785 self._addParentSchemaKeys(schema)
786 self._addChildSchemaKeys(schema)
787 self._addSharedSchemaKeys(schema)
788 self.schema = schema
789 self.toCopyFromParent = [
790 name
791 for item in self.schema
792 if (name := item.field.getName()).startswith("merge_footprint")
793 ]
795 def _addParentSchemaKeys(self, schema: afwTable.Schema):
796 """Add parent specific keys to the schema"""
797 # Parent (blend) fields
798 schema.addField(
799 "deblend_runtime", type=np.float32, doc="runtime in ms"
800 )
801 schema.addField(
802 "deblend_iterations", type=np.int32, doc="iterations to converge"
803 )
804 schema.addField(
805 "deblend_nChild",
806 type=np.int32,
807 doc="Number of children this object has (defaults to 0)",
808 )
809 schema.addField(
810 "deblend_nPeaks",
811 type=np.int32,
812 doc="Number of initial peaks in the blend. "
813 "This includes peaks that may have been culled "
814 "during deblending or failed to deblend",
815 )
816 schema.addField(
817 "deblend_spectrumInitFlag",
818 type="Flag",
819 doc="True when scarlet initializes sources "
820 "in the blend with a more accurate spectrum. "
821 "The algorithm uses a lot of memory, "
822 "so large dense blends will use "
823 "a less accurate initialization.",
824 )
825 # Skipped flags
826 schema.addField(
827 "deblend_skipped", type="Flag", doc="Deblender skipped this source"
828 )
829 schema.addField(
830 "deblend_isolatedParent",
831 type="Flag",
832 doc="The source has only a single peak " "and was not deblended",
833 )
834 schema.addField(
835 "deblend_isPseudo",
836 type="Flag",
837 doc='The source is identified as a "pseudo" source and '
838 "was not deblended",
839 )
840 schema.addField(
841 "deblend_tooManyPeaks",
842 type="Flag",
843 doc="Source had too many peaks; " "only the brightest were included",
844 )
845 schema.addField(
846 "deblend_parentTooBig",
847 type="Flag",
848 doc="Parent footprint covered too many pixels",
849 )
850 schema.addField(
851 "deblend_masked",
852 type="Flag",
853 doc="Parent footprint had too many masked pixels",
854 )
855 # Convergence flags
856 schema.addField(
857 "deblend_blendConvergenceFailedFlag",
858 type="Flag",
859 doc="at least one source in the blend" "failed to converge",
860 )
861 # Error flags
862 schema.addField(
863 "deblend_failed", type="Flag", doc="Deblending failed on source"
864 )
865 schema.addField(
866 "deblend_childFailed",
867 type="Flag",
868 doc="Deblending failed on at least one child blend. "
869 " This is set in a parent when at least one of its children "
870 "is a blend that failed to deblend.",
871 )
872 schema.addField(
873 "deblend_error",
874 type="String",
875 size=25,
876 doc="Name of error if the blend failed",
877 )
878 schema.addField(
879 "deblend_incompleteData",
880 type="Flag",
881 doc="True when a blend has at least one band "
882 "that could not generate a PSF and was "
883 "not included in the model.",
884 )
886 def _addChildSchemaKeys(self, schema: afwTable.Schema):
887 """Add deblender specific keys to the schema"""
888 afwTable.Point2IKey.addFields(
889 schema,
890 name="deblend_peak_center",
891 doc="Center used to apply constraints in scarlet",
892 unit="pixel",
893 )
894 schema.addField(
895 "deblend_peakId",
896 type=np.int32,
897 doc="ID of the peak in the parent footprint. "
898 "This is not unique, but the combination of 'parent'"
899 "and 'peakId' should be for all child sources. "
900 "Top level blends with no parents have 'peakId=0'",
901 )
902 schema.addField(
903 "deblend_peak_instFlux",
904 type=float,
905 units="count",
906 doc="The instFlux at the peak position of deblended mode",
907 )
908 schema.addField(
909 "deblend_scarletFlux", type=np.float32, doc="Flux measurement from scarlet"
910 )
911 schema.addField(
912 "deblend_edgePixels",
913 type="Flag",
914 doc="Source had flux on the edge of the parent footprint",
915 )
916 schema.addField(
917 "deblend_dataCoverage",
918 type=np.float32,
919 doc="Fraction of pixels with data. "
920 "In other words, 1 - fraction of pixels with NO_DATA set.",
921 )
922 schema.addField(
923 "deblend_zeroFlux", type="Flag", doc="Source has zero flux."
924 )
925 # Blendedness/classification metrics
926 schema.addField(
927 "deblend_maxOverlap",
928 type=np.float32,
929 doc="Maximum overlap with all of the other neighbors flux "
930 "combined."
931 "This is useful as a metric for determining how blended a "
932 "source is because if it only overlaps with other sources "
933 "at or below the noise level, it is likely to be a mostly "
934 "isolated source in the deconvolved model frame.",
935 )
936 schema.addField(
937 "deblend_fluxOverlap",
938 type=np.float32,
939 doc="This is the total flux from neighboring objects that "
940 "overlaps with this source.",
941 )
942 schema.addField(
943 "deblend_fluxOverlapFraction",
944 type=np.float32,
945 doc="This is the fraction of "
946 "`flux from neighbors/source flux` "
947 "for a given source within the source's"
948 "footprint.",
949 )
950 schema.addField(
951 "deblend_blendedness",
952 type=np.float32,
953 doc="The Bosch et al. 2018 metric for 'blendedness.' ",
954 )
955 schema.addField(
956 "deblend_blendId",
957 type=np.int64,
958 doc="Parents in the catalog may be subdivided by deblending "
959 "into multiple deconvolved blends that are each "
960 "deblended separately. This is the ID of the "
961 "deconvolved blend in the catalog."
962 )
963 schema.addField(
964 "deblend_blendNChild",
965 type=np.int32,
966 doc="The number of children in the deconvolved blend."
967 )
969 def _addSharedSchemaKeys(self, schema: afwTable.Schema):
970 """Add parent and child specific keys to the schema"""
971 # Measurement keys
972 schema.addField(
973 "deblend_logL",
974 type=np.float32,
975 doc="Final logL, used to identify regressions in scarlet.",
976 )
977 schema.addField(
978 "deblend_chi2",
979 type=np.float32,
980 doc="Final reduced chi2 (per pixel), used to identify goodness of fit.",
981 )
982 schema.addField(
983 "deblend_parentNPeaks",
984 type=np.int32,
985 doc="deblend_nPeaks from this records parent.",
986 )
987 schema.addField(
988 "deblend_parentNChild",
989 type=np.int32,
990 doc="deblend_nChild from this records parent.",
991 )
992 schema.addField(
993 "deblend_nComponents",
994 type=np.int32,
995 doc="Number of components in a ScarletLiteSource. "
996 "If `config.version != 'lite'`then "
997 "this column is set to zero.",
998 )
1000 @timeMethod
1001 def run(
1002 self,
1003 mExposure: afwImage.MultibandExposure,
1004 mDeconvolved: afwImage.MultibandExposure,
1005 mergedSources: afwTable.SourceCatalog,
1006 ) -> pipeBase.Struct:
1007 """Get the psf from each exposure and then run deblend().
1009 Parameters
1010 ----------
1011 mExposure :
1012 The exposures should be co-added images of the same
1013 shape and region of the sky.
1014 mDeconvolved :
1015 The deconvolved images of the same shape and region of the sky.
1016 mergedSources :
1017 The merged `SourceCatalog` that contains parent footprints
1018 to (potentially) deblend.
1020 Returns
1021 -------
1022 templateCatalogs: dict
1023 Keys are the names of the bands and the values are
1024 `lsst.afw.table.source.source.SourceCatalog`'s.
1025 These are catalogs with heavy footprints that are the templates
1026 created by the multiband templates.
1027 """
1028 # Create a table to hold parent record information
1029 table = afwTable.SourceTable.make(self.blendSchema)
1030 blendCatalog = afwTable.SourceCatalog(table)
1031 return self.deblend(mExposure, mDeconvolved, mergedSources, blendCatalog)
1033 @timeMethod
1034 def deblend(
1035 self,
1036 mExposure: afwImage.MultibandExposure,
1037 mDeconvolved: afwImage.MultibandExposure,
1038 catalog: afwTable.SourceCatalog,
1039 blendCatalog: afwTable.SourceCatalog,
1040 ) -> pipeBase.Struct:
1041 """Deblend a data cube of multiband images
1043 Deblending iterates over sources from the input catalog,
1044 which are blends of peaks with overlapping PSFs (depth 0 parents).
1045 In many cases those footprints can be subdived into multiple
1046 deconvolved footprints, which have an intermediate
1047 parent record added to the catalog and are be deblended separately.
1048 All deblended peaks have a source record added to the catalog,
1049 each of which has a depth one greater than the parent.
1051 Parameters
1052 ----------
1053 mExposure :
1054 The exposures should be co-added images of the same
1055 shape and region of the sky.
1056 mDeconvolved :
1057 The deconvolved images of the same shape and region of the sky.
1058 catalog :
1059 The merged `SourceCatalog` that contains parent footprints
1060 to (potentially) deblend. The new deblended sources are
1061 appended to this catalog in place.
1063 Returns
1064 -------
1065 catalog :
1066 The ``deblendedCatalog`` with parents and child sources.
1067 modelData :
1068 The persistable data model for the deblender.
1069 """
1071 # Cull footprints if required by ci
1072 if self.config.useCiLimits:
1073 self.log.info(
1074 "Using CI catalog limits, the original number of sources to deblend was %d.",
1075 len(catalog),
1076 )
1077 # Select parents with a number of children in the range
1078 # config.ciDeblendChildRange
1079 minChildren, maxChildren = self.config.ciDeblendChildRange
1080 nPeaks = np.array([len(src.getFootprint().peaks) for src in catalog])
1081 childrenInRange = np.where(
1082 (nPeaks >= minChildren) & (nPeaks <= maxChildren)
1083 )[0]
1084 if len(childrenInRange) < self.config.ciNumParentsToDeblend:
1085 raise ValueError(
1086 "Fewer than ciNumParentsToDeblend children were contained in the range "
1087 "indicated by ciDeblendChildRange. Adjust this range to include more "
1088 "parents."
1089 )
1090 # Keep all of the isolated parents and the first
1091 # `ciNumParentsToDeblend` children
1092 parents = nPeaks == 1
1093 children = np.zeros((len(catalog),), dtype=bool)
1094 children[childrenInRange[: self.config.ciNumParentsToDeblend]] = True
1095 catalog = catalog[parents | children]
1096 # We need to update the IdFactory, otherwise the the source ids
1097 # will not be sequential
1098 idFactory = catalog.getIdFactory()
1099 maxId = np.max(catalog["id"])
1100 idFactory.notify(maxId)
1101 del children
1103 self.log.info(
1104 "Deblending %d sources in %d exposure bands", len(catalog), len(mExposure),
1105 )
1106 periodicLog = PeriodicLogger(self.log)
1108 # Add the NOT_DEBLENDED mask to the mask plane in each band
1109 if self.config.notDeblendedMask:
1110 for mask in mExposure.mask:
1111 mask.addMaskPlane(self.config.notDeblendedMask)
1113 # Create the context for the entire coadd
1114 context = ScarletDeblendContext.build(
1115 mExposure=mExposure,
1116 mDeconvolved=mDeconvolved,
1117 config=self.config,
1118 catalog=catalog,
1119 )
1120 nBands = len(context.observation.bands)
1122 # Initialize the persistable ScarletModelData object
1123 modelData = scl.io.ScarletModelData(metadata={
1124 "model_psf": context.observation.model_psf[0],
1125 "psf": context.observation.psfs,
1126 "bands": context.observation.bands,
1127 })
1129 # Attach full image objects to the task to simplify the API
1130 # and use for debugging.
1131 self.catalog = catalog
1132 self.blendCatalog = blendCatalog
1133 self.context = context
1134 self.modelData = modelData
1135 self.mExposure = mExposure
1137 # Subdivide the psf blended parents into deconvolved parents
1138 # using the deconvolved footprints stored in the context.
1139 nParents = len(catalog)
1140 self._buildBlendCatalog(
1141 catalog,
1142 blendCatalog,
1143 context,
1144 )
1145 nBlends = len(blendCatalog)
1147 self.log.info(
1148 "Subdivided %d top level parents to create %d deconvolved parents.",
1149 nParents,
1150 nBlends,
1151 )
1153 # Deblend sources
1154 for parentIndex in range(nParents):
1155 # Log a message if it has been a while since the last log.
1156 periodicLog.log(
1157 "Deblended %d out of %d parents",
1158 parentIndex,
1159 nParents,
1160 )
1162 parentRecord = catalog[parentIndex]
1163 blendRecords = blendCatalog[blendCatalog["parent"] == parentRecord.getId()]
1165 if len(blendRecords) == 0:
1166 # There are no children so we must not be processing singles,
1167 # but we still need to update the parent record
1168 self._updateParentRecord(
1169 parentRecord=parentRecord,
1170 nPeaks=len(parentRecord.getFootprint().peaks),
1171 nChild=0,
1172 nComponents=0,
1173 runtime=0.0,
1174 iterations=0,
1175 logL=np.nan,
1176 chi2=np.nan,
1177 spectrumInit=False,
1178 converged=True, # No children, so no convergence issues
1179 )
1180 continue
1182 self.log.trace(
1183 "Split parent %d into %d deconvolved parents",
1184 parentRecord.getId(),
1185 len(blendRecords),
1186 )
1187 # Create an image to keep track of the cumulative model
1188 # for all sub blends in the parent footprint.
1189 parentFootprint = parentRecord.getFootprint()
1190 bbox = parentFootprint.getBBox()
1191 width, height = bbox.getDimensions()
1192 x0, y0 = bbox.getMin()
1193 emptyModel = np.zeros(
1194 (nBands, height, width),
1195 dtype=mExposure.image.array.dtype,
1196 )
1197 parentModel = scl.Image(
1198 emptyModel,
1199 bands=context.observation.images.bands,
1200 yx0=(y0, x0),
1201 )
1203 sourceRecords = []
1204 parentBlends = {}
1205 for blendRecord in blendRecords:
1206 try:
1207 blend, blendModel, chi2 = self._deblendParent(blendRecord)
1208 except DeblenderSkippedError as e:
1209 self._skipBlend(blendRecord, e.skipKey, e.message)
1210 parentRecord.set("deblend_skipped", True)
1211 parentRecord.set(e.skipKey, True)
1212 continue
1213 except DeblenderError as e:
1214 blendRecord.set("deblend_error", e.errorName)
1215 blendRecord.set("deblend_failed", True)
1216 self._skipBlend(blendRecord, "deblend_failed", e.message)
1217 parentRecord.set("deblend_childFailed", True)
1218 continue
1220 # Update the parent model
1221 parentModel.insert(blendModel)
1223 # Add each deblended source to the catalog
1224 for scarletSource in blend.sources:
1225 # Add all fields except the HeavyFootprint to the
1226 # source record
1227 scarletSource.peak_id = scarletSource.detectedPeak.getId()
1228 sourceRecord = self._addDeblendedSource(
1229 parent=parentRecord,
1230 blendRecord=blendRecord,
1231 peak=scarletSource.detectedPeak,
1232 catalog=self.catalog,
1233 scarletSource=scarletSource,
1234 chi2=chi2,
1235 )
1236 scarletSource.record_id = sourceRecord.getId()
1237 sourceRecords.append(sourceRecord)
1239 # Store the blend information so that it can be persisted
1240 blendData = scl.io.ScarletBlendData.from_blend(blend)
1241 parentBlends[blendRecord.getId()] = blendData
1243 # Calculate the reduced chi2 for the PSF parent
1244 parentFootprintImage = parentModel.data > 0
1245 chi2 = utils.calcChi2(parentModel, context.observation, parentFootprintImage)
1247 # Update the parent record with the deblending results
1248 self._updateParentRecord(
1249 parentRecord=parentRecord,
1250 nPeaks=len(parentFootprint.peaks),
1251 nChild=np.sum([child["deblend_nChild"] for child in blendRecords]),
1252 nComponents=np.sum([child["deblend_nComponents"] for child in blendRecords]),
1253 runtime=np.sum([child["deblend_runtime"] for child in blendRecords]),
1254 iterations=np.sum([child["deblend_iterations"] for child in blendRecords]),
1255 logL=np.nan,
1256 chi2=np.sum(chi2.data)/np.sum(parentFootprintImage),
1257 spectrumInit=np.all([
1258 child["deblend_spectrumInitFlag"]
1259 for child in blendRecords
1260 ]), # type: ignore
1261 converged=np.all([
1262 child["deblend_blendConvergenceFailedFlag"]
1263 for child in blendRecords
1264 ]), # type: ignore
1265 )
1266 # Persist parent columns to the children
1267 for child in sourceRecords:
1268 for key in self.toCopyFromParent:
1269 child.set(key, parentRecord.get(key))
1270 for parentColumn, childColumn in self.config.columnInheritance.items():
1271 child.set(childColumn, parentRecord.get(parentColumn))
1273 # Persist the blend data
1274 modelData.blends[parentRecord.getId()] = scl.io.HierarchicalBlendData(
1275 children=parentBlends,
1276 )
1278 nDeblendedSources = np.sum(catalog["parent"] != 0)
1279 self.log.info(
1280 "Deblender results: %d parent sources were "
1281 "split into %d deconvolved parents,"
1282 "resulting in %d deblended sources, "
1283 "for a total catalog size of %d sources",
1284 nParents,
1285 nBlends,
1286 nDeblendedSources,
1287 len(catalog),
1288 )
1290 table = afwTable.SourceTable.make(self.schema)
1291 sortedCatalog = afwTable.SourceCatalog(table)
1292 sortedCatalog.extend(catalog, deep=True)
1293 table = afwTable.SourceTable.make(self.blendSchema)
1294 sortedBlendCatalog = afwTable.SourceCatalog(table)
1295 sortedBlendCatalog.extend(blendCatalog, deep=True)
1297 return pipeBase.Struct(
1298 deblendedCatalog=sortedCatalog,
1299 scarletModelData=modelData,
1300 objectParents=sortedBlendCatalog,
1301 )
1303 def _deblendParent(
1304 self,
1305 blendRecord: afwTable.SourceRecord,
1306 ) -> tuple[scl.Blend, scl.Image, scl.Image]:
1307 """Deblend a parent source record
1309 Parameters
1310 ----------
1311 blendRecord :
1312 The parent source record to deblend.
1314 Returns
1315 -------
1316 blend :
1317 The `scl.Blend` object that contains the deblended sources.
1318 blendModel :
1319 The `scl.Image` model of the blend.
1320 chi2 :
1321 The reduced chi2 of the blend model.
1322 """
1323 footprint = blendRecord.getFootprint()
1324 bbox = footprint.getBBox()
1325 peaks = footprint.getPeaks()
1327 # Since we use the first peak for the parent object, we should
1328 # propagate its flags to the parent source.
1329 blendRecord.assign(peaks[0], self.blendPeakSchemaMapper)
1331 # Skip the source if it meets the skipping criteria
1332 isSkipped = self._checkSkipped(blendRecord, self.mExposure)
1333 if isSkipped is not None:
1334 skipKey, skipMessage = isSkipped
1335 raise DeblenderSkippedError(
1336 skipMessage,
1337 blendRecord.getId(),
1338 skipKey,
1339 )
1341 self.log.trace(
1342 "Blend %d: deblending {%d} peaks",
1343 blendRecord.getId(),
1344 len(peaks),
1345 )
1346 # Choose whether or not to use improved spectral initialization.
1347 # This significantly cuts down on the number of iterations
1348 # that the optimizer needs and usually results in a better
1349 # fit.
1350 # But using least squares on a very large blend causes memory
1351 # issues, so it is not done for large blends
1352 if self.config.setSpectra:
1353 if self.config.maxSpectrumCutoff <= 0:
1354 spectrumInit = True
1355 else:
1356 spectrumInit = (
1357 len(footprint.peaks) * bbox.getArea() < self.config.maxSpectrumCutoff
1358 )
1359 else:
1360 spectrumInit = False
1362 try:
1363 t0 = time.monotonic()
1364 # Build the parameter lists with the same ordering
1365 blend = deblend(self.context, footprint, self.config, spectrumInit)
1366 tf = time.monotonic()
1367 runtime = (tf - t0) * 1000
1368 converged = _checkBlendConvergence(blend, self.config.relativeError)
1369 # Store the number of components in the blend
1370 nComponents = len(blend.components)
1371 nChild = len(blend.sources)
1372 # Catch all errors and filter out the ones that we know about
1373 except Exception as e:
1374 blendError = type(e).__name__
1375 if self.config.catchFailures:
1376 # Make it easy to find UnknownErrors in the log file
1377 self.log.warn("UnknownError")
1378 import traceback
1379 traceback.print_exc()
1380 else:
1381 raise
1383 raise DeblenderError(
1384 f"Unable to deblend parent {blendRecord.getId()}: {blendError}",
1385 blendRecord.getId(),
1386 blendError,
1387 )
1389 # Calculate the reduced chi2
1390 blendModel = blend.get_model(convolve=False)
1391 blendFootprintImage = blendModel.data > 0
1392 chi2 = utils.calcChi2(blendModel, self.context.observation, blendFootprintImage)
1394 # Update the blend record with the deblending results
1395 self._updateParentRecord(
1396 parentRecord=blendRecord,
1397 nPeaks=len(peaks),
1398 nChild=nChild,
1399 nComponents=nComponents,
1400 runtime=runtime,
1401 iterations=len(blend.loss),
1402 logL=blend.loss[-1],
1403 chi2=np.sum(chi2.data)/np.sum(blendFootprintImage),
1404 spectrumInit=spectrumInit,
1405 converged=converged,
1406 )
1408 return blend, blendModel, chi2
1410 def _isLargeFootprint(self, footprint: afwDet.Footprint) -> bool:
1411 """Returns whether a Footprint is large
1413 'Large' is defined by thresholds on the area, size and axis ratio,
1414 and total area of the bounding box multiplied by
1415 the number of children.
1416 These may be disabled independently by configuring them to be
1417 non-positive.
1418 """
1419 if (
1420 self.config.maxFootprintArea > 0
1421 and footprint.getArea() > self.config.maxFootprintArea
1422 ):
1423 return True
1424 if self.config.maxFootprintSize > 0:
1425 bbox = footprint.getBBox()
1426 if max(bbox.getWidth(), bbox.getHeight()) > self.config.maxFootprintSize:
1427 return True
1428 if self.config.minFootprintAxisRatio > 0:
1429 axes = afwEll.Axes(footprint.getShape())
1430 if axes.getB() < self.config.minFootprintAxisRatio * axes.getA():
1431 return True
1432 if self.config.maxAreaTimesPeaks > 0:
1433 if (
1434 footprint.getBBox().getArea() * len(footprint.peaks)
1435 > self.config.maxAreaTimesPeaks
1436 ):
1437 return True
1438 return False
1440 def _isMasked(self, footprint: afwDet.Footprint, mExposure: afwImage.MultibandExposure) -> bool:
1441 """Returns whether the footprint violates the mask limits
1443 Parameters
1444 ----------
1445 footprint :
1446 The footprint to check for masked pixels
1447 mExposure :
1448 The multiband exposure to check for masked pixels.
1450 Returns
1451 -------
1452 isMasked : `bool`
1453 `True` if `self.config.maskPlaneLimits` is less than the
1454 fraction of pixels for a given mask in
1455 `self.config.maskLimits`.
1456 """
1457 bbox = footprint.getBBox()
1458 mask = np.bitwise_or.reduce(mExposure.mask[:, bbox].array, axis=0)
1459 size = float(footprint.getArea())
1460 for maskName, limit in self.config.maskLimits.items():
1461 maskVal = mExposure.mask.getPlaneBitMask(maskName)
1462 _mask = afwImage.MaskX(mask & maskVal, xy0=bbox.getMin())
1463 # spanset of masked pixels
1464 maskedSpan = footprint.spans.intersect(_mask, maskVal)
1465 if (maskedSpan.getArea()) / size > limit:
1466 return True
1467 return False
1469 def _skipBlend(
1470 self,
1471 blendRecord: afwTable.SourceRecord,
1472 skipKey: str,
1473 logMessage: str | None,
1474 ):
1475 """Update a parent record that is not being deblended.
1477 This is a fairly trivial function but is implemented to ensure
1478 that a skipped blend updates the appropriate columns
1479 consistently, and always has a flag to mark the reason that
1480 it is being skipped.
1482 Parameters
1483 ----------
1484 blendRecord :
1485 The blend record to flag as skipped.
1486 skipKey :
1487 The name of the flag to mark the reason for skipping.
1488 logMessage :
1489 The message to display in a log.trace when a source
1490 is skipped.
1491 """
1492 if logMessage is not None:
1493 self.log.trace(logMessage)
1494 footprint = blendRecord.getFootprint()
1495 self._updateParentRecord(
1496 parentRecord=blendRecord,
1497 nPeaks=len(footprint.peaks),
1498 nChild=0,
1499 nComponents=0,
1500 runtime=np.nan,
1501 iterations=0,
1502 logL=np.nan,
1503 chi2=np.nan,
1504 spectrumInit=False,
1505 converged=False,
1506 )
1508 # Mark the source as skipped by the deblender and
1509 # flag the reason why.
1510 blendRecord.set("deblend_skipped", True)
1511 blendRecord.set(skipKey, True)
1513 # Add the NOT_DEBLENDED mask to the mask plane in each band
1514 if self.config.notDeblendedMask:
1515 for mask in self.mExposure.mask:
1516 footprint.spans.setMask(
1517 mask, mask.getPlaneBitMask(self.config.notDeblendedMask)
1518 )
1520 def _checkSkipped(
1521 self,
1522 parent: afwTable.SourceRecord,
1523 mExposure: afwImage.MultibandExposure
1524 ) -> tuple[afwTable.Key, str] | None:
1525 """Update a parent record that is not being deblended.
1527 This is a fairly trivial function but is implemented to ensure
1528 that a skipped parent updates the appropriate columns
1529 consistently, and always has a flag to mark the reason that
1530 it is being skipped.
1532 Parameters
1533 ----------
1534 parent :
1535 The parent record to flag as skipped.
1536 mExposure :
1537 The exposures should be co-added images of the same
1538 shape and region of the sky.
1540 Returns
1541 -------
1542 skip :
1543 `True` if the deblender will skip the parent
1544 """
1545 skipKey = None
1546 skipMessage = ""
1547 footprint = parent.getFootprint()
1548 if isPseudoSource(parent, self.config.pseudoColumns):
1549 # We also skip pseudo sources, like sky objects, which
1550 # are intended to be skipped.
1551 skipKey = "deblend_isPseudo"
1552 elif self._isLargeFootprint(footprint):
1553 # The footprint is above the maximum footprint size limit
1554 skipKey = "deblend_parentTooBig"
1555 skipMessage = f"Parent {parent.getId()}: skipping large footprint"
1556 elif self._isMasked(footprint, mExposure):
1557 # The footprint exceeds the maximum number of masked pixels
1558 skipKey = "deblend_masked"
1559 skipMessage = f"Parent {parent.getId()}: skipping masked footprint"
1560 elif (
1561 self.config.maxNumberOfPeaks > 0
1562 and len(footprint.peaks) > self.config.maxNumberOfPeaks
1563 ):
1564 # Unlike meas_deblender, in scarlet we skip the entire blend
1565 # if the number of peaks exceeds max peaks, since neglecting
1566 # to model any peaks often results in catastrophic failure
1567 # of scarlet to generate models for the brighter sources.
1568 skipKey = "deblend_tooManyPeaks"
1569 skipMessage = f"Parent {parent.getId()}: skipping blend with too many peaks"
1570 if skipKey is not None:
1571 return (cast(afwTable.Key, skipKey), skipMessage)
1572 return None
1574 def _updateParentRecord(
1575 self,
1576 parentRecord: afwTable.SourceRecord,
1577 nPeaks: int,
1578 nChild: int,
1579 nComponents: int,
1580 runtime: float,
1581 iterations: int,
1582 logL: float,
1583 chi2: float,
1584 spectrumInit: bool,
1585 converged: bool,
1586 ):
1587 """Update a parent record in all of the single band catalogs.
1589 Ensure that all locations that update a blend record,
1590 whether it is skipped or updated after deblending,
1591 update all of the appropriate columns.
1593 Parameters
1594 ----------
1595 blendRecord :
1596 The catalog record to update.
1597 nPeaks :
1598 Number of peaks in the parent footprint.
1599 nChild :
1600 Number of children deblended from the parent.
1601 This may differ from `nPeaks` if some of the peaks
1602 were culled and have no deblended model.
1603 nComponents :
1604 Total number of components in the parent.
1605 This is usually different than the number of children,
1606 since it is common for a single source to have multiple
1607 components.
1608 runtime :
1609 Total runtime for deblending.
1610 iterations :
1611 Total number of iterations in scarlet before convergence.
1612 logL :
1613 Final log likelihood of the blend.
1614 chi2 :
1615 Final reduced chi2 of the blend.
1616 spectrumInit :
1617 True when scarlet used `set_spectra` to initialize all
1618 sources with better initial intensities.
1619 converged :
1620 True when the optimizer reached convergence before
1621 reaching the maximum number of iterations.
1622 """
1623 parentRecord.set("deblend_nPeaks", nPeaks)
1624 parentRecord.set("deblend_nChild", nChild)
1625 parentRecord.set("deblend_nComponents", nComponents)
1626 parentRecord.set("deblend_runtime", runtime)
1627 parentRecord.set("deblend_iterations", iterations)
1628 parentRecord.set("deblend_logL", logL)
1629 parentRecord.set("deblend_spectrumInitFlag", spectrumInit)
1630 parentRecord.set("deblend_blendConvergenceFailedFlag", converged)
1631 parentRecord.set("deblend_chi2", chi2)
1633 def _buildBlendCatalog(
1634 self,
1635 catalog: afwTable.SourceCatalog,
1636 blendCatalog: afwTable.SourceCatalog,
1637 context: ScarletDeblendContext,
1638 ) -> None:
1639 """Create footprints for deconvolved parents
1641 Each parent may be subdivided into multiple blends that are
1642 isolated in deconvolved space but still blended in the image.
1643 This method finds all of the deconvolved footprints that overlap
1644 with a single parent footprint from the input catalog and
1645 returns a dictionary to map the parent ids to a list of
1646 deconvolved footprints.
1648 Parameters
1649 ----------
1650 catalog :
1651 The merged `SourceCatalog` that contains parent footprints.
1652 blendCatalog :
1653 The catalog of deconvolved parents used for deblending.
1654 context :
1655 The context for the entire coadd.
1656 """
1657 nParents = len(catalog)
1659 for n in range(nParents):
1660 parent = catalog[n]
1661 parentFoot = parent.getFootprint()
1662 # Since we use the first peak for the parent object, we should
1663 # propagate its flags to the parent source.
1664 # For example, this propagates `merge_peak_sky` to the parent
1665 parent.assign(parent.getFootprint().peaks[0], self.peakSchemaMapper)
1667 isPseudo = isPseudoSource(parent, self.config.pseudoColumns)
1668 skipIsolated = len(parentFoot.peaks) < 2 and not self.config.processSingles
1670 if isPseudo or skipIsolated:
1671 # Skip pseudo sources and if processSingles is turned off,
1672 # Skip isolated sources.
1673 # Note: this does not flag isolated sources as skipped or
1674 # set the NOT_DEBLENDED mask in the exposure,
1675 # since these aren't really any skipped blends.
1676 continue
1678 # Find deconvolved footprints that intersect with the parent
1679 # and add them to the blend catalog.
1680 parentId = parent.getId()
1681 self._buildIntersectingFootprints(
1682 parentId,
1683 parentFoot,
1684 blendCatalog,
1685 context.footprints,
1686 context.footprintImage
1687 )
1689 def _buildIntersectingFootprints(
1690 self,
1691 parentId: int,
1692 afwFootprint: afwDet.Footprint,
1693 blendCatalog: afwTable.SourceCatalog,
1694 sclFootprints: list[scl.detect.Footprint],
1695 footprintImage: scl.Image,
1696 ) -> None:
1697 """Get the intersection of two footprints
1699 Parameters
1700 ----------
1701 parentId :
1702 The parent id containing the footprints.
1703 afwFootprint :
1704 The afw footprint
1705 blendCatalog :
1706 The catalog of deconvolved parents to add the new
1707 deconvolved parent to.
1708 sclFootprints :
1709 List of scarlet lite Footprints.
1710 footprintImage :
1711 An indexed image of the scarlet footprints so that the value
1712 of a pixel gives the index + 1 of the footprints that
1713 contain that pixel.
1715 Returns
1716 -------
1717 intersection :
1718 The intersection of the two footprints
1719 """
1720 footprintIndices = set()
1721 ymin, xmin = footprintImage.bbox.origin
1723 # Get the index of the deconvolved footprint at the peak location
1724 for peak in afwFootprint.peaks:
1725 x = peak["i_x"] - xmin
1726 y = peak["i_y"] - ymin
1727 try:
1728 footprintIndex = footprintImage.data[y, x] - 1
1729 except IndexError:
1730 raise RuntimeError(f"no footprint at ({y}, {x})")
1731 if footprintIndex >= 0:
1732 footprintIndices.add(footprintIndex)
1734 # Get the intersection of each deconvolved footprint with
1735 # the parent footprint.
1736 for index in footprintIndices:
1737 _sclFootprint = scarletFootprintToAfw(sclFootprints[index])
1738 intersection = afwFootprint.intersect(_sclFootprint, copyPeaks=True)
1739 if len(intersection.peaks) > 0:
1740 self._addBlendRecord(
1741 parentId=parentId,
1742 blendCatalog=blendCatalog,
1743 footprint=intersection,
1744 )
1746 def _addBlendRecord(
1747 self,
1748 parentId: int,
1749 blendCatalog: afwTable.SourceCatalog,
1750 footprint: afwDet.Footprint,
1751 ) -> None:
1752 """Add deconvolved parents to the parent catalog
1754 Each parent may be subdivided into multiple blends that are
1755 isolated in deconvolved space but still blended in the image.
1756 This function adds the sub-parents to the catalog.
1758 Parameters
1759 ----------
1760 parent :
1761 The parent of the sub-parents.
1762 blendCatalog :
1763 The catalog of deconvolved parents to add the new
1764 deconvolved parent to.
1765 footprint :
1766 The footprint of the deconvolved parent.
1767 """
1768 blendRecord = blendCatalog.addNew()
1769 blendRecord.setParent(parentId)
1770 blendRecord.setFootprint(footprint)
1772 def _addDeblendedSource(
1773 self,
1774 parent: afwTable.SourceRecord,
1775 blendRecord: afwTable.SourceRecord,
1776 peak: afwDet.PeakRecord,
1777 catalog: afwTable.SourceCatalog,
1778 scarletSource: scl.Source,
1779 chi2: scl.Image,
1780 ):
1781 """Add a deblended source to a catalog.
1783 This creates a new child in the source catalog,
1784 assigning it a parent id, and adding all columns
1785 that are independent across all filter bands and not
1786 updated after deblending.
1788 Parameters
1789 ----------
1790 parent :
1791 The parent of the new child record.
1792 blendRecord :
1793 The deconvolved parent of the new child record.
1794 peak :
1795 The peak record for the peak from the parent peak catalog.
1796 catalog :
1797 The source catalog that the child is added to.
1798 scarletSource :
1799 The scarlet model for the new source record.
1800 chi2 :
1801 The chi2 for each pixel.
1803 Returns
1804 -------
1805 src :
1806 The new child source record.
1807 """
1808 src = catalog.addNew()
1809 # The peak catalog is the same for all bands,
1810 # so we just use the first peak catalog
1811 src.assign(peak, self.peakSchemaMapper)
1812 src.setParent(parent.getId())
1813 src.set("deblend_blendId", blendRecord.getId())
1814 src.set("deblend_nPeaks", 1)
1815 src.set("deblend_nChild", 0)
1816 src.set("deblend_blendNChild", len(blendRecord.getFootprint().peaks))
1817 # We set the runtime to zero so that summing up the
1818 # runtime column will give the total time spent
1819 # running the deblender for the catalog.
1820 src.set("deblend_runtime", 0)
1822 # Set the position of the peak from the parent footprint
1823 # This will make it easier to match the same source across
1824 # deblenders and across observations, where the peak
1825 # position is unlikely to change unless enough time passes
1826 # for a source to move on the sky.
1827 src.set("deblend_peak_center_x", peak["i_x"])
1828 src.set("deblend_peak_center_y", peak["i_y"])
1829 src.set("deblend_peakId", peak["id"])
1831 # Store the number of components for the source
1832 src.set("deblend_nComponents", len(scarletSource.components))
1834 # Flag sources missing one or more bands
1835 src.set("deblend_incompleteData", blendRecord.get("deblend_incompleteData"))
1837 # Calculate the reduced chi2 for the source
1838 area = np.sum(scarletSource.get_model().data > 0)
1839 src.set("deblend_chi2", np.sum(chi2[:, scarletSource.bbox].data/area))
1841 return src