Coverage for python/lsst/images/fields/_chebyshev.py: 81%
223 statements
« prev ^ index » next coverage.py v7.14.3, created at 2026-07-01 09:12 +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__ = ("ChebyshevField", "ChebyshevFieldSerializationModel")
16from collections.abc import Iterator
17from typing import TYPE_CHECKING, Any, ClassVar, Literal, final
19import astropy.units
20import numpy as np
21import pydantic
23from .._concrete_bounds import BoundsSerializationModel
24from .._geom import YX, Bounds, Box
25from .._image import Image
26from ..serialization import ArchiveTree, InlineArray, InputArchive, InvalidParameterError, OutputArchive, Unit
27from ._base import BaseField
29if TYPE_CHECKING:
30 try:
31 from lsst.afw.math import BackgroundMI as LegacyBackground
32 from lsst.afw.math import Chebyshev1Function2D as LegacyChebyshev1Function2D
33 from lsst.afw.math import ChebyshevBoundedField as LegacyChebyshevBoundedField
34 except ImportError:
35 type LegacyBackground = Any # type: ignore[no-redef]
36 type LegacyChebyshevBoundedField = Any # type: ignore[no-redef]
37 type LegacyChebyshev1Function2D = Any # type: ignore[no-redef]
40@final
41class ChebyshevField(BaseField):
42 """A 2-d Chebyshev polynomial over a rectangular region.
44 Parameters
45 ----------
46 bounds
47 The region where this field can be evaluated. The ``bbox`` of this
48 region is grown by half a pixel on all sides and then used to remap
49 coordinates to ``[-1, 1]x[-1, 1]``, which is the natural domain of a
50 2-d Chebyshev polynomial.
51 coefficients
52 Coefficients for the 2-d Chebyshev polynomial of the first kind, as a
53 2-d matrix in which element ``[p, q]`` corresponds to the coefficient
54 of ``T_p(y) T_q(x)``. Will be set to read-only in place.
55 unit
56 Units of the field.
57 """
59 def __init__(
60 self, bounds: Bounds, coefficients: np.ndarray, *, unit: astropy.units.UnitBase | None = None
61 ) -> None:
62 self._bounds = bounds
63 self._coefficients = coefficients
64 self._coefficients.flags.writeable = False
65 self._unit = unit
66 # Compute the scaling and translation that map points in the bbox
67 # (including an extra 0.5 on all sides, since the bbox is int-based)
68 # to [-1, 1].
69 bbox = bounds.bbox
70 self._xs = 2.0 / bbox.x.size
71 self._xt = bbox.x.min + 0.5 * bbox.x.size - 0.5
72 self._ys = 2.0 / bbox.y.size
73 self._yt = bbox.y.min + 0.5 * bbox.y.size - 0.5
75 def __eq__(self, other: object) -> bool:
76 if type(other) is not ChebyshevField: 76 ↛ 77line 76 didn't jump to line 77 because the condition on line 76 was never true
77 return NotImplemented
78 return (
79 self._bounds == other._bounds
80 and self._unit == other._unit
81 and np.array_equal(self._coefficients, other._coefficients, equal_nan=True)
82 )
84 __hash__ = None # type: ignore[assignment]
86 @staticmethod
87 def fit(
88 bounds: Bounds,
89 data: np.ndarray | astropy.units.Quantity,
90 order: int | None = None,
91 *,
92 y: np.ndarray,
93 x: np.ndarray,
94 weight: np.ndarray | None = None,
95 y_order: int | None = None,
96 x_order: int | None = None,
97 triangular: bool = True,
98 unit: astropy.units.UnitBase | None = None,
99 ) -> ChebyshevField:
100 """Fit a Chebyshev field to data points using linear least squares.
102 Parameters
103 ----------
104 bounds
105 Bounding box over which the Chebyshev field is defined.
106 data
107 Data points to fit. If this is an `astropy.units.Quantity`,
108 this sets the units of the field and the ``unit`` argument cannot
109 also be provided.
110 order
111 Maximum order for the Chebyshev polynomial in both dimensions.
112 y
113 Y coordinates of the data points. Must have either the same
114 shape as ``data`` (providing the coordinates for all points
115 directly), or be a 1-d array with the same size as
116 ``data.shape[0]`` (when ``data`` is a 2-d image and ``y`` provides
117 the coordinates of the rows).
118 x
119 X coordinates of the data points. Must have either the same
120 shape as ``data`` (providing the coordinates for all points
121 directly), or be a 1-d array with the same size as
122 ``data.shape[1]`` (when ``data`` is a 2-d image and ``x`` provides
123 the coordinates of the columns).
124 weight
125 Weights to apply to the data points. Must have the same shape as
126 ``data``.
127 y_order
128 Maximum order for the Chebyshev polynomial in ``y``. Requires
129 ``x_order`` to also be provided. Incompatible with ``order``.
130 x_order
131 Maximum order for the Chebyshev polynomial in ``x``. Requires
132 ``y_order`` to also be provided. Incompatible with ``order``.
133 triangular
134 If `True`, only fit for coefficients of ``T_p(y) T_q(x)`` where
135 ``p + q <= max(y_order, x_order)``.
136 unit
137 Units of the returned field.
138 """
139 match (order, x_order, y_order):
140 case (int(), None, None):
141 x_order = order
142 y_order = order
143 case (None, int(), int()): 143 ↛ 145line 143 didn't jump to line 145 because the pattern on line 143 always matched
144 pass
145 case _:
146 raise TypeError("Either 'order' (only) or both 'x_order' and 'y_order' must be provided.")
147 if weight is not None and weight.shape != data.shape: 147 ↛ 148line 147 didn't jump to line 148 because the condition on line 147 was never true
148 raise ValueError(f"Shape of 'data' {data.shape} does not match 'weight' {weight.shape}.")
149 if isinstance(data, astropy.units.Quantity):
150 if unit is not None: 150 ↛ 151line 150 didn't jump to line 151 because the condition on line 150 was never true
151 raise TypeError("If 'data' is a Quantity, 'unit' cannot be provided separately.")
152 unit = data.unit
153 data = data.to_value()
154 result = ChebyshevField(bounds, np.zeros((y_order + 1, x_order + 1), dtype=np.float64), unit=unit)
155 if len(data.shape) == 2 and len(x.shape) == 1 and len(y.shape) == 1:
156 if data.shape != y.shape + x.shape: 156 ↛ 157line 156 didn't jump to line 157 because the condition on line 156 was never true
157 raise ValueError(
158 f"Shape of 2-d 'data' {data.shape} does not match 1-d 'y' {y.shape} and/or 'x' {x.shape}."
159 )
160 matrix = result._make_grid_matrix(x=x, y=y, triangular=triangular)
161 else:
162 if data.shape != y.shape: 162 ↛ 163line 162 didn't jump to line 163 because the condition on line 162 was never true
163 raise ValueError(f"Shape of 'data' {data.shape} does not match 'y' {y.shape}.")
164 if data.shape != x.shape: 164 ↛ 165line 164 didn't jump to line 165 because the condition on line 164 was never true
165 raise ValueError(f"Shape of 'data' {data.shape} does not match 'x' {x.shape}.")
166 matrix = result._make_general_matrix(x=x, y=y, triangular=triangular)
167 if weight is not None:
168 weight = weight.ravel() # copies only if needed
169 matrix *= weight[:, np.newaxis]
170 data = data.flatten() # always copies
171 data *= weight
172 mask = np.logical_and(weight > 0, np.isfinite(data))
173 else:
174 data = data.ravel()
175 mask = np.isfinite(data)
176 n_good = mask.sum()
177 if n_good == 0: 177 ↛ 178line 177 didn't jump to line 178 because the condition on line 177 was never true
178 raise ValueError("No good data points.")
179 if n_good < data.size:
180 data = data[mask]
181 matrix = matrix[mask, :]
182 packed_coefficients, *_ = np.linalg.lstsq(matrix, data)
183 result._coefficients.flags.writeable = True
184 for i, pq in result._packing_indices(triangular):
185 result._coefficients[pq.y, pq.x] = packed_coefficients[i]
186 result._coefficients.flags.writeable = False
187 return result
189 @property
190 def bounds(self) -> Bounds:
191 return self._bounds
193 @property
194 def unit(self) -> astropy.units.UnitBase | None:
195 return self._unit
197 @property
198 def x_order(self) -> int:
199 """Maximum polynomial order in the column dimension (`int`)."""
200 return self._coefficients.shape[1] - 1
202 @property
203 def y_order(self) -> int:
204 """Maximum polynomial order in the row dimension (`int`)."""
205 return self._coefficients.shape[0] - 1
207 @property
208 def order(self) -> int:
209 """Maximum polynomial order in either dimension (`int`)."""
210 return max(self.x_order, self.y_order)
212 @property
213 def coefficients(self) -> np.ndarray:
214 """Coefficients for the 2-d Chebyshev polynomial of the first kind,
215 as a 2-d matrix in which element ``[p, q]`` corresponds to the
216 coefficient of ``T_p(y) T_q(x)``.
217 """
218 return self._coefficients
220 @property
221 def is_constant(self) -> bool:
222 return self.x_order == 0 and self.y_order == 0
224 def evaluate(
225 self, *, x: np.ndarray, y: np.ndarray, quantity: bool
226 ) -> np.ndarray | astropy.units.Quantity:
227 m = self._remap(x=x.copy(), y=y.copy())
228 # We swap x and y relative to Numpy's docs because that's how our
229 # coefficients are ordered.
230 v = np.polynomial.chebyshev.chebval2d(m.y, m.x, self._coefficients)
231 if quantity:
232 return astropy.units.Quantity(v, self.unit)
233 return v
235 def render(self, bbox: Box | None = None, *, dtype: np.typing.DTypeLike | None = None) -> Image:
236 if bbox is None:
237 bbox = self.bounds.bbox
238 m = self._remap(
239 x=bbox.x.arange.astype(np.float64),
240 y=bbox.y.arange.astype(np.float64),
241 )
242 # We swap x and y relative to Numpy's docs because that's how our
243 # coefficients and images are ordered.
244 v = np.polynomial.chebyshev.chebgrid2d(m.y, m.x, self._coefficients)
245 return Image(v, bbox=bbox, unit=self.unit, dtype=dtype)
247 def multiply_constant(
248 self, factor: float | astropy.units.Quantity | astropy.units.UnitBase
249 ) -> ChebyshevField:
250 factor, unit = self._handle_factor_units(factor)
251 return ChebyshevField(self.bounds, self.coefficients * factor, unit=unit)
253 def serialize(self, archive: OutputArchive[Any]) -> ChebyshevFieldSerializationModel:
254 """Serialize the Chebyshev field to an output archive.
256 Parameters
257 ----------
258 archive
259 Archive to write to.
260 """
261 return ChebyshevFieldSerializationModel(
262 bounds=self.bounds.serialize(),
263 coefficients=self.coefficients,
264 unit=self.unit,
265 )
267 @staticmethod
268 def _get_archive_tree_type(
269 pointer_type: type[Any],
270 ) -> type[ChebyshevFieldSerializationModel]:
271 """Return the serialization model type for this object for an archive
272 type that uses the given pointer type.
273 """
274 return ChebyshevFieldSerializationModel
276 @staticmethod
277 def from_legacy(
278 legacy: LegacyChebyshevBoundedField,
279 unit: astropy.units.UnitBase | None = None,
280 bounds: Bounds | None = None,
281 ) -> ChebyshevField:
282 """Convert from a legacy `lsst.afw.math.ChebyshevBoundedField`.
284 Parameters
285 ----------
286 legacy
287 Legacy field to convert.
288 unit
289 The units of the returned field (`lsst.afw.math.BoundedField`
290 objects do not know their units).
291 bounds
292 The bounds of the returned field, if they should be different from
293 the bounding box of ``legacy``.
294 """
295 bbox = Box.from_legacy(legacy.getBBox())
296 if bounds is not None: 296 ↛ 297line 296 didn't jump to line 297 because the condition on line 296 was never true
297 if bounds.bbox != bbox:
298 raise ValueError(
299 "Custom bounds when converting a ChebyshevBoundedField must not change the bbox."
300 )
301 else:
302 bounds = bbox
303 return ChebyshevField(bounds=bounds, coefficients=legacy.getCoefficients(), unit=unit)
305 def to_legacy(self) -> LegacyChebyshevBoundedField:
306 """Convert to a legacy `lsst.afw.math.ChebyshevBoundedField`."""
307 from lsst.afw.math import ChebyshevBoundedField as LegacyChebyshevBoundedField
309 return LegacyChebyshevBoundedField(self.bounds.bbox.to_legacy(), self.coefficients)
311 @staticmethod
312 def from_legacy_background(
313 legacy_background: LegacyBackground,
314 bounds: Bounds | None = None,
315 unit: astropy.units.UnitBase | None = None,
316 ) -> ChebyshevField:
317 """Convert from a legacy `lsst.afw.math.BackgroundMI` instance.
319 Parameters
320 ----------
321 legacy_background
322 Legacy background object to convert.
323 bounds
324 The bounds of the returned field, if they should be different from
325 the bounding box of ``legacy_background``.
326 unit
327 The units of the returned field (`lsst.afw.math.Background`
328 objects do not know their units).
329 """
330 from lsst.afw.math import ApproximateControl
332 approx_control = legacy_background.getBackgroundControl().getApproximateControl()
333 stats_image = legacy_background.getStatsImage()
334 if approx_control.getStyle() != ApproximateControl.CHEBYSHEV:
335 raise TypeError("Legacy background does not use Chebyshev approximation.")
336 if approx_control.getWeighting():
337 weight = stats_image.variance.array ** (-0.5)
338 else:
339 weight = None
340 x = legacy_background.getBinCentersX()
341 y = legacy_background.getBinCentersY()
342 bbox = Box.from_legacy(legacy_background.getImageBBox())
343 if bounds is not None:
344 if bounds.bbox != bbox:
345 raise ValueError(
346 "Custom bounds when converting a Chebyshev background must not change the bbox."
347 )
348 else:
349 bounds = bbox
350 return ChebyshevField.fit(
351 bounds,
352 stats_image.image.array,
353 x=x,
354 y=y,
355 x_order=approx_control.getOrderX(),
356 y_order=approx_control.getOrderY(),
357 weight=weight,
358 unit=unit,
359 )
361 @staticmethod
362 def from_legacy_function2(
363 legacy_function2: LegacyChebyshev1Function2D,
364 bounds: Bounds | None = None,
365 unit: astropy.units.Unit | None = None,
366 ) -> ChebyshevField:
367 """Convert from a legacy `lsst.afw.math.Chebyshev1Function2D`.
369 Parameters
370 ----------
371 legacy_function2
372 Legacy function object to convert.
373 bounds
374 The bounds of the returned field, if they should be different from
375 the bounding box of ``legacy_background``.
376 unit
377 The units of the returned field.
378 """
379 xy_range = legacy_function2.getXYRange()
380 bbox = Box.factory[
381 _require_int(xy_range.y.min + 0.5) : _require_int(xy_range.y.max + 0.5),
382 _require_int(xy_range.x.min + 0.5) : _require_int(xy_range.x.max + 0.5),
383 ]
384 if bounds is not None: 384 ↛ 385line 384 didn't jump to line 385 because the condition on line 384 was never true
385 if bounds.bbox != bbox:
386 raise ValueError(
387 "Custom bounds when converting a Chebyshev background must not change the bbox."
388 )
389 else:
390 bounds = bbox
391 order = legacy_function2.getOrder()
392 coefficients = np.zeros((order + 1, order + 1), dtype=np.float64)
393 for i, pq in ChebyshevField._legacy_function2_indices(order):
394 coefficients[pq.y, pq.x] = legacy_function2.getParameter(i)
395 return ChebyshevField(bbox, coefficients, unit=unit)
397 def to_legacy_function2(self) -> LegacyChebyshev1Function2D:
398 """Convert to a legacy `lsst.afw.math.Chebyshev1Function2D`."""
399 from lsst.afw.math import Chebyshev1Function2D as LegacyChebyshev1Function2D
400 from lsst.geom import Box2D as LegacyBox2D
402 order = max(self.y_order, self.x_order)
403 result = LegacyChebyshev1Function2D(order, LegacyBox2D(self.bounds.bbox.to_legacy()))
404 for i, pq in self._legacy_function2_indices(order):
405 result.setParameter(
406 i,
407 (
408 self._coefficients[pq.y, pq.x]
409 if pq.y < self._coefficients.shape[0] and pq.x < self._coefficients.shape[1]
410 else 0.0
411 ),
412 )
413 return result
415 @staticmethod
416 def _legacy_function2_indices(order: int) -> Iterator[tuple[int, YX[int]]]:
417 i = 0
418 for n in range(order + 1):
419 for p in range(0, n + 1):
420 q = n - p
421 yield i, YX(y=p, x=q)
422 i += 1
424 def _remap(self, *, x: np.ndarray, y: np.ndarray) -> YX[np.ndarray]:
425 x -= self._xt
426 x *= self._xs
427 y -= self._yt
428 y *= self._ys
429 return YX(y=y, x=x)
431 def _packing_indices(self, triangular: bool) -> Iterator[tuple[int, YX[int]]]:
432 i = 0
433 for p in range(self.y_order + 1):
434 for q in range(self.x_order + 1):
435 if not triangular or p + q <= self.order:
436 yield i, YX(y=p, x=q)
437 i += 1
439 def _make_grid_matrix(self, *, x: np.ndarray, y: np.ndarray, triangular: bool) -> np.ndarray:
440 r = self._remap(
441 x=np.asarray(x, dtype=np.float64, copy=True),
442 y=np.asarray(y, dtype=np.float64, copy=True),
443 )
444 yv = np.polynomial.chebyshev.chebvander(r.y, self.y_order)
445 xv = np.polynomial.chebyshev.chebvander(r.x, self.x_order)
446 indices = list(self._packing_indices(triangular))
447 tensor = np.zeros(r.y.shape + r.x.shape + (len(indices),), dtype=np.float64)
448 for i, pq in indices:
449 tensor[:, :, i] = np.multiply.outer(yv[:, pq.y], xv[:, pq.x])
450 return tensor.reshape(y.shape[0] * x.shape[0], len(indices))
452 def _make_general_matrix(self, *, x: np.ndarray, y: np.ndarray, triangular: bool) -> np.ndarray:
453 r = self._remap(
454 x=np.asarray(x, dtype=np.float64, copy=True).ravel(),
455 y=np.asarray(y, dtype=np.float64, copy=True).ravel(),
456 )
457 yv = np.polynomial.chebyshev.chebvander(r.y, self.y_order)
458 xv = np.polynomial.chebyshev.chebvander(r.x, self.x_order)
459 indices = list(self._packing_indices(triangular))
460 matrix = np.zeros(r.y.shape + (len(indices),), dtype=np.float64)
461 for i, pq in indices:
462 matrix[:, i] = yv[:, pq.y] * xv[:, pq.x]
463 return matrix
466class ChebyshevFieldSerializationModel(ArchiveTree):
467 """Serialization model for `ChebyshevField`."""
469 SCHEMA_NAME: ClassVar[str] = "chebyshev_field"
470 SCHEMA_VERSION: ClassVar[str] = "1.0.0"
471 MIN_READ_VERSION: ClassVar[int] = 1
472 PUBLIC_TYPE: ClassVar[type] = ChebyshevField
474 bounds: BoundsSerializationModel = pydantic.Field(
475 description=(
476 "The region where this field can be evaluated. "
477 "The bbox of this region is grown by half a pixel on all sides and then used to remap "
478 "coordinates to [-1, 1]x[-1, 1], which is the natural domain of a 2-d Chebyshev polynomial."
479 )
480 )
482 coefficients: InlineArray = pydantic.Field(
483 description=(
484 "Coefficients for a 2-d Chebyshev polynomial of the first kind, as a 2-d matrix in which "
485 "element [p, q] corresponds to the coefficient of T_p(y) T_q(x)."
486 )
487 )
489 unit: Unit | None = pydantic.Field(default=None, description="Units of the field.")
491 field_type: Literal["CHEBYSHEV"] = "CHEBYSHEV"
493 def deserialize(self, archive: InputArchive, **kwargs: Any) -> ChebyshevField:
494 """Deserialize the Chebyshev field from an input archive.
496 Parameters
497 ----------
498 archive
499 Archive to read from.
500 **kwargs
501 Unsupported keyword arguments are accepted only to provide
502 better error messages (raising
503 `.serialization.InvalidParameterError`).
504 """
505 if kwargs: 505 ↛ 506line 505 didn't jump to line 506 because the condition on line 505 was never true
506 raise InvalidParameterError(f"Unrecognized parameters for ChebyshevField: {set(kwargs.keys())}.")
507 return ChebyshevField(self.bounds.deserialize(), self.coefficients, unit=self.unit)
510def _require_int(v: float) -> int:
511 if (z := int(v)) == v: 511 ↛ 513line 511 didn't jump to line 513 because the condition on line 511 was always true
512 return z
513 raise ValueError("Legacy Chebyshev1Function2 XY range must be at half-integer positions.")