Coverage for tests/test_polygon.py: 99%
190 statements
« prev ^ index » next coverage.py v7.17.2, created at 2026-11-10 09:42 +0000
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
16import astropy.units as u
17import numpy as np
18import pydantic
19import pytest
20import shapely
22from lsst.images import (
23 XY,
24 YX,
25 Box,
26 GeneralFrame,
27 NoOverlapError,
28 Polygon,
29 Region,
30 RegionDomainError,
31 RegionSerializationModel,
32 Transform,
33)
34from lsst.images.tests import assert_values_equal, check_bounds_contains_broadcasting
36try:
37 import lsst.afw.geom # noqa: F401
39 have_legacy = True
40except ImportError:
41 have_legacy = False
43skip_no_legacy = pytest.mark.skipif(not have_legacy, reason="lsst legacy packages could not be imported.")
46class _PolygonHolder(pydantic.BaseModel):
47 polygon: Polygon
50class _RegionHolder(pydantic.BaseModel):
51 region: Region
54def _make_polygon() -> Polygon:
55 """Return a near-box quadrilateral that is easy to reason about."""
56 x_vertices = [32.0, 31.0, 50.0, 53.0]
57 y_vertices = [-5.0, 7.0, 7.2, -4.8]
58 return Polygon(x_vertices=x_vertices, y_vertices=y_vertices)
61@dataclasses.dataclass
62class _TestRegions:
63 """Four overlapping box-polygons for region operation tests.
65 Rough layout (y increasing upwards)::
67 ┌─────┐
68 ┌───┼─┐B┌─┼────┐
69 │ A└─┼─┼─┘ ┌─┐│
70 └─────┘ │ C │D││
71 │ └─┘│
72 └──────┘
73 """
75 a: Polygon
76 b: Polygon
77 c: Polygon
78 d: Polygon
80 def __init__(self) -> None:
81 self.a = Polygon.from_box(Box.factory[3:6, 0:5])
82 self.b = Polygon.from_box(Box.factory[4:7, 3:8])
83 self.c = Polygon.from_box(Box.factory[0:6, 6:12])
84 self.d = Polygon.from_box(Box.factory[1:4, 9:10])
87def test_vertices() -> None:
88 """Test the vertices accessors."""
89 polygon = _make_polygon()
90 x_vertices = [32.0, 31.0, 50.0, 53.0]
91 y_vertices = [-5.0, 7.0, 7.2, -4.8]
92 assert polygon.n_vertices == 4
93 np.testing.assert_array_equal(polygon.x_vertices, np.asarray(x_vertices))
94 np.testing.assert_array_equal(polygon.y_vertices, np.asarray(y_vertices))
95 with pytest.raises(ValueError, match="read-only"):
96 polygon.x_vertices[0] = 0.0
97 with pytest.raises(ValueError, match="read-only"):
98 polygon.y_vertices[0] = 0.0
101def test_boxes() -> None:
102 """Test 'from_box', the ``area`` property, and the 'contains' method
103 with polygon arguments.
104 """
105 polygon = _make_polygon()
106 small = Polygon.from_box(Box.factory[-3:3, 40:45])
107 assert small.area == 30.0
108 assert small.bbox == Box.factory[-3:3, 40:45]
109 assert polygon.contains(small)
110 assert not small.contains(polygon)
111 assert small.centroid == XY(x=42.0, y=-0.5)
112 big = Polygon.from_box(Box.factory[-10:10, 20:60])
113 assert big.area == 800.0
114 assert not polygon.contains(big)
115 assert big.contains(polygon)
116 medium = Polygon.from_box(Box.factory[-4:8, 31:52])
117 assert medium.area == 252.0
118 assert not polygon.contains(medium)
119 assert not medium.contains(polygon)
120 assert polygon.contains(polygon)
123def test_transform() -> None:
124 """Test applying a coordinate transform to a polygon."""
125 polygon = _make_polygon()
126 matrix = np.array([[0.4, 0.25], [-0.20, 0.6]])
127 t = Transform.affine(GeneralFrame(unit=u.pix), GeneralFrame(unit=u.pix), matrix)
128 tp = polygon.transform(t)
129 assert isinstance(tp, Polygon)
130 assert_values_equal(tp.area, polygon.area * np.linalg.det(matrix), rtol=1e-5)
131 xyt = t.apply_forward(x=polygon.x_vertices, y=polygon.y_vertices)
132 # Slicing below is because shapely sometimes adds a duplicate closing
133 # vertex.
134 assert_values_equal(tp.x_vertices[: len(xyt.x)], xyt.x, rtol=1e-14)
135 assert_values_equal(tp.y_vertices[: len(xyt.y)], xyt.y, rtol=1e-14)
138def test_transform_outside_domain() -> None:
139 """Test that a transform yielding non-finite vertices raises a
140 `RegionDomainError` instead of producing a broken region.
141 """
142 polygon = _make_polygon()
143 matrix = np.array([[1e308, 0.0], [0.0, 1e308]])
144 t = Transform.affine(GeneralFrame(unit=u.pix), GeneralFrame(unit=u.pix), matrix)
145 assert not np.isfinite(t.apply_forward(x=polygon.x_vertices, y=polygon.y_vertices).x).all()
146 with pytest.raises(RegionDomainError):
147 polygon.transform(t)
150def test_contains_points() -> None:
151 """Test the 'contains' method with points."""
152 polygon = _make_polygon()
153 assert polygon.contains(x=40.0, y=0.0)
154 assert not polygon.contains(x=0.0, y=0.0)
155 assert not polygon.contains(x=40.0, y=10.0)
156 np.testing.assert_array_equal(
157 polygon.contains(x=np.array([40.0, 0.0, 40.0]), y=np.array([0.0, 0.0, 10.0])),
158 np.array([True, False, False]),
159 )
162def test_contains_points_xy_yx() -> None:
163 """Verify that Region.contains accepts XY and YX positional arguments."""
164 polygon = _make_polygon()
165 # Scalar XY and YX — results must match the keyword form.
166 assert polygon.contains(XY(x=40.0, y=0.0)) == polygon.contains(x=40.0, y=0.0)
167 assert polygon.contains(YX(y=0.0, x=40.0)) == polygon.contains(x=40.0, y=0.0)
168 assert not polygon.contains(XY(x=0.0, y=0.0))
169 assert not polygon.contains(YX(y=0.0, x=0.0))
170 # Array XY and YX.
171 xv = np.array([40.0, 0.0, 40.0])
172 yv = np.array([0.0, 0.0, 10.0])
173 np.testing.assert_array_equal(polygon.contains(XY(xv, yv)), polygon.contains(x=xv, y=yv))
174 np.testing.assert_array_equal(polygon.contains(YX(yv, xv)), polygon.contains(x=xv, y=yv))
175 # Mixing positional point with keyword x= or y= must raise TypeError.
176 with pytest.raises(TypeError):
177 polygon.contains(XY(x=40.0, y=0.0), x=40.0)
178 with pytest.raises(TypeError):
179 polygon.contains(YX(y=0.0, x=40.0), y=0.0)
182def test_region_contains_broadcasting() -> None:
183 """Test that Region.contains broadcasts like a numpy ufunc."""
184 check_bounds_contains_broadcasting(_make_polygon())
187def test_polygon_io() -> None:
188 """Test serialization and stringification."""
189 polygon = _make_polygon()
190 assert (
191 RegionSerializationModel.model_validate_json(polygon.serialize().model_dump_json()).deserialize()
192 == polygon
193 )
194 assert Polygon.from_wkt(polygon.wkt) == polygon
195 assert Polygon.from_wkt(str(polygon)) == polygon
196 assert eval(repr(polygon), {"array": np.array, "Polygon": Polygon}) == polygon
199def test_polygon_model_field() -> None:
200 """Test that we can use a Polygon directly as a Pydantic model field."""
201 polygon = _make_polygon()
202 holder = _PolygonHolder(polygon=polygon)
203 assert polygon == holder.model_validate_json(holder.model_dump_json()).polygon
204 assert (
205 _PolygonHolder.model_json_schema()["properties"]["polygon"]
206 == RegionSerializationModel.model_json_schema()
207 )
210@skip_no_legacy
211def test_polygon_legacy() -> None:
212 """Test conversion to/from lsst.afw.geom.Polygon."""
213 polygon = _make_polygon()
214 legacy_polygon = polygon.to_legacy()
215 assert legacy_polygon.calculateArea() == polygon.area
216 assert Polygon.from_legacy(legacy_polygon) == polygon
219def test_intersection() -> None:
220 """Test region intersection."""
221 regions = _TestRegions()
222 # Usual case:
223 assert regions.a.intersection(regions.b) == Polygon.from_box(Box.factory[4:6, 3:5])
224 # No-overlap case:
225 with pytest.raises(NoOverlapError):
226 regions.a.intersection(regions.c)
227 # LHS fully contains RHS:
228 assert regions.c.intersection(regions.d) == regions.d
229 # Intersections with the boxes themselves should return boxes:
230 assert regions.a.intersection(regions.b.bbox) == Box.factory[4:6, 3:5]
231 assert regions.a.bbox.intersection(regions.b) == Box.factory[4:6, 3:5]
232 assert (
233 # A Box is not possible when the result is not simple.
234 regions.a.union(regions.c).intersection(regions.b.bbox)
235 == regions.a.union(regions.c).intersection(regions.b)
236 )
239def test_union() -> None:
240 """Test region union."""
241 regions = _TestRegions()
242 # Usual case:
243 assert regions.a.union(regions.b).bbox == Box.factory[3:7, 0:8]
244 assert regions.a.union(regions.b).area == 15 + 15 - 4
245 # Operands are disjoint, so union is not a single Polygon:
246 assert not isinstance(regions.a.union(regions.c), Polygon)
247 assert regions.a.union(regions.c).area == regions.a.area + regions.c.area
248 # LHS fully contains RHS:
249 assert regions.c.union(regions.d) == regions.c
252def test_difference() -> None:
253 """Test region difference."""
254 regions = _TestRegions()
255 # Usual case:
256 assert regions.a.difference(regions.b).bbox == regions.a.bbox
257 assert regions.a.difference(regions.b).area == 15 - 4
258 # Operands are disjoint, so difference is just the LHS.
259 assert regions.a.difference(regions.c) == regions.a
260 # LHS fully contains RHS -> polygon with hole is not a Polygon:
261 assert not isinstance(regions.c.difference(regions.d), Polygon)
262 assert regions.c.difference(regions.d).bbox == regions.c.bbox
263 assert regions.c.difference(regions.d).area == regions.c.area - regions.d.area
264 # RHS fully contains LHS -> region is empty.
265 assert regions.d.difference(regions.d).area == 0
268def test_dilated_by() -> None:
269 """Test padding a region."""
270 regions = _TestRegions()
271 # The polygon of a box pads to the polygon of the padded box, both ways.
272 box = Box.factory[0:6, 6:12]
273 assert regions.c.dilated_by(2.0) == Polygon.from_box(box.dilated_by(2))
274 assert regions.c.dilated_by(-1.0) == Polygon.from_box(box.dilated_by(-1))
275 # A tilted quadrilateral keeps its sharp corners, so padding it inward
276 # again gives it back.
277 polygon = _make_polygon()
278 padded = polygon.dilated_by(0.5)
279 assert isinstance(padded, Polygon)
280 assert padded.contains(polygon)
281 assert shapely.hausdorff_distance(padded.dilated_by(-0.5).to_shapely(), polygon.to_shapely()) < 1e-12
282 # Padding closes the gap between disjoint operands.
283 assert not isinstance(regions.a.union(regions.c), Polygon)
284 assert isinstance(regions.a.union(regions.c).dilated_by(1.0), Polygon)
285 # Padding inward by half the width or more leaves nothing.
286 assert regions.d.dilated_by(-0.5).area == 0
289def test_convex_hull() -> None:
290 """Test the convex hull of a region."""
291 regions = _TestRegions()
292 # A convex polygon is its own hull.
293 polygon = _make_polygon()
294 assert polygon.convex_hull() == polygon
295 # Disjoint pieces have a single polygon as their hull.
296 union = regions.a.union(regions.c)
297 hull = union.convex_hull()
298 assert isinstance(hull, Polygon)
299 assert hull.contains(regions.a)
300 assert hull.contains(regions.c)
301 assert hull.bbox == union.bbox
302 # The hull of a region with a hole fills the hole.
303 assert regions.c.difference(regions.d).convex_hull() == regions.c
304 # An empty region has no hull.
305 with pytest.raises(ValueError, match="no area"):
306 regions.d.difference(regions.d).convex_hull()
309def test_region_io() -> None:
310 """Test serialization and stringification of non-polygon regions."""
311 regions = _TestRegions()
312 # A two-polygon region with a hole:
313 region = regions.a.union(regions.c).difference(regions.d)
314 assert (
315 RegionSerializationModel.model_validate_json(region.serialize().model_dump_json()).deserialize()
316 == region
317 )
318 assert Region.from_wkt(region.wkt) == region
319 assert Region.from_wkt(str(region)) == region
320 assert eval(repr(region), {"Region": Region}) == region
323def test_region_model_field() -> None:
324 """Test that we can use a Region directly as a Pydantic model field."""
325 regions = _TestRegions()
326 region = regions.a.union(regions.c).difference(regions.d)
327 holder = _RegionHolder(region=region)
328 assert region == holder.model_validate_json(holder.model_dump_json()).region
329 assert (
330 _RegionHolder.model_json_schema()["properties"]["region"]
331 == RegionSerializationModel.model_json_schema()
332 )
335def test_region_polygon_repr_str_pinned() -> None:
336 """Region and Polygon str/repr match their documented forms."""
337 poly = Polygon(x_vertices=[0.0, 4.0, 4.0, 0.0], y_vertices=[0.0, 0.0, 5.0, 5.0])
338 assert repr(poly) == f"Polygon(x_vertices={poly.x_vertices!r}, y_vertices={poly.y_vertices!r})"
339 # str is inherited from Region and uses WKT.
340 assert str(poly) == poly.wkt
342 # Region.from_wkt on a simple polygon returns a Polygon via try_to_polygon.
343 # To test Region.__str__/__repr__ directly, construct a true Region
344 # (multi-polygon), which is not coerced to Polygon.
345 multi = shapely.MultiPolygon(
346 [
347 shapely.Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]),
348 shapely.Polygon([(2, 2), (3, 2), (3, 3), (2, 3)]),
349 ]
350 )
351 region = Region(multi)
352 assert str(region) == region.wkt
353 assert repr(region) == f"Region.from_wkt({region.wkt!r})"