Coverage for tests/test_masked_image.py: 90%
180 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-16 15:24 -0700
1# This file is part of lsst-images.
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# Use of this source code is governed by a 3-clause BSD-style
10# license that can be found in the LICENSE file.
12from __future__ import annotations
14import dataclasses
15import os
16from pathlib import Path
17from typing import Any
19import astropy.io.fits
20import astropy.units as u
21import numpy as np
22import pytest
24from lsst.images import Box, Image, MaskedImage, MaskPlane, MaskSchema, get_legacy_visit_image_mask_planes
25from lsst.images.fits import FitsCompressionOptions
26from lsst.images.tests import (
27 RoundtripFits,
28 RoundtripJson,
29 RoundtripNdf,
30 assert_masked_images_equal,
31 compare_masked_image_to_legacy,
32)
34try:
35 import h5py # noqa: F401
37 HAVE_H5PY = True
38except ImportError:
39 HAVE_H5PY = False
41try:
42 from lsst.afw.image import MaskedImageReader as LegacyMaskedImageReader
44except ImportError:
45 type LegacyMaskedImageReader = Any # type: ignore[no-redef]
47EXTERNAL_DATA_DIR = os.environ.get("TESTDATA_IMAGES_DIR", None)
49skip_no_h5py = pytest.mark.skipif(not HAVE_H5PY, reason="h5py is not installed")
52@dataclasses.dataclass
53class _LegacyTestData:
54 masked_image: MaskedImage
55 reader: LegacyMaskedImageReader
56 plane_map: dict[str, MaskPlane]
59@pytest.fixture(scope="session")
60def legacy_test_data() -> _LegacyTestData:
61 """Return a Mask read directly from the legacy test dataset and a legacy
62 reader for that image.
64 Skips if TESTDATA_IMAGES_DIR is unset or lsst.afw.image is unavailable.
65 """
66 if EXTERNAL_DATA_DIR is None: 66 ↛ 68line 66 didn't jump to line 68 because the condition on line 66 was always true
67 pytest.skip("TESTDATA_IMAGES_DIR is not in the environment.")
68 try:
69 from lsst.afw.image import MaskedImageFitsReader
70 except ImportError:
71 pytest.skip("'lsst.afw.image' could not be imported.")
72 filename = os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", "visit_image.fits")
73 plane_map = get_legacy_visit_image_mask_planes()
74 masked_image = MaskedImage.read_legacy(filename, plane_map=plane_map)
75 reader = MaskedImageFitsReader(filename)
76 return _LegacyTestData(masked_image=masked_image, reader=reader, plane_map=plane_map)
79def make_masked_image() -> MaskedImage:
80 """Return a freshly-constructed MaskedImage with BAD and HUNGRY mask
81 planes set.
82 """
83 rng = np.random.default_rng(500)
84 masked_image = MaskedImage(
85 Image(rng.normal(100.0, 8.0, size=(200, 251)), dtype=np.float64, unit=u.nJy, yx0=(5, 8)),
86 mask_schema=MaskSchema(
87 [
88 MaskPlane("BAD", "Pixel is very bad, possibly downright evil."),
89 MaskPlane("HUNGRY", "Pixel hasn't had enough to eat today."),
90 ]
91 ),
92 metadata={"fifty": "5 * 10"},
93 )
94 masked_image.mask.array |= np.multiply.outer(
95 masked_image.image.array < 102.0,
96 masked_image.mask.schema.bitmask("BAD"),
97 )
98 masked_image.mask.array |= np.multiply.outer(
99 masked_image.image.array > 98.0,
100 masked_image.mask.schema.bitmask("HUNGRY"),
101 )
102 masked_image.variance.array = rng.normal(64.0, 0.5, size=masked_image.bbox.shape)
103 return masked_image
106def test_construction() -> None:
107 """Verify the MaskedImage constructed by make_masked_image has the
108 expected attributes.
109 """
110 mi = make_masked_image()
111 assert mi.bbox == Box.factory[5:205, 8:259]
112 assert mi.mask.bbox == mi.bbox
113 assert mi.variance.bbox == mi.bbox
114 assert mi.image.array.shape == mi.bbox.shape
115 assert mi.mask.array.shape == mi.bbox.shape + (1,)
116 assert mi.variance.array.shape == mi.bbox.shape
117 assert mi.unit == u.nJy
118 assert mi.variance.unit == u.nJy**2
119 assert mi.metadata == {"fifty": "5 * 10"}
120 # The checks below are subject to the vagaries of the RNG, but we want
121 # the seed to be such that they all pass, or other tests will be weaker.
122 assert np.sum(mi.mask.array == mi.mask.schema.bitmask("BAD")) > 0
123 assert np.sum(mi.mask.array == mi.mask.schema.bitmask("HUNGRY")) > 0
124 assert np.sum(mi.mask.array == mi.mask.schema.bitmask("BAD", "HUNGRY")) > 0
126 assert mi[...] is not mi
127 assert str(mi) == "MaskedImage(Image([y=5:205, x=8:259], float64), ['BAD', 'HUNGRY'])"
128 assert (
129 repr(mi)
130 == "MaskedImage(Image(..., bbox=Box(y=Interval(start=5, stop=205), x=Interval(start=8, stop=259)), "
131 "dtype=dtype('float64')), mask_schema=MaskSchema([MaskPlane(name='BAD', description='Pixel is "
132 "very bad, possibly downright evil.'), MaskPlane(name='HUNGRY', description=\"Pixel hasn't had "
133 "enough to eat today.\")], dtype=dtype('uint8')))"
134 )
135 copy = mi.copy()
136 original = mi.image.array[0, 0]
137 copy.image.array[0, 0] = 38.0
138 assert mi.image.array[0, 0] == original
139 assert copy.image.array[0, 0] == 38.0
141 # Test error conditions.
142 with pytest.raises(ValueError):
143 # Disagreement over mask bbox.
144 MaskedImage(Image(42.0, shape=(5, 6)), mask=mi.mask)
145 with pytest.raises(TypeError):
146 # No mask definition.
147 MaskedImage(mi.image, variance=mi.variance)
148 with pytest.raises(TypeError):
149 # Can not provide mask and mask schema.
150 MaskedImage(
151 Image(42.0, shape=(5, 5)),
152 mask=mi.mask,
153 mask_schema=mi.mask.schema,
154 )
155 with pytest.raises(ValueError):
156 # image and variance bbox disagreement.
157 MaskedImage(
158 Image(42.0, shape=(5, 5)),
159 mask_schema=mi.mask.schema,
160 variance=mi.variance,
161 )
162 with pytest.raises(ValueError):
163 # no image unit but there is variance unit.
164 MaskedImage(
165 Image(42.0, shape=(5, 5)),
166 mask_schema=mi.mask.schema,
167 variance=Image(1.0, shape=(5, 5), unit=u.nJy),
168 )
169 with pytest.raises(ValueError):
170 # image and variance units disagree.
171 MaskedImage(
172 Image(42.0, shape=(5, 5), unit=u.nJy),
173 mask_schema=mi.mask.schema,
174 variance=Image(1.0, shape=(5, 5), unit=u.nJy),
175 )
178def test_subset() -> None:
179 """Verify assignment of a subset into a MaskedImage copy."""
180 mi = make_masked_image()
181 copy = mi.copy()
182 subset = copy.local[0:10, 20:30].copy()
183 subset.image[...] = Image(42.0, shape=(10, 10), unit=u.nJy)
184 copy[subset.bbox] = subset
185 assert copy.image.array[0, 20] == 42.0
186 assert copy.image.array[0, 0] == mi.image.array[0, 0]
189def test_mask_setter() -> None:
190 """Verify the mask plane can be replaced with one grown by add_plane."""
191 mi = make_masked_image()
192 bad = mi.mask.get("BAD")
193 mi.mask = mi.mask.add_plane("OUTSIDE_STENCIL", "Pixel lies outside the stencil.")
194 assert "OUTSIDE_STENCIL" in mi.mask.schema.names
195 assert mi.mask.bbox == mi.image.bbox
196 np.testing.assert_array_equal(mi.mask.get("BAD"), bad)
197 assert not mi.mask.get("OUTSIDE_STENCIL").any()
198 # A mask whose bounding box disagrees with the image is rejected.
199 with pytest.raises(ValueError):
200 mi.mask = mi.mask[Box.factory[10:20, 12:22]]
203def test_fits_roundtrip() -> None:
204 """Verify MaskedImage round-trips correctly through FITS, including
205 subimage reads.
206 """
207 mi = make_masked_image()
208 subbox = Box.factory[11:20, 25:30]
209 subslices = (slice(6, 15), slice(17, 22))
210 np.testing.assert_array_equal(mi.image.array[subslices], mi.image[subbox].array)
211 with RoundtripFits(mi, "MaskedImageV2") as roundtrip:
212 subimage = roundtrip.get(bbox=subbox)
213 # Check that we used lossless compression (the default).
214 fits = roundtrip.inspect()
215 assert fits[1].header["ZCMPTYPE"] == "GZIP_2"
216 assert fits[2].header["ZCMPTYPE"] == "GZIP_2"
217 assert fits[3].header["ZCMPTYPE"] == "GZIP_2"
218 assert_masked_images_equal(roundtrip.result, mi, expect_view=False)
219 assert_masked_images_equal(subimage, roundtrip.result[subbox], expect_view=False)
222def test_fits_roundtrip_legacy_read() -> None:
223 """Verify a round-tripped MaskedImageV2 can be read back as a legacy afw
224 MaskedImage.
225 """
226 try:
227 import lsst.afw.image
228 except ImportError:
229 pytest.skip("afw could not be imported")
230 mi = make_masked_image()
231 with RoundtripFits(mi, "MaskedImageV2") as roundtrip:
232 legacy_masked_image = roundtrip.get(storageClass="MaskedImage")
233 assert isinstance(legacy_masked_image, lsst.afw.image.MaskedImage)
234 compare_masked_image_to_legacy(mi, legacy_masked_image, expect_view=False)
237def test_fits_roundtrip_lossy(tmp_path: Path) -> None:
238 """Verify MaskedImage round-trips correctly through FITS with lossy
239 compression.
240 """
241 mi = make_masked_image()
242 subbox = Box.factory[11:20, 25:30]
243 subslices = (slice(6, 15), slice(17, 22))
244 np.testing.assert_array_equal(mi.image.array[subslices], mi.image[subbox].array)
245 path = tmp_path / "lossy.fits"
246 mi.write(
247 path,
248 compression_options={
249 "image": FitsCompressionOptions.LOSSY,
250 "variance": FitsCompressionOptions.LOSSY,
251 },
252 compression_seed=50,
253 )
254 roundtripped = MaskedImage.read(path)
255 subimage = MaskedImage.read(path, bbox=subbox)
256 with astropy.io.fits.open(path, disable_image_compression=True) as fits:
257 assert fits[1].header["ZCMPTYPE"] == "RICE_1"
258 assert fits[2].header["ZCMPTYPE"] == "GZIP_2"
259 assert fits[3].header["ZCMPTYPE"] == "RICE_1"
260 assert_masked_images_equal(roundtripped, mi, expect_view=False, rtol=0.01)
261 assert_masked_images_equal(subimage, roundtripped[subbox], expect_view=False)
264@skip_no_h5py
265def test_round_trip_ndf_compatible_mask() -> None:
266 """Verify NDF round-trip for a MaskedImage with ≤8 mask planes."""
267 mi = make_masked_image()
268 with RoundtripNdf(mi, "MaskedImageV2") as roundtrip:
269 assert_masked_images_equal(roundtrip.result, mi, expect_view=False)
272@skip_no_h5py
273def test_round_trip_ndf_incompatible_mask() -> None:
274 """Verify NDF round-trip for a MaskedImage with more than 8 mask planes."""
275 rng = np.random.default_rng(7)
276 planes = [MaskPlane(f"P{i}", f"plane {i}") for i in range(12)]
277 wide = MaskedImage(
278 Image(
279 rng.normal(100.0, 8.0, size=(50, 60)),
280 dtype=np.float64,
281 unit=u.nJy,
282 yx0=(0, 0),
283 ),
284 mask_schema=MaskSchema(planes),
285 )
286 wide.variance.array = rng.normal(64.0, 0.5, size=wide.bbox.shape)
287 with RoundtripNdf(wide, "MaskedImageV2") as roundtrip:
288 assert_masked_images_equal(roundtrip.result, wide, expect_view=False)
291@skip_no_h5py
292def test_round_trip_ndf_many_plane_mask() -> None:
293 """Verify NDF round-trip for a mask that needs more than one int32
294 chunk.
295 """
296 rng = np.random.default_rng(11)
297 planes = [MaskPlane(f"P{i}", f"plane {i}") for i in range(40)]
298 wide = MaskedImage(
299 Image(
300 rng.normal(100.0, 8.0, size=(10, 12)),
301 dtype=np.float64,
302 unit=u.nJy,
303 yx0=(0, 0),
304 ),
305 mask_schema=MaskSchema(planes),
306 )
307 wide.mask.set("P0", wide.image.array > 100.0)
308 wide.mask.set("P17", wide.image.array < 95.0)
309 wide.mask.set("P39", wide.image.array > 110.0)
310 wide.variance.array = rng.normal(64.0, 0.5, size=wide.bbox.shape)
311 with RoundtripNdf(wide, "MaskedImageV2") as roundtrip:
312 assert_masked_images_equal(roundtrip.result, wide, expect_view=False)
315@skip_no_h5py
316def test_fits_ndf_consistency() -> None:
317 """Verify FITS and NDF backends produce equal MaskedImages on
318 round-trip.
319 """
320 mi = make_masked_image()
321 with (
322 RoundtripFits(mi) as fits_rt,
323 RoundtripNdf(mi) as ndf_rt,
324 ):
325 assert_masked_images_equal(mi, fits_rt.result, expect_view=False)
326 assert_masked_images_equal(mi, ndf_rt.result, expect_view=False)
327 assert_masked_images_equal(fits_rt.result, ndf_rt.result, expect_view=False)
330def test_fits_json_consistency() -> None:
331 """Verify FITS and JSON backends produce equal MaskedImages on
332 round-trip.
333 """
334 mi = make_masked_image()
335 with (
336 RoundtripFits(mi) as fits_rt,
337 RoundtripJson(mi) as json_rt,
338 ):
339 assert_masked_images_equal(mi, fits_rt.result, expect_view=False)
340 assert_masked_images_equal(mi, json_rt.result, expect_view=False)
341 assert_masked_images_equal(fits_rt.result, json_rt.result, expect_view=False)
344def test_legacy(legacy_test_data: _LegacyTestData) -> None:
345 """Test MaskedImage.read_legacy, MaskedImage.to_legacy, and
346 MaskedImage.from_legacy.
347 """
348 legacy_masked_image = legacy_test_data.reader.read()
349 compare_masked_image_to_legacy(
350 legacy_test_data.masked_image,
351 legacy_masked_image,
352 plane_map=legacy_test_data.plane_map,
353 expect_view=False,
354 )
355 compare_masked_image_to_legacy(
356 legacy_test_data.masked_image,
357 legacy_test_data.masked_image.to_legacy(plane_map=legacy_test_data.plane_map),
358 plane_map=legacy_test_data.plane_map,
359 expect_view=True,
360 )
361 compare_masked_image_to_legacy(
362 MaskedImage.from_legacy(legacy_masked_image, plane_map=legacy_test_data.plane_map),
363 legacy_masked_image,
364 expect_view=True,
365 plane_map=legacy_test_data.plane_map,
366 )