Coverage for python/lsst/images/

1# This file is part of lsst-images.

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.

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

14__all__ = ("Polygon", "Region", "RegionDomainError", "RegionSerializationModel")

16from typing import TYPE_CHECKING, Any, Literal, assert_type, overload

18import numpy as np

19import numpy.typing as npt

20import pydantic

21import pydantic_core.core_schema as pcs

22import shapely

23from pydantic.json_schema import JsonSchemaValue

25from ._geom import XY, YX, Bounds, Box

26from .describe import DescribableMixin, DescribeOptions, Report, ReportField

27from .utils import round_half_down, round_half_up

29if TYPE_CHECKING:

30 from ._geom import Bounds

31 from ._transforms import Transform

33 try:

34 from lsst.afw.geom import Polygon as LegacyPolygon

35 except ImportError:

36 type LegacyPolygon = Any # type: ignore[no-redef]

39class RegionDomainError(ValueError):

40 """Exception raised when a region is transformed to non-finite

41 coordinates, i.e. it lies outside the domain of the transform.

42 """

45class Region(DescribableMixin):

46 """A 2-d Euclidean region represented as one or more polygons with

47 optional holes.

49 Parameters

50 ----------

51 geometry

52 A polygon or multi-polygon from the Shapely library.

53 """

55 def __init__(self, geometry: shapely.Polygon | shapely.MultiPolygon) -> None:

56 self._impl = geometry

58 @property

59 def area(self) -> float:

60 """The area of the region (`float`)."""

61 return self._impl.area

63 @property

64 def bbox(self) -> Box:

65 """The integer-coordinate bounding box of the region (`Box`).

67 Because a `Box` logically contains the entirety of the pixels on its

68 boundary, but the centers of those pixels are the numerical values of

69 its bounds, the region may contain vertices that are up to 0.5 beyond

70 the integer box coordinates in either dimension.

71 """

72 x_min, y_min, x_max, y_max = self._impl.bounds

73 return Box.factory[

74 round_half_up(y_min) : round_half_down(y_max) + 1,

75 round_half_up(x_min) : round_half_down(x_max) + 1,

76 ]

78 @property

79 def wkt(self) -> str:

80 """The 'Well-Known Text' representation of this region (`str`)."""

81 return self._impl.wkt

83 @staticmethod

84 def from_wkt(wkt: str) -> Region:

85 """Construct from a 'Well-Known Text' string.

87 Parameters

88 ----------

89 wkt

90 Well-Known Text representation of the region.

91 """

92 impl = shapely.from_wkt(wkt)

93 if not isinstance(impl, shapely.Polygon | shapely.MultiPolygon): 93 ↛ 94line 93 didn't jump to line 94 because the condition on line 93 was never true

94 raise ValueError("Only Polygon and MultiPolygon geometries can be converted to Regions.")

95 return Region(impl).try_to_polygon()

97 def _describe(self, options: DescribeOptions = DescribeOptions(), /) -> Report:

98 """Return a `Report` describing this region.

100 Parameters

101 ----------

102 options : `DescribeOptions`, optional

103 Unused; accepted for interface compatibility.

104 """

105 return Report(

106 type_name="Region",

107 summary=self._impl.wkt,

108 fields=[

109 ReportField(

110 label="wkt",

111 value=self._impl.wkt,

112 repr_value=repr(self._impl.wkt),

113 positional=True,

114 )

115 ],

116 )

117

118 def __repr__(self) -> str:

119 return f"Region.from_wkt({self._impl.wkt!r})"

120

121 def __eq__(self, other: object) -> bool:

122 if isinstance(other, Region):

123 return bool(shapely.equals(self._impl, other._impl))

124 return NotImplemented

125

126 @overload

127 def contains(self, point: XY[int | float] | YX[int | float], /) -> bool: ...

128

129 @overload

130 def contains(self, point: XY[npt.ArrayLike] | YX[npt.ArrayLike], /) -> np.ndarray: ...

131

132 @overload

133 def contains(self, /, *, x: int | float, y: int | float) -> bool: ...

134

135 @overload

136 def contains(self, /, *, x: npt.ArrayLike, y: npt.ArrayLike) -> np.ndarray: ...

137

138 @overload

139 def contains(self, other: Polygon, /) -> bool: ...

140

141 def contains(

142 self,

143 other: Region | XY[Any] | YX[Any] | None = None,

144 /,

145 *,

146 x: Any = None,

147 y: Any = None,

148 ) -> bool | np.ndarray:

149 """Test whether the geometry contains the given points or another

150 geometry.

151

152 Parameters

153 ----------

154 other

155 Another geometry, or an `XY` or `YX` coordinate pair, to compare

156 to. Mutually exclusive with ``x`` and ``y``.

157 x

158 One or more X coordinates to test for containment, as a scalar or

159 any array-like. Results are broadcast against ``y``.

160 Mutually exclusive with ``other``.

161 y

162 One or more Y coordinates to test for containment, as a scalar or

163 any array-like. Results are broadcast against ``x``.

164 Mutually exclusive with ``other``.

165 """

166 match other:

167 case None:

168 if x is None or y is None: 168 ↛ 169line 168 didn't jump to line 169 because the condition on line 168 was never true

169 raise TypeError("Pass either a point or both x= and y= to 'Region.contains'.")

170 case Region():

171 if x is not None or y is not None: 171 ↛ 172line 171 didn't jump to line 172 because the condition on line 171 was never true

172 raise TypeError("'Region.contains' point argument is mutually exclusive with x= and y=.")

173 return self._impl.contains(other._impl)

174 case XY() | YX(): 174 ↛ 178line 174 didn't jump to line 178 because the pattern on line 174 always matched

175 if x is not None or y is not None:

176 raise TypeError("'Region.contains' point argument is mutually exclusive with x= and y=.")

177 x, y = other.x, other.y

178 case _:

179 raise TypeError(f"Unexpected positional argument type: {type(other)!r}.")

180 return shapely.contains_xy(self._impl, x=x, y=y)

181

182 @overload

183 def intersection(self, other: Region) -> Region: ...

184

185 @overload

186 def intersection(self, other: Box) -> Region | Box: ...

187

188 @overload

189 def intersection(self, other: Bounds) -> Bounds: ...

190

191 def intersection(self, other: Bounds) -> Bounds:

192 """Compute the intersection of this region with a `Bounds` object.

193

194 Notes

195 -----

196 Because `Region` implements the `Bounds` interface, its intersections

197 need to support all other `Bounds` objects. This is not true of other

198 `Region` point-set operations like `union` and `difference`.

199

200 Parameters

201 ----------

202 other

203 Bounds to intersect with this region.

204 """

205 from ._concrete_bounds import _intersect_region

206

207 return _intersect_region(self, other)

208

209 def union(self, other: Region) -> Region:

210 """Compute the point-set union of this region with another.

211

212 Parameters

213 ----------

214 other

215 Region to union with this one.

216 """

217 impl = shapely.union(self._impl, other._impl)

218 assert isinstance(impl, shapely.Polygon | shapely.MultiPolygon), (

219 "A union of Polygons and MultiPolygons should be one of those."

220 )

221 return Region(impl).try_to_polygon()

222

223 def difference(self, other: Region) -> Region:

224 """Compute the point-set difference of this region with another.

225

226 Parameters

227 ----------

228 other

229 Region to subtract from this one.

230 """

231 impl = shapely.difference(self._impl, other._impl)

232 assert isinstance(impl, shapely.Polygon | shapely.MultiPolygon), (

233 "A difference of Polygons and MultiPolygons should be one of those."

234 )

235 return Region(impl).try_to_polygon()

236

237 def dilated_by(self, padding: float) -> Region:

238 """Return a new region padded by the given distance on all sides.

239

240 Edges move parallel to themselves, so padding the polygon of a `Box`

241 gives the polygon of the padded box.

242

243 Parameters

244 ----------

245 padding

246 Distance to move the edges outward; a negative value moves them

247 inward.

248

249 Returns

250 -------

251 `Region`

252 The padded region, as a `Polygon` if it is a single polygon with

253 no holes.

254 """

255 impl = shapely.buffer(self._impl, padding, join_style="mitre")

256 assert isinstance(impl, shapely.Polygon | shapely.MultiPolygon), (

257 "Padding Polygons and MultiPolygons should give one of those."

258 )

259 return Region(impl).try_to_polygon()

260

261 def convex_hull(self) -> Polygon:

262 """Return the smallest convex polygon that contains the region.

263

264 Returns

265 -------

266 `Polygon`

267 The convex hull, a single polygon with no holes even when the

268 region has several pieces or holes.

269

270 Raises

271 ------

272 ValueError

273 Raised if the region has no area.

274 """

275 impl = shapely.convex_hull(self._impl)

276 if not isinstance(impl, shapely.Polygon):

277 raise ValueError(f"Region {self} has no area, so it has no convex hull.")

278 result = Region(impl).try_to_polygon()

279 assert isinstance(result, Polygon), "A convex hull should have no holes."

280 return result

281

282 def try_to_polygon(self) -> Region:

283 """If the underlying geometry is a single polygon with no holes,

284 return a `Polygon` instance holding it.

285

286 In all other cases ``self`` is returned.

287 """

288 impl = self._impl

289 if isinstance(impl, shapely.MultiPolygon) and len(impl.geoms) == 1: 289 ↛ 290line 289 didn't jump to line 290 because the condition on line 289 was never true

290 impl = impl.geoms[0]

291 if isinstance(impl, shapely.Polygon) and not impl.interiors:

292 vertices = np.array(impl.exterior.coords)

293 return Polygon(x_vertices=vertices[:, 0], y_vertices=vertices[:, 1])

294 return self

295

296 def try_to_box(self) -> Region | Box:

297 """If the underlying geometry is a rectangle that fully covers integer

298 pixels (i.e. has all vertices at half-integer positions), return the

299 equivalent `Box`.

300

301 In all other cases ``self`` is returned.

302 """

303 if Polygon.from_box(self.bbox) == self:

304 return self.bbox

305 return self

306

307 def to_shapely(self) -> shapely.Polygon | shapely.MultiPolygon:

308 """Convert to a `shapely.Polygon` or `shapely.MultiPolygon` object."""

309 return self._impl

310

311 def transform(self, transform: Transform[Any, Any]) -> Region:

312 """Transform the coordinates of the region, returning a new one.

313

314 Parameters

315 ----------

316 transform

317 Coordinate transform to apply (in the forward direction).

318

319 Returns

320 -------

321 `Region`

322 The transformed region, as a `Polygon` if it is a single polygon

323 with no holes.

324

325 Raises

326 ------

327 RegionDomainError

328 Raised if any transformed vertex is not finite, as happens when

329 the region lies outside the domain of a projection.

330

331 Notes

332 -----

333 This applies the transform to all vertices, assuming that the

334 transform is close enough to affine that the topology of the geometry

335 does not change and straight-line edges do not need to be subsampled.

336 """

337

338 def wrapper(x: np.ndarray, y: np.ndarray) -> XY[np.ndarray]:

339 result = transform.apply_forward(x=x, y=y)

340 if not (np.isfinite(result.x).all() and np.isfinite(result.y).all()):

341 raise RegionDomainError(f"Region {self} transforms to non-finite coordinates.")

342 return result

343

344 return Region(

345 # Shapely overloads don't seem to have been annotated rigorously

346 shapely.transform(self._impl, wrapper, interleaved=False) # type: ignore[arg-type]

347 ).try_to_polygon()

348

349 def serialize(self) -> RegionSerializationModel:

350 """Serialize the region to a Pydantic model.

351

352 Region serialization uses a subset of the GeoJSON specification (IETF

353 RFC 7946).

354 """

355 return RegionSerializationModel.model_validate_json(shapely.to_geojson(self._impl))

356

357 @classmethod

358 def __get_pydantic_core_schema__(

359 cls, source_type: Any, handler: pydantic.GetCoreSchemaHandler

360 ) -> pcs.CoreSchema:

361 from_model_schema = pcs.chain_schema(

362 [

363 handler(RegionSerializationModel),

364 pcs.no_info_plain_validator_function(RegionSerializationModel.deserialize),

365 ]

366 )

367 return pcs.json_or_python_schema(

368 json_schema=from_model_schema,

369 python_schema=pcs.union_schema([pcs.is_instance_schema(cls), from_model_schema]),

370 serialization=pcs.plain_serializer_function_ser_schema(cls.serialize),

371 )

372

373 @classmethod

374 def __get_pydantic_json_schema__(

375 cls, schema: pcs.CoreSchema, handler: pydantic.GetJsonSchemaHandler

376 ) -> JsonSchemaValue:

377 return handler(RegionSerializationModel.__pydantic_core_schema__)

378

379

380class Polygon(Region):

381 """A simple 2-d polygon in Euclidean coordinates, with no holes.

382

383 Parameters

384 ----------

385 x_vertices

386 The x coordinates of the vertices of the polygon.

387 y_vertices

388 The y coordinate of the vertices of the polygon.

389 """

390

391 def __init__(self, *, x_vertices: npt.ArrayLike, y_vertices: npt.ArrayLike) -> None:

392 self._vertices = np.stack(

393 [np.asarray(x_vertices).flat, np.asarray(y_vertices).flat], dtype=np.float64

394 ).transpose()

395 self._vertices.flags.writeable = False

396 super().__init__(shapely.Polygon(self._vertices))

397

398 @staticmethod

399 def from_box(box: Box) -> Polygon:

400 """Construct from an integer-coordinate box.

401

402 Parameters

403 ----------

404 box

405 Integer-coordinate box to convert to a polygon.

406

407 Notes

408 -----

409 Because the integer min and max coordinates of the box are

410 interpreted as pixel centers, these are expanded by 0.5 on all sides

411 before using them to form the polygon vertices.

412 """

413 return Polygon(

414 x_vertices=[box.x.min - 0.5, box.x.min - 0.5, box.x.max + 0.5, box.x.max + 0.5],

415 y_vertices=[box.y.min - 0.5, box.y.max + 0.5, box.y.max + 0.5, box.y.min - 0.5],

416 )

417

418 @property

419 def n_vertices(self) -> int:

420 """The number of vertices in the polygon."""

421 return self._vertices.shape[0]

422

423 @property

424 def x_vertices(self) -> np.ndarray:

425 """The x coordinates of the vertices of the polygon.

426

427 This is a read-only array; polygons are immutable.

428 """

429 return self._vertices[:, 0]

430

431 @property

432 def y_vertices(self) -> np.ndarray:

433 """The y coordinates of the vertices of the polygon.

434

435 This is a read-only array; polygons are immutable.

436 """

437 return self._vertices[:, 1]

438

439 @property

440 def centroid(self) -> XY[float]:

441 """The centroid of the polygon (`XY` [`float`])."""

442 c = self._impl.centroid

443 return XY(x=c.x, y=c.y)

444

445 def _describe(self, options: DescribeOptions = DescribeOptions(), /) -> Report:

446 """Return a `Report` describing this polygon.

447

448 Parameters

449 ----------

450 options : `DescribeOptions`, optional

451 Unused; accepted for interface compatibility.

452 """

453 return Report(

454 type_name="Polygon",

455 summary=self._impl.wkt,

456 fields=[

457 ReportField(label="x_vertices", value=self.x_vertices, repr_value=repr(self.x_vertices)),

458 ReportField(label="y_vertices", value=self.y_vertices, repr_value=repr(self.y_vertices)),

459 ],

460 )

461

462 def __repr__(self) -> str:

463 return f"Polygon(x_vertices={self.x_vertices!r}, y_vertices={self.y_vertices!r})"

464

465 def transform(self, transform: Transform[Any, Any]) -> Polygon:

466 # Docstring inherited.

467 result = super().transform(transform)

468 assert isinstance(result, Polygon), "Transforming a polygon should not change its topology."

469 return result

470

471 @staticmethod

472 def from_legacy(legacy: LegacyPolygon) -> Polygon:

473 """Convert from a legacy `lsst.afw.geom.Polygon` instance.

474

475 Parameters

476 ----------

477 legacy

478 Legacy `lsst.afw.geom.Polygon` to convert.

479 """

480 vertices = legacy.getVertices()

481 x_vertices = np.zeros(len(vertices), dtype=np.float64)

482 y_vertices = np.zeros(len(vertices), dtype=np.float64)

483 for n, point in enumerate(vertices):

484 x_vertices[n] = point.x

485 y_vertices[n] = point.y

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

487

488 def to_legacy(self) -> LegacyPolygon:

489 """Convert to a legacy `lsst.afw.geom.Polygon` instance."""

490 from lsst.afw.geom import Polygon as LegacyPolygon

491 from lsst.geom import Point2D

492

493 return LegacyPolygon([Point2D(x, y) for x, y in zip(self.x_vertices, self.y_vertices)])

494

495

496class RegionSerializationModel(pydantic.BaseModel):

497 """Serialization model for `Region` and `Polygon`.

498

499 This model is a subset of the GeoJSON specification (IETF RFC 7946).

500 """

501

502 type: Literal["Polygon", "MultiPolygon"] = pydantic.Field(description="Geometry type.")

503

504 coordinates: list[list[tuple[float, float] | list[tuple[float, float]]]] = pydantic.Field(

505 description="Vertices of the polygon or polygons."

506 )

507

508 def deserialize(self) -> Region:

509 """Deserialize into a `Region` (a `Polygon`, if possible)."""

510 region_impl = shapely.from_geojson(self.model_dump_json())

511 assert isinstance(region_impl, shapely.Polygon | shapely.MultiPolygon), (

512 "Other geometry types are not used."

513 )

514 return Region(region_impl).try_to_polygon()

515

516

517if TYPE_CHECKING:

518

519 def _test_types() -> None:

520 region = Region(shapely.Point(0, 0).buffer(1.0))

521 arr = np.zeros(3)

522

523 # Region satisfies the Bounds Protocol.

524 bounds: Bounds = region

525

526 # Region.contains: XY/YX, scalar, array-like

527 assert_type(region.contains(x=1.0, y=2.0), bool)

528 assert_type(region.contains(x=arr, y=arr), np.ndarray)

529 assert_type(region.contains(XY(1.0, 2.0)), bool)

530 assert_type(region.contains(YX(2.0, 1.0)), bool)

531 assert_type(region.contains(XY(arr, arr)), np.ndarray)

532 assert_type(region.contains(YX(arr, arr)), np.ndarray)

533

534 # Via the Bounds Protocol view, same signatures hold.

535 assert_type(bounds.contains(x=1, y=1), bool)

536 assert_type(bounds.contains(x=arr, y=arr), np.ndarray)

537 assert_type(bounds.contains(XY(1, 1)), bool)

538 assert_type(bounds.contains(YX(1, 1)), bool)

539 assert_type(bounds.contains(XY(arr, arr)), np.ndarray)

540 assert_type(bounds.contains(YX(arr, arr)), np.ndarray)

Coverage for python/lsst/images/_polygon.py: 51%

173 statements