Coverage for tests/test_transforms.py: 54%
196 statements
« prev ^ index » next coverage.py v7.14.3, created at 2026-07-04 10:14 -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 functools
16import os
17import unittest
18from typing import Any, ClassVar
20import astropy.units as u
21import numpy as np
22import pydantic
24from lsst.images import (
25 ICRS,
26 XY,
27 Box,
28 CameraFrameSet,
29 CameraFrameSetSerializationModel,
30 DetectorFrame,
31 FocalPlaneFrame,
32 GeneralFrame,
33 SkyProjection,
34 Transform,
35 TransformSerializationModel,
36)
37from lsst.images._transforms import _ast as astshim
38from lsst.images.fits import PointerModel
39from lsst.images.serialization import ArchiveTree, InputArchive, JsonRef, OutputArchive
40from lsst.images.tests import (
41 DP2_VISIT_DETECTOR_DATA_ID,
42 RoundtripFits,
43 RoundtripJson,
44 check_transform,
45 compare_sky_projection_to_legacy_wcs,
46 legacy_points_to_xy_array,
47)
49DATA_DIR = os.environ.get("TESTDATA_IMAGES_DIR", None)
52class TransformTestCase(unittest.TestCase):
53 """Tests for the Transform, SkyProjection, and FrameSet classes."""
55 def test_identity(self) -> None:
56 """Test an identity transform."""
57 frame = DetectorFrame(**DP2_VISIT_DETECTOR_DATA_ID, bbox=Box.factory[:5, :4])
58 xy = frame.bbox.meshgrid().map(np.ravel)
59 identity = Transform.identity(frame)
60 check_transform(self, identity, xy, xy, frame, frame)
61 self.assertEqual(identity.decompose(), [])
62 with RoundtripJson(self, identity) as roundtrip:
63 pass
64 check_transform(self, roundtrip.result, xy, xy, frame, frame)
66 def test_transform_equality(self) -> None:
67 """Test ``Transform.__eq__`` across all of its comparison branches."""
68 pixel_frame = DetectorFrame(**DP2_VISIT_DETECTOR_DATA_ID, bbox=Box.factory[:5, :4])
69 focal_plane = FocalPlaneFrame(instrument="LSSTCam", visit=1, unit=u.mm)
70 # A distinct frame for the in-frame and out-frame branches.
71 alt_frame = DetectorFrame(instrument="LSSTCam", visit=1, detector=12, bbox=Box.factory[:5, :4])
72 in_bounds = Box.factory[:5, :4]
73 out_bounds = Box.factory[:10, :8]
75 def make(
76 *,
77 in_frame: Any = pixel_frame,
78 out_frame: Any = focal_plane,
79 ast_mapping: astshim.Mapping | None = None,
80 in_bounds_: Box | None = in_bounds,
81 out_bounds_: Box | None = out_bounds,
82 components: Any = (),
83 ) -> Transform[Any, Any]:
84 return Transform(
85 in_frame,
86 out_frame,
87 ast_mapping if ast_mapping is not None else astshim.UnitMap(2),
88 in_bounds=in_bounds_,
89 out_bounds=out_bounds_,
90 components=components,
91 )
93 base = make()
95 # Identity short-circuit: an object is always equal to itself.
96 self.assertEqual(base, base)
98 # Two independently constructed but equivalent transforms are equal,
99 # and equality is symmetric.
100 self.assertEqual(base, make())
101 self.assertEqual(make(), base)
103 # Comparison against a non-Transform yields NotImplemented, so Python
104 # falls back to identity: the objects are unequal and ``!=`` is True.
105 self.assertFalse(base == "not a transform")
106 self.assertTrue(base != "not a transform")
107 self.assertNotEqual(base, None)
108 self.assertNotEqual(base, 42)
110 # Each remaining branch differs from ``base`` in exactly one attribute.
111 self.assertNotEqual(base, make(ast_mapping=astshim.ShiftMap([1.0, 2.0])))
112 self.assertNotEqual(base, make(in_bounds_=out_bounds))
113 self.assertNotEqual(base, make(out_bounds_=in_bounds))
114 self.assertNotEqual(base, make(in_frame=alt_frame))
115 self.assertNotEqual(base, make(out_frame=alt_frame))
116 self.assertNotEqual(base, make(components=[Transform.identity(alt_frame)]))
118 def test_sky_projection_equality(self) -> None:
119 """Test ``SkyProjection.__eq__`` across all of its comparison
120 branches.
121 """
122 pixel_frame = DetectorFrame(**DP2_VISIT_DETECTOR_DATA_ID, bbox=Box.factory[:5, :4])
124 # Check the two failure modes.
125 with self.assertRaises(ValueError):
126 SkyProjection(Transform(ICRS, ICRS, astshim.UnitMap(2)))
128 with self.assertRaises(ValueError):
129 SkyProjection(Transform(pixel_frame, pixel_frame, astshim.UnitMap(2)))
131 def make_pixel_to_sky(ast_mapping: astshim.Mapping | None = None) -> Transform[Any, Any]:
132 mapping = ast_mapping if ast_mapping is not None else astshim.UnitMap(2)
133 return Transform(pixel_frame, ICRS, mapping)
135 base = SkyProjection(make_pixel_to_sky())
137 # Identity short-circuit: an object is always equal to itself.
138 self.assertEqual(base, base)
140 # Two independently constructed but equivalent projections are equal.
141 self.assertEqual(base, SkyProjection(make_pixel_to_sky()))
143 # Comparison against a non-SkyProjection yields NotImplemented, so
144 # Python falls back to identity.
145 self.assertFalse(base == "not a projection")
146 self.assertTrue(base != "not a projection")
147 self.assertNotEqual(base, None)
149 # Differ only in the pixel-to-sky transform.
150 self.assertNotEqual(base, SkyProjection(make_pixel_to_sky(astshim.ShiftMap([1.0, 2.0]))))
152 # The fits_approximation branch: absent on ``base`` but present here.
153 with_approx = SkyProjection(
154 make_pixel_to_sky(), fits_approximation=make_pixel_to_sky(astshim.ShiftMap([0.1, 0.2]))
155 )
156 self.assertNotEqual(base, with_approx)
158 # Equal pixel-to-sky and equal fits_approximations are equal.
159 with_approx_again = SkyProjection(
160 make_pixel_to_sky(), fits_approximation=make_pixel_to_sky(astshim.ShiftMap([0.1, 0.2]))
161 )
162 self.assertEqual(with_approx, with_approx_again)
164 # Same pixel-to-sky transform but a different fits_approximation.
165 other_approx = SkyProjection(
166 make_pixel_to_sky(), fits_approximation=make_pixel_to_sky(astshim.ShiftMap([0.3, 0.4]))
167 )
168 self.assertNotEqual(with_approx, other_approx)
170 def test_affine_2x2(self) -> None:
171 """Test an affine transform constructed from a 2x2 matrix."""
172 in_frame = DetectorFrame(**DP2_VISIT_DETECTOR_DATA_ID, bbox=Box.factory[:5, :4])
173 out_frame = GeneralFrame(unit=u.pix)
174 transform_matrix = np.array([[2.0, 0.25], [-0.75, 0.8]])
175 in_xy = in_frame.bbox.meshgrid().map(np.ravel)
176 in_matrix = np.array([in_xy.x, in_xy.y])
177 out_matrix = np.dot(transform_matrix, in_matrix)
178 check_transform(
179 self,
180 Transform.affine(in_frame, out_frame, transform_matrix),
181 in_xy,
182 XY(x=out_matrix[0, :], y=out_matrix[1, :]),
183 in_frame,
184 out_frame,
185 in_atol=1e-15 * u.pix,
186 out_atol=1e-15 * u.pix,
187 )
189 def test_affine_3x3(self) -> None:
190 """Test an affine transform constructed from a 3x3 matrix."""
191 in_frame = DetectorFrame(**DP2_VISIT_DETECTOR_DATA_ID, bbox=Box.factory[:5, :4])
192 out_frame = GeneralFrame(unit=u.pix)
193 transform_matrix = np.array([[2.0, 0.25, -0.5], [-0.75, 0.8, 0.4], [0.0, 0.0, 1.0]])
194 in_xy = in_frame.bbox.meshgrid().map(np.ravel)
195 in_matrix = np.array([in_xy.x, in_xy.y, np.ones(in_xy.x.shape)])
196 out_matrix = np.dot(transform_matrix, in_matrix)
197 check_transform(
198 self,
199 Transform.affine(in_frame, out_frame, transform_matrix),
200 in_xy,
201 XY(x=out_matrix[0, :], y=out_matrix[1, :]),
202 in_frame,
203 out_frame,
204 in_atol=1e-15 * u.pix,
205 out_atol=1e-15 * u.pix,
206 )
208 @unittest.skipUnless(DATA_DIR is not None, "TESTDATA_IMAGES_DIR is not in the environment.")
209 def test_camera(self) -> None:
210 """Test that we can:
212 - make a CameraFrameSet from the AST representation returned by afw;
213 - transform points and get the same result as afw;
214 - round-trip the CameraFrameSet through FITS serialization and still
215 do all of that;
216 - also roundtrip a Transform that can be obtained from the
217 CameraFrameSet, by referencing the mappings in the frame set.
219 This test is skipped if legacy modules cannot be imported.
221 This test provides coverage for the archive system's pointer and
222 frame-set reference machinery.
223 """
224 try:
225 from lsst.afw.cameraGeom import Camera
226 except ImportError:
227 raise unittest.SkipTest("'lsst.afw.cameraGeom' could not be imported.") from None
228 assert DATA_DIR is not None, "Guaranteed by decorator."
229 filename = os.path.join(DATA_DIR, "dp2", "legacy", "camera.fits")
230 legacy_camera = Camera.readFits(filename)
231 frame_set = CameraFrameSet.from_legacy(legacy_camera)
232 detector_id: int = DP2_VISIT_DETECTOR_DATA_ID["detector"]
233 self.compare_to_legacy_camera(legacy_camera, frame_set)
234 test_holder = FrameSetTestHolder(
235 frames=frame_set,
236 pixels_to_fp=frame_set[frame_set.detector(detector_id), frame_set.focal_plane()],
237 )
238 with RoundtripFits(self, test_holder) as roundtrip1:
239 self.assertEqual(len(roundtrip1.serialized.pixels_to_fp.frames), 2)
240 self.assertEqual(len(roundtrip1.serialized.pixels_to_fp.bounds), 2)
241 self.assertEqual(len(roundtrip1.serialized.pixels_to_fp.mappings), 1)
242 # Instead of storing the AST mapping directly, we should have
243 # stored a reference to the frame set:
244 self.assertIsInstance(roundtrip1.serialized.pixels_to_fp.mappings[0], PointerModel)
245 self.compare_to_legacy_camera(legacy_camera, roundtrip1.result.frames)
246 self.assertEqual(roundtrip1.result.pixels_to_fp.in_frame, frame_set.detector(detector_id))
247 self.assertEqual(roundtrip1.result.pixels_to_fp.out_frame, frame_set.focal_plane())
248 self.assertEqual(
249 roundtrip1.result.pixels_to_fp._ast_mapping.simplified().show(),
250 test_holder.pixels_to_fp._ast_mapping.simplified().show(),
251 )
252 with RoundtripJson(self, test_holder) as roundtrip2:
253 self.assertEqual(len(roundtrip2.serialized.pixels_to_fp.frames), 2)
254 self.assertEqual(len(roundtrip2.serialized.pixels_to_fp.bounds), 2)
255 self.assertEqual(len(roundtrip2.serialized.pixels_to_fp.mappings), 1)
256 # Instead of storing the AST mapping directly, we should have
257 # stored a reference to the frame set:
258 self.assertIsInstance(roundtrip2.serialized.pixels_to_fp.mappings[0], JsonRef)
259 raw_data = roundtrip2.inspect()
260 self.assertEqual(len(raw_data["indirect"]), 1)
261 self.assertEqual(raw_data["frames"], {"$ref": "#/indirect/0"})
262 self.compare_to_legacy_camera(legacy_camera, roundtrip2.result.frames)
263 self.assertEqual(roundtrip2.result.pixels_to_fp.in_frame, frame_set.detector(detector_id))
264 self.assertEqual(roundtrip2.result.pixels_to_fp.out_frame, frame_set.focal_plane())
265 self.assertEqual(
266 roundtrip2.result.pixels_to_fp._ast_mapping.simplified().show(),
267 test_holder.pixels_to_fp._ast_mapping.simplified().show(),
268 )
270 def compare_to_legacy_camera(self, legacy_camera: Any, frame_set: CameraFrameSet) -> None:
271 """Test the transforms extracted from a CameraFrameSet against the
272 legacy lsst.afw.cameraGeom implementations.
273 """
274 from lsst.afw.cameraGeom import FIELD_ANGLE, FOCAL_PLANE, PIXELS
275 from lsst.geom import Point2D
277 legacy_detector = legacy_camera[16]
278 pixel_legacy_points = [Point2D(50.0, 60.0), Point2D(801.2, 322.8), Point2D(33.5, 22.1)]
279 fp_legacy_points = [legacy_detector.transform(p, PIXELS, FOCAL_PLANE) for p in pixel_legacy_points]
280 fa_legacy_points = [legacy_detector.transform(p, PIXELS, FIELD_ANGLE) for p in pixel_legacy_points]
281 pixel_xy_array = legacy_points_to_xy_array(pixel_legacy_points)
282 fp_xy_array = legacy_points_to_xy_array(fp_legacy_points)
283 fa_xy_array = legacy_points_to_xy_array(fa_legacy_points)
284 # Test transforms extracted directly from the frame set.
285 pixel_to_fp = frame_set[frame_set.detector(16), frame_set.focal_plane()]
286 check_transform(
287 self, pixel_to_fp, pixel_xy_array, fp_xy_array, frame_set.detector(16), frame_set.focal_plane()
288 )
289 pixel_to_fa = frame_set[frame_set.detector(16), frame_set.field_angle()]
290 check_transform(
291 self, pixel_to_fa, pixel_xy_array, fa_xy_array, frame_set.detector(16), frame_set.field_angle()
292 )
293 fp_to_fa = frame_set[frame_set.focal_plane(), frame_set.field_angle()]
294 check_transform(
295 self, fp_to_fa, fp_xy_array, fa_xy_array, frame_set.focal_plane(), frame_set.field_angle()
296 )
297 # Test a composition.
298 pixel_to_fa_indirect = pixel_to_fp.then(fp_to_fa)
299 check_transform(
300 self,
301 pixel_to_fa_indirect,
302 pixel_xy_array,
303 fa_xy_array,
304 frame_set.detector(16),
305 frame_set.field_angle(),
306 )
307 pixel_to_fp_d, fp_to_fa_d = pixel_to_fa_indirect.decompose()
308 check_transform(
309 self, pixel_to_fp_d, pixel_xy_array, fp_xy_array, frame_set.detector(16), frame_set.focal_plane()
310 )
311 check_transform(
312 self, fp_to_fa_d, fp_xy_array, fa_xy_array, frame_set.focal_plane(), frame_set.field_angle()
313 )
314 fa_to_fp_d, fp_to_pixel_d = pixel_to_fa_indirect.inverted().decompose()
315 check_transform(
316 self, fa_to_fp_d, fa_xy_array, fp_xy_array, frame_set.field_angle(), frame_set.focal_plane()
317 )
318 check_transform(
319 self, fp_to_pixel_d, fp_xy_array, pixel_xy_array, frame_set.focal_plane(), frame_set.detector(16)
320 )
322 @unittest.skipUnless(DATA_DIR is not None, "TESTDATA_IMAGES_DIR is not in the environment.")
323 def test_detector_wcs(self) -> None:
324 """Test the Transform/SkyProjection representation of a detector
325 WCS.
326 """
327 try:
328 from lsst.afw.image import ExposureFitsReader
329 except ImportError:
330 raise unittest.SkipTest("'lsst.afw.image' could not be imported.") from None
331 assert DATA_DIR is not None, "Guaranteed by decorator."
332 filename = os.path.join(DATA_DIR, "dp2", "legacy", "visit_image.fits")
333 reader = ExposureFitsReader(filename)
334 legacy_wcs = reader.readWcs()
335 wcs_bbox = Box.from_legacy(reader.readDetector().getBBox())
336 subimage_bbox = Box.from_legacy(reader.readBBox())
337 detector_frame = DetectorFrame(**DP2_VISIT_DETECTOR_DATA_ID, bbox=wcs_bbox)
338 sky_projection = SkyProjection.from_legacy(legacy_wcs, detector_frame)
339 assert sky_projection.fits_approximation is not None
340 compare_sky_projection_to_legacy_wcs(self, sky_projection, legacy_wcs, detector_frame, subimage_bbox)
341 # When we convert from a legacy SkyWcs, the internal AST Mapping needs
342 # to really be an AST FrameSet in order to be able to convert back.
343 self.assertIn("Begin FrameSet", sky_projection.show())
344 compare_sky_projection_to_legacy_wcs(
345 self, sky_projection, sky_projection.to_legacy(), detector_frame, subimage_bbox
346 )
347 self.assertIn("Begin FrameSet", sky_projection.fits_approximation.show())
348 compare_sky_projection_to_legacy_wcs(
349 self,
350 sky_projection.fits_approximation,
351 sky_projection.fits_approximation.to_legacy(),
352 detector_frame,
353 subimage_bbox,
354 is_fits=True,
355 )
356 with RoundtripJson(self, sky_projection, "SkyProjection") as roundtrip:
357 pass
358 compare_sky_projection_to_legacy_wcs(
359 self, roundtrip.result, legacy_wcs, detector_frame, subimage_bbox
360 )
361 # The AST FrameSet-ness needs to propagate through serialization.
362 self.assertIn("Begin FrameSet", roundtrip.result.show())
363 compare_sky_projection_to_legacy_wcs(
364 self, sky_projection, roundtrip.result.to_legacy(), detector_frame, subimage_bbox
365 )
366 with RoundtripJson(self, sky_projection.fits_approximation, "SkyProjection") as roundtrip:
367 pass
368 compare_sky_projection_to_legacy_wcs(
369 self,
370 roundtrip.result,
371 legacy_wcs.getFitsApproximation(),
372 detector_frame,
373 subimage_bbox,
374 is_fits=True,
375 )
376 self.assertIn("Begin FrameSet", roundtrip.result.show())
377 compare_sky_projection_to_legacy_wcs(
378 self,
379 sky_projection.fits_approximation,
380 roundtrip.result.to_legacy(),
381 detector_frame,
382 subimage_bbox,
383 is_fits=True,
384 )
387@dataclasses.dataclass
388class FrameSetTestHolder:
389 """A top-level object that holds a CameraFrameSet and a transform
390 extracted from it, for testing archive pointers and frame set references.
391 """
393 frames: CameraFrameSet
394 pixels_to_fp: Transform[DetectorFrame, FocalPlaneFrame]
396 def serialize[P: pydantic.BaseModel](self, archive: OutputArchive[P]) -> FrameSetTestHolderModel[P]:
397 frames_model = archive.serialize_frame_set(
398 "frames", self.frames, self.frames.serialize, key=id(self.frames)
399 )
400 pixels_to_fp_model = archive.serialize_direct(
401 "pixels_to_fp", functools.partial(self.pixels_to_fp.serialize, use_frame_sets=True)
402 )
403 return FrameSetTestHolderModel[P](frames=frames_model, pixels_to_fp=pixels_to_fp_model)
405 @staticmethod
406 def _get_archive_tree_type[P: pydantic.BaseModel](
407 pointer_type: type[P],
408 ) -> type[FrameSetTestHolderModel[P]]:
409 return FrameSetTestHolderModel[pointer_type] # type: ignore
412class FrameSetTestHolderModel[P: pydantic.BaseModel](ArchiveTree):
413 """The serialization model for FrameSetTestHolder."""
415 SCHEMA_NAME: ClassVar[str] = "_test_frame_set_holder"
416 SCHEMA_VERSION: ClassVar[str] = "1.0.0"
417 MIN_READ_VERSION: ClassVar[int] = 1
418 PUBLIC_TYPE: ClassVar[type] = FrameSetTestHolder
420 frames: CameraFrameSetSerializationModel | P
421 pixels_to_fp: TransformSerializationModel[P]
423 def deserialize(self, archive: InputArchive[Any]) -> FrameSetTestHolder:
424 assert not isinstance(self.frames, CameraFrameSetSerializationModel), "Archive pointer expected."
425 frames = archive.deserialize_pointer(
426 self.frames, CameraFrameSetSerializationModel, CameraFrameSetSerializationModel.deserialize
427 )
428 pixels_to_fp = self.pixels_to_fp.deserialize(archive)
429 return FrameSetTestHolder(frames, pixels_to_fp)
432if __name__ == "__main__":
433 unittest.main()