Coverage for python/lsst/images/fields/_spline.py: 26%
150 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-06 08:34 +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
14__all__ = ("SplineField", "SplineFieldSerializationModel")
16from typing import TYPE_CHECKING, Any, ClassVar, Literal, final
18import astropy.units
19import numpy as np
20import pydantic
21from scipy.interpolate import Akima1DInterpolator
23from .._concrete_bounds import SerializableBounds
24from .._geom import Bounds, Box
25from .._image import Image
26from ..serialization import (
27 ArchiveTree,
28 ArrayReferenceModel,
29 InlineArray,
30 InputArchive,
31 InvalidParameterError,
32 OutputArchive,
33 Unit,
34)
35from ._base import BaseField
37if TYPE_CHECKING:
38 try:
39 from lsst.afw.math import BackgroundMI as LegacyBackground
40 except ImportError:
41 type LegacyBackground = Any # type: ignore[no-redef]
44@final
45class SplineField(BaseField):
46 """A 2-d Akima spline interpolation of data on a regular grid.
48 Parameters
49 ----------
50 bounds
51 The region where this field can be evaluated.
52 data
53 The data points to be interpolated. Missing values (indicated by NaN)
54 are allowed. Will be set to read-only in place.
55 y
56 Coordinates for the first dimension of ``data``. Will be set to
57 read-only in place.
58 x
59 Coordinates for the second dimension of ``data``. Will be set to
60 read-only in place.
61 unit
62 Units of the field.
64 Notes
65 -----
66 This field is much faster to evaluate on a grid via `render` than at
67 arbitrary points via the function-call operator.
68 """
70 def __init__(
71 self,
72 bounds: Bounds,
73 data: np.ndarray,
74 *,
75 y: np.ndarray,
76 x: np.ndarray,
77 unit: astropy.units.UnitBase | None = None,
78 ):
79 if isinstance(data, astropy.units.Quantity):
80 if unit is not None:
81 raise TypeError("If 'data' is a Quantity, 'unit' cannot be provided separately.")
82 unit = data.unit
83 data = data.to_value()
84 if data.ndim != 2:
85 raise ValueError("'data' must be 2-d.")
86 if y.ndim != 1:
87 raise ValueError("'y' must be 1-d.")
88 if x.ndim != 1:
89 raise ValueError("'x' must be 1-d.")
90 if data.shape != y.shape + x.shape:
91 raise ValueError(
92 f"Shape of 2-d 'data' {data.shape} does not match "
93 f"expected 1-d 'y' {y.shape} and/or 'x' {x.shape}."
94 )
95 self._bounds = bounds
96 self._data = data
97 self._data.flags.writeable = False
98 self._x = x
99 self._x.flags.writeable = False
100 self._y = y
101 self._y.flags.writeable = False
102 self._unit = unit
104 def __eq__(self, other: object) -> bool:
105 if type(other) is not SplineField:
106 return NotImplemented
107 return (
108 self._bounds == other._bounds
109 and self._unit == other._unit
110 and np.array_equal(self._data, other._data, equal_nan=True)
111 and np.array_equal(self._x, other._x, equal_nan=True)
112 and np.array_equal(self._y, other._y, equal_nan=True)
113 )
115 __hash__ = None # type: ignore[assignment]
117 @property
118 def bounds(self) -> Bounds:
119 return self._bounds
121 @property
122 def unit(self) -> astropy.units.UnitBase | None:
123 return self._unit
125 @property
126 def data(self) -> np.ndarray:
127 """The data points to be interpolated (`numpy.ndarray`).
129 May have missing values indicated by NaNs.
130 """
131 return self._data
133 @property
134 def x(self) -> np.ndarray:
135 """Coordinates for the second dimension of `data` (`numpy.ndarray`)."""
136 return self._x
138 @property
139 def y(self) -> np.ndarray:
140 """Coordinates for the first dimension of `data` (`numpy.ndarray`)."""
141 return self._y
143 @property
144 def is_constant(self) -> bool:
145 # We really do want an exact floating-point comparison here.
146 return (self._data == self._data[0, 0]).all()
148 def evaluate(
149 self, *, x: np.ndarray, y: np.ndarray, quantity: bool = False
150 ) -> np.ndarray | astropy.units.Quantity:
151 y, x = np.broadcast_arrays(y, x)
152 xg = self._x
153 y_render = np.zeros(xg.shape + y.shape, dtype=np.float64)
154 mask = np.zeros(xg.size, dtype=bool)
155 for j in range(xg.size):
156 if (y_interpolator := self._make_y_interpolator(j)) is not None:
157 y_render[j, ...] = y_interpolator(y)
158 mask[j] = True
159 if not np.all(mask):
160 y_render = y_render[mask, ...]
161 xg = xg[mask]
162 result = np.zeros(y.shape, dtype=np.float64)
163 # There doesn't seem to be a way to avoid looping in Python here;
164 # maybe someday we'll push this down to a compiled language.
165 for i, xv in np.ndenumerate(x):
166 if (x_interpolator := self._make_1d_interpolator(xg, y_render[:, *i])) is None:
167 raise ValueError("No valid data points.")
168 v = x_interpolator(xv)
169 result[*i] = v
170 if quantity:
171 return astropy.units.Quantity(result, self._unit)
172 return result
174 def render(self, bbox: Box | None = None, *, dtype: np.typing.DTypeLike | None = None) -> Image:
175 if bbox is None:
176 bbox = self.bounds.bbox
177 xg = self._x
178 y_render = np.zeros((xg.size, bbox.y.size), dtype=dtype)
179 mask = np.zeros(xg.size, dtype=bool)
180 for j in range(xg.size): # we have to loop, but only over bins, not evaluation points.
181 if (y_interpolator := self._make_y_interpolator(j)) is not None:
182 y_render[j, :] = y_interpolator(bbox.y.arange)
183 mask[j] = True
184 if not np.all(mask):
185 y_render = y_render[mask, :]
186 xg = xg[mask]
187 if (x_interpolator := self._make_1d_interpolator(xg, y_render)) is None:
188 raise ValueError("No valid data points.")
189 rendered_array = x_interpolator(bbox.x.arange)
190 return Image(rendered_array.transpose().copy(), bbox=bbox, unit=self._unit, dtype=dtype)
192 def multiply_constant(
193 self, factor: float | astropy.units.Quantity | astropy.units.UnitBase
194 ) -> SplineField:
195 factor, unit = self._handle_factor_units(factor)
196 return SplineField(self._bounds, self._data * factor, y=self._y, x=self._x, unit=unit)
198 def serialize(self, archive: OutputArchive[Any]) -> SplineFieldSerializationModel:
199 """Serialize the spline field to an output archive."""
200 return SplineFieldSerializationModel(
201 bounds=self.bounds.serialize(),
202 data=archive.add_array(self._data, name="data"),
203 y=self._y,
204 x=self._x,
205 unit=self._unit,
206 )
208 @staticmethod
209 def _get_archive_tree_type(
210 pointer_type: type[Any],
211 ) -> type[SplineFieldSerializationModel]:
212 """Return the serialization model type for this object for an archive
213 type that uses the given pointer type.
214 """
215 return SplineFieldSerializationModel
217 @staticmethod
218 def from_legacy_background(
219 legacy_background: LegacyBackground,
220 bounds: Bounds | None = None,
221 unit: astropy.units.UnitBase | None = None,
222 ) -> SplineField:
223 """Convert from a legacy `lsst.afw.math.BackgroundMI` instance.
225 Parameters
226 ----------
227 legacy
228 Legacy background object to convert.
229 bounds
230 The bounds of the returned field, if they should be different from
231 the bounding box of ``legacy_background``.
232 unit
233 The units of the returned field (`lsst.afw.math.Background`
234 objects do not know their units).
236 Notes
237 -----
238 `SplineField.render` and the `lsst.afw` background interpolator both
239 use Akima splines, but with slightly different boundary conditions.
240 They should produce equivalent to single-precision round-off error
241 when evaluated within the region enclosed by bin centers (i.e. where
242 no extrapolation is necessary) and when there are five or more
243 points to be interpolated in each row and column.
244 """
245 from lsst.afw.math import ApproximateControl, Interpolate
247 bg_control = legacy_background.getBackgroundControl()
248 approx_control = bg_control.getApproximateControl()
249 stats_image = legacy_background.getStatsImage()
250 # In the afw background system, "approximate" is the opposite of
251 # "interpolate", but it also implied Chebyshev since that's the only
252 # approximation algorithm we every implemented. All of the
253 # interpolation options are similarly splines, and non-Akima splines
254 # are *mostly* only used when there aren't enough control points for
255 # Akima splines. Since SciPy automatically falls back to non-Akima
256 # splines in those cases (or maybe they're formally a limit of Akima
257 # splines, I don't know), we just always assume what we get can be
258 # Akima-spline interpolated by SciPy to good enough approximation with
259 # what afw would do.
260 if approx_control.getStyle() != ApproximateControl.UNKNOWN:
261 raise TypeError("Legacy background uses Chebyshev approximation, not splines.")
262 if bg_control.getInterpStyle() == Interpolate.UNKNOWN:
263 raise TypeError("Legacy background does not use spline interpolation.")
264 x = legacy_background.getBinCentersX()
265 y = legacy_background.getBinCentersY()
266 return SplineField(
267 Box.from_legacy(legacy_background.getImageBBox()) if bounds is None else bounds,
268 stats_image.image.array,
269 x=x,
270 y=y,
271 unit=unit,
272 )
274 def _make_1d_interpolator(self, loc: np.ndarray, val: np.ndarray) -> Akima1DInterpolator | None:
275 match len(loc):
276 case 0:
277 return None
278 case 1:
279 # SciPy can handle only two points by downgrading to linear
280 # interpolation, but it raises if given only one. Mock up
281 # two for the nearest-neighbor fallback.
282 return Akima1DInterpolator(np.array([loc[0], loc[0]]), np.array([val[0], val[0]]))
283 case _:
284 return Akima1DInterpolator(loc, val, extrapolate=True)
286 def _make_y_interpolator(self, j: int) -> Akima1DInterpolator | None:
287 y = self._y
288 z = self._data[:, j]
289 mask = np.isfinite(z)
290 if not np.all(mask):
291 y = y[mask]
292 z = z[mask]
293 del mask
294 return self._make_1d_interpolator(y, z)
297class SplineFieldSerializationModel(ArchiveTree):
298 """Serialization model for `SplineField`."""
300 SCHEMA_NAME: ClassVar[str] = "spline_field"
301 SCHEMA_VERSION: ClassVar[str] = "1.0.0"
302 MIN_READ_VERSION: ClassVar[int] = 1
303 PUBLIC_TYPE: ClassVar[type] = SplineField
305 bounds: SerializableBounds = pydantic.Field(description=("The region where this field can be evaluated."))
307 data: ArrayReferenceModel = pydantic.Field(
308 description="2-d data to interpolate. NaNs indicate missing values."
309 )
311 y: InlineArray = pydantic.Field(description="Row positions of the data points.")
313 x: InlineArray = pydantic.Field(description="Column positions of the data points.")
315 unit: Unit | None = pydantic.Field(default=None, description="Units of the field.")
317 field_type: Literal["SPLINE"] = "SPLINE"
319 def deserialize(self, archive: InputArchive, **kwargs: Any) -> SplineField:
320 """Deserialize the spline field from an input archive."""
321 if kwargs:
322 raise InvalidParameterError(f"Unrecognized parameters for SplineField: {set(kwargs.keys())}.")
323 return SplineField(
324 self.bounds.deserialize(),
325 archive.get_array(self.data),
326 y=self.y,
327 x=self.x,
328 unit=self.unit,
329 )