Coverage for python / lsst / scarlet / lite / initialization.py: 11%
220 statements
« prev ^ index » next coverage.py v7.14.0, created at 2026-05-22 07:47 +0000
1# This file is part of scarlet_lite.
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/>.
22import logging
23from typing import Sequence, cast
25import numpy as np
26from deprecated.sphinx import deprecated # type: ignore
28from .bbox import Box
29from .component import Component, FactorizedComponent
30from .detect import bounds_to_bbox, get_detect_wavelets
31from .image import Image
32from .measure import calculate_snr
33from .observation import Observation
34from .operators import Monotonicity, prox_monotonic_mask, prox_uncentered_symmetry
35from .source import Source
37logger = logging.getLogger("scarlet.lite.initialization")
40def trim_morphology(
41 morph: np.ndarray,
42 threshold: float = 0,
43 bg_thresh: float | None = None,
44) -> tuple[np.ndarray, Box]:
45 """Trim the morphology up to pixels above a threshold
47 Parameters
48 ----------
49 morph:
50 The morphology to be trimmed.
51 threshold:
52 The morphology is trimmed to pixels above the threshold.
53 bg_thresh:
54 Deprecated in favor of ``threshold``.
56 Returns
57 -------
58 morph:
59 The trimmed morphology.
60 box:
61 A tight bounding box around the non-zero pixels of the trimmed
62 morphology. The caller is responsible for any padding.
63 """
64 # Temporarily support bg_thresh
65 if bg_thresh is not None:
66 logger.warning("bg_thresh is deprecated and will be after v29.0. " "Use threshold instead.")
67 threshold = bg_thresh
69 # trim morph to pixels above threshold
70 trimmed = np.where(morph > threshold, morph, 0)
71 bbox = Box.from_data(trimmed, threshold=0)
72 return trimmed, bbox
75def init_monotonic_morph(
76 detect: np.ndarray,
77 center: tuple[int, int],
78 full_box: Box,
79 padding: int = 5,
80 normalize: bool = True,
81 monotonicity: Monotonicity | None = None,
82 threshold: float = 0,
83 max_iter: int = 0,
84 center_radius: int = 1,
85 variance_factor: float = 0,
86) -> tuple[Box, np.ndarray | None]:
87 """Initialize a morphology for a monotonic source
89 Parameters
90 ----------
91 detect:
92 The 2D detection image contained in `full_box`.
93 center:
94 The center of the monotonic source.
95 full_box:
96 The bounding box of `detect`.
97 padding:
98 The number of pixels to grow the morphology in each direction.
99 This can be useful if initializing a source with a kernel that
100 is known to be narrower than the expected value of the source.
101 normalize:
102 Whether or not to normalize the morphology.
103 monotonicity:
104 When `monotonicity` is `None`,
105 the component is initialized with only the
106 monotonic pixels, otherwise the monotonicity operator is used to
107 project the morphology to a monotonic solution.
108 threshold:
109 The minimum value to use for trimming the
110 morphology.
111 max_iter:
112 The maximum number of iterations to use in the monotonicity operator.
113 Only used if `monotonicity` is `None`.
114 center_radius:
115 The amount that the center can be shifted to a local maximum.
116 Only used if `monotonicity` is `None`.
117 variance_factor:
118 The average variance in the image.
119 This is used to allow pixels to be non-monotonic up to
120 `variance` * `noise_rms`, so setting `variance = 0`
121 will force strict monotonicity in the mask.
122 Only used if `monotonicity` is `None`.
124 Returns
125 -------
126 bbox:
127 The bounding box of the morphology.
128 morph:
129 The initialized morphology.
130 """
131 center: tuple[int, int] = tuple(center[i] - full_box.origin[i] for i in range(2)) # type: ignore
133 if monotonicity is None:
134 _, morph, bounds = prox_monotonic_mask(
135 x=detect,
136 center=center,
137 center_radius=center_radius,
138 variance=variance_factor,
139 max_iter=max_iter,
140 )
141 bbox = bounds_to_bbox(bounds)
142 if bbox.shape == (1, 1) and morph[bbox.slices][0, 0] == 0:
143 return Box((0, 0)), None
145 if threshold > 0:
146 morph, bbox = trim_morphology(morph, threshold=threshold)
148 else:
149 morph = monotonicity(detect, center)
151 # truncate morph at thresh * bg_rms
152 morph, bbox = trim_morphology(morph, threshold=threshold)
154 # Shift the bounding box to account for the non-zero origin
155 bbox += full_box.origin
157 if np.max(morph) == 0:
158 return Box((0, 0), origin=full_box.origin), None
160 if normalize:
161 morph /= np.max(morph)
163 if padding is not None and padding > 0:
164 # Pad the morphology to allow it to grow
165 bbox = bbox.grow(padding)
167 # Ensure that the bounding box is inside the full box,
168 # even after padding.
169 bbox = bbox & full_box
170 return bbox, morph
173def multifit_spectra(
174 observation: Observation,
175 morphs: Sequence[Image],
176 model: Image | None = None,
177) -> np.ndarray:
178 """Fit the spectra of multiple components simultaneously
180 Parameters
181 ----------
182 observation:
183 The class containing the observation data.
184 morphs:
185 The morphology of each component.
186 model:
187 An optional model for sources that are not factorized,
188 and thus will not have their spectra fit.
189 This model is subtracted from the data before fitting the other
190 spectra.
192 Returns
193 -------
194 spectra:
195 The spectrum for each component, in the same order as `morphs`.
196 """
197 _bands = observation.bands
198 n_bands = len(_bands)
199 dtype = observation.images.dtype
201 if model is not None:
202 image = observation.images - model
203 else:
204 image = observation.images.copy()
206 morph_images = np.zeros((n_bands, len(morphs), image.data[0].size), dtype=dtype)
207 for idx, morph in enumerate(morphs):
208 _image = morph.repeat(observation.bands)
209 _image = Image.from_box(image.bbox, bands=image.bands).insert(_image)
210 morph_images[:, idx] = observation.convolve(_image).data.reshape(n_bands, -1)
212 spectra = np.zeros((len(morphs), n_bands), dtype=dtype)
214 for b in range(n_bands):
215 a = morph_images[b].T
216 spectra[:, b] = np.linalg.lstsq(a, image[observation.bands[b]].data.flatten(), rcond=None)[0]
217 spectra[spectra < 0] = 0
218 return spectra
221class FactorizedInitialization:
222 """Common variables and methods for both Factorized Component schemes
224 There are a large number of parameters that are universal for all of the
225 sources being initialized from the same set of observed images.
226 To simplify the API those parameters are all initialized by this class
227 and passed to `init_source` for each source.
228 It also creates temporary objects that only need to be created once for
229 all of the sources in a blend.
231 Parameters
232 ----------
233 observation:
234 The observation containing the blend
235 centers:
236 The center of each source to initialize.
237 detect:
238 The array that contains a 2D image used for detection.
239 min_snr:
240 The minimum SNR required per component.
241 So a 2-component source requires at least `2*min_snr` while sources
242 with SNR < `min_snr` will be initialized with the PSF.
243 monotonicity:
244 When `monotonicity` is `None`,
245 the component is initialized with only the
246 monotonic pixels, otherwise the monotonicity operator is used to
247 project the morphology to a monotonic solution.
248 disk_percentile:
249 The percentage of the overall flux to attribute to the disk.
250 bg_thresh:
251 The fraction of the background RMS to use as a threshold for the
252 morphology (in other words the threshold is set at
253 `bg_thresh * noise_rms`).
254 initital_bg_thresh:
255 The same as bg_thresh, but this allows a different background
256 threshold to be used for initialization.
257 thresh:
258 Deprecated. Use `initial_bg_thresh` instead.
259 padding:
260 The amount to pad the morphology to allow for extra flux
261 in the first few iterations before resizing.
262 use_sparse_init:
263 Use a monotonic mask to prevent initial source models from growing
264 too large.
265 max_components:
266 The maximum number of components in the source.
267 This should be one or two.
268 convolved:
269 Deprecated. This is now calculated in __init__, but the
270 old API is supported until v29.0.
271 is_symmetric:
272 Whether or not the sources are symmetric.
273 This is used to determine whether to use the symmetry operator
274 for initialization.
275 """
277 def __init__(
278 self,
279 observation: Observation,
280 centers: Sequence[tuple[int, int]],
281 detect: np.ndarray | None = None,
282 min_snr: float = 50,
283 monotonicity: Monotonicity | None = None,
284 disk_percentile: float = 25,
285 initial_bg_thresh: float = 0.5,
286 bg_thresh: float = 0.25,
287 thresh: float | None = None,
288 padding: int = 2,
289 use_sparse_init: bool = True,
290 max_components: int = 2,
291 convolved: Image | None = None,
292 is_symmetric: bool = False,
293 ):
294 if detect is None:
295 # Build the morphology detection image
296 detect = np.sum(
297 observation.images.data / (observation.noise_rms**2)[:, None, None],
298 axis=0,
299 )
300 self.detect = detect
302 if convolved is None:
303 _detect = Image(detect)
304 convolved = observation.convolve(_detect.repeat(observation.bands), mode="real")
305 else:
306 logger.warning(
307 "convolved is deprecated and will be removed after v29.0. "
308 "The convolved image is now calculated in __init__ and does "
309 "not need to be specified."
310 )
312 if thresh is not None:
313 initial_bg_thresh = thresh
315 self.observation = observation
316 self.convolved = convolved
317 # Ensure that each center is a tuple of integers
318 centers = [(int(round(c[0])), int(round(c[1]))) for c in centers]
319 self.centers = centers
320 self.min_snr = min_snr
321 self.monotonicity = monotonicity
322 self.use_sparse_init = use_sparse_init
323 self.is_symmetric = is_symmetric
325 # Get the model PSF
326 # Convolve the PSF in order to set the spectrum
327 # of a point source correctly.
328 model_psf = Image(cast(np.ndarray, observation.model_psf)[0])
329 convolved_psf = model_psf.repeat(observation.bands)
330 self.convolved_psf = observation.convolve(convolved_psf, mode="real").data
331 # Get the "spectrum" of the PSF
332 self.py = model_psf.shape[0] // 2
333 self.px = model_psf.shape[1] // 2
334 self.psf_spectrum = self.convolved_psf[:, self.py, self.px]
336 # Set the maximum number of components for any source in the blend
337 if max_components < 0 or max_components > 2:
338 raise ValueError(f"max_components must be 0, 1 or 2, got {max_components}")
339 self.max_components = max_components
340 self.initial_bg_thresh = initial_bg_thresh
341 self.bg_thresh = bg_thresh
342 self.padding = padding
343 self.disk_percentile = disk_percentile
345 # Initalize all of the sources
346 sources = []
347 for center in centers:
348 if max_components == 0:
349 source = Source([self.get_psf_component(center)])
350 else:
351 source = self.init_source(center)
352 sources.append(source)
353 self.sources = sources
355 @property
356 def thresh(self):
357 logger.warning(
358 "thresh is deprecated and will be removed after v29.0. " "Use initial_bg_thresh instead."
359 )
360 return self.initial_bg_thresh
362 def get_snr(self, center: tuple[int, int]) -> float:
363 """Get the SNR at the center of a component
365 Parameters
366 ----------
367 center:
368 The location of the center of the source.
370 Returns
371 -------
372 result:
373 The SNR at the center of the component.
374 """
375 snr = np.floor(
376 calculate_snr(
377 self.observation.images,
378 self.observation.variance,
379 self.observation.psfs,
380 center,
381 )
382 )
383 return snr / self.min_snr
385 def get_psf_component(self, center: tuple[int, int]) -> FactorizedComponent:
386 """Create a factorized component with a PSF morphology
388 Parameters
389 ----------
390 center:
391 The center of the component.
393 Returns
394 -------
395 component:
396 A `FactorizedComponent` with a PSF-like morphology.
397 """
398 if not self.observation.bbox.contains(center):
399 raise ValueError(f"Source center {center} is outside the observation {self.observation.bbox}")
400 local_center = (
401 center[0] - self.observation.bbox.origin[0],
402 center[1] - self.observation.bbox.origin[1],
403 )
404 # There wasn't sufficient flux for an extended source,
405 # so create a PSF source.
406 spectrum_center = (slice(None), local_center[0], local_center[1])
407 img_center = self.observation.images.data[spectrum_center]
408 spectrum = np.divide(
409 img_center,
410 self.psf_spectrum,
411 out=np.zeros_like(img_center),
412 where=self.psf_spectrum > 0,
413 )
414 spectrum[spectrum < 0] = 0
416 psf = cast(np.ndarray, self.observation.model_psf)[0].copy()
417 py = psf.shape[0] // 2
418 px = psf.shape[1] // 2
419 psf_bbox = Box(psf.shape, origin=(-py + center[0], -px + center[1]))
420 bbox = self.observation.bbox & psf_bbox
421 morph = Image(psf, yx0=cast(tuple[int, int], psf_bbox.origin))[bbox].data
422 component = FactorizedComponent(
423 self.observation.bands,
424 spectrum,
425 morph,
426 bbox,
427 center,
428 self.observation.noise_rms,
429 monotonicity=self.monotonicity,
430 is_symmetric=self.is_symmetric,
431 )
432 return component
434 def get_single_component(
435 self,
436 center: tuple[int, int],
437 detect: np.ndarray,
438 thresh: float,
439 padding: int,
440 ) -> FactorizedComponent | None:
441 """Initialize parameters for a `FactorizedComponent`
443 Parameters
444 ----------
445 center:
446 The location of the center of the source to detect in the
447 full image.
448 detect:
449 The image used for detection of the morphology.
450 thresh:
451 The lower cutoff threshold to use for the morphology.
452 padding:
453 The amount to pad the morphology to allow for extra flux
454 in the first few iterations before resizing.
456 Returns
457 -------
458 component:
459 A `FactorizedComponent` created from the detection image.
461 """
462 local_center = (
463 center[0] - self.observation.bbox.origin[0],
464 center[1] - self.observation.bbox.origin[1],
465 )
467 if self.use_sparse_init:
468 monotonicity = None
469 else:
470 monotonicity = self.monotonicity
471 bbox, morph = init_monotonic_morph(
472 detect,
473 center,
474 self.observation.bbox,
475 padding=padding,
476 normalize=False,
477 monotonicity=monotonicity,
478 threshold=thresh,
479 )
481 if morph is None:
482 return None
483 morph = morph[(bbox - self.observation.bbox.origin).slices]
485 spectrum_center = (slice(None), local_center[0], local_center[1])
486 images = self.observation.images
488 convolved = self.convolved
489 img_center = images.data[spectrum_center]
490 conv_center = convolved.data[spectrum_center]
491 spectrum = np.divide(
492 img_center,
493 conv_center,
494 out=np.zeros_like(img_center),
495 where=conv_center > 0,
496 )
497 spectrum[spectrum < 0] = 0
498 morph_max = np.max(morph)
499 spectrum *= morph_max
500 morph /= morph_max
502 return FactorizedComponent(
503 bands=self.observation.bands,
504 spectrum=spectrum,
505 morph=morph,
506 bbox=bbox,
507 peak=center,
508 bg_rms=self.observation.noise_rms,
509 bg_thresh=self.bg_thresh,
510 monotonicity=self.monotonicity,
511 is_symmetric=self.is_symmetric,
512 )
514 def init_source(self, center: tuple[int, int]) -> Source:
515 """Initialize a source from a chi^2 detection.
517 Parameters
518 ----------
519 center:
520 The center of the source.
521 """
522 # Some operators need the local center, not center in the full image
523 local_center = (
524 center[0] - self.observation.bbox.origin[0],
525 center[1] - self.observation.bbox.origin[1],
526 )
528 # Calculate the signal to noise at the center of this source
529 component_snr = self.get_snr(center)
531 # Initialize the bbox, morph, and spectrum
532 # for a single component source
533 detect = prox_uncentered_symmetry(self.detect.copy(), local_center, fill=0)
534 thresh = np.mean(self.observation.noise_rms) * self.initial_bg_thresh
535 component = self.get_single_component(center, detect, thresh, self.padding)
537 if component is None:
538 # There wasn't enough flux to initialize the source as
539 # as single component, so initialize it with the model PSF.
540 components = [self.get_psf_component(center)]
541 elif component_snr < 2 or self.max_components < 2:
542 # There isn't sufficient flux to add a second component.
543 components = [component]
544 else:
545 # There was enough flux for a 2-component source,
546 # so split the single component model into two components,
547 # using the same algorithm as scarlet main.
548 bulge_morph = component.morph.copy()
549 disk_morph = component.morph.copy()
550 # Set the threshold for the bulge.
551 # Since the morphology is monotonic, this selects the inner
552 # of the single component morphology and assigns it to the bulge.
553 flux_thresh = self.disk_percentile / 100
554 mask = disk_morph > flux_thresh
555 # Remove the flux above the threshold so that the disk will have
556 # a flat center.
557 disk_morph[mask] = flux_thresh
558 # Subtract off the thresholded flux (since we're normalizing the
559 # morphology anyway) so that it does not have a sharp
560 # discontinuity at the edge.
561 bulge_morph -= flux_thresh
562 bulge_morph[bulge_morph < 0] = 0
564 bulge_morph /= np.max(bulge_morph)
565 disk_morph /= np.max(disk_morph)
567 # Fit the spectra assuming that all of the flux in the image
568 # is due to both components. This is not true, but for the
569 # vast majority of sources this is a good approximation.
570 try:
571 bulge_spectrum, disk_spectrum = multifit_spectra(
572 self.observation,
573 [
574 Image(bulge_morph, yx0=cast(tuple[int, int], component.bbox.origin)),
575 Image(disk_morph, yx0=cast(tuple[int, int], component.bbox.origin)),
576 ],
577 )
579 components = [
580 FactorizedComponent(
581 bands=self.observation.bands,
582 spectrum=bulge_spectrum,
583 morph=bulge_morph,
584 bbox=component.bbox.copy(),
585 peak=center,
586 bg_rms=self.observation.noise_rms,
587 bg_thresh=self.bg_thresh,
588 monotonicity=self.monotonicity,
589 is_symmetric=self.is_symmetric,
590 ),
591 FactorizedComponent(
592 bands=self.observation.bands,
593 spectrum=disk_spectrum,
594 morph=disk_morph,
595 bbox=component.bbox.copy(),
596 peak=center,
597 bg_rms=self.observation.noise_rms,
598 bg_thresh=self.bg_thresh,
599 monotonicity=self.monotonicity,
600 is_symmetric=self.is_symmetric,
601 ),
602 ]
603 except np.linalg.LinAlgError:
604 components = [component]
606 return Source(components) # type: ignore
609@deprecated(
610 reason="This class is replaced by FactorizedInitialization and will be removed after v29.0",
611 version="v29.0",
612 category=FutureWarning,
613)
614class FactorizedChi2Initialization(FactorizedInitialization):
615 """Initialize all sources with chi^2 detections
617 There are a large number of parameters that are universal for all of the
618 sources being initialized from the same set of observed images.
619 To simplify the API those parameters are all initialized by this class
620 and passed to `init_main_source` for each source.
621 It also creates temporary objects that only need to be created once for
622 all of the sources in a blend.
624 Parameters
625 ----------
626 observation:
627 The observation containing the blend
628 centers:
629 The center of each source to initialize.
630 detect:
631 The array that contains a 2D image used for detection.
632 min_snr:
633 The minimum SNR required per component.
634 So a 2-component source requires at least `2*min_snr` while sources
635 with SNR < `min_snr` will be initialized with the PSF.
636 monotonicity:
637 When `monotonicity` is `None`,
638 the component is initialized with only the
639 monotonic pixels, otherwise the monotonicity operator is used to
640 project the morphology to a monotonic solution.
641 disk_percentile:
642 The percentage of the overall flux to attribute to the disk.
643 thresh:
644 The threshold used to trim the morphology,
645 so all pixels below `thresh * bg_rms` are set to zero.
646 padding:
647 The amount to pad the morphology to allow for extra flux
648 in the first few iterations before resizing.
649 """
651 pass
654@deprecated(
655 reason="FactorizedWaveletInitialization will be removed after v29.0 "
656 "since it does not appear to offer any advantages over "
657 "FactorizedChi2Initialization. Consider switching to "
658 "FactorizedInitialization now.",
659 version="v29.0",
660 category=FutureWarning,
661)
662class FactorizedWaveletInitialization(FactorizedInitialization):
663 """Parameters used to initialize all sources with wavelet detections
665 There are a large number of parameters that are universal for all of the
666 sources being initialized from the same set of wavelet coefficients.
667 To simplify the API those parameters are all initialized by this class
668 and passed to `init_wavelet_source` for each source.
670 Parameters
671 ----------
672 observation:
673 The multiband observation of the blend.
674 centers:
675 The center of each source to initialize.
676 bulge_slice, disk_slice:
677 The slice used to select the wavelet scales used for the
678 bulge/disk.
679 bulge_padding, disk_padding:
680 The number of pixels to grow the bounding box of the bulge/disk
681 to leave extra room for growth in the first few iterations.
682 use_psf:
683 Whether or not to use the PSF for single component sources.
684 If `use_psf` is `False` then only sources with low signal
685 at all scales are initialized with the PSF morphology.
686 scales:
687 Number of wavelet scales to use.
688 wavelets:
689 The array of wavelet coefficients `(scale, y, x)`
690 used for detection.
691 monotonicity:
692 When `monotonicity` is `None`,
693 the component is initialized with only the
694 monotonic pixels, otherwise the monotonicity operator is used to
695 project the morphology to a monotonic solution.
696 min_snr:
697 The minimum SNR required per component.
698 So a 2-component source requires at least `2*min_snr` while sources
699 with SNR < `min_snr` will be initialized with the PSF.
700 """
702 def __init__(
703 self,
704 observation: Observation,
705 centers: Sequence[tuple[int, int]],
706 bulge_slice: slice = slice(None, 2),
707 disk_slice: slice = slice(2, -1),
708 bulge_padding: int = 5,
709 disk_padding: int = 5,
710 use_psf: bool = True,
711 scales: int = 5,
712 wavelets: np.ndarray | None = None,
713 monotonicity: Monotonicity | None = None,
714 min_snr: float = 50,
715 use_sparse_init: bool = True,
716 ):
717 if wavelets is None:
718 wavelets = get_detect_wavelets(
719 observation.images.data,
720 observation.variance.data,
721 scales=scales,
722 )
723 wavelets[wavelets < 0] = 0
724 # The detection coadd for single component sources
725 detectlets = np.sum(wavelets[:-1], axis=0)
726 # The detection coadd for the bulge
727 bulgelets = np.sum(wavelets[bulge_slice], axis=0)
728 # The detection coadd for the disk
729 disklets = np.sum(wavelets[disk_slice], axis=0)
731 self.detectlets = detectlets
732 self.bulgelets = bulgelets
733 self.disklets = disklets
734 self.bulge_grow = bulge_padding
735 self.disk_grow = disk_padding
736 self.use_psf = use_psf
738 # Initialize the sources
739 super().__init__(
740 observation=observation,
741 centers=centers,
742 detect=detectlets,
743 min_snr=min_snr,
744 monotonicity=monotonicity,
745 use_sparse_init=use_sparse_init,
746 )
748 def init_source(self, center: tuple[int, int]) -> Source:
749 """Initialize a source from a chi^2 detection.
751 Parameters
752 ----------
753 center:
754 The center of the source.
755 """
756 local_center = (
757 center[0] - self.observation.bbox.origin[0],
758 center[1] - self.observation.bbox.origin[1],
759 )
760 nbr_components = self.get_snr(center)
761 observation = self.observation
763 components: list[Component] | None = None
765 if (nbr_components < 1 and self.use_psf) or self.detectlets[local_center[0], local_center[1]] <= 0:
766 # Initialize the source as an PSF source
767 components = [self.get_psf_component(center)]
768 elif nbr_components < 2:
769 # Inititialize with a single component
770 component = self.get_single_component(center, self.detectlets, 0, self.disk_grow)
771 if component is not None:
772 components = [component]
773 else:
774 # Initialize with a 2 component model
775 bulge_box, bulge_morph = init_monotonic_morph(
776 self.bulgelets, center, observation.bbox, self.bulge_grow
777 )
778 disk_box, disk_morph = init_monotonic_morph(
779 self.disklets, center, observation.bbox, self.disk_grow
780 )
781 if bulge_morph is None or disk_morph is None:
782 if bulge_morph is None:
783 if disk_morph is None:
784 raise RuntimeError("Both components are None")
785 # One of the components was null,
786 # so initialize as a single component
787 component = self.get_single_component(center, self.detectlets, 0, self.disk_grow)
788 if component is not None:
789 components = [component]
790 else:
791 local_bulge_box = bulge_box - self.observation.bbox.origin
792 local_disk_box = disk_box - self.observation.bbox.origin
793 bulge_morph = bulge_morph[local_bulge_box.slices]
794 disk_morph = disk_morph[local_disk_box.slices]
796 bulge_spectrum, disk_spectrum = multifit_spectra(
797 observation,
798 [
799 Image(bulge_morph, yx0=cast(tuple[int, int], bulge_box.origin)),
800 Image(disk_morph, yx0=cast(tuple[int, int], disk_box.origin)),
801 ],
802 )
804 components = []
805 if np.any(bulge_spectrum != 0):
806 components.append(
807 FactorizedComponent(
808 observation.bands,
809 bulge_spectrum,
810 bulge_morph,
811 bulge_box,
812 center,
813 monotonicity=self.monotonicity,
814 )
815 )
816 else:
817 logger.debug("cut bulge")
818 if np.any(disk_spectrum != 0):
819 components.append(
820 FactorizedComponent(
821 observation.bands,
822 disk_spectrum,
823 disk_morph,
824 disk_box,
825 center,
826 monotonicity=self.monotonicity,
827 )
828 )
829 else:
830 logger.debug("cut disk")
832 # If every init path above either failed (left components as
833 # None) or produced no components (empty list when bulge and
834 # disk spectra were both cut), fall back to a PSF source -- a
835 # point-source model is the most conservative non-empty model
836 # we can return.
837 if not components:
838 logger.debug("fall back to PSF source")
839 components = [self.get_psf_component(center)]
840 return Source(components)