Coverage for tests/test_stencils.py: 96%
220 statements
« prev ^ index » next coverage.py v7.14.3, created at 2026-06-27 08:47 +0000
1# This file is part of dax_images_cutout.
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/>.
22from __future__ import annotations
24import unittest
26import astropy.coordinates
27import astropy.io.fits
28import astropy.units as u
29import astropy.wcs
30import numpy as np
32import lsst.sphgeom
33from lsst.dax.images.cutout.stencils import (
34 MaskBackend,
35 SkyCircle,
36 SkyPolygon,
37 SkyStencil,
38 StencilNotContainedError,
39)
40from lsst.images import Box, GeneralFrame, Mask, MaskPlane, MaskSchema, SkyProjection
41from lsst.sphgeom import Angle, LonLat, UnitVector3d # noqa: F401 (used by eval(repr))
43# Bounding box for the cutout tests, in [y, x] (stop exclusive). Slightly
44# bigger in x to catch x<->y transposition bugs.
45TEST_BOX = Box.factory[-13:28, -16:27]
48def _arcsec(value: float) -> Angle:
49 """Return a `lsst.sphgeom.Angle` for ``value`` arcseconds."""
50 return Angle((value * u.arcsec).to_value(u.rad))
53def _make_wcs() -> astropy.wcs.WCS:
54 """Build a gnomonic FITS WCS with 0.1 arcsec pixels at (12, 13) deg.
56 The reference pixel is placed at pixel (5, 7) so the stencils land at an
57 arbitrary nonzero offset within `TEST_BOX`.
58 """
59 wcs = astropy.wcs.WCS(naxis=2)
60 # FITS CRPIX is 1-based, so 0-based pixel (5, 7) is CRPIX (6, 8).
61 wcs.wcs.crpix = [6.0, 8.0]
62 wcs.wcs.crval = [12.0, 13.0]
63 scale = 0.1 / 3600.0
64 wcs.wcs.cd = [[-scale, 0.0], [0.0, scale]]
65 wcs.wcs.ctype = ["RA---TAN", "DEC--TAN"]
66 return wcs
69def _make_car_wcs() -> astropy.wcs.WCS:
70 """Build a plate-carree (CAR) WCS referenced on the equator.
72 CAR is non-gnomonic, so great circles do not map to straight lines in
73 pixel space. The equatorial reference keeps pixel coordinates a simple
74 (lon, lat) grid; a polygon placed at high declination then has edges that
75 bow well away from the straight chords joining its projected vertices.
76 """
77 wcs = astropy.wcs.WCS(naxis=2)
78 wcs.wcs.crpix = [1.0, 1.0]
79 wcs.wcs.crval = [0.0, 0.0]
80 scale = 0.02
81 wcs.wcs.cd = [[-scale, 0.0], [0.0, scale]]
82 wcs.wcs.ctype = ["RA---CAR", "DEC--CAR"]
83 return wcs
86class ModuleHelpersTestCase(unittest.TestCase):
87 """Tests for module-level helpers that survive the rewrite."""
89 def test_mask_backend_members(self) -> None:
90 self.assertEqual({b.name for b in MaskBackend}, {"AST", "SPHGEOM"})
93class SkyCircleTestCase(unittest.TestCase):
94 """Tests for `SkyCircle`."""
96 def setUp(self) -> None:
97 self.center = LonLat.fromDegrees(12.0, 13.0)
98 self.instance = SkyCircle(self.center, _arcsec(1.0))
100 def test_from_astropy(self) -> None:
101 other = SkyCircle.from_astropy(
102 astropy.coordinates.SkyCoord(
103 frame="icrs", ra=12.0 * astropy.units.deg, dec=13.0 * astropy.units.deg
104 ),
105 astropy.coordinates.Angle(1.0 * astropy.units.arcsec),
106 )
107 self.assertEqual(self.instance.region, other.region)
109 def test_repr(self) -> None:
110 self.assertEqual(eval(repr(self.instance)).region, self.instance.region)
112 def test_to_pixel(self) -> None:
113 _check_to_pixel(self, self.instance, _make_wcs(), backend=MaskBackend.AST, max_missing=2, max_extra=2)
115 def test_to_polygon(self) -> None:
116 polygon_stencil = self.instance.to_polygon()
117 self.assertNotEqual(
118 self.instance.region.relate(polygon_stencil.region.getBoundingCircle()), lsst.sphgeom.DISJOINT
119 )
120 _check_to_pixel(
121 self, polygon_stencil, _make_wcs(), backend=MaskBackend.AST, max_missing=6, max_extra=6
122 )
124 def test_ast_sky_region_circle_contains_center(self) -> None:
125 region = self.instance._ast_sky_region()
126 self.assertTrue(
127 region.pointinregion([self.center.getLon().asRadians(), self.center.getLat().asRadians()])
128 )
130 def test_to_pixel_sphgeom(self) -> None:
131 _check_to_pixel(
132 self, self.instance, _make_wcs(), backend=MaskBackend.SPHGEOM, max_missing=0, max_extra=0
133 )
135 def test_to_pixel_sphgeom_polygon(self) -> None:
136 polygon_stencil = self.instance.to_polygon()
137 _check_to_pixel(
138 self, polygon_stencil, _make_wcs(), backend=MaskBackend.SPHGEOM, max_missing=0, max_extra=0
139 )
142class SkyPolygonTestCase(unittest.TestCase):
143 """Tests for `SkyPolygon` orientation handling."""
145 def setUp(self) -> None:
146 self.instance = SkyCircle(LonLat.fromDegrees(12.0, 13.0), _arcsec(2.0)).to_polygon(n_vertices=8)
148 def test_ast_sky_region_polygon_contains_centroid(self) -> None:
149 region = self.instance._ast_sky_region()
150 lonlat = lsst.sphgeom.LonLat(self.instance.region.getCentroid())
151 self.assertTrue(region.pointinregion([lonlat.getLon().asRadians(), lonlat.getLat().asRadians()]))
154class BackendComparisonTestCase(unittest.TestCase):
155 """Assert the AST and sphgeom backends agree on bbox and masked pixels."""
157 def setUp(self) -> None:
158 self.center = LonLat.fromDegrees(12.0, 13.0)
159 self.projection = SkyProjection.from_fits_wcs(_make_wcs(), GeneralFrame(unit=u.pix))
160 self.box = TEST_BOX
162 def _masked_array(self, stencil: SkyStencil, backend: MaskBackend) -> tuple[np.ndarray, Box]:
163 pixel_stencil = stencil.to_pixels(self.projection, self.box, backend=backend)
164 mask = Mask(schema=MaskSchema([MaskPlane("STENCIL", "stencil coverage")]), bbox=self.box)
165 pixel_stencil.set_mask(mask, "STENCIL")
166 return mask.get("STENCIL"), pixel_stencil.bbox
168 def test_backends_agree_circle(self) -> None:
169 circle = SkyCircle(self.center, _arcsec(1.0))
170 ast_mask, ast_box = self._masked_array(circle, MaskBackend.AST)
171 sph_mask, sph_box = self._masked_array(circle, MaskBackend.SPHGEOM)
172 self.assertEqual(ast_box, sph_box)
173 self.assertEqual(int(np.sum(ast_mask != sph_mask)), 0)
175 def test_backends_agree_polygon(self) -> None:
176 polygon = SkyCircle(self.center, _arcsec(1.0)).to_polygon()
177 ast_mask, ast_box = self._masked_array(polygon, MaskBackend.AST)
178 sph_mask, sph_box = self._masked_array(polygon, MaskBackend.SPHGEOM)
179 self.assertEqual(ast_box, sph_box)
180 self.assertLessEqual(int(np.sum(ast_mask != sph_mask)), 12)
182 def test_set_mask_covered_false_marks_outside(self) -> None:
183 """``set_mask(covered=False)`` flags exactly the pixels the stencil
184 does not cover, including the region of the mask outside the stencil's
185 bounding box.
186 """
187 circle = SkyCircle(self.center, _arcsec(1.0))
188 pixel_stencil = circle.to_pixels(self.projection, self.box)
190 inside = Mask(schema=MaskSchema([MaskPlane("STENCIL", "stencil coverage")]), bbox=self.box)
191 pixel_stencil.set_mask(inside, "STENCIL")
193 outside = Mask(schema=MaskSchema([MaskPlane("STENCIL", "stencil coverage")]), bbox=self.box)
194 pixel_stencil.set_mask(outside, "STENCIL", covered=False)
196 inside_arr = inside.get("STENCIL")
197 outside_arr = outside.get("STENCIL")
198 # The two planes partition the mask: every pixel is flagged in exactly
199 # one of them.
200 np.testing.assert_array_equal(outside_arr, np.logical_not(inside_arr))
201 # The stencil covers some pixels but not the whole box, so neither
202 # plane is empty.
203 self.assertTrue(inside_arr.any())
204 self.assertTrue(outside_arr.any())
207class GreatCircleCurvatureTestCase(unittest.TestCase):
208 """Polygon stencils whose great-circle edges curve in pixel space.
210 The other tests use a gnomonic (TAN) projection, which maps great circles
211 to exactly straight lines and so cannot exercise edge curvature. These
212 tests use a plate-carree (CAR) projection referenced on the equator with a
213 polygon at high declination, where the great-circle edges bow well away
214 from the straight pixel-space chords joining the projected vertices. Both
215 mask backends must follow the true great circle rather than the chord.
217 The vertices land exactly on pixel centers in this geometry (lon ``+/-4``
218 and ``0`` degrees, dec ``70`` and ``66`` degrees map to integer pixels at
219 this reference and scale), so the handful of pixels of residual
220 disagreement allowed by the tolerances below are the vertex pixels
221 themselves: their centers sit exactly on the polygon boundary, where
222 containment is a tie that each backend's edge test resolves differently.
223 Vertices at generic sub-pixel positions would typically agree exactly.
224 """
226 def setUp(self) -> None:
227 self.wcs = _make_car_wcs()
228 self.projection = SkyProjection.from_fits_wcs(self.wcs, GeneralFrame(unit=u.pix))
229 self.polygon = SkyPolygon(
230 [
231 LonLat.fromDegrees(-4.0, 70.0),
232 LonLat.fromDegrees(4.0, 70.0),
233 LonLat.fromDegrees(0.0, 66.0),
234 ]
235 )
236 # The tight pixel bounding box is backend-independent, so any backend
237 # may be used to obtain it from a generous reference box.
238 self.box = self.polygon.to_pixels(self.projection, Box.factory[-10000:10000, -10000:10000]).bbox
240 def _coverage(self, backend: MaskBackend) -> np.ndarray:
241 pixel_stencil = self.polygon.to_pixels(self.projection, self.box, backend=backend)
242 mask = Mask(schema=MaskSchema([MaskPlane("STENCIL", "stencil coverage")]), bbox=self.box)
243 pixel_stencil.set_mask(mask, "STENCIL")
244 return mask.get("STENCIL")
246 def test_scenario_exercises_curvature(self) -> None:
247 """The true spherical coverage differs substantially from a straight-
248 edged pixel-space approximation, so the backend checks below are a
249 meaningful test of great-circle handling rather than vacuously true.
250 """
251 truth = _brute_force_stencil_array(self.polygon, self.wcs, self.box)
252 cartesian = _cartesian_pixel_coverage(self.polygon, self.wcs, self.box)
253 self.assertGreater(int(np.sum(truth != cartesian)), 500)
255 def test_ast_backend_follows_great_circle(self) -> None:
256 # Only the three vertex pixels may disagree (see class docstring); a
257 # larger count would mean the edges were rasterized as straight pixel
258 # chords, as happens under AST's default ``SimpVertices=1``.
259 truth = _brute_force_stencil_array(self.polygon, self.wcs, self.box)
260 got = self._coverage(MaskBackend.AST)
261 self.assertLessEqual(int(np.sum(got != truth)), 3)
263 def test_sphgeom_backend_follows_great_circle(self) -> None:
264 truth = _brute_force_stencil_array(self.polygon, self.wcs, self.box)
265 got = self._coverage(MaskBackend.SPHGEOM)
266 self.assertLessEqual(int(np.sum(got != truth)), 3)
269def _brute_force_stencil_array(sky_stencil: SkyStencil, wcs: astropy.wcs.WCS, box: Box) -> np.ndarray:
270 """Make a boolean ``(ny, nx)`` array, `True` where a center is inside.
272 The pixel grid is transformed to the sky with the FITS WCS (independent of
273 the `SkyProjection` under test) and tested against the stencil's sphgeom
274 region.
275 """
276 grid = box.meshgrid()
277 sky = wcs.pixel_to_world(grid.x.ravel(), grid.y.ravel())
278 contained = sky_stencil.region.contains(sky.ra.rad, sky.dec.rad)
279 return contained.reshape(box.shape)
282def _cartesian_pixel_coverage(polygon: SkyPolygon, wcs: astropy.wcs.WCS, box: Box) -> np.ndarray:
283 """Make a boolean ``(ny, nx)`` array for the polygon with straight edges.
285 The vertices are projected to pixels and joined by straight chords (a
286 convex point-in-polygon test). This models the cartesian rasterization
287 that ignores great-circle curvature, so it can be compared against the
288 true spherical coverage to show the curved scenario is non-trivial.
289 """
290 vertices = polygon._boundary_skycoord()
291 vx, vy = wcs.world_to_pixel_values(vertices.ra.deg, vertices.dec.deg)
292 grid = box.meshgrid()
293 px = grid.x.ravel().astype(float)
294 py = grid.y.ravel().astype(float)
295 n = len(vx)
296 cross = np.array(
297 [
298 (vx[(i + 1) % n] - vx[i]) * (py - vy[i]) - (vy[(i + 1) % n] - vy[i]) * (px - vx[i])
299 for i in range(n)
300 ]
301 )
302 inside = np.all(cross >= 0.0, axis=0) | np.all(cross <= 0.0, axis=0)
303 return inside.reshape(box.shape)
306def _check_to_pixel(
307 test_case: unittest.TestCase,
308 sky_stencil: SkyStencil,
309 wcs: astropy.wcs.WCS,
310 *,
311 box: Box = TEST_BOX,
312 expected_bbox: Box | None = None,
313 backend: MaskBackend = MaskBackend.AST,
314 max_missing: int = 0,
315 max_extra: int = 0,
316 plot: bool = False,
317) -> None:
318 """Check a `SkyStencil.to_pixels` result against brute force.
320 ``box`` is the reference bounding box passed to `to_pixels`; when it does
321 not fully contain the stencil the result is clipped to it. Brute force is
322 evaluated over ``box`` too, which yields the correct expected coverage for
323 a clipped stencil: a pixel that is inside the region and inside ``box`` is
324 necessarily inside the clipped bounding box, since the region is contained
325 by its own tight bounding box. ``expected_bbox``, if given, is asserted to
326 equal the clipped result bounding box.
327 """
328 projection = SkyProjection.from_fits_wcs(wcs, GeneralFrame(unit=u.pix))
329 pixel_stencil = sky_stencil.to_pixels(projection, box, backend=backend)
330 test_case.assertTrue(box.contains(pixel_stencil.bbox))
331 if expected_bbox is not None:
332 test_case.assertEqual(pixel_stencil.bbox, expected_bbox)
333 mask = Mask(schema=MaskSchema([MaskPlane("STENCIL", "stencil coverage")]), bbox=box)
334 pixel_stencil.set_mask(mask, "STENCIL")
335 got = mask.get("STENCIL")
336 check_array = _brute_force_stencil_array(sky_stencil, wcs, box)
337 missing = np.logical_and(check_array, np.logical_not(got))
338 extra = np.logical_and(got, np.logical_not(check_array))
339 if plot: 339 ↛ 340line 339 didn't jump to line 340 because the condition on line 339 was never true
340 from matplotlib import pyplot
342 display_array = np.zeros((box.shape.y, box.shape.x, 3), dtype=np.uint8)
343 display_array[:, :, 0] = 255 * check_array
344 display_array[:, :, 1] = 255 * got
345 pyplot.imshow(display_array, origin="lower", interpolation="nearest")
346 pyplot.title("red=check, green=SkyStencil.to_pixel, yellow=both")
347 pyplot.show()
348 test_case.assertLessEqual(int(missing.sum()), max_missing)
349 test_case.assertLessEqual(int(extra.sum()), max_extra)
352class StencilContainmentTestCase(unittest.TestCase):
353 """Clipping and raising when a stencil only partially overlaps, or does not
354 overlap at all, the reference bounding box passed to `to_pixels`.
356 The 1 arcsec circle used throughout has the fixed tight pixel bounding box
357 ``Box.factory[-3:18, -5:16]`` under `_make_wcs`, so the reference boxes
358 below produce exactly predictable intersections.
359 """
361 # Reference boxes relative to the circle's tight pixel bbox
362 # [y=-3:18, x=-5:16].
363 PARTIAL_BOX = Box.factory[5:30, 5:30]
364 PARTIAL_CLIPPED = Box.factory[5:18, 5:16]
365 INSIDE_STENCIL_BOX = Box.factory[12:18, 12:16]
366 TOUCHING_BOX = Box.factory[-3:18, 16:30]
367 DISJOINT_BOX = Box.factory[100:120, 100:120]
369 # Per-backend rasterization tolerance, matching the existing circle tests.
370 BACKEND_TOLERANCE = {MaskBackend.AST: 2, MaskBackend.SPHGEOM: 0}
372 def setUp(self) -> None:
373 self.center = LonLat.fromDegrees(12.0, 13.0)
374 self.wcs = _make_wcs()
375 self.projection = SkyProjection.from_fits_wcs(self.wcs, GeneralFrame(unit=u.pix))
377 def _circle(self, *, clip: bool) -> SkyCircle:
378 return SkyCircle(self.center, _arcsec(1.0), clip=clip)
380 # Box resolution happens in `to_pixels` before any mask backend is
381 # selected, so the raising behavior is backend-independent; the default
382 # backend is sufficient for the raising tests below.
384 def test_clip_false_raises_on_partial_overlap(self) -> None:
385 with self.assertRaises(StencilNotContainedError):
386 self._circle(clip=False).to_pixels(self.projection, self.PARTIAL_BOX)
388 def test_clip_false_raises_when_box_inside_stencil(self) -> None:
389 with self.assertRaises(StencilNotContainedError):
390 self._circle(clip=False).to_pixels(self.projection, self.INSIDE_STENCIL_BOX)
392 def test_clip_false_raises_when_disjoint(self) -> None:
393 with self.assertRaises(StencilNotContainedError):
394 self._circle(clip=False).to_pixels(self.projection, self.DISJOINT_BOX)
396 def test_clip_true_raises_when_touching(self) -> None:
397 # The box starts one pixel beyond the tight bbox's max x, so the two
398 # share no pixel and clipping cannot produce an overlap.
399 with self.assertRaises(StencilNotContainedError):
400 self._circle(clip=True).to_pixels(self.projection, self.TOUCHING_BOX)
402 def test_clip_true_raises_when_disjoint(self) -> None:
403 with self.assertRaises(StencilNotContainedError):
404 self._circle(clip=True).to_pixels(self.projection, self.DISJOINT_BOX)
406 def test_clip_true_unchanged_when_contained(self) -> None:
407 # A fully contained stencil keeps its tight bbox even when clipping.
408 for backend, tolerance in self.BACKEND_TOLERANCE.items():
409 with self.subTest(backend=str(backend)):
410 _check_to_pixel(
411 self,
412 self._circle(clip=True),
413 self.wcs,
414 box=TEST_BOX,
415 expected_bbox=Box.factory[-3:18, -5:16],
416 backend=backend,
417 max_missing=tolerance,
418 max_extra=tolerance,
419 )
421 def test_clip_true_clips_to_intersection_on_partial_overlap(self) -> None:
422 for backend, tolerance in self.BACKEND_TOLERANCE.items():
423 with self.subTest(backend=str(backend)):
424 _check_to_pixel(
425 self,
426 self._circle(clip=True),
427 self.wcs,
428 box=self.PARTIAL_BOX,
429 expected_bbox=self.PARTIAL_CLIPPED,
430 backend=backend,
431 max_missing=tolerance,
432 max_extra=tolerance,
433 )
435 def test_clip_true_clips_to_box_when_box_inside_stencil(self) -> None:
436 for backend, tolerance in self.BACKEND_TOLERANCE.items():
437 with self.subTest(backend=str(backend)):
438 _check_to_pixel(
439 self,
440 self._circle(clip=True),
441 self.wcs,
442 box=self.INSIDE_STENCIL_BOX,
443 expected_bbox=self.INSIDE_STENCIL_BOX,
444 backend=backend,
445 max_missing=tolerance,
446 max_extra=tolerance,
447 )
450class StencilFitsMetadataTestCase(unittest.TestCase):
451 """`SkyStencil.to_fits_metadata` returns an `astropy.io.fits.Header` whose
452 cards carry the descriptive comments.
453 """
455 def test_circle(self) -> None:
456 circle = SkyCircle(LonLat.fromDegrees(12.0, 13.0), _arcsec(1.0))
457 header = circle.to_fits_metadata()
458 self.assertIsInstance(header, astropy.io.fits.Header)
459 self.assertEqual(header["ST_TYPE"], "CIRCLE")
460 self.assertEqual(header.comments["ST_TYPE"], "Type of stencil used to create this cutout")
461 self.assertAlmostEqual(header["ST_RA"], 12.0)
462 self.assertAlmostEqual(header["ST_DEC"], 13.0)
463 self.assertAlmostEqual(header["ST_RAD"], (1.0 * u.arcsec).to_value(u.deg))
464 self.assertEqual(header.comments["ST_RAD"], "[deg] Circle radius")
466 def test_polygon(self) -> None:
467 polygon = SkyCircle(LonLat.fromDegrees(12.0, 13.0), _arcsec(2.0)).to_polygon(n_vertices=4)
468 header = polygon.to_fits_metadata()
469 self.assertIsInstance(header, astropy.io.fits.Header)
470 self.assertEqual(header["ST_TYPE"], "POLYGON")
471 self.assertEqual(header.comments["ST_TYPE"], "Type of stencil used to create this cutout")
472 self.assertIn("ST_RA00", header)
473 self.assertIn("ST_DEC00", header)
474 self.assertEqual(header.comments["ST_RA00"], "[deg] Vertex 0 Right Ascension")
475 self.assertEqual(header.comments["ST_DEC00"], "[deg] Vertex 0 Declination")
478if __name__ == "__main__": 478 ↛ 479line 478 didn't jump to line 479 because the condition on line 478 was never true
479 unittest.main()