From 5dcdee44ac22a75bdf8b0ea40e93ccca6b52cdcd Mon Sep 17 00:00:00 2001 From: Kieran Leschinski Date: Wed, 19 Aug 2026 16:03:35 +0200 Subject: [PATCH] Add fast bilinear path to make_image_interpolations Rectification evaluates the detector-image interpolators at every pixel of the rectified grid via grid=False, which goes through FITPACK's slow scattered-point path. For kx=ky=1 (the case used by SpectralTrace.rectify and MetisLMS rectify_cube) return a direct bilinear interpolator with identical results instead; other spline degrees are unchanged and keep using RectBivariateSpline. The interpolator keeps the (j, i) call convention and the clamping behaviour at the image edges. Adds an equivalence test against RectBivariateSpline (including out-of-image coordinates) and a fallback-type test. --- scopesim/effects/spectral_trace_list_utils.py | 43 ++++++++++++++++++- .../test_SpectralTraceListUtils.py | 25 +++++++++++ 2 files changed, 67 insertions(+), 1 deletion(-) diff --git a/scopesim/effects/spectral_trace_list_utils.py b/scopesim/effects/spectral_trace_list_utils.py index da583fdc..8cc9aefc 100644 --- a/scopesim/effects/spectral_trace_list_utils.py +++ b/scopesim/effects/spectral_trace_list_utils.py @@ -1021,16 +1021,57 @@ def _xiy2xlam_fit(layout, params): return xiy2x, xiy2lam +class BilinearImageInterpolation: + """ + Bilinear interpolation of an image on its integer pixel grid. + + Drop-in replacement for + ``RectBivariateSpline(arange(ny), arange(nx), image, kx=1, ky=1)``: + called with (j, i) pixel coordinates and ``grid=False``, it returns the + same values (coordinates outside the image are clamped to the edges), + but avoids FITPACK's slow scattered-point evaluation. + """ + + def __init__(self, image): + self.image = np.asarray(image) + + def __call__(self, jarr, iarr, grid=False): + if grid: + raise NotImplementedError( + "BilinearImageInterpolation only supports grid=False") + n_j, n_i = self.image.shape + jarr = np.clip(np.asarray(jarr, dtype=float), 0, n_j - 1) + iarr = np.clip(np.asarray(iarr, dtype=float), 0, n_i - 1) + j_0 = np.clip(jarr.astype(int), 0, n_j - 2) + i_0 = np.clip(iarr.astype(int), 0, n_i - 2) + w_j = jarr - j_0 + w_i = iarr - i_0 + img = self.image + return ((1 - w_j) * (1 - w_i) * img[j_0, i_0] + + (1 - w_j) * w_i * img[j_0, i_0 + 1] + + w_j * (1 - w_i) * img[j_0 + 1, i_0] + + w_j * w_i * img[j_0 + 1, i_0 + 1]) + + def make_image_interpolations(hdulist, **kwargs): """Create 2D interpolation functions for images. The interpolation functions are called with (j, i) pixel coordinates, i.e. the first argument corresponds to the row (NAXIS2) axis of the image, the second to the column (NAXIS1) axis. + + For ``kx=1, ky=1`` (the case used for rectification) a direct bilinear + interpolator is returned, which evaluates large scattered coordinate + arrays orders of magnitude faster than FITPACK. Other spline degrees + fall back to ``RectBivariateSpline``. """ interps = [] for hdu in hdulist: - if isinstance(hdu, fits.ImageHDU): + if not isinstance(hdu, fits.ImageHDU): + continue + if kwargs.get("kx") == 1 and kwargs.get("ky") == 1: + interps.append(BilinearImageInterpolation(hdu.data)) + else: interps.append( RectBivariateSpline(np.arange(hdu.header["NAXIS2"]), np.arange(hdu.header["NAXIS1"]), diff --git a/scopesim/tests/tests_effects/test_SpectralTraceListUtils.py b/scopesim/tests/tests_effects/test_SpectralTraceListUtils.py index e1bcbb46..2e8404c8 100644 --- a/scopesim/tests/tests_effects/test_SpectralTraceListUtils.py +++ b/scopesim/tests/tests_effects/test_SpectralTraceListUtils.py @@ -203,6 +203,31 @@ def test_interpolation_is_accurate(self): imginterp = interps[0](yy, xx, grid=False) assert np.allclose(imginterp, img) + def test_bilinear_matches_spline_evaluation(self): + """The kx=ky=1 fast path must reproduce RectBivariateSpline.""" + from scipy.interpolate import RectBivariateSpline + + rng = np.random.default_rng(3) + img = rng.random((50, 80)) + hdul = fits.HDUList([fits.PrimaryHDU(), fits.ImageHDU(data=img)]) + fast = make_image_interpolations(hdul, kx=1, ky=1)[0] + spline = RectBivariateSpline(np.arange(50), np.arange(80), img, + kx=1, ky=1) + + # scattered points, deliberately extending beyond the image + jarr = rng.uniform(-2, 52, 1000) + iarr = rng.uniform(-2, 82, 1000) + np.testing.assert_allclose(fast(jarr, iarr, grid=False), + spline(jarr, iarr, grid=False), + rtol=1e-10) + + def test_other_spline_degrees_fall_back_to_fitpack(self): + from scipy.interpolate import RectBivariateSpline + img = np.random.rand(20, 20) + hdul = fits.HDUList([fits.PrimaryHDU(), fits.ImageHDU(data=img)]) + interps = make_image_interpolations(hdul) # default: cubic + assert isinstance(interps[0], RectBivariateSpline) + def test_works_with_non_square_image(self): # The axes passed to RectBivariateSpline were transposed, which # raised ValueError for any image with NAXIS1 != NAXIS2