Coverage for tests / test_stencils.py: 31%
61 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-30 09:26 +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.units
28import numpy as np
30import lsst.sphgeom
31from lsst.afw.geom import SkyWcs, makeCdMatrix, makeSkyWcs
32from lsst.afw.image import Mask
33from lsst.dax.images.cutout.stencils import SkyCircle, SkyStencil
34from lsst.geom import Angle, Box2I, Point2D, Point2I, SpherePoint, arcseconds, degrees
37class SkyCircleTestCase(unittest.TestCase):
38 """Tests for `SkyCircle."""
40 def setUp(self) -> None:
41 self.center = SpherePoint(12.0, 13.0, degrees)
42 self.instance = SkyCircle(self.center, Angle(1.0, arcseconds))
44 def test_from_astropy(self) -> None:
45 """Test that the from_astropy factory function is equivalent to
46 passing the same values to the constructor, by comparing `region`
47 result.
49 This is a pretty circular test, but it's at least runs a lot of code
50 to make sure it's not completely broken.
51 """
52 other = SkyCircle.from_astropy(
53 astropy.coordinates.SkyCoord(
54 frame="icrs", ra=12.0 * astropy.units.deg, dec=13.0 * astropy.units.deg
55 ),
56 astropy.coordinates.Angle(1.0 * astropy.units.arcsec),
57 )
58 self.assertEqual(self.instance.region, other.region)
60 def test_repr(self) -> None:
61 """Test that eval(repr(...)) round-trips."""
62 self.assertEqual(eval(repr(self.instance)).region, self.instance.region)
64 def test_to_pixel(self) -> None:
65 """Test `SkyCircle.to_pixel`, implicitly testing the
66 `PixelStencil` implementation it returns.
67 """
68 # WCS is gnomonic, with projection point at the circle center and
69 # 0.1" pixels, so 1" radius circle will have be about 20 pixels across
70 # the diameter. Make the offset arbitray but nonzero.
71 wcs = makeSkyWcs(Point2D(5.0, 7.0), self.center, makeCdMatrix(0.1 * arcseconds))
72 # Bounding box should be roughly twice the size of the circle.
73 # Make it slightly bigger in x to catch x<->y transposition bugs.
74 bbox = Box2I(Point2I(-16, -13), Point2I(26, 27))
75 # Check to_pixel(...).set_mask(...) against a brute-force
76 # implementation. The max_missing threshold (and implicitly, the
77 # max_extra=0 threshold) was set after inspection with plot=True to
78 # prevent accidental regressions.
79 _check_to_pixel(self, self.instance, wcs, bbox, max_missing=12)
81 def test_to_polygon(self) -> None:
82 """Test conversion to a `SkyPolygon`.
84 This provides a lot of test coverage for `SkyPolygon` as well,
85 but it doesn't check the non-convex case (hard, since sphgeom doesn't
86 have a non-convex polygon), and it doesn't check conversion from
87 astropy or the documented orientation conventions.
88 """
89 # Same WCS and bbox as used in test_to_pixel.
90 wcs = makeSkyWcs(Point2D(5.0, 7.0), self.center, makeCdMatrix(0.1 * arcseconds))
91 bbox = Box2I(Point2I(-16, -13), Point2I(26, 27))
92 # Convert to polygon.
93 polygon_stencil = self.instance.to_polygon()
94 # Make sure sphgeom regions at least aren't disjoint; in exact math
95 # the true circle would contain the polygon, but round-off error makes
96 # that not guaranteed with floats.
97 self.assertNotEqual(
98 self.instance.region.relate(polygon_stencil.region.getBoundingCircle()), lsst.sphgeom.DISJOINT
99 )
100 # Check the polygon's to_pixel implementation. The max_missing
101 # threshold (and implicitly, the max_extra=0 threshold) was set after
102 # inspection with plot=True to prevent accidental regressions.
103 _check_to_pixel(self, polygon_stencil, wcs, bbox, max_missing=4)
106def _brute_force_stencil_array(sky_stencil: SkyStencil, wcs: SkyWcs, bbox: Box2I) -> np.ndarray:
107 """Create a boolean Numpy array of a `SkyStencil` by transforming and
108 checking every pixel within it.
110 Parameters
111 ----------
112 sky_stencil : `SkyStencil`
113 Stencil to create an image of.
114 wcs : `SkyWcs`
115 WCS that transforms sky coordinates to pixel coordinates.
116 bbox : `Box2I`
117 Bounding box that must contain the pixel-coordinate stencil.
119 Returns
120 -------
121 array : `np.ndarray`
122 2-d boolean array of shape ``(bbox.getWidth(), bbox.getHeight())``.
123 `True` pixels are those whose WCS-transformed centers are within
124 ``sky_stencil``.
125 """
126 # Make Numpy arrays of the pixel coordinates.
127 x1 = np.arange(bbox.getBeginX(), bbox.getEndX())
128 y1 = np.arange(bbox.getBeginY(), bbox.getEndY())
129 x2, y2 = np.meshgrid(x1, y1)
130 # Stuff those into one (2, nPoints) array so we can pass it to the AST
131 # mapping inside a WCS.
132 pixels = np.zeros((2, bbox.getArea()), dtype=float)
133 pixels[0, :] = x2.flatten()
134 pixels[1, :] = y2.flatten()
135 # Transform all those points to sky coordinates, yielding another
136 # (2, nPoint) array.
137 sky = wcs.getTransform().getMapping().applyForward(pixels)
138 # Test which of those points are inside the sky stencil's sphgeom region,
139 # and reshape back to image coordinates.
140 contained = sky_stencil.region.contains(sky[0, :], sky[1, :])
141 # Reshape the boolean 'contained' array back to an image and return it.
142 return contained.reshape(bbox.getHeight(), bbox.getWidth())
145def _check_to_pixel(
146 test_case: unittest.TestCase,
147 sky_stencil: SkyStencil,
148 wcs: SkyWcs,
149 bbox: Box2I,
150 *,
151 max_missing: int = 0,
152 max_extra: int = 0,
153 plot: bool = False,
154) -> None:
155 """Test helper function that checks a `SkyStencil.to_pixel` implementation.
157 Parameters
158 ----------
159 test_case : `unittest.TestCase`
160 Test case object that provides assertion methods.
161 sky_stencil : `SkyStencil`
162 Stencil to test.
163 wcs : `SkyWcs`
164 WCS that transforms sky coordinates to pixel coordinates.
165 bbox : `Box2I`
166 Bounding box that must contain the pixel-coordinate stencil.
167 max_missing : `int`, optional
168 The number of WCS-transformed pixels that are in the stencil but not
169 masked by `PixelStencil.set_mask` must not exceed this number (defaults
170 to 0).
171 max_extra : `int`, optional
172 The number of WCS-transformed pixels that are not in the stencil but
173 are masked by `PixelStencil.set_mask` must not exceed this number
174 (defaults to 0).
175 plot : `bool`, optional
176 If `True` (`False` is default), create an interactive matplotlib image
177 of the comparison for human inspection. This should never be set
178 except when debugging or actively developing a new test, and even then
179 it probably requires that the test not be run by pytest.
180 """
181 pixel_stencil = sky_stencil.to_pixels(wcs, bbox)
182 test_case.assertTrue(bbox.contains(pixel_stencil.bbox))
183 mask = Mask(bbox)
184 mask.addMaskPlane("STENCIL")
185 bits = mask.getPlaneBitMask("STENCIL")
186 pixel_stencil.set_mask(mask, bits)
187 check_array = _brute_force_stencil_array(sky_stencil, wcs, bbox)
188 missing = np.logical_and(check_array, np.logical_not(mask.array & bits))
189 extra = np.logical_and(mask.array & bits, np.logical_not(check_array))
190 if plot:
191 from matplotlib import pyplot
193 display_array = np.zeros((bbox.getHeight(), bbox.getWidth(), 3), dtype=np.uint8)
194 display_array[:, :, 0] = 255 * check_array
195 display_array[:, :, 1] = 255 * (mask.array & bits).astype(bool)
196 pyplot.imshow(display_array, origin="lower", interpolation="nearest")
197 pyplot.title("red=check, green=SkyStencil.to_pixel, yellow=both")
198 pyplot.show()
199 test_case.assertLessEqual(sum(missing.flatten()), max_missing)
200 test_case.assertLessEqual(sum(extra.flatten()), max_extra)
203if __name__ == "__main__": 203 ↛ 204line 203 didn't jump to line 204 because the condition on line 203 was never true
204 unittest.main()