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