Coverage for tests/test_imageUtils.py: 96%
243 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-24 10:06 +0000
1# This file is part of analysis_ap.
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 of reading AP images as `lsst.images` types, and of the functions
23that use them.
24"""
26import collections
27import unittest
28import uuid
29from unittest import mock
31import astropy.time
32import astropy.units as u
33import numpy as np
34import PIL
36import lsst.afw.cameraGeom.testUtils as cameraGeomTestUtils
37import lsst.afw.detection as afwDetection
38import lsst.afw.geom as afwGeom
39import lsst.afw.image as afwImage
40import lsst.daf.base as dafBase
41import lsst.geom
42import lsst.utils.tests
43from lsst.analysis.ap import imageQA, nb_utils, plotImageSubtractionCutouts, ppdb
44from lsst.analysis.ap.imageUtils import (
45 _pixel_to_sky, _sky_to_pixel, as_masked_image, as_visit_image, cutout, cutout_box, get_visit_image,
46 pixel_scale,
47)
48from lsst.daf.butler import DatasetNotFoundError, DimensionUniverse, Timespan
49from lsst.images import (
50 Box, DifferenceImage, DifferenceImageTemplateInfo, MaskedImage, NoOverlapError, Polygon, VisitImage,
51)
53SKY_CENTER = lsst.geom.SpherePoint(45.0, 30.0, lsst.geom.degrees)
56def make_exposure(seed=1):
57 """Make a small legacy exposure that `lsst.images` can convert.
59 The detector comes from `~lsst.afw.cameraGeom.testUtils.CameraWrapper`,
60 since conversion needs amplifier raw geometry and a field-angle
61 transform.
62 """
63 detector = list(cameraGeomTestUtils.CameraWrapper().camera)[0]
64 bbox = detector.getBBox()
65 exposure = afwImage.ExposureF(bbox)
66 rng = np.random.default_rng(seed)
67 exposure.image.array[:] = rng.normal(scale=10.0, size=exposure.image.array.shape)
68 exposure.variance.array[:] = 100.0
69 exposure.mask.array[10:20, 10:30] = exposure.mask.getPlaneBitMask("BAD")
70 exposure.mask.array[50:55, 60:62] = exposure.mask.getPlaneBitMask("DETECTED")
71 exposure.setDetector(detector)
72 exposure.setWcs(afwGeom.makeSkyWcs(lsst.geom.Point2D(bbox.getCenter()), SKY_CENTER,
73 afwGeom.makeCdMatrix(0.2*lsst.geom.arcseconds)))
74 exposure.setPsf(afwDetection.GaussianPsf(21, 21, 2.0))
75 exposure.setFilter(afwImage.FilterLabel(band="g", physical="g_test"))
76 exposure.setPhotoCalib(afwImage.PhotoCalib(1.0))
77 exposure.info.setVisitInfo(afwImage.VisitInfo(
78 id=12345, exposureTime=30.0,
79 date=dafBase.DateTime(60815.000174, dafBase.DateTime.MJD, dafBase.DateTime.TAI)))
80 return exposure
83def make_exposure_record():
84 """Make the ``exposure`` dimension record of the test exposure."""
85 begin = astropy.time.Time(60815.0, format="mjd", scale="tai")
86 return DimensionUniverse()["exposure"].RecordClass(
87 instrument="testCam", id=12345, obs_id="test_12345", group="12345", day_obs=20250520,
88 physical_filter="g_test", exposure_time=30.0, seq_num=1, seq_start=1, seq_end=1,
89 can_see_sky=True, azimuth=120.0, zenith_angle=30.0,
90 timespan=Timespan(begin=begin, end=begin + astropy.time.TimeDelta(30.0, format="sec")),
91 )
94FakeRef = collections.namedtuple("FakeRef", ["name", "dataId"])
97class FakeButler:
98 """Butler stand-in that returns fixed objects by dataset name.
100 Component reads (``"<dataset>.<component>"``) return that attribute of
101 the legacy exposure or `lsst.images` image stored under ``<dataset>``.
102 """
104 def __init__(self, datasets, record=None):
105 self.datasets = datasets
106 self.record = record
107 self.record_queries = []
109 def find_dataset(self, name, data_id=None, **kwargs):
110 return FakeRef(name, data_id) if name in self.datasets else None
112 def get(self, target, *args, **kwargs):
113 name = target.name if isinstance(target, FakeRef) else target
114 if "." in name: 114 ↛ 115line 114 didn't jump to line 115 because the condition on line 114 was never true
115 parent, component = name.split(".")
116 dataset = self.datasets[parent]
117 if isinstance(dataset, afwImage.Exposure):
118 return {"detector": dataset.getDetector, "bbox": dataset.getBBox}[component]()
119 return getattr(dataset, component)
120 return self.datasets[name]
122 def query_dimension_records(self, element, data_id=None, explain=True):
123 self.record_queries.append(data_id)
124 return [self.record] if self.record is not None else []
127class AsVisitImageTestCase(lsst.utils.tests.TestCase):
128 def setUp(self):
129 self.exposure = make_exposure()
131 def testNativePassesThrough(self):
132 image = as_visit_image(self.exposure, image_type="difference")
133 self.assertIs(as_visit_image(image), image)
134 self.assertIs(as_masked_image(image), image)
136 def testRoles(self):
137 for image_type, image_class in [("science", VisitImage), ("template", DifferenceImage),
138 ("difference", DifferenceImage)]:
139 with self.subTest(image_type=image_type):
140 image = as_visit_image(self.exposure, image_type=image_type,
141 exposure_record=make_exposure_record())
142 self.assertIs(type(image), image_class)
143 self.assertEqual(image.bbox, Box.from_legacy(self.exposure.getBBox()))
144 np.testing.assert_array_equal(image.image.array, self.exposure.image.array)
145 np.testing.assert_array_equal(image.variance.array, self.exposure.variance.array)
146 bad = self.exposure.mask.array == self.exposure.mask.getPlaneBitMask("BAD")
147 np.testing.assert_array_equal(image.mask.get("BAD"), bad)
148 self.assertEqual(image.band, "g")
150 def testUnits(self):
151 self.assertEqual(as_visit_image(self.exposure).unit, u.nJy)
152 self.exposure.setPhotoCalib(afwImage.PhotoCalib(57.5))
153 image = as_visit_image(self.exposure)
154 self.assertEqual(image.unit, u.electron)
155 self.assertIsNotNone(image.photometric_scaling)
156 self.exposure.metadata["BUNIT"] = "nJy"
157 self.assertEqual(as_visit_image(self.exposure).unit, u.nJy)
158 # The caller's metadata is left alone, so a second conversion works.
159 self.assertEqual(self.exposure.metadata["BUNIT"], "nJy")
160 self.assertEqual(as_visit_image(self.exposure).unit, u.nJy)
162 def testTemplateDetector(self):
163 detector = self.exposure.getDetector()
164 self.exposure.setDetector(None)
165 with self.assertRaisesRegex(ValueError, "no detector"):
166 as_visit_image(self.exposure, image_type="template")
167 native_detector = as_visit_image(make_exposure()).detector
168 for given in (detector, native_detector):
169 image = as_visit_image(self.exposure, image_type="template", detector=given)
170 self.assertEqual(image.detector.to_legacy().getId(), detector.getId())
171 self.assertIsNone(self.exposure.getDetector())
173 def testUnsupported(self):
174 with self.assertRaises(TypeError):
175 as_visit_image(self.exposure.maskedImage)
176 with self.assertRaises(TypeError):
177 as_masked_image(self.exposure.maskedImage)
179 def testMaskedImage(self):
180 # A post-ISR image has no PSF or WCS.
181 self.exposure.setPsf(None)
182 self.exposure.setWcs(None)
183 image = as_masked_image(self.exposure)
184 self.assertIs(type(image), MaskedImage)
185 np.testing.assert_array_equal(image.image.array, self.exposure.image.array)
186 self.assertEqual(int(image.mask.get("BAD").sum()), 200)
189class GetVisitImageTestCase(lsst.utils.tests.TestCase):
190 def setUp(self):
191 self.exposure = make_exposure()
192 self.record = make_exposure_record()
194 def testLegacy(self):
195 butler = FakeButler({"difference_image": self.exposure}, record=self.record)
196 image = get_visit_image(butler, "difference_image", {"visit": 12345}, image_type="difference")
197 self.assertIsInstance(image, DifferenceImage)
198 self.assertEqual(butler.record_queries, [{"visit": 12345}])
199 self.assertEqual(image.obs_info.exposure_time, 30.0*u.s)
201 def testNative(self):
202 native = as_visit_image(self.exposure, image_type="difference")
203 butler = FakeButler({"difference_image": native})
204 self.assertIs(get_visit_image(butler, "difference_image", {"visit": 1}), native)
205 self.assertEqual(butler.record_queries, [])
207 def testMissing(self):
208 with self.assertRaises(DatasetNotFoundError):
209 get_visit_image(FakeButler({}), "difference_image", {"visit": 1})
212class CutoutTestCase(lsst.utils.tests.TestCase):
213 def setUp(self):
214 self.exposure = make_exposure()
215 self.image = as_visit_image(self.exposure)
217 def testMatchesLegacyCutout(self):
218 # Inside the image, and across a corner, where the cutout is padded.
219 for x, y in [(100.3, 80.7), (101.0, 60.0), (1.2, 2.6)]:
220 with self.subTest(x=x, y=y):
221 center = self.exposure.wcs.pixelToSky(x, y)
222 legacy = self.exposure.getCutout(center, lsst.geom.Extent2I(9, 9))
223 stamp = cutout(self.image, center, 9)
224 self.assertEqual(stamp.bbox, Box.from_legacy(legacy.getBBox()))
225 np.testing.assert_array_equal(stamp.array, legacy.image.array)
227 def testOffImage(self):
228 far = self.exposure.wcs.pixelToSky(-500.0, -500.0)
229 with self.assertRaises(NoOverlapError):
230 cutout(self.image, far, 9)
232 def testCutoutBox(self):
233 box = Box.from_legacy(lsst.geom.Box2I(lsst.geom.Point2I(5, 7), lsst.geom.Extent2I(10, 12)))
234 np.testing.assert_array_equal(cutout_box(self.image, box).array,
235 self.exposure.image[box.to_legacy()].array)
237 def testPixelScale(self):
238 expected = self.exposure.wcs.getPixelScale(self.exposure.getBBox().getCenter()).asArcseconds()
239 self.assertAlmostEqual(pixel_scale(self.image), expected, places=9)
241 def testSkyPixelTransforms(self):
242 x, y = np.array([3.0, 150.5]), np.array([9.0, 40.25])
243 ra, dec = _pixel_to_sky(self.image.sky_projection, x, y)
244 np.testing.assert_allclose((ra, dec), _pixel_to_sky(self.exposure.wcs, x, y), rtol=0, atol=1e-10)
245 for wcs in (self.image.sky_projection, self.exposure.wcs):
246 np.testing.assert_allclose(_sky_to_pixel(wcs, ra, dec), (x, y), rtol=0, atol=1e-6)
249class NativeImageFunctionsTestCase(lsst.utils.tests.TestCase):
250 """Test that functions give the same result for a legacy image and its
251 `lsst.images` form.
252 """
254 def setUp(self):
255 self.legacy = make_exposure()
256 self.native = as_visit_image(self.legacy, image_type="difference")
258 def testImageDiffStats(self):
259 legacy = imageQA.image_diff_stats(FakeButler({"difference_image": self.legacy}), 1, 2)
260 native = imageQA.image_diff_stats(FakeButler({"difference_image": self.native}), 1, 2)
261 self.assertEqual(native, legacy)
262 n_pixels = self.legacy.image.array.size
263 self.assertEqual(legacy["frac_BAD"], 200/n_pixels)
264 self.assertEqual(legacy["frac_DETECTED"], 10/n_pixels)
266 def testPixelCompare(self):
267 result = imageQA.pixel_compare(FakeButler({"difference_image": self.legacy}),
268 FakeButler({"difference_image": self.native}), 1, 2)
269 self.assertEqual(result.summary["diff_stddev"], 0.0)
270 self.assertEqual(result.summary["n_mask_pixels_changed"], 0)
272 def testMaskDifference(self):
273 other = make_exposure()
274 other.mask.array[0, 0] |= other.mask.getPlaneBitMask("SAT")
275 other = as_visit_image(other)
276 differ = imageQA._mask_difference(self.native.mask, other.mask)
277 self.assertEqual(int(differ.sum()), 1)
278 self.assertTrue(differ[0, 0])
280 def testRegionFromExposure(self):
281 bbox = self.legacy.getBBox()
282 center = self.legacy.wcs.pixelToSky(bbox.getCenter())
283 radius = max(center.separation(self.legacy.wcs.pixelToSky(lsst.geom.Point2D(corner)))
284 for corner in bbox.getCorners()) + 5*lsst.geom.arcseconds
285 expected = (center.getRa().asDegrees(), center.getDec().asDegrees(), radius.asDegrees())
286 for image in (self.native, self.legacy):
287 np.testing.assert_allclose(ppdb.region_from_exposure(image, padding=5.0), expected,
288 rtol=0, atol=1e-10)
290 def testGenerateImage(self):
291 task = plotImageSubtractionCutouts.PlotImageSubtractionCutoutsTask(output_path="")
292 legacy_template = make_exposure(seed=2)
293 legacy_template.setDetector(None)
294 # A legacy template is given the science image's detector.
295 stamp = task.generate_image(self.legacy, legacy_template, self.native, SKY_CENTER, 0.2)
296 with PIL.Image.open(stamp) as image:
297 self.assertGreater(image.width, 0)
299 def testButlerCache(self):
300 config = plotImageSubtractionCutouts.PlotImageSubtractionCutoutsConfig()
301 config.diff_image_type = None
302 legacy_template = make_exposure(seed=2)
303 legacy_template.setDetector(None)
304 butler = FakeButler({"preliminary_visit_image": self.legacy,
305 "template_detector": legacy_template,
306 "difference_image": self.native})
307 plotImageSubtractionCutouts.butler_cache.set(butler, config)
308 science, template, difference = plotImageSubtractionCutouts.butler_cache.get_images("testCam", 2, 1)
309 self.assertIs(type(science), VisitImage)
310 self.assertIs(type(template), DifferenceImage)
311 self.assertIs(difference, self.native)
313 def testLoadCutout(self):
314 butler = FakeButler({"preliminary_visit_image": self.legacy, "difference_image": self.native})
315 row = mock.Mock(visit=1, detector=2, ra=SKY_CENTER.getRa().asDegrees(),
316 dec=SKY_CENTER.getDec().asDegrees())
317 data = nb_utils._load_cutout(butler, row, size=15, image_type="difference",
318 image_datasets=nb_utils._IMAGE_DATASETS)
319 legacy = self.legacy.getCutout(SKY_CENTER, lsst.geom.Extent2I(15, 15))
320 np.testing.assert_array_equal(data["data"], legacy.image.array)
321 self.assertEqual((data["x0"], data["y0"]), (legacy.getX0(), legacy.getY0()))
324class TemplatePatchOutlinesTestCase(lsst.utils.tests.TestCase):
325 """Test finding the coadds of a template stored as an
326 `lsst.images.DifferenceImage`.
327 """
329 def setUp(self):
330 exposure = make_exposure()
331 # The coadd of patch 0 shares the difference image's pixel grid.
332 self.coadd_wcs = exposure.wcs
333 self.sky_projection = as_visit_image(exposure).sky_projection
334 self.ids = [uuid.uuid4(), uuid.uuid4()]
335 # Listed out of order, to check that the outlines are sorted.
336 self.templates = [
337 DifferenceImageTemplateInfo(
338 skymap="test", tract=1, patch=patch, dataset_id=dataset_id, dataset_run="run",
339 bounds=Polygon.from_box(Box.from_shape((100, 60), start=(0, x0))),
340 psf_shape_xx=1.0, psf_shape_yy=1.0, psf_shape_xy=0.0, psf_shape_flag=False,
341 )
342 for patch, dataset_id, x0 in [(1, self.ids[1], 40), (0, self.ids[0], 0)]
343 ]
344 self.coadd_bbox = lsst.geom.Box2I(lsst.geom.Point2I(-10, -20), lsst.geom.Extent2I(70, 140))
345 self.ref = mock.MagicMock()
346 self.ref.datasetType.storageClass.allComponents.return_value = {"wcs": None, "bbox": None}
347 self.ref.makeComponentRef.side_effect = lambda component: component
349 def get(target, *args, storageClass=None):
350 return {"template_detector.templates": self.templates, "wcs": self.coadd_wcs,
351 "bbox": self.coadd_bbox}[target]
353 self.butler = mock.MagicMock()
354 self.butler.get_dataset_type.return_value.storageClass.allComponents.return_value = {
355 "templates": None}
356 self.butler.get.side_effect = get
357 # The coadd of patch 1 is missing from the butler.
358 self.butler.get_dataset.side_effect = {self.ids[0]: self.ref, self.ids[1]: None}.get
360 def testOutlines(self):
361 outlines, refs = nb_utils._template_patch_outlines(
362 self.butler, "template_detector", {"visit": 1, "detector": 2}, self.sky_projection,
363 coadd_dataset="unused", band="g")
364 self.assertEqual(list(outlines), [(1, 0), (1, 1)])
365 self.assertEqual(refs, {(1, 0): self.ref})
366 # A found coadd is outlined by its own bounding box.
367 expected = lsst.geom.Box2D(self.coadd_bbox).getCorners()
368 np.testing.assert_allclose(outlines[(1, 0)], [(p.getX(), p.getY()) for p in expected], atol=1e-6)
369 # A missing coadd is outlined by its recorded overlap.
370 corners = outlines[(1, 1)]
371 self.assertEqual(len(corners), 4)
372 self.assertEqual({x for x, _ in corners}, {39.5, 99.5})
373 self.assertEqual({y for _, y in corners}, {-0.5, 99.5})
376class MemoryTester(lsst.utils.tests.MemoryTestCase):
377 pass
380def setup_module(module):
381 lsst.utils.tests.init()
384if __name__ == "__main__": 384 ↛ 385line 384 didn't jump to line 385 because the condition on line 384 was never true
385 lsst.utils.tests.init()
386 unittest.main()