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