Coverage for tests/test_deconvolve_task.py: 98%
143 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-15 03:23 +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 ``DeconvolveExposureTask``.
24Compares the per-band deconvolved exposures produced by the task
25against the truth ``deconvolved`` image (the model rendered with the
26narrow model PSF only). The two existing tests cover the default
27``useFootprints=True`` path and the catalog-free
28``useFootprints=False`` path; both consume the cached
29``multi-blend`` scene from ``pipeline.py``.
30"""
32import unittest
33import warnings
35import lsst.afw.image as afwImage
36import lsst.geom as geom
37import lsst.meas.extensions.scarlet as mes
38import lsst.scarlet.lite as scl
39import lsst.utils.tests
40import numpy as np
41from lsst.meas.extensions.scarlet.deconvolveExposureTask import (
42 DeconvolveExposureTask,
43 calculateUpdateStep,
44 calculate_update_step,
45)
46from lsst.meas.extensions.scarlet.scarletDeblendTask import ScarletDeblendTask
48import pipeline
49from scenes import SCENES
52class TestDeconvolveTask(lsst.utils.tests.TestCase):
53 """Tests for ``DeconvolveExposureTask`` in
54 ``lsst.meas.extensions.scarlet.deconvolveExposureTask``.
56 Both tests run the deconvolve task on the cached ``multi-blend``
57 scene and compare its output to the truth ``deconvolved`` image
58 (the model convolved with the narrow model PSF only). The
59 assertions ignore a 3×3 region around each source center because
60 Sersic models have sharp peaks that the deconvolver does not
61 recover bit-exactly.
62 """
64 def test_default_deconvolve(self):
65 image = pipeline.build_image(SCENES["multi-blend"])
66 detection = pipeline.detect(image)
67 deconv = pipeline.deconvolve(detection)
69 diff = image.deconvolved.data - deconv.mDeconvolved.image.array
70 # Due to peakiness of Sersic models the center has a sharp peak,
71 # so we ignore a 3x3 region around each source center
72 for model in SCENES["multi-blend"].models:
73 yc, xc = model.center
74 for x in (-1, 0, 1):
75 for y in (-1, 0, 1):
76 diff[:, yc+y, xc+x] = 0
77 self.assertTrue(np.max(diff[:2]) < 10*np.std(image.noise))
78 self.assertTrue(np.max(diff[2]) < 20*np.std(image.noise))
80 context = mes.scarletDeblendTask.ScarletDeblendContext.build(
81 image.mCoadd,
82 deconv.mDeconvolved,
83 detection.catalog,
84 ScarletDeblendTask.ConfigClass()
85 )
87 self.assertEqual(len(context.footprints), 4)
89 def test_catalog_free_deconvolve(self):
90 config = DeconvolveExposureTask.ConfigClass()
91 config.useFootprints = False
92 image = pipeline.build_image(SCENES["multi-blend"])
93 detection = pipeline.detect(image)
94 deconv = pipeline.deconvolve(detection, config=config)
96 diff = image.deconvolved.data - deconv.mDeconvolved.image.array
97 # Due to peakiness of Sersic models the center has a sharp peak,
98 # so we ignore a 3x3 region around each source center
99 for model in SCENES["multi-blend"].models:
100 yc, xc = model.center
101 for x in (-1, 0, 1):
102 for y in (-1, 0, 1):
103 diff[:, yc+y, xc+x] = 0
104 self.assertTrue(np.max(diff[:2]) < 10*np.std(image.noise))
105 self.assertTrue(np.max(diff[2]) < 20*np.std(image.noise))
107 def test_deconvolve_one_isolated_psf(self):
108 """Deconvolving a single isolated PSF source conserves flux.
110 Integrated flux in a 5×5 box around the source center matches
111 the truth within 2σ of expected sum-of-noise (≈ √25·σ) in
112 every band. The deconvolver redistributes flux between
113 adjacent pixels — per-pixel peak amplitude can drift by
114 several σ — but the total flux over the source's support is
115 preserved. See audit finding DC-1 for the per-pixel
116 amplitude residual.
117 """
118 scene = SCENES["one_isolated_psf"]
119 image = pipeline.build_image(scene)
120 detection = pipeline.detect(image)
121 deconv = pipeline.deconvolve(detection)
123 # The truth image's ``yx0`` is the offset of the local array
124 # within the scene's global coordinate frame, so model centers
125 # must be translated before indexing into the array.
126 yx0_y, yx0_x = image.deconvolved.yx0
127 expected_noise = np.sqrt(25) * np.std(image.noise)
128 for model in scene.models:
129 yc, xc = model.center
130 lyc, lxc = yc - yx0_y, xc - yx0_x
131 box = (slice(None), slice(lyc - 2, lyc + 3), slice(lxc - 2, lxc + 3))
132 truth_flux = image.deconvolved.data[box].sum(axis=(1, 2))
133 recovered_flux = deconv.mDeconvolved.image.array[box].sum(axis=(1, 2))
134 for b, band in enumerate(image.bands):
135 self.assertLess(
136 abs(recovered_flux[b] - truth_flux[b]),
137 2 * expected_noise,
138 f"{band} band: flux not conserved "
139 f"(diff {recovered_flux[b] - truth_flux[b]:.4f}, "
140 f"limit {2 * expected_noise:.4f})",
141 )
143 def test_deconvolve_preserves_image_metadata(self):
144 """Deconvolved output preserves bbox, PSF, and WCS from input.
146 For each band, the deconvolved exposure has the same bounding
147 box as its input coadd, the same WCS (deconvolution is per-pixel,
148 no geometric change), and a PSF whose kernel image matches the
149 input PSF's (the task does not synthesize a new PSF).
150 """
151 image = pipeline.build_image(SCENES["multi-blend"])
152 detection = pipeline.detect(image)
153 deconv = pipeline.deconvolve(detection)
155 for band in image.bands:
156 in_exp = image.mCoadd[band]
157 out_exp = deconv.mDeconvolved[band]
158 self.assertEqual(out_exp.getBBox(), in_exp.getBBox())
159 self.assertEqual(out_exp.getWcs(), in_exp.getWcs())
161 out_psf = out_exp.getPsf()
162 self.assertIsNotNone(out_psf)
163 in_psf = in_exp.getPsf()
164 np.testing.assert_array_equal(
165 out_psf.computeImage(out_psf.getAveragePosition()).array,
166 in_psf.computeImage(in_psf.getAveragePosition()).array,
167 )
169 def test_deconvolve_breaks_on_nonfinite_residual(self):
170 """The deconvolution loop stops early when every residual
171 pixel is non-finite rather than running every iteration to
172 ``maxIter`` on meaningless data.
174 The original loop computed ``loss = -0.5 * np.sum(residual**2)``
175 with no NaN guard, so a single NaN in ``residual`` poisoned
176 every subsequent ``loss`` entry; both the convergence test
177 ``np.abs(loss[-1] - loss[-2]) < eRel * np.abs(loss[-1])`` and
178 the divergence test ``loss[-1] < loss[-2]`` return ``False``
179 for NaN, so neither convergence nor step-halving triggered
180 and the loop ran to ``maxIter`` on garbage. The fix moves to
181 ``np.nansum`` for partial-NaN robustness and breaks the loop
182 when ``residual`` is entirely non-finite.
184 The production ``_buildObservation`` sanitizes both
185 ``images`` and ``weights``, so this test bypasses it and
186 constructs an ``scl.Observation`` with all-NaN images
187 directly. It stands in for any mid-iteration scenario where
188 ``convolve`` produces NaN everywhere (numerical artifacts,
189 degenerate PSF). Regression test for finding C-11 of the
190 ``audits/audit-2026-05-05.md`` audit.
191 """
192 config = DeconvolveExposureTask.ConfigClass()
193 config.maxIter = 20
194 config.minIter = 0
195 task = DeconvolveExposureTask(config=config)
197 shape = (1, 8, 8)
198 psf = scl.utils.integrated_circular_gaussian(sigma=0.8).astype(
199 np.float32
200 )
201 observation = scl.Observation(
202 images=np.full(shape, np.nan, dtype=np.float32),
203 variance=np.ones(shape, dtype=np.float32),
204 weights=np.ones(shape, dtype=np.float32),
205 psfs=psf[None],
206 model_psf=psf[None],
207 bands=("dummy",),
208 convolution_mode="fft",
209 )
211 _, loss = task._deconvolve(observation)
213 self.assertLess(len(loss), config.maxIter)
214 self.assertFalse(np.isfinite(loss[-1]))
216 def test_calculate_update_step_excludes_masked_pixels(self):
217 """``calculateUpdateStep`` divides by the count of unmasked
218 pixels rather than the full image size.
220 The previous implementation computed ``sparsity =
221 np.sum(signal_mask) / image.size``; the denominator counted
222 every pixel in the array even when many of them carried zero
223 weight (border, NO_DATA, BAD). On heavily masked inputs such
224 as tract edges this artificially shrinks ``sparsity`` and in
225 turn the update step. The fix restricts both numerator and
226 denominator to pixels with non-zero weight, so the sparsity
227 reflects the fraction of *valid* pixels carrying signal.
229 Two observations are built that differ only in their weight
230 plane: ``full`` has weights ``1`` everywhere; ``half`` masks
231 the bottom half of the image (which contains no signal). A
232 signal-amplitude/noise pair is chosen so the resulting scale
233 does not saturate at the ``1.0`` cap. Under the previous
234 formula both observations yielded the same step (the denominator
235 ignored the mask); under the fix the masked observation yields
236 a step that is roughly twice as large because the denominator
237 halves while the signal count is preserved.
239 Regression test for finding DC-8 of the
240 ``audits/audit-2026-05-05.md`` audit.
241 """
242 shape = (1, 32, 32)
243 noise = 1.0
244 image = np.zeros(shape, dtype=np.float32)
245 # 16-pixel signal block in the top-left; amplitude tuned so
246 # ``scale = sparsity * sqrt(snr) / 0.1`` lands well below 1.0
247 # in the unmasked case.
248 image[0, :4, :4] = 5.0
249 variance = np.full(shape, noise**2, dtype=np.float32)
250 psf = scl.utils.integrated_circular_gaussian(sigma=0.8).astype(np.float32)
252 full_weights = np.ones(shape, dtype=np.float32)
253 half_weights = np.ones(shape, dtype=np.float32)
254 half_weights[:, 16:, :] = 0
256 def _make_obs(weights):
257 return scl.Observation(
258 images=image,
259 variance=variance,
260 weights=weights,
261 psfs=psf[None],
262 model_psf=psf[None],
263 bands=("dummy",),
264 convolution_mode="fft",
265 )
267 step_full = calculateUpdateStep(_make_obs(full_weights))
268 step_half = calculateUpdateStep(_make_obs(half_weights))
270 self.assertLess(step_full, 1.0)
271 self.assertGreater(step_half, step_full)
272 # The signal count is preserved across both observations and
273 # the masked denominator is exactly half the full denominator,
274 # so the masked step should be ~2× larger when neither caps.
275 self.assertAlmostEqual(step_half / step_full, 2.0, places=5)
277 def test_calculate_update_step_deprecation_wrapper(self):
278 """The snake_case ``calculate_update_step`` shim emits a
279 ``FutureWarning`` and forwards to ``calculateUpdateStep``.
281 The function was renamed to match the surrounding LSST
282 camelCase style; the legacy name is retained as a thin
283 deprecation wrapper so external callers continue to work for
284 one release.
286 Regression test for finding DC-10 of the
287 ``audits/audit-2026-05-05.md`` audit.
288 """
289 shape = (1, 8, 8)
290 psf = scl.utils.integrated_circular_gaussian(sigma=0.8).astype(np.float32)
291 observation = scl.Observation(
292 images=np.ones(shape, dtype=np.float32),
293 variance=np.ones(shape, dtype=np.float32),
294 weights=np.ones(shape, dtype=np.float32),
295 psfs=psf[None],
296 model_psf=psf[None],
297 bands=("dummy",),
298 convolution_mode="fft",
299 )
301 expected = calculateUpdateStep(observation)
302 with warnings.catch_warnings(record=True) as caught:
303 warnings.simplefilter("always")
304 actual = calculate_update_step(observation)
306 self.assertEqual(actual, expected)
307 deprecation_warnings = [
308 w for w in caught if issubclass(w.category, FutureWarning)
309 ]
310 self.assertEqual(len(deprecation_warnings), 1)
311 self.assertIn("calculateUpdateStep", str(deprecation_warnings[0].message))
313 def test_model_to_exposure_decouples_mask_and_variance(self):
314 """``_modelToExposure`` detaches the output mask/variance from
315 the input coadd and invalidates the variance plane.
317 Convolution-then-deconvolution changes the per-pixel noise
318 covariance, so the input coadd's variance plane no longer
319 corresponds to the pixel values of the deconvolved model;
320 propagating it unchanged would advertise an incorrect variance
321 as if it were valid. The previous implementation also aliased
322 the output's mask and variance to the input coadd's by
323 reference, so any downstream mutation of the deconvolved
324 exposure's mask/variance would silently leak back into the
325 input.
327 The output exposure now carries (a) a deep-copied mask and
328 (b) a fresh ``inf``-filled variance plane signalling "no
329 information about the noise here". Regression test for finding
330 C-12 of the ``audits/audit-2026-05-05.md`` audit.
331 """
332 bbox = geom.Box2I(geom.Point2I(0, 0), geom.Extent2I(16, 16))
333 coadd = afwImage.ExposureF(bbox)
334 coadd.image.array[:] = 1.0
335 coadd.variance.array[:] = 5.0
336 edge_bit = coadd.mask.getPlaneBitMask("EDGE")
337 coadd.mask.array[0, 0] = edge_bit
339 task = DeconvolveExposureTask()
340 model = np.full((16, 16), 2.0, dtype=coadd.image.array.dtype)
341 out = task._modelToExposure(model, coadd)
343 # Pre-existing mask bits survive the copy.
344 self.assertTrue(out.mask.array[0, 0] & edge_bit != 0)
345 # Variance plane is invalidated by filling with inf.
346 np.testing.assert_array_equal(out.variance.array, np.inf)
348 # Mutating the output mask/variance does not affect the input.
349 out.mask.array[5, 5] |= edge_bit
350 out.variance.array[5, 5] = 999.0
351 self.assertEqual(coadd.mask.array[5, 5], 0)
352 self.assertEqual(coadd.variance.array[5, 5], 5.0)
354 def test_deconvolve_with_nan_input(self):
355 """A NaN pixel in the input does not propagate to the
356 deconvolved output.
358 The NaN is inserted at the center of the first source so it
359 falls inside a detected footprint and is actually processed by
360 the deconvolution (rather than being skipped as outside all
361 footprints). Audit area C-11 lives here; if the contract
362 changes, this test moves with it.
363 """
364 image = pipeline.build_image(SCENES["multi-blend"])
365 detection = pipeline.detect(image)
367 # Clone the multiband exposure so the NaN insertion does not
368 # pollute the cached input shared with other tests.
369 cloned = afwImage.MultibandExposure.fromExposures(
370 image.bands,
371 [image.mCoadd[band].clone() for band in image.bands],
372 )
373 yc, xc = SCENES["multi-blend"].models[0].center
374 for band in image.bands:
375 cloned[band].image.array[yc, xc] = np.nan
377 # Run the deconvolve task directly; ``pipeline.deconvolve`` is
378 # cached and would otherwise reuse the un-mutated input.
379 task = DeconvolveExposureTask()
380 for band in image.bands:
381 result = task.run(cloned[band], detection.catalog)
382 self.assertFalse(
383 np.any(np.isnan(result.deconvolved.image.array)),
384 f"NaN propagated to deconvolved output in {band} band",
385 )
388def setup_module(module):
389 lsst.utils.tests.init()
392class MemoryTester(lsst.utils.tests.MemoryTestCase):
393 pass
396if __name__ == "__main__": 396 ↛ 397line 396 didn't jump to line 397 because the condition on line 396 was never true
397 lsst.utils.tests.init()
398 unittest.main()