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