Coverage for tests/test_polygon.py: 99%

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. 

11 

12from __future__ import annotations 

13 

14import dataclasses 

15 

16import astropy.units as u 

17import numpy as np 

18import pydantic 

19import pytest 

20 

21from lsst.images import ( 

22 XY, 

23 YX, 

24 Box, 

25 GeneralFrame, 

26 NoOverlapError, 

27 Polygon, 

28 Region, 

29 RegionSerializationModel, 

30 Transform, 

31) 

32from lsst.images.tests import assert_close, check_bounds_contains_broadcasting 

33 

34try: 

35 import lsst.afw.geom # noqa: F401 

36 

37 have_legacy = True 

38except ImportError: 

39 have_legacy = False 

40 

41skip_no_legacy = pytest.mark.skipif(not have_legacy, reason="lsst legacy packages could not be imported.") 

42 

43 

44class _PolygonHolder(pydantic.BaseModel): 

45 polygon: Polygon 

46 

47 

48class _RegionHolder(pydantic.BaseModel): 

49 region: Region 

50 

51 

52def _make_polygon() -> Polygon: 

53 """Return a near-box quadrilateral that is easy to reason about.""" 

54 x_vertices = [32.0, 31.0, 50.0, 53.0] 

55 y_vertices = [-5.0, 7.0, 7.2, -4.8] 

56 return Polygon(x_vertices=x_vertices, y_vertices=y_vertices) 

57 

58 

59@dataclasses.dataclass 

60class _TestRegions: 

61 """Four overlapping box-polygons for region operation tests. 

62 

63 Rough layout (y increasing upwards):: 

64 

65 ┌─────┐ 

66 ┌───┼─┐B┌─┼────┐ 

67 │ A└─┼─┼─┘ ┌─┐│ 

68 └─────┘ │ C │D││ 

69 │ └─┘│ 

70 └──────┘ 

71 """ 

72 

73 a: Polygon 

74 b: Polygon 

75 c: Polygon 

76 d: Polygon 

77 

78 def __init__(self) -> None: 

79 self.a = Polygon.from_box(Box.factory[3:6, 0:5]) 

80 self.b = Polygon.from_box(Box.factory[4:7, 3:8]) 

81 self.c = Polygon.from_box(Box.factory[0:6, 6:12]) 

82 self.d = Polygon.from_box(Box.factory[1:4, 9:10]) 

83 

84 

85def test_vertices() -> None: 

86 """Test the vertices accessors.""" 

87 polygon = _make_polygon() 

88 x_vertices = [32.0, 31.0, 50.0, 53.0] 

89 y_vertices = [-5.0, 7.0, 7.2, -4.8] 

90 assert polygon.n_vertices == 4 

91 np.testing.assert_array_equal(polygon.x_vertices, np.asarray(x_vertices)) 

92 np.testing.assert_array_equal(polygon.y_vertices, np.asarray(y_vertices)) 

93 with pytest.raises(ValueError): 

94 polygon.x_vertices[0] = 0.0 

95 with pytest.raises(ValueError): 

96 polygon.y_vertices[0] = 0.0 

97 

98 

99def test_boxes() -> None: 

100 """Test 'from_box', the ``area`` property, and the 'contains' method 

101 with polygon arguments. 

102 """ 

103 polygon = _make_polygon() 

104 small = Polygon.from_box(Box.factory[-3:3, 40:45]) 

105 assert small.area == 30.0 

106 assert small.bbox == Box.factory[-3:3, 40:45] 

107 assert polygon.contains(small) 

108 assert not small.contains(polygon) 

109 assert small.centroid == XY(x=42.0, y=-0.5) 

110 big = Polygon.from_box(Box.factory[-10:10, 20:60]) 

111 assert big.area == 800.0 

112 assert not polygon.contains(big) 

113 assert big.contains(polygon) 

114 medium = Polygon.from_box(Box.factory[-4:8, 31:52]) 

115 assert medium.area == 252.0 

116 assert not polygon.contains(medium) 

117 assert not medium.contains(polygon) 

118 assert polygon.contains(polygon) 

119 

120 

121def test_transform() -> None: 

122 """Test applying a coordinate transform to a polygon.""" 

123 polygon = _make_polygon() 

124 matrix = np.array([[0.4, 0.25], [-0.20, 0.6]]) 

125 t = Transform.affine(GeneralFrame(unit=u.pix), GeneralFrame(unit=u.pix), matrix) 

126 tp = polygon.transform(t) 

127 assert isinstance(tp, Polygon) 

128 assert_close(tp.area, polygon.area * np.linalg.det(matrix)) 

129 xyt = t.apply_forward(x=polygon.x_vertices, y=polygon.y_vertices) 

130 # Slicing below is because shapely sometimes adds a duplicate closing 

131 # vertex. 

132 assert_close(tp.x_vertices[: len(xyt.x)], xyt.x) 

133 assert_close(tp.y_vertices[: len(xyt.y)], xyt.y) 

134 

135 

136def test_contains_points() -> None: 

137 """Test the 'contains' method with points.""" 

138 polygon = _make_polygon() 

139 assert polygon.contains(x=40.0, y=0.0) 

140 assert not polygon.contains(x=0.0, y=0.0) 

141 assert not polygon.contains(x=40.0, y=10.0) 

142 np.testing.assert_array_equal( 

143 polygon.contains(x=np.array([40.0, 0.0, 40.0]), y=np.array([0.0, 0.0, 10.0])), 

144 np.array([True, False, False]), 

145 ) 

146 

147 

148def test_contains_points_xy_yx() -> None: 

149 """Verify that Region.contains accepts XY and YX positional arguments.""" 

150 polygon = _make_polygon() 

151 # Scalar XY and YX — results must match the keyword form. 

152 assert polygon.contains(XY(x=40.0, y=0.0)) == polygon.contains(x=40.0, y=0.0) 

153 assert polygon.contains(YX(y=0.0, x=40.0)) == polygon.contains(x=40.0, y=0.0) 

154 assert not polygon.contains(XY(x=0.0, y=0.0)) 

155 assert not polygon.contains(YX(y=0.0, x=0.0)) 

156 # Array XY and YX. 

157 xv = np.array([40.0, 0.0, 40.0]) 

158 yv = np.array([0.0, 0.0, 10.0]) 

159 np.testing.assert_array_equal(polygon.contains(XY(xv, yv)), polygon.contains(x=xv, y=yv)) 

160 np.testing.assert_array_equal(polygon.contains(YX(yv, xv)), polygon.contains(x=xv, y=yv)) 

161 # Mixing positional point with keyword x= or y= must raise TypeError. 

162 with pytest.raises(TypeError): 

163 polygon.contains(XY(x=40.0, y=0.0), x=40.0) 

164 with pytest.raises(TypeError): 

165 polygon.contains(YX(y=0.0, x=40.0), y=0.0) 

166 

167 

168def test_region_contains_broadcasting() -> None: 

169 """Test that Region.contains broadcasts like a numpy ufunc.""" 

170 check_bounds_contains_broadcasting(_make_polygon()) 

171 

172 

173def test_polygon_io() -> None: 

174 """Test serialization and stringification.""" 

175 polygon = _make_polygon() 

176 assert ( 

177 RegionSerializationModel.model_validate_json(polygon.serialize().model_dump_json()).deserialize() 

178 == polygon 

179 ) 

180 assert Polygon.from_wkt(polygon.wkt) == polygon 

181 assert Polygon.from_wkt(str(polygon)) == polygon 

182 assert eval(repr(polygon), {"array": np.array, "Polygon": Polygon}) == polygon 

183 

184 

185def test_polygon_model_field() -> None: 

186 """Test that we can use a Polygon directly as a Pydantic model field.""" 

187 polygon = _make_polygon() 

188 holder = _PolygonHolder(polygon=polygon) 

189 assert polygon == holder.model_validate_json(holder.model_dump_json()).polygon 

190 assert ( 

191 _PolygonHolder.model_json_schema()["properties"]["polygon"] 

192 == RegionSerializationModel.model_json_schema() 

193 ) 

194 

195 

196@skip_no_legacy 

197def test_polygon_legacy() -> None: 

198 """Test conversion to/from lsst.afw.geom.Polygon.""" 

199 polygon = _make_polygon() 

200 legacy_polygon = polygon.to_legacy() 

201 assert legacy_polygon.calculateArea() == polygon.area 

202 assert Polygon.from_legacy(legacy_polygon) == polygon 

203 

204 

205def test_intersection() -> None: 

206 """Test region intersection.""" 

207 regions = _TestRegions() 

208 # Usual case: 

209 assert regions.a.intersection(regions.b) == Polygon.from_box(Box.factory[4:6, 3:5]) 

210 # No-overlap case: 

211 with pytest.raises(NoOverlapError): 

212 regions.a.intersection(regions.c) 

213 # LHS fully contains RHS: 

214 assert regions.c.intersection(regions.d) == regions.d 

215 # Intersections with the boxes themselves should return boxes: 

216 assert regions.a.intersection(regions.b.bbox) == Box.factory[4:6, 3:5] 

217 assert regions.a.bbox.intersection(regions.b) == Box.factory[4:6, 3:5] 

218 assert ( 

219 # A Box is not possible when the result is not simple. 

220 regions.a.union(regions.c).intersection(regions.b.bbox) 

221 == regions.a.union(regions.c).intersection(regions.b) 

222 ) 

223 

224 

225def test_union() -> None: 

226 """Test region union.""" 

227 regions = _TestRegions() 

228 # Usual case: 

229 assert regions.a.union(regions.b).bbox == Box.factory[3:7, 0:8] 

230 assert regions.a.union(regions.b).area == 15 + 15 - 4 

231 # Operands are disjoint, so union is not a single Polygon: 

232 assert not isinstance(regions.a.union(regions.c), Polygon) 

233 assert regions.a.union(regions.c).area == regions.a.area + regions.c.area 

234 # LHS fully contains RHS: 

235 assert regions.c.union(regions.d) == regions.c 

236 

237 

238def test_difference() -> None: 

239 """Test region difference.""" 

240 regions = _TestRegions() 

241 # Usual case: 

242 assert regions.a.difference(regions.b).bbox == regions.a.bbox 

243 assert regions.a.difference(regions.b).area == 15 - 4 

244 # Operands are disjoint, so difference is just the LHS. 

245 assert regions.a.difference(regions.c) == regions.a 

246 # LHS fully contains RHS -> polygon with hole is not a Polygon: 

247 assert not isinstance(regions.c.difference(regions.d), Polygon) 

248 assert regions.c.difference(regions.d).bbox == regions.c.bbox 

249 assert regions.c.difference(regions.d).area == regions.c.area - regions.d.area 

250 # RHS fully contains LHS -> region is empty. 

251 assert regions.d.difference(regions.d).area == 0 

252 

253 

254def test_region_io() -> None: 

255 """Test serialization and stringification of non-polygon regions.""" 

256 regions = _TestRegions() 

257 # A two-polygon region with a hole: 

258 region = regions.a.union(regions.c).difference(regions.d) 

259 assert ( 

260 RegionSerializationModel.model_validate_json(region.serialize().model_dump_json()).deserialize() 

261 == region 

262 ) 

263 assert Region.from_wkt(region.wkt) == region 

264 assert Region.from_wkt(str(region)) == region 

265 assert eval(repr(region), {"Region": Region}) == region 

266 

267 

268def test_region_model_field() -> None: 

269 """Test that we can use a Region directly as a Pydantic model field.""" 

270 regions = _TestRegions() 

271 region = regions.a.union(regions.c).difference(regions.d) 

272 holder = _RegionHolder(region=region) 

273 assert region == holder.model_validate_json(holder.model_dump_json()).region 

274 assert ( 

275 _RegionHolder.model_json_schema()["properties"]["region"] 

276 == RegionSerializationModel.model_json_schema() 

277 )