Coverage for python / lsst / dax / images / cutout / stencils.py: 44%
136 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-18 09:19 +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
24__all__ = (
25 "PixelPolygon",
26 "PixelStencil",
27 "SkyCircle",
28 "SkyPolygon",
29 "SkyStencil",
30 "StencilNotContainedError",
31)
33import struct
34from abc import ABC, abstractmethod
35from collections.abc import Iterable
36from hashlib import blake2b
38import astropy.coordinates
39import numpy as np
41import lsst.afw.detection
42import lsst.afw.geom.ellipses
43import lsst.afw.geom.polygon
44import lsst.sphgeom
45from lsst.afw.geom import SkyWcs, makeCdMatrix, makeSkyWcs
46from lsst.afw.image import Mask
47from lsst.daf.base import PropertyList
48from lsst.geom import Angle, Box2D, Box2I, Point2D, SpherePoint, radians
51class StencilNotContainedError(RuntimeError):
52 """Exception that may be raised when a stencil is not with a desired
53 bounding box.
54 """
57class PixelStencil(ABC):
58 """An image cutout stencil defined in pixel coordinates."""
60 @property
61 @abstractmethod
62 def bbox(self) -> Box2I:
63 """Bounding box of this stencil in pixel coordinates."""
64 raise NotImplementedError()
66 @abstractmethod
67 def set_mask(self, mask: Mask, bits: int) -> None:
68 """Set mask planes for the pixels covered by this stencil.
70 Parameters
71 ----------
72 mask : `Mask`
73 Mask to modify in-place.
74 bits : `int`
75 Integer bitmask to bitwise-OR into pixels covered by the stencil.
77 Notes
78 -----
79 The what "pixels covered by" means in detail is implementation-defined,
80 at least at present.
81 """
82 raise NotImplementedError()
85class PixelPolygon(PixelStencil):
86 """A pixel-coordinate stencil backed by a polygon.
88 Parameters
89 ----------
90 polygon : `lsst.afw.geom.polygon.Polygon`
91 Backing polygon.
92 bbox : `Box2I`, optional
93 If provided (`None` is default), a pixel bounding box to clip
94 the polygon to.
96 Notes
97 -----
98 This class is backed by a polygon with straight-line edges in the image
99 pixel coordinate system; this corresponds exactly to a great-circle polygon
100 on the sky only for gnomonic projections.
102 The `set_mask` implementation for this class currently masks pixels that
103 are more than 50% covered by the polygon, which may not be the same as
104 pixels whose centers are contained by the polygon when a vertex lies within
105 a pixel. This may change in the future.
106 """
108 def __init__(self, polygon: lsst.afw.geom.polygon.Polygon, bbox: Box2I | None = None):
109 self._polygon = polygon
110 self._bbox = Box2I(polygon.getBBox())
111 if bbox is not None:
112 self._bbox.clip(bbox)
113 self._polygon = self._polygon.intersectionSingle(Box2D(self._bbox))
115 @property
116 def bbox(self) -> Box2I:
117 # Docstring inherited.
118 return self._bbox
120 def set_mask(self, mask: Mask, bits: int) -> None:
121 # Docstring inherited.
122 image = self._polygon.createImage(self.bbox)
123 submask = mask[self.bbox]
124 submask.array[:, :] |= bits * (image.array > 0.5)
127class SkyStencil(ABC):
128 """An image cutout stencil defined in sky (ICRS) coordinates."""
130 @abstractmethod
131 def to_pixels(self, wcs: SkyWcs, bbox: Box2I) -> PixelStencil:
132 """Transform to a pixel-coordinate set of spans.
134 Parameters
135 ----------
136 wcs : `SkyWcs`
137 Mapping from sky coordinates to pixel coordinates.
138 bbox : `Box2I`
139 Bounds that the returned stencil must lie within.
141 Returns
142 -------
143 pixels : `PixelStencil`
144 Pixel-coordinate stencil object. `PixelStencil.bbox` is guaranteed
145 to be contained by the given ``bbox`` for this object.
147 Raises
148 ------
149 StencilNotContainedError
150 May be raised if the pixel-coordinate stencil does not lie within
151 ``bbox``. Implementations may also clip instead.
153 Notes
154 -----
155 This operation may be an approximation; see concrete class
156 documentation for additional information.
157 """
158 raise NotImplementedError()
160 @property
161 @abstractmethod
162 def region(self) -> lsst.sphgeom.Region:
163 """A `lsst.sphgeom.Region` that bounds this stencil on the sky."""
164 raise NotImplementedError()
166 @abstractmethod
167 def to_fits_metadata(self, metadata: PropertyList) -> None:
168 """Write FITS header entries that describe the stencil.
170 Parameters
171 ----------
172 metadata : `PropertyList`
173 Metadata, to be modified in place.
174 """
175 raise NotImplementedError()
177 @property
178 @abstractmethod
179 def fingerprint(self) -> bytes:
180 """A 16-byte blob that is unique to this stencil.
182 This may be a hash or a reversible encoding.
183 """
184 raise NotImplementedError()
187class SkyCircle(SkyStencil):
188 """A sky-coordinate circular stencil.
190 Parameters
191 ----------
192 center : `SpherePoint`
193 The center of the circle, in ICRS (ra, dec).
194 radius : `Angle`
195 Radius of the circle.
196 clip : `bool`, optional
197 If `True` (`False` is default), clip pixel stencils returned by
198 `to_pixels` instead of raising `StencilNotContainedError`.
200 Notes
201 -----
202 The `to_pixels` implementation for this class converts the sky region to a
203 polygon with approximately arcsecond vertices (but a minimum of 12
204 vertices), and then converts that to pixel coordinates.
205 """
207 MAX_POLYGON_VERTICES = 64
209 def __init__(self, center: SpherePoint, radius: Angle, clip: bool = False):
210 self._center = center
211 self._radius = radius
212 self._clip = clip
214 def __repr__(self) -> str:
215 return f"SkyCircle({self._center!r}, {self._radius!r}, clip={self._clip!r})"
217 @classmethod
218 def from_astropy(
219 cls, center: astropy.coordinates.SkyCoord, radius: astropy.coordinates.Angle, clip: bool = False
220 ) -> SkyCircle:
221 """Construct from `astropy.coordinates` arguments.
223 Parameters
224 ----------
225 center : `astropy.coordinates.SkyCoord`
226 The center of the circle, in ICRS (ra, dec).
227 radius : `Angle`
228 Radius of the circle.
229 clip : `bool`, optional
230 If `True` (`False` is default), clip pixel stencils returned by
231 `to_pixels` instead of raising `StencilNotContainedError`.
233 Returns
234 -------
235 stencil : `SkyCircle`
236 Circular stencil.
238 Raises
239 ------
240 ValueError
241 Raised if ``center`` or ``radius`` is not scalar-valued.
242 """
243 return cls(
244 center=_spherepoint_from_astropy(center),
245 radius=_angle_from_astropy(radius),
246 clip=clip,
247 )
249 @classmethod
250 def from_sphgeom(cls, circle: lsst.sphgeom.Circle, clip: bool = False) -> SkyCircle:
251 """Construct from a `lsst.sphgeom.Circle` instance."""
252 return cls(SpherePoint(circle.getCenter()), Angle(circle.getOpeningAngle()), clip=clip)
254 def to_polygon(self, n_vertices: int = 16) -> SkyPolygon:
255 """Return a polygon sky stencil that approximates this circle.
257 Parameters
258 ----------
259 n_vertices : `int`, optional
260 Number of polygon vertices in the approximation.
262 Returns
263 -------
264 polygon : `SkyPolygon`
265 Polygon approximation.
267 Notes
268 -----
269 For large circles and/or highly nonlinear projections, this polygon
270 approximation can be mapped much more accurately to pixel coordinates.
271 """
272 factor = (2 * np.pi / n_vertices) * radians
273 return SkyPolygon(
274 (self._center.offset(b * factor, self._radius) for b in range(n_vertices)), clip=self._clip
275 )
277 def to_pixels(self, wcs: SkyWcs, bbox: Box2I) -> PixelStencil:
278 # Docstring inherited.
279 # convert to a polygon with ~arcsecond vertices
280 circumference = 2 * np.pi * np.sin(self._radius.asRadians()) * radians
281 n_vertices = min(max(16, int(np.round(circumference.asArcseconds()))), self.MAX_POLYGON_VERTICES)
282 return self.to_polygon(n_vertices).to_pixels(wcs, bbox)
284 @property
285 def region(self) -> lsst.sphgeom.Region:
286 # Docstring inherited.
287 return lsst.sphgeom.Circle(self._center.getVector(), self._radius)
289 def to_fits_metadata(self, metadata: PropertyList) -> None:
290 # Docstring inherited.
291 metadata.set("ST_TYPE", "CIRCLE", "Type of stencil used to create this cutout.")
292 metadata.set("ST_RA", self._center.getRa().asDegrees(), "Circle center right ascension in degrees.")
293 metadata.set("ST_DEC", self._center.getDec().asDegrees(), "Circle center declination in degrees.")
294 metadata.set("ST_RAD", self._radius.asDegrees(), "Circle radius in degrees.")
296 @property
297 def fingerprint(self) -> bytes:
298 # Docstring inherited.
299 hasher = blake2b(digest_size=16)
300 hasher.update(b"CIRCLE")
301 hasher.update(struct.pack("!d", self._center.getRa().asRadians()))
302 hasher.update(struct.pack("!d", self._center.getDec().asRadians()))
303 hasher.update(struct.pack("!d", self._radius.asRadians()))
304 return hasher.digest()
307class SkyPolygon(SkyStencil):
308 """A sky-coordinate stencil in the shape of a great-circle polygon.
310 Parameters
311 ----------
312 vertices : `Iterable` [ `SpherePoint` ]
313 Vertices of the polygon, CCW when looking out from the origin.
314 Implicitly closed (the first vertex should not be duplicated as the
315 last).
316 clip : `bool`, optional
317 If `True` (`False` is default), clip pixel stencils returned by
318 `to_pixels` instead of raising `StencilNotContainedError`.
320 Notes
321 -----
322 Vertex orientation is not checked at construction, and incorrect
323 orientation may result in unspecified failures in `to_pixels`.
324 """
326 def __init__(self, vertices: Iterable[SpherePoint], clip: bool = False):
327 self._vertices = tuple(vertices)
328 self._clip = clip
330 @classmethod
331 def from_astropy(cls, vertices: astropy.coordinates.SkyCoord, clip: bool = False) -> SkyPolygon:
332 """Construct from `astropy.coordinates` arguments.
334 Parameters
335 ----------
336 vertices : `astropy.coordinates.SkyCoord`
337 Array of vertices. CCW when looking out from the origin.
338 Implicitly closed (the first vertex should not be duplicated as the
339 last).
340 clip : `bool`, optional
341 If `True` (`False` is default), clip pixel stencils returned by
342 `to_pixels` instead of raising `StencilNotContainedError`.
344 Returns
345 -------
346 stencil : `SkyPolygon`
347 Polygon stencil.
348 """
349 return cls((_spherepoint_from_astropy(v) for v in vertices), clip=clip)
351 def to_pixels(self, wcs: SkyWcs, bbox: Box2I) -> PixelStencil:
352 # Docstring inherited.
353 pixel_vertices = wcs.skyToPixel(self._vertices)
354 pixel_vertices.append(pixel_vertices[0]) # afw.polygon expects to be explicitly closed
355 # Input sky coordinates should be CCW looking out, and afw.polygon
356 # expects pixel coordinates to be CW. First question is whether the
357 # WCS has a parity flip, which is frustratingly not something we can
358 # ask it (or the underlying AST object) directly, so we compute the
359 # determinant of a linear approximation. It doesn't matter where,
360 # since the polarity can't actually change with position
361 affine = wcs.linearizeSkyToPixel(self._vertices[0], radians)
362 if affine.getLinear().computeDeterminant() > 1:
363 # No parity flip, so we have to reverse the vertices ourselves.
364 pixel_vertices = reversed(pixel_vertices)
365 result = PixelPolygon(
366 lsst.afw.geom.polygon.Polygon(pixel_vertices), bbox=(bbox if self._clip else None)
367 )
368 if not self._clip and not bbox.contains(result.bbox):
369 raise StencilNotContainedError(
370 f"{self} has bbox {result.bbox} in pixel coordinates, which is not within {bbox}."
371 )
372 return result
374 @property
375 def region(self) -> lsst.sphgeom.Region:
376 # Docstring inherited.
377 return lsst.sphgeom.ConvexPolygon([v.getVector() for v in self._vertices])
379 def to_fits_metadata(self, metadata: PropertyList) -> None:
380 # Docstring inherited.
381 metadata.set("ST_TYPE", "POLYGON", "Type of stencil used to create this cutout.")
382 if len(self._vertices) > 100:
383 raise NotImplementedError(
384 "TODO: FITS limitations make it difficult to serialize big stencils to the header."
385 )
386 for n, v in enumerate(self._vertices):
387 metadata.set(f"ST_RA{n:02d}", v.getRa().asDegrees(), f"Vertex {n} right ascension in degrees.")
388 metadata.set(f"ST_DEC{n:02d}", v.getDec().asDegrees(), f"Vertex {n} declination in degrees.")
390 @property
391 def fingerprint(self) -> bytes:
392 # Docstring inherited.
393 hasher = blake2b(digest_size=16)
394 hasher.update(b"POLYGON")
395 for v in self._vertices:
396 hasher.update(struct.pack("!d!d", v.getRa().asRadians(), v.getDecx().asRadians()))
397 return hasher.digest()
400def _angle_from_astropy(angle: astropy.coordinates.Angle) -> Angle:
401 """Convert an `astropy.coordinates.Angle` to a `lsst.geom.Angle`.
403 Parameters
404 ----------
405 angle : `astropy.coordinates.Angle`
406 Astropy Angle to convert. Must be a scalar.
408 Returns
409 -------
410 angle : `lsst.geom.Angle`
411 LSST angle.
412 """
413 if not angle.isscalar:
414 raise ValueError("Only scalar angles are supported.")
415 return Angle(angle.rad * radians)
418def _spherepoint_from_astropy(skycoord: astropy.coordinates.SkyCoord) -> SpherePoint:
419 """Convert an `astropy.coordinates.SkyCoord` to a `lsst.geom.SpherePoint`.
421 Parameters
422 ----------
423 skycoord : `astropy.coordinates.SkyCoord`
424 Astropy coordinates to convert. Must be a scalar.
426 Returns
427 -------
428 angle : `lsst.geom.SpherePoint`
429 LSST spherical point.
430 """
431 if not skycoord.isscalar:
432 raise ValueError("Only scalar coordinates are supported.")
433 icrs = skycoord.transform_to("icrs")
434 return SpherePoint(icrs.ra.rad * radians, icrs.dec.rad * radians)
437def _make_local_gnomonic_wcs(position: SpherePoint) -> SkyWcs:
438 """Construct a WCS that represents a local gnomonic transform at a point.
440 Parameters
441 ----------
442 position : `SpherePoint`
443 Sky coordinate of the center of the projection.
445 Returns
446 -------
447 wcs : `SkyWcs`
448 A gnomonic (TAN) WCS with ``CRVAL=position``, ``CRPIX=(0,0)``, a pixel
449 scale of ``1/rad``, aligned such that ``(x, y)`` correspond to
450 ``(ra, dec)``.
451 """
452 return makeSkyWcs(Point2D(0.0, 0.0), position, makeCdMatrix(1.0 * radians))