Coverage for tests / test_fields.py: 23%
190 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-25 08:35 +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
14import os
15import unittest
17import astropy.units as u
18import numpy as np
20from lsst.images import Box, Image
21from lsst.images.fields import (
22 BaseField,
23 ChebyshevField,
24 ProductField,
25 SplineField,
26 SumField,
27 field_from_legacy,
28 field_from_legacy_background,
29)
30from lsst.images.tests import (
31 RoundtripFits,
32 assert_images_equal,
33 compare_field_to_legacy,
34)
36try:
37 from lsst.afw.math import BackgroundList as LegacyBackgroundList
38 from lsst.afw.math import ChebyshevBoundedField as LegacyChebyshevBoundedField
39 from lsst.afw.math import ProductBoundedField as LegacyProductBoundedField
40except ImportError:
41 HAVE_LEGACY = False
42else:
43 HAVE_LEGACY = True
45DATA_DIR = os.environ.get("TESTDATA_IMAGES_DIR", None)
48class FieldTestCase(unittest.TestCase):
49 """Tests for the Field classes that do not require legacy code."""
51 def setUp(self) -> None:
52 self.rng = np.random.default_rng(10)
53 self.box = Box.factory[6:32, -7:26]
54 self.cheby = ChebyshevField(self.box, np.array([[0.5, -0.25], [0.40, 0.0]]))
55 self.spline = SplineField(
56 self.box,
57 self.rng.standard_normal(size=(6, 7)),
58 y=self.box.y.linspace(6),
59 x=self.box.x.linspace(7),
60 )
61 self.product = ProductField([self.spline, self.cheby])
62 self.sum = SumField([self.spline, self.cheby])
64 def test_chebyshev_call_limits(self) -> None:
65 """Test that ChebyshevField.__call__ evaluates correctly at low order
66 at the corners of its box.
68 Testing at first order is enough to catch (y, x) swap errors, and is
69 easy because ``T_0(x) = 1`` and ``T_1(x) = x``. Higher-order
70 evaluation is covered in the legacy comparison tests.
71 """
72 result = self.cheby(x=np.array([-7.5, 25.5, 25.5, -7.5]), y=np.array([5.5, 5.5, 31.5, 31.5]))
73 self.assertEqual(result[0], 0.5 + 0.25 - 0.4) # [x=-1, y=-1] after remap
74 self.assertEqual(result[1], 0.5 - 0.25 - 0.4) # [x=1, y=-1] after remap
75 self.assertEqual(result[2], 0.5 - 0.25 + 0.4) # [x=1, y=1] after remap
76 self.assertEqual(result[3], 0.5 + 0.25 + 0.4) # [x=-1, y=1] after remap
78 def test_chebyshev_attributes(self) -> None:
79 """Test the basic properties of a Chebyshev field."""
80 self.assertEqual(self.cheby.bounds, self.box)
81 self.assertIsNone(self.cheby.unit)
82 self.assertEqual(self.cheby.x_order, 1)
83 self.assertEqual(self.cheby.y_order, 1)
84 self.assertEqual(self.cheby.order, 1)
85 np.testing.assert_array_equal(self.cheby.coefficients, np.array([[0.5, -0.25], [0.40, 0.0]]))
87 def test_chebyshev_fit(self) -> None:
88 """Test that we can fit a ChebyshevField to gridded data that should
89 have zero residuals.
90 """
91 data_image = self.cheby.render()
92 cheby2 = ChebyshevField.fit(
93 self.box, data_image.array, order=1, y=self.box.y.arange, x=self.box.x.arange
94 )
95 self.assertEqual(cheby2.bounds, self.box)
96 np.testing.assert_array_almost_equal(cheby2.coefficients, self.cheby.coefficients)
97 # Repeat the experiment in a few different ways that don't affect the
98 # result (at least in this scenario where we have an exact fit).
99 # Fit to order 2 in x (will get us extra zero-valued coefficients):
100 cheby3 = ChebyshevField.fit(
101 self.box, data_image.array, x_order=2, y_order=1, y=self.box.y.arange, x=self.box.x.arange
102 )
103 self.assertEqual(cheby3.bounds, self.box)
104 np.testing.assert_array_almost_equal(
105 cheby3.coefficients,
106 np.array([[0.5, -0.25, 0.0], [0.40, 0.0, 0.0]], dtype=np.float64),
107 )
108 # Fit with triangular=False (would allow the (1, 1) term to be nonzero,
109 # but it will still be zero here):
110 cheby4 = ChebyshevField.fit(
111 self.box, data_image.array, order=1, y=self.box.y.arange, x=self.box.x.arange, triangular=False
112 )
113 self.assertEqual(cheby4.bounds, self.box)
114 np.testing.assert_array_almost_equal(cheby4.coefficients, self.cheby.coefficients)
115 # Fit with weights.
116 cheby5 = ChebyshevField.fit(
117 self.box,
118 data_image.array,
119 order=1,
120 y=self.box.y.arange,
121 x=self.box.x.arange,
122 weight=self.rng.uniform(0.8, 1.2, size=data_image.array.shape),
123 )
124 self.assertEqual(cheby5.bounds, self.box)
125 np.testing.assert_array_almost_equal(cheby5.coefficients, self.cheby.coefficients)
126 # Fit to a Quantity.
127 cheby6 = ChebyshevField.fit(
128 self.box, data_image.array * u.nJy, order=1, y=self.box.y.arange, x=self.box.x.arange
129 )
130 self.assertEqual(cheby6.bounds, self.box)
131 self.assertEqual(cheby6.unit, u.nJy)
132 np.testing.assert_array_almost_equal(cheby6.coefficients, self.cheby.coefficients)
133 # Fit with units provided separately.
134 cheby7 = ChebyshevField.fit(
135 self.box, data_image.array, order=1, y=self.box.y.arange, x=self.box.x.arange, unit=u.nJy
136 )
137 self.assertEqual(cheby7.bounds, self.box)
138 self.assertEqual(cheby7.unit, u.nJy)
139 np.testing.assert_array_almost_equal(cheby7.coefficients, self.cheby.coefficients)
140 # Fit with x and y labeling every data point, not the grid.
141 m = self.box.meshgrid()
142 cheby8 = ChebyshevField.fit(self.box, data_image.array, order=1, y=m.y, x=m.x)
143 self.assertEqual(cheby8.bounds, self.box)
144 np.testing.assert_array_almost_equal(cheby8.coefficients, self.cheby.coefficients)
145 # Fit with x and y labeling every data point, not the grid, as well as
146 # weights.
147 cheby9 = ChebyshevField.fit(
148 self.box,
149 data_image.array,
150 order=1,
151 y=m.y,
152 x=m.x,
153 weight=self.rng.uniform(0.8, 1.2, size=data_image.array.shape),
154 )
155 self.assertEqual(cheby9.bounds, self.box)
156 np.testing.assert_array_almost_equal(cheby9.coefficients, self.cheby.coefficients)
157 # Fit with one data points replaced by NaN, which should be ignored
158 # (the fit is still overconstrained).
159 new_data = data_image.array.copy()
160 new_data[5, 7] = np.nan
161 cheby10 = ChebyshevField.fit(self.box, new_data, order=1, y=self.box.y.arange, x=self.box.x.arange)
162 self.assertEqual(cheby10.bounds, self.box)
163 np.testing.assert_array_almost_equal(cheby10.coefficients, self.cheby.coefficients)
165 def test_chebyshev_evaluation_consistency(self) -> None:
166 self.check_evaluation_consistency(self.cheby)
168 def test_chebyshev_units(self) -> None:
169 self.check_units(self.cheby)
171 def test_spline_knot_evaluation(self) -> None:
172 """Test that SplineField.__call__ evaluates to the input data points
173 at their positions.
174 """
175 xv, yv = np.meshgrid(self.spline.x, self.spline.y)
176 result = self.spline(x=xv, y=yv)
177 np.testing.assert_array_almost_equal(result, self.spline.data)
179 def test_spline_evaluation_consistency(self) -> None:
180 self.check_evaluation_consistency(self.spline)
182 def test_spline_units(self) -> None:
183 self.check_units(self.spline)
185 def test_product_evaluation(self) -> None:
186 """Test ProductField.__call__ against direct calls to its operands."""
187 xv, yv = self.box.meshgrid(n=3)
188 z = self.product(x=xv, y=yv)
189 np.testing.assert_array_equal(z, self.cheby(x=xv, y=yv) * self.spline(x=xv, y=yv))
191 def test_product_evaluation_consistency(self) -> None:
192 self.check_evaluation_consistency(self.product)
194 def test_product_units(self) -> None:
195 self.check_units(self.product)
196 self.assertEqual(ProductField([self.cheby, self.spline * u.nJy]).unit, u.nJy)
197 self.assertEqual(ProductField([self.cheby * u.nJy, self.spline]).unit, u.nJy)
198 self.assertEqual(
199 ProductField([self.cheby * u.nJy, self.spline / u.arcsec**2]).unit, u.nJy / u.arcsec**2
200 )
202 def test_sum_evaluation(self) -> None:
203 """Test sumField.__call__ against direct calls to its operands."""
204 xv, yv = self.box.meshgrid(n=3)
205 z = self.sum(x=xv, y=yv)
206 np.testing.assert_array_equal(z, self.cheby(x=xv, y=yv) + self.spline(x=xv, y=yv))
208 def test_sum_evaluation_consistency(self) -> None:
209 self.check_evaluation_consistency(self.sum)
211 def test_sum_units(self) -> None:
212 self.check_units(self.sum)
213 with self.assertRaises(u.UnitConversionError):
214 SumField([self.cheby, self.spline * u.nJy])
215 with self.assertRaises(u.UnitConversionError):
216 SumField([self.cheby * u.nJy, self.spline])
217 # Test a SumField where the operands have different but compatible
218 # units.
219 mixed = SumField([self.cheby * u.rad, self.spline * u.deg])
220 self.assertEqual(mixed.unit, u.rad)
221 small_box = Box.factory[10:12, 2:5]
222 cheby_render = self.cheby.render(small_box)
223 spline_render = self.spline.render(small_box)
224 mixed_render = mixed.render(small_box)
225 self.assertEqual(mixed_render.unit, u.rad)
226 np.testing.assert_array_almost_equal(
227 mixed_render.array, cheby_render.array + (spline_render.array * np.pi / 180.0)
228 )
229 self.check_evaluation_consistency(mixed)
231 def check_evaluation_consistency(self, field: BaseField) -> None:
232 """Test that `BaseField.__call__` and `BaseField.render` agree on a
233 concrete field.
234 """
235 image_1 = field.render()
236 p = field.bounds.bbox.meshgrid()
237 image_2 = Image(field(x=p.x, y=p.y), bbox=field.bounds.bbox, unit=field.unit)
238 assert_images_equal(self, image_1, image_2)
239 scaled_field = field * 2.0
240 image_3 = scaled_field.render()
241 image_3.array *= 0.5
242 assert_images_equal(self, image_1, image_3)
243 image_4 = field.render(Box.factory[9:11, -3:1])
244 assert_images_equal(self, image_4, image_1[image_4.bbox])
246 def check_units(self, field: BaseField) -> None:
247 """Test that the methods of a `BaseField` implementation correctly
248 propogate and check units.
249 """
250 self.assertIsNone(field.unit)
251 with_units_1 = field * u.nJy
252 self.assertEqual(with_units_1.unit, u.nJy)
253 self.assertEqual(with_units_1(x=np.array([0.0]), y=np.array([10.0]), quantity=True).unit, u.nJy)
254 image_1 = with_units_1.render(bbox=Box.factory[10:12, 0:3])
255 self.assertEqual(image_1.unit, u.nJy)
256 self.assertEqual((with_units_1 * 2.0).unit, u.nJy)
257 self.assertEqual((with_units_1 / u.arcsec**2).unit, u.nJy / u.arcsec**2)
260@unittest.skipUnless(HAVE_LEGACY, "This test requires lsst.afw.math to be importable.")
261class FieldLegacyTestCase(unittest.TestCase):
262 """Test the Field classes against legacy implementations.
264 This includes many tests for correct evaluation, since the legacy types
265 serve as our reference implementation.
266 """
268 def setUp(self) -> None:
269 self.rng = np.random.default_rng(10)
270 self.box = Box.factory[6:32, -7:26]
271 # This Chebyshev coefficient matrix is unusual in that it has nonzero
272 # entries for the whole matrix, not just the p + q < order triangle,
273 # and it has different orders for y and x.
274 # But we want to make sure those round-trip, too.
275 self.cheby_coeffs = self.rng.random((6, 3))
276 self.legacy_cheby = LegacyChebyshevBoundedField(self.box.to_legacy(), self.cheby_coeffs)
277 cheby2 = LegacyChebyshevBoundedField(self.box.to_legacy(), self.rng.standard_normal(size=(2, 2)))
278 self.legacy_product = LegacyProductBoundedField([self.legacy_cheby, cheby2])
280 def test_chebyshev_roundtrip(self) -> None:
281 """Test converting ChebyshevField from and to legacy, and serializing
282 it in between.
283 """
284 cheby = field_from_legacy(self.legacy_cheby)
285 assert isinstance(cheby, ChebyshevField)
286 compare_field_to_legacy(self, cheby, self.legacy_cheby, subimage_bbox=Box.factory[8:12, -3:2])
287 with RoundtripFits(self, cheby) as roundtrip:
288 pass
289 compare_field_to_legacy(
290 self, roundtrip.result, self.legacy_cheby, subimage_bbox=Box.factory[8:12, -3:2]
291 )
292 compare_field_to_legacy(
293 self, roundtrip.result, cheby.to_legacy(), subimage_bbox=Box.factory[8:12, -3:2]
294 )
296 def test_product_roundtrip(self) -> None:
297 """Test converting ProductField from and to legacy, and serializing
298 it in between.
299 """
300 product = field_from_legacy(self.legacy_product)
301 assert isinstance(product, ProductField)
302 compare_field_to_legacy(self, product, self.legacy_product, subimage_bbox=Box.factory[8:12, -3:2])
303 with RoundtripFits(self, product) as roundtrip:
304 pass
305 compare_field_to_legacy(
306 self, roundtrip.result, self.legacy_product, subimage_bbox=Box.factory[8:12, -3:2]
307 )
308 compare_field_to_legacy(
309 self, roundtrip.result, product.to_legacy(), subimage_bbox=Box.factory[8:12, -3:2]
310 )
312 def test_spline_simple(self) -> None:
313 """Test SplineField against `lsst.afw.math.BackgroundMI`, when there
314 is no missing data.
315 """
316 from lsst.afw.image import MaskedImageF
317 from lsst.afw.math import BackgroundMI
319 bins = MaskedImageF(Box.factory[0:5, 0:6].to_legacy())
320 bins.image.array[:, :] = self.rng.standard_normal(bins.image.array.shape)
321 bins.variance.array[::] = 1.0
322 legacy_bg = BackgroundMI(self.box.to_legacy(), bins)
323 spline = field_from_legacy_background(legacy_bg)
324 render_bbox = self.box.padded(-3)
325 assert_images_equal(
326 self,
327 spline.render(render_bbox),
328 Image.from_legacy(
329 legacy_bg.getImageF(
330 render_bbox.to_legacy(), legacy_bg.getBackgroundControl().getInterpStyle()
331 )
332 ),
333 rtol=1e-7,
334 )
336 def test_spline_one_nan(self) -> None:
337 """Test SplineField against `lsst.afw.math.BackgroundMI`, when there
338 is missing data.
339 """
340 from lsst.afw.image import MaskedImageF
341 from lsst.afw.math import BackgroundMI
343 bins = MaskedImageF(Box.factory[0:7, 0:6].to_legacy())
344 bins.image.array[:, :] = self.rng.standard_normal(bins.image.array.shape)
345 bins.image.array[3, 2] = np.nan
346 bins.variance.array[::] = 1.0
347 legacy_bg = BackgroundMI(self.box.to_legacy(), bins)
348 spline = field_from_legacy_background(legacy_bg)
349 render_bbox = self.box.padded(-3)
350 assert_images_equal(
351 self,
352 spline.render(render_bbox),
353 Image.from_legacy(
354 legacy_bg.getImageF(
355 render_bbox.to_legacy(),
356 legacy_bg.getBackgroundControl().getInterpStyle(),
357 )
358 ),
359 rtol=1e-7,
360 )
363@unittest.skipIf(DATA_DIR is None, "This test requires the TESTDATA_IMAGES_DIR envvar to be set.")
364@unittest.skipUnless(HAVE_LEGACY, "This test requires lsst.afw.math to be importable.")
365class FieldLegacyDataTestCase(unittest.TestCase):
366 """Tests for using Field classes using legacy datasets."""
368 def test_visit_background(self) -> None:
369 assert DATA_DIR is not None, "Guaranteed by decorator."
370 filename = os.path.join(DATA_DIR, "dp2", "legacy", "visit_image_background.fits")
371 legacy_bg_list = LegacyBackgroundList.readFits(filename)
372 bg_field = field_from_legacy_background(legacy_bg_list)
373 assert_images_equal(self, bg_field.render(), Image.from_legacy(legacy_bg_list.getImage()), rtol=1e-6)
376if __name__ == "__main__":
377 unittest.main()