Coverage for tests/test_utils.py: 99%
259 statements
« prev ^ index » next coverage.py v7.14.2, created at 2026-06-22 08:44 +0000
1# This file is part of meas_extensions_scarlet.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
22"""Tests for ``lsst.meas.extensions.scarlet.utils``."""
24import unittest
25import warnings
27import lsst.afw.image as afwImage
28import lsst.meas.extensions.scarlet as mes
29import lsst.scarlet.lite as scl
30import lsst.utils.tests
31import numpy as np
32import scipy.signal
33from lsst.afw.detection import Footprint, GaussianPsf, InvalidPsfError, PeakTable, Psf
34from lsst.afw.geom import SpanSet
35from lsst.afw.table import SourceCatalog, SourceTable
36from lsst.geom import Extent2I, Point2D, Point2I
37from lsst.pipe.base import NoWorkFound
40class BadPsf(Psf):
41 def __init__(self, validPoint: Point2D, psf: GaussianPsf):
42 self.validPoint = validPoint
43 self.psf = psf
44 super().__init__()
46 def computeKernelImage(self, location: Point2D):
47 if location == self.validPoint:
48 return self.psf.computeKernelImage(location)
49 raise InvalidPsfError(f"Invalid PSF at location {location}")
52class MultiPointBadPsf(Psf):
53 """A PSF valid at multiple discrete integer locations, with a
54 potentially different underlying ``GaussianPsf`` per location.
56 Used to drive the multiband-PSF fallback search through scenarios
57 where the same band would have multiple acceptable fallback
58 locations with distinguishable kernels.
59 """
61 def __init__(self, validPsfs: dict[tuple[int, int], GaussianPsf]):
62 self.validPsfs = validPsfs
63 super().__init__()
65 def computeKernelImage(self, location: Point2D):
66 key = (int(location.getX()), int(location.getY()))
67 if key in self.validPsfs:
68 return self.validPsfs[key].computeKernelImage(location)
69 raise InvalidPsfError(f"Invalid PSF at location {location}")
72class TestUtils(lsst.utils.tests.TestCase):
73 def setUp(self):
74 self.bands = tuple("gri")
76 def test_computeNearestPsfGood(self):
77 # Test that using a valid PSF works normally
78 psf, psfImage = self._generateGoodPsf()
79 coadd = self._generateCoadd(psf)
81 # Test that computing the PSF works
82 derivedPsf, center, dist = mes.utils.computeNearestPsf(coadd, None, "g", Point2D(25, 25))
83 np.testing.assert_array_equal(derivedPsf.array, psfImage)
84 self.assertEqual(center, Point2D(25, 25))
85 self.assertEqual(dist, 0)
87 def test_computeNearestPsfRecoverable(self):
88 # Test that using a PSF not defined at the initial location
89 # will fallback to a valid location.
90 psf, psfImage = self._generateGoodPsf()
91 coadd = self._generateCoadd(BadPsf(Point2D(1, 1), psf))
92 catalog = self._generateCatalog(self.bands, [[(1, 1, 10)]])
94 # Test that computing the PSF works after finding a new location.
95 # Since the PSF above is *only* defined at (1, 1) it will fail to
96 # compute a PSF image at (4, 5) but should fall back to (1, 1).
97 # Per finding U-6 of the ``audits/audit-2026-05-05.md`` audit,
98 # the fallback path previously returned a ``Point2I`` while the
99 # direct-success path returned a ``Point2D``; both now uniformly
100 # return ``Point2D`` and ``Point2D != Point2I`` even at the same
101 # coordinates, so the type check below is also a regression guard
102 # on the unified return type.
103 derivedPsf, center, dist = mes.utils.computeNearestPsf(coadd, catalog, None, Point2D(4, 5))
104 np.testing.assert_array_equal(derivedPsf.array, psfImage)
105 self.assertIsInstance(center, Point2D)
106 self.assertEqual(center, Point2D(1, 1))
107 self.assertEqual(dist, 5)
109 def test_computeNearestPsfBad(self):
110 # Test that a PSF that cannot find a matching location returns None
111 psf = self._generateGoodPsf()
112 coadd = self._generateCoadd(BadPsf(Point2D(1, 1), psf))
113 catalog = self._generateCatalog(self.bands)
115 # Test that computing the PSF cannot generate a PSF
116 derivedPsf, center, dist = mes.utils.computeNearestPsf(coadd, catalog, None, Point2D(4, 5))
117 self.assertIsNone(derivedPsf)
118 self.assertIsNone(center)
119 self.assertIsNone(dist)
121 def test_computeNearestPsfMultiBandGood(self):
122 # Test that a valid PSF in every band works normally
123 bands = tuple("gri")
124 psfs, psfImage = self._generateMultibandPsf([1.0, 1.2, 1.4])
125 mCoadd = self._generateMultibandCoadd(psfs, bands)
127 # Test that computing the PSF works
128 psfArray, newCoadd = mes.utils.computeNearestPsfMultiBand(mCoadd, Point2D(25, 25), None)
129 np.testing.assert_array_equal(psfArray, psfImage)
130 self.assertTupleEqual(newCoadd.bands, bands)
132 def test_computeNearestPsfMultiBandRecoverable(self):
133 # Test that a Psf at a different location is still recoverable
134 bands = tuple("gri")
135 psfs, psfImage = self._generateMultibandPsf([1.0, 1.2, 1.4])
136 psfs[1] = BadPsf(Point2D(1, 1), psfs[1])
137 mCoadd = self._generateMultibandCoadd(psfs, bands)
138 catalog = self._generateCatalog(self.bands, [[(1, 1, 10)]])
140 # Test that computing the PSF works because the catalog has a peak
141 # at the location of the BadPsf.
142 psfArray, newCoadd = mes.utils.computeNearestPsfMultiBand(mCoadd, Point2D(25, 25), catalog)
143 np.testing.assert_array_equal(psfArray, psfImage)
144 self.assertTupleEqual(newCoadd.bands, bands)
146 def test_computeNearestPsfMultiBandIncomplete(self):
147 # Test that missing a PSF in one band returns a PSF and
148 # an exposure that are missing bands.
149 bands = tuple("gri")
150 psfs, psfImage = self._generateMultibandPsf([1.0, 1.2, 1.4])
151 psfs[1] = BadPsf(Point2D(1, 1), psfs[1])
152 mCoadd = self._generateMultibandCoadd(psfs, bands)
153 catalog = self._generateCatalog(self.bands)
155 # Test that computing the PSF works for the g- and i-band PSFs that
156 # are not BadPsf.
157 psfArray, newCoadd = mes.utils.computeNearestPsfMultiBand(mCoadd, Point2D(25, 25), catalog)
158 np.testing.assert_array_equal(psfArray, np.delete(psfImage, 1, axis=0))
159 self.assertTupleEqual(newCoadd.bands, tuple("gi"))
161 def test_computeNearestPsfMultiBandBad(self):
162 # Test that None is returned if none of the PSFs can be computed
163 bands = tuple("gri")
164 psfs, psfImage = self._generateMultibandPsf([1.0, 1.2, 1.4])
165 psfs = [BadPsf(Point2D(1, 1), psfs) for psf in psfs]
166 mCoadd = self._generateMultibandCoadd(psfs, bands)
167 catalog = self._generateCatalog(self.bands)
169 psfArray, newCoadd = mes.utils.computeNearestPsfMultiBand(mCoadd, Point2D(25, 25), catalog)
170 self.assertIsNone(psfArray)
171 self.assertIsNone(newCoadd)
173 def test_computeNearestPsfMultiBand_search_center_does_not_drift(self):
174 """When one band falls back, the search for subsequent bands
175 stays anchored at the requested center.
177 Requested center is ``(25, 25)``. Band g's PSF is valid only at
178 ``(1, 1)`` (≈34 px from center). Band r's PSF is valid at both
179 ``(10, 10)`` with σ=2.0 and ``(30, 30)`` with σ=1.2. Sorting the
180 catalog peaks by distance from the *requested center* puts
181 ``(30, 30)`` first for band r; sorting by distance from band g's
182 fallback ``(1, 1)`` puts ``(10, 10)`` first. The bug carried band
183 g's fallback into r's search, so r returned the σ=2.0 kernel at
184 ``(10, 10)``. The fix re-anchors r's search at ``(25, 25)``, so r
185 returns the σ=1.2 kernel at ``(30, 30)``. Regression test for
186 finding C-6 of the ``audits/audit-2026-05-05.md`` audit.
187 """
188 bands = ("g", "r")
189 g_inner = GaussianPsf(41, 41, 1.0)
190 r_at_10 = GaussianPsf(41, 41, 2.0)
191 r_at_30 = GaussianPsf(41, 41, 1.2)
192 psfs = [
193 BadPsf(Point2D(1, 1), g_inner),
194 MultiPointBadPsf({(10, 10): r_at_10, (30, 30): r_at_30}),
195 ]
196 mCoadd = self._generateMultibandCoadd(psfs, bands)
197 catalog = self._generateCatalog(
198 bands, [[(1, 1, 10), (10, 10, 10), (30, 30, 10)]]
199 )
201 psfArray, newCoadd = mes.utils.computeNearestPsfMultiBand(
202 mCoadd, Point2D(25, 25), catalog
203 )
205 self.assertTupleEqual(newCoadd.bands, bands)
206 # GaussianPsf kernels are stationary so their bboxes don't depend
207 # on the position they were computed at; the peak amplitude is
208 # 1/(2πσ²), so each σ value has a distinct kernel max we can
209 # discriminate on. Band g should land at (1, 1) (its only valid
210 # location, σ=1.0); band r should land at (30, 30) (σ=1.2 — the
211 # closest valid r position to the requested center). Under the
212 # bug, band r would land at (10, 10) and return the σ=2.0 kernel.
213 arr = np.asarray(psfArray)
214 expected_g = g_inner.computeKernelImage(Point2D(1, 1)).array
215 expected_r = r_at_30.computeKernelImage(Point2D(30, 30)).array
216 self.assertAlmostEqual(arr[0].max(), expected_g.max(), places=4)
217 self.assertAlmostEqual(arr[1].max(), expected_r.max(), places=4)
219 def test_computeNearestPsfMultiBand_upgrades_all_bands_to_common(self):
220 """When the fallback location found for a failing band is also
221 valid for the bands that succeeded at the center, every band is
222 re-sampled at that common location — including bands that
223 already had a PSF at the center. The previously-successful
224 center PSFs are discarded.
226 Band g is invalid at the requested center ``(25, 25)`` but valid
227 at ``(40, 40)`` with σ=1.0. Band r is valid at *both* ``(25, 25)``
228 with σ=1.2 *and* ``(40, 40)`` with σ=1.5. Because the common
229 fallback ``(40, 40)`` is also valid for r, the upgrade fires and
230 r's returned kernel is the σ=1.5 one (re-sampled at (40, 40)),
231 not the σ=1.2 kernel r had at the center.
232 """
233 bands = ("g", "r")
234 g_at_40 = GaussianPsf(41, 41, 1.0)
235 r_at_center = GaussianPsf(41, 41, 1.2)
236 r_at_40 = GaussianPsf(41, 41, 1.5)
237 psfs = [
238 BadPsf(Point2D(40, 40), g_at_40),
239 MultiPointBadPsf({(25, 25): r_at_center, (40, 40): r_at_40}),
240 ]
241 mCoadd = self._generateMultibandCoadd(psfs, bands)
242 catalog = self._generateCatalog(bands, [[(40, 40, 10)]])
244 psfArray, newCoadd = mes.utils.computeNearestPsfMultiBand(
245 mCoadd, Point2D(25, 25), catalog
246 )
248 self.assertTupleEqual(newCoadd.bands, bands)
249 # Both bands re-sampled at (40, 40): g matches σ=1.0, and r
250 # matches σ=1.5 (the (40, 40) kernel), not σ=1.2 (the kernel r
251 # held at the center). Failing this assertion would mean either
252 # the upgrade did not run or r kept its center PSF.
253 arr = np.asarray(psfArray)
254 expected_g = g_at_40.computeKernelImage(Point2D(40, 40)).array
255 expected_r = r_at_40.computeKernelImage(Point2D(40, 40)).array
256 self.assertAlmostEqual(arr[0].max(), expected_g.max(), places=4)
257 self.assertAlmostEqual(arr[1].max(), expected_r.max(), places=4)
259 def test_computeNearestPsfMultiBand_falls_back_at_two_locations(self):
260 """When a band needs to fall back but the fallback location is
261 invalid for a band that succeeded at the center, the successful
262 band keeps its center PSF — it is not dropped.
264 Requested center is ``(25, 25)``. Band g's PSF is valid only at
265 ``(40, 40)``; band r's PSF is valid only at ``(25, 25)``. The
266 catalog has a single peak at ``(40, 40)``, so g falls back there.
267 ``(40, 40)`` is invalid for r, so r stays at the center — both
268 bands are kept. Under the bug, band g's fallback shifted r's
269 search center to ``(40, 40)``; r's direct compute at
270 ``(40, 40)`` then failed and the only catalog peak also failed
271 for r, so r was silently dropped from the multiband PSF.
272 """
273 bands = ("g", "r")
274 g_inner = GaussianPsf(41, 41, 1.0)
275 r_inner = GaussianPsf(41, 41, 1.2)
276 psfs = [
277 BadPsf(Point2D(40, 40), g_inner),
278 BadPsf(Point2D(25, 25), r_inner),
279 ]
280 mCoadd = self._generateMultibandCoadd(psfs, bands)
281 catalog = self._generateCatalog(bands, [[(40, 40, 10)]])
283 psfArray, newCoadd = mes.utils.computeNearestPsfMultiBand(
284 mCoadd, Point2D(25, 25), catalog
285 )
287 self.assertTupleEqual(newCoadd.bands, bands)
288 # g lands at (40, 40) (σ=1.0); r stays at the requested center
289 # (σ=1.2). The (40, 40) fallback is *not* valid for r, so the
290 # upgrade-to-common path is correctly skipped.
291 arr = np.asarray(psfArray)
292 expected_g = g_inner.computeKernelImage(Point2D(40, 40)).array
293 expected_r = r_inner.computeKernelImage(Point2D(25, 25)).array
294 self.assertAlmostEqual(arr[0].max(), expected_g.max(), places=4)
295 self.assertAlmostEqual(arr[1].max(), expected_r.max(), places=4)
297 def test_computePsfKernelImage_catalog_emits_future_warning(self):
298 """Passing the deprecated ``catalog`` argument to
299 ``computePsfKernelImage`` emits a ``FutureWarning``.
301 Per finding U-7 of the ``audits/audit-2026-05-05.md`` audit,
302 ``catalog`` is a dead argument -- the body never references it.
303 Rather than remove the parameter (which would silently break
304 any external caller passing it positionally or by keyword) the
305 fix deprecates it and steers callers toward
306 ``computeNearestPsfMultiBand`` for nearest-PSF fallback. The
307 warning lets users find and remove the dead-arg call site
308 before the parameter is dropped after v31.
310 The test pins the warning channel (``FutureWarning``) and
311 confirms the function still returns a valid result -- the
312 ``catalog`` argument is ignored, so the output is identical to
313 the ``catalog=None`` path.
314 """
315 bands = tuple("gri")
316 psfs, psfImage = self._generateMultibandPsf([1.0, 1.2, 1.4])
317 mCoadd = self._generateMultibandCoadd(psfs, bands)
318 catalog = self._generateCatalog(bands)
320 with self.assertWarns(FutureWarning):
321 psfArray, newCoadd = mes.utils.computePsfKernelImage(
322 mCoadd, Point2D(25, 25), catalog=catalog,
323 )
325 # The catalog argument is ignored, so the output matches the
326 # catalog=None / catalog-absent path exactly.
327 np.testing.assert_array_equal(psfArray, psfImage)
328 self.assertTupleEqual(newCoadd.bands, bands)
330 def test_buildObservation_no_divide_warning_on_zero_variance(self):
331 """``buildObservation`` does not emit numpy ``RuntimeWarning``
332 when the input variance plane contains zeros.
334 Per finding U-5 of the ``audits/audit-2026-05-05.md`` audit,
335 the inverse-variance weights were computed as
336 ``weights = 1 / mExposure.variance.array`` without an
337 ``errstate`` guard. Any zero pixel in the variance plane
338 produced ``RuntimeWarning: divide by zero encountered in
339 divide``, and any non-finite pixel produced
340 ``RuntimeWarning: invalid value encountered in divide``. The
341 warnings are spurious -- the immediately following
342 ``weights[~np.isfinite(weights)] = 0`` line replaces every
343 offending value with the intended sentinel -- but they pollute
344 production logs and look like real numerical problems. The
345 fix suppresses the spurious warnings via ``np.errstate``.
347 The fixture coadd's variance plane defaults to all zeros,
348 which under the bug fires the warning on every pixel; the
349 modelPsf and a valid per-band PSF satisfy
350 ``buildObservation``'s preconditions so the function runs all
351 the way through and the test exercises both the divide site
352 and the downstream weight-zeroing.
353 """
354 modelPsf = scl.utils.integrated_circular_gaussian(sigma=0.8).astype(np.float32)
355 bands = tuple("gri")
356 psfs, _ = self._generateMultibandPsf([1.0, 1.2, 1.4])
357 mCoadd = self._generateMultibandCoadd(psfs, bands)
359 with warnings.catch_warnings():
360 warnings.simplefilter("error", RuntimeWarning)
361 observation = mes.utils.buildObservation(
362 modelPsf, Point2I(25, 25), mCoadd
363 )
365 # Sanity check that the call actually went through the
366 # divide-by-zero path: every weight should have been zeroed.
367 np.testing.assert_array_equal(
368 observation.weights, np.zeros_like(observation.weights)
369 )
371 def test_buildObservationBadPsfs(self):
372 # Test that creating an observation with all bad PSFs
373 # raises NoWorkFound
374 modelPsf = scl.utils.integrated_circular_gaussian(sigma=0.8).astype(np.float32)
375 bands = tuple("gri")
376 psfs, psfImage = self._generateMultibandPsf([1.0, 1.2, 1.4])
377 psfs = [BadPsf(Point2D(1, 1), psf) for psf in psfs]
378 mCoadd = self._generateMultibandCoadd(psfs, bands)
379 catalog = self._generateCatalog(self.bands)
381 # Test that building the observation fails without a catalog
382 with self.assertRaises(NoWorkFound):
383 mes.utils.buildObservation(modelPsf, Point2I(25, 25), mCoadd)
385 # Test that building the observation fails even with a catalog
386 with self.assertRaises(NoWorkFound):
387 mes.utils.buildObservation(modelPsf, Point2I(25, 25), mCoadd, catalog=catalog)
389 def _generateGoodPsf(self, sigma: float = 1.0):
390 # Generate a PSF and Image of the PSF
391 psfRadius = 20
392 psfShape = (2 * psfRadius + 1, 2 * psfRadius + 1)
393 psf = GaussianPsf(psfShape[1], psfShape[0], sigma)
394 psfImage = psf.computeImage(psf.getAveragePosition()).array
395 return psf, psfImage
397 def _generateMultibandPsf(self, sigmas: list[float]):
398 # Generate a multiband PSF with a BadPsf for each None value in sigmas
399 psfs = []
400 psfImages = []
401 for sigma in sigmas:
402 psf, psfImage = self._generateGoodPsf(sigma)
403 psfs.append(psf)
404 psfImages.append(psfImage)
405 return psfs, np.asarray(psfImages)
407 def _generateCoadd(self, psf: Psf):
408 # Create an empty exposure
409 masked_image = afwImage.MaskedImage(Extent2I(50, 50), dtype=np.float32)
410 coadd = afwImage.Exposure(masked_image, dtype=np.float32)
411 coadd.setPsf(psf)
412 return coadd
414 def _generateMultibandCoadd(self, psfs: Psf, bands: list[str]):
415 # Create an empty multi-band exposure
416 coadds = []
417 for psf in psfs:
418 coadds.append(self._generateCoadd(psf))
419 return afwImage.MultibandExposure.fromExposures(bands, coadds)
421 def _generateCatalog(self, bands, footprints: list[list[tuple[int, int, int]]] | None = None):
422 # Generate a catalog with a source for each footprint
423 if footprints is None:
424 footprints = []
425 schema = SourceTable.makeMinimalSchema()
426 peakSchema = PeakTable.makeMinimalSchema()
427 for band in bands:
428 schema.addField(f"merge_footprint_{band}", type="Flag")
429 peakSchema.addField(f"merge_peak_{band}", type="Flag")
431 table = SourceTable.make(schema)
432 catalog = SourceCatalog(table)
434 for peaks in footprints:
435 src = catalog.addNew()
436 footprint = Footprint(SpanSet(), peakSchema)
437 for peak in peaks:
438 footprint.addPeak(*peak)
439 src.setFootprint(footprint)
441 for band in bands:
442 src[f"merge_footprint_{band}"] = True
443 footprint.peaks[f"merge_peak_{band}"] = True
444 return catalog
447class TestNonzeroBandSupport(lsst.utils.tests.TestCase):
448 """Tests for ``nonzeroBandSupport`` in
449 ``lsst.meas.extensions.scarlet.utils``.
451 The helper consolidates three previously inconsistent idioms for
452 "this pixel is in the source's support across bands" (``> 0``,
453 ``np.max != 0``, ``np.any != 0``) into a single canonical
454 ``np.any(data != 0, axis=0)``. The discriminator between the
455 canonical form and the historical idioms is a pixel whose band
456 values are all zero except for a negative entry, or a mix of
457 negative and zero (which ``np.max != 0`` excludes when the
458 largest value is exactly zero). Regression coverage for
459 finding U-2 of the ``audits/audit-2026-05-05.md`` audit.
460 """
462 def test_nonzeroBandSupport_includes_negative_only_pixels(self):
463 """A pixel that is negative in some bands and zero in others
464 counts as in the support.
466 Layout of the 2-band ``(2, 3, 3)`` input:
468 - ``(0, 0)``: ``[+1, 0]`` — positive, in support.
469 - ``(0, 1)``: ``[-1, 0]`` — negative-and-zero mix (``max == 0``
470 excludes this; the canonical helper includes it).
471 - ``(0, 2)``: ``[0, 0]`` — all-zero, not in support.
472 - ``(1, 0)``: ``[-1, -1]`` — uniformly negative; ``> 0``
473 excludes, the canonical helper includes.
474 - All other pixels zero.
475 """
476 data = np.zeros((2, 3, 3), dtype=np.float32)
477 data[0, 0, 0] = 1.0
478 data[0, 0, 1] = -1.0
479 data[:, 1, 0] = -1.0
481 result = mes.utils.nonzeroBandSupport(data)
483 expected = np.array(
484 [
485 [True, True, False],
486 [True, False, False],
487 [False, False, False],
488 ]
489 )
490 np.testing.assert_array_equal(result, expected)
492 def test_nonzeroBandSupport_all_zero(self):
493 """An all-zero cube returns an all-False support mask."""
494 data = np.zeros((3, 4, 4), dtype=np.float32)
495 result = mes.utils.nonzeroBandSupport(data)
496 np.testing.assert_array_equal(result, np.zeros((4, 4), dtype=bool))
498 def test_nonzeroBandSupport_all_positive(self):
499 """A strictly-positive cube returns an all-True support mask."""
500 data = np.ones((3, 2, 2), dtype=np.float32)
501 result = mes.utils.nonzeroBandSupport(data)
502 np.testing.assert_array_equal(result, np.ones((2, 2), dtype=bool))
504 def test_nonzeroBandSupport_single_band(self):
505 """A single-band cube reduces along the band axis cleanly.
507 The historical ``np.max != 0`` form silently failed for a
508 ``(1, h, w)`` slice whose only non-zero pixel was negative —
509 the test pixel at ``(0, 1)`` distinguishes ``!= 0`` from
510 ``> 0`` even with only one band.
511 """
512 data = np.array(
513 [[[0.0, -1.0], [2.0, 0.0]]], dtype=np.float32
514 )
515 result = mes.utils.nonzeroBandSupport(data)
516 np.testing.assert_array_equal(
517 result, np.array([[False, True], [True, False]])
518 )
521class TestMultibandConvolve(lsst.utils.tests.TestCase):
522 """Tests for ``multiband_convolve`` in
523 ``lsst.meas.extensions.scarlet.utils``.
525 ``multiband_convolve`` iterates over ``zip(images, psfs, strict=True)``
526 and calls ``scipy.signal.convolve(..., mode="same")`` per band. Both
527 arguments must be 3-D ``(bands, h, w)`` — the function does *not*
528 broadcast a 2-D PSF across bands; the caller is responsible for
529 that (see ``tests/utils.py::DeblenderTestModel.render``).
530 """
532 def test_multiband_convolve_per_band_psf(self):
533 """Each band is convolved with its own PSF.
535 Passes three distinct Gaussian PSFs (sigma = 0.8, 1.2, 1.6) and
536 verifies that ``result[b]`` equals
537 ``scipy.signal.convolve(images[b], psfs[b], mode="same")``
538 computed independently for each band. A regression that
539 cross-routed bands (e.g. always using ``psfs[0]``) would fail.
540 """
541 rng = np.random.RandomState(0)
542 images = rng.rand(3, 21, 21).astype(np.float32)
543 psfs = np.stack([
544 scl.utils.integrated_circular_gaussian(sigma=s).astype(np.float32)
545 for s in (0.8, 1.2, 1.6)
546 ])
548 result = mes.utils.multiband_convolve(images, psfs)
550 self.assertEqual(result.shape, images.shape)
551 for b in range(3):
552 expected = scipy.signal.convolve(images[b], psfs[b], mode="same")
553 np.testing.assert_allclose(result[b], expected, atol=1e-6)
555 def test_multiband_convolve_identity_psf(self):
556 """A centered delta PSF returns the input unchanged.
558 Pins the ``mode="same"`` contract: with a 3×3 PSF that is zero
559 everywhere except a 1 at the center, the per-band convolution
560 is an identity transformation. Any shape or centering bug in
561 the wrapper would shift or truncate the output.
562 """
563 rng = np.random.RandomState(1)
564 images = rng.rand(3, 11, 11).astype(np.float32)
565 psfs = np.zeros((3, 3, 3), dtype=np.float32)
566 psfs[:, 1, 1] = 1.0
568 result = mes.utils.multiband_convolve(images, psfs)
570 np.testing.assert_allclose(result, images, atol=1e-6)
572 def test_multiband_convolve_shape_mismatch_raises(self):
573 """Mismatched band counts raise ``ValueError``.
575 Pins the ``zip(images, psfs, strict=True)`` contract; a
576 regression that drops ``strict=True`` would silently broadcast
577 or truncate.
578 """
579 images = np.zeros((3, 11, 11), dtype=np.float32)
580 psfs = np.zeros((2, 5, 5), dtype=np.float32)
582 with self.assertRaises(ValueError):
583 mes.utils.multiband_convolve(images, psfs)
586def setup_module(module):
587 lsst.utils.tests.init()
590class MemoryTester(lsst.utils.tests.MemoryTestCase):
591 pass
594if __name__ == "__main__": 594 ↛ 595line 594 didn't jump to line 595 because the condition on line 594 was never true
595 lsst.utils.tests.init()
596 unittest.main()