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