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