From 873f21f4ee83135876f62c88c82cec5bfb8f7cb3 Mon Sep 17 00:00:00 2001 From: Kieran Leschinski Date: Wed, 19 Aug 2026 15:58:48 +0200 Subject: [PATCH] Vectorise the wavelength resampling in XiLamImage The per-eta loop built a full 2D RectBivariateSpline over each cube plane and evaluated it exactly at the xi knot positions, i.e. the spline degenerated to independent 1D linear interpolation along the wavelength axis. Replace it with a direct searchsorted-based linear interpolation over all xi rows at once. Out-of-range wavelengths are clamped to the boundary values, exactly as the degree-1 spline did. For a METIS LSS field of view this removes ~1.7 s of spline construction per FieldOfView; for METIS LMS the loop runs 28 times per observation. Adds an equivalence test against the previous spline-based computation on random cube data (including the clamped edge region). The full basic_instrument spectroscopy/ifu suite passes unchanged. --- scopesim/effects/spectral_trace_list_utils.py | 17 +++++++-- .../test_SpectralTraceListUtils.py | 37 +++++++++++++++++-- 2 files changed, 47 insertions(+), 7 deletions(-) diff --git a/scopesim/effects/spectral_trace_list_utils.py b/scopesim/effects/spectral_trace_list_utils.py index d523aee7..a01d6cb4 100644 --- a/scopesim/effects/spectral_trace_list_utils.py +++ b/scopesim/effects/spectral_trace_list_utils.py @@ -741,10 +741,19 @@ def __init__(self, fov, dlam_per_pix): # lam0 is the target wavelength. We need to check that this # overlaps with the wavelength range covered by the cube if lam0.min() < cube_lam.max() and lam0.max() > cube_lam.min(): - plane = fov.cube.data[:, i, :].T - plane_interp = RectBivariateSpline(cube_xi, cube_lam, plane, - kx=1, ky=1) - self.image += plane_interp(cube_xi, lam0) + # Linear interpolation of the cube plane onto lam0 along + # the wavelength axis (the xi axis is not resampled). + # Values outside the cube's wavelength range are clamped + # to the boundary values, matching the behaviour of a + # degree-1 RectBivariateSpline. + plane = fov.cube.data[:, i, :] # (n_lam, n_xi) + jlo = np.clip(np.searchsorted(cube_lam, lam0) - 1, + 0, n_lam - 2) + weight = np.clip((lam0 - cube_lam[jlo]) + / (cube_lam[jlo + 1] - cube_lam[jlo]), + 0.0, 1.0) + self.image += (plane[jlo, :] * (1 - weight[:, None]) + + plane[jlo + 1, :] * weight[:, None]).T self.image *= d_eta # ph/s/um/arcsec2 --> ph/s/um/arcsec diff --git a/scopesim/tests/tests_effects/test_SpectralTraceListUtils.py b/scopesim/tests/tests_effects/test_SpectralTraceListUtils.py index ca2f19d2..4b0429d8 100644 --- a/scopesim/tests/tests_effects/test_SpectralTraceListUtils.py +++ b/scopesim/tests/tests_effects/test_SpectralTraceListUtils.py @@ -161,7 +161,7 @@ def test_grid_false_shape_is_preserved(self, tf2d): class MockCubeFov: """Minimal stand-in for a FieldOfView carrying a spectral cube.""" - def __init__(self, n_lam=20, n_eta=3, n_xi=11): + def __init__(self, n_lam=20, n_eta=3, n_xi=11, data=None): hdr = fits.Header() hdr["NAXIS"] = 3 hdr["NAXIS1"], hdr["NAXIS2"], hdr["NAXIS3"] = n_xi, n_eta, n_lam @@ -170,8 +170,9 @@ def __init__(self, n_lam=20, n_eta=3, n_xi=11): hdr["CRVAL1"], hdr["CRVAL2"], hdr["CRVAL3"] = 0., 0., 2.0 hdr["CDELT1"], hdr["CDELT2"], hdr["CDELT3"] = 0.1, 0.1, 0.001 hdr["CUNIT1"], hdr["CUNIT2"], hdr["CUNIT3"] = "arcsec", "arcsec", "um" - self.cube = fits.ImageHDU( - data=np.ones((n_lam, n_eta, n_xi)), header=hdr) + if data is None: + data = np.ones((n_lam, n_eta, n_xi)) + self.cube = fits.ImageHDU(data=data, header=hdr) self.meta = {"xi_min": -0.5 * u.arcsec, "xi_max": 0.5 * u.arcsec} @@ -182,6 +183,36 @@ def test_primary_wcs_keeps_arcsec_cunit(self): assert list(xilam.wcs.wcs.cunit) == [u.um, u.arcsec] assert list(xilam.wcsa.wcs.cunit) == [u.um, u.dimensionless_unscaled] + def test_image_matches_rect_bivariate_spline_reference(self): + """The vectorised wavelength interpolation must reproduce the + previous per-plane RectBivariateSpline(kx=1, ky=1) computation.""" + from scipy.interpolate import RectBivariateSpline + + n_lam, n_eta, n_xi = 40, 5, 13 + rng = np.random.default_rng(42) + data = rng.random((n_lam, n_eta, n_xi)) + fov = MockCubeFov(n_lam, n_eta, n_xi, data=data) + dlam_per_pix = 0.0015 + + xilam = XiLamImage(fov, dlam_per_pix) + + # reference: the removed spline-per-plane implementation + hdr = fov.cube.header + d_xi, d_eta, d_lam = hdr["CDELT1"], hdr["CDELT2"], hdr["CDELT3"] + cube_xi = d_xi * np.arange(n_xi) + fov.meta["xi_min"].value + cube_eta = d_eta * (np.arange(n_eta) - (n_eta - 1) / 2) + cube_lam = hdr["CRVAL3"] + d_lam * np.arange(n_lam) + reference = np.zeros((n_xi, n_lam)) + for i, eta in enumerate(cube_eta): + lam0 = cube_lam + dlam_per_pix * eta / d_eta + plane = data[:, i, :].T + spline = RectBivariateSpline(cube_xi, cube_lam, plane, + kx=1, ky=1) + reference += spline(cube_xi, lam0) + reference *= d_eta + + np.testing.assert_allclose(xilam.image, reference, rtol=1e-6) + class TestImageInterpolations: """Tests for function make_image_interpolations"""