PSF interpolation order

scopesim/defaults.yaml

@@ -32,6 +32,10 @@ properties :
     preload_field_of_views : False
     nan_fill_value: 0.
 
+  psf:
+    inter_psf: True
+    interp_order: 1
+
   file :
     local_packages_path : "./inst_pkgs/"
     server_base_url : "https://scopesim.univie.ac.at/InstPkgSvr/"

scopesim/effects/psfs/discrete.py

@@ -8,7 +8,7 @@
 from tqdm.auto import tqdm
 from scipy.signal import convolve
 from scipy.ndimage import zoom
-from scipy.interpolate import griddata, RegularGridInterpolator
+from scipy.interpolate import griddata, RectBivariateSpline, RegularGridInterpolator
 
 from astropy import units as u
 from astropy.io import fits
@@ -166,7 +166,7 @@ class FieldConstantPSF(DiscretePSF):
     .. versionchanged:: 0.10.0
 
        PSF interpolation can now be limited to the central wavelength by
-       setting the "!OBS.interp_psf" keyword to False.
+       setting the "!SIM.psf.interp_psf" keyword to False.
 
     .. versionchanged:: 0.11.2
 
@@ -197,9 +197,11 @@ def get_kernel(self, fov):
         if ext == self.current_layer_id:
             return self.kernel
 
+        spline_order = from_currsys(self.meta["interp_order"], cmds=self.cmds)
+
         if fov.hdu.header["NAXIS"] == 3:
             self.current_layer_id = ext
-            self.make_psf_cube(fov)
+            self.make_psf_cube(fov, spline_order)
             return self.kernel
 
         self.kernel = self._file[ext].data
@@ -220,9 +222,8 @@ def get_kernel(self, fov):
         # rescaling kept inside loop to avoid rescaling for every fov
         pix_ratio = kernel_pixel_scale / fov_pixel_scale
         if abs(pix_ratio - 1) > self.meta["flux_accuracy"]:
-            spline_order = from_currsys(
-                "!SIM.computing.spline_order", cmds=self.cmds)
             self.kernel = _rescale_kernel(self.kernel, pix_ratio,
+                                          spline_order=spline_order,
                                           image_header=hdr)
 
         if ((fov.header["NAXIS1"] < hdr["NAXIS1"]) or
@@ -232,8 +233,19 @@ def get_kernel(self, fov):
 
         return np.atleast_2d(self.kernel.squeeze())
 
-    def make_psf_cube(self, fov):
-        """Create a wavelength-dependent psf cube."""
+    def make_psf_cube(self, fov, spline_order=1):
+        """Create a wavelength-dependent psf cube.
+
+        Parameters
+        ----------
+        fov : FieldOfView
+
+        spline_order : int
+            spline order for interpolation by ``RectBivariateSpline``.
+            Default is 1 (linear interpolation). The method uses the same order in both
+            directions (i.e. `kx = ky = spline_order`).
+
+        """
         # Some data from the fov
         nxfov, nyfov = fov.hdu.header["NAXIS1"], fov.hdu.header["NAXIS2"]
         fov_pixel_scale = fov.hdu.header["CDELT1"]
@@ -245,7 +257,7 @@ def make_psf_cube(self, fov):
             + fov.hdu.header["CRVAL3"]
         )
 
-        if not self.cmds.get("!OBS.interp_psf", True):
+        if not self.cmds.get("!SIM.psf.interp_psf", True):
             lam = np.array([lam[len(lam)//2]])
 
         # Some data from the psf file
@@ -265,13 +277,8 @@ def make_psf_cube(self, fov):
         psf = psf / psf.sum()         # normalisation of the input psf
         nxin, nyin = psf.shape
 
-        # We need linear interpolation to preserve positivity. Might think of
-        # more elaborate positivity-preserving schemes.
-        # An alternative would be to use method 'pchip', which uses monotonic
-        # splines to prevent overshooting (including to negative values)
-        ipsf = RegularGridInterpolator((np.arange(nyin), np.arange(nxin)),
-                                       psf, method='linear', bounds_error=False,
-                                       fill_value=0)
+        ipsf = RectBivariateSpline(
+            np.arange(nxin), np.arange(nyin), psf, kx=spline_order, ky=spline_order)
 
         # adapt the size of the output cube to the FOV's spatial shape
         psf_maxsize = self.cmds.get("!SIM.computing.psf_maxsize", np.nan)
@@ -290,19 +297,20 @@ def make_psf_cube(self, fov):
         xworld, yworld = cubewcs.all_pix2world(xcube, ycube, 1)
         outcube = np.zeros((lam.shape[0], nypsf, nxpsf), dtype=np.float32)
         logger.info("Interpolating PSF onto %s cube", outcube.shape)
-
+        logger.info("Interpolation order %d", spline_order)
         nlam = len(lam)
         for i, wave in tqdm(enumerate(lam), desc=" PSF slices", position=2,
                             total=nlam, disable=(nlam <= 1)):
             psf_wave_pixscale = ref_pixel_scale * wave / refwave
             psfwcs.wcs.cdelt = [psf_wave_pixscale,
                                 psf_wave_pixscale]
             xpsf, ypsf = psfwcs.all_world2pix(xworld, yworld, 0)
-            outcube[i,] = (ipsf( (ypsf.ravel(), xpsf.ravel()) ).reshape(outcube[i,].shape)
+            outcube[i,] = (ipsf( xpsf.ravel(), ypsf.ravel(), grid=False ).reshape(outcube[i,].shape)
                            * fov_pixel_scale**2 / psf_wave_pixscale**2)
 
         # .squeeze() gets rid of any axes with length one
         self.kernel = outcube.reshape((lam.shape[0], nypsf, nxpsf)).squeeze()
+
         # fits.writeto("test_psfcube.fits", data=self.kernel, overwrite=True)
 
     def plot(self):
@@ -510,7 +518,7 @@ def _make_strehl_map_from_table(tbl, pixel_scale=1*u.arcsec):
     return map_hdu
 
 
-def _rescale_kernel(image, scale_factor, method="linear",
+def _rescale_kernel(image, scale_factor, spline_order=1,
                     image_header=None):
     """Rescale `image` by `scale_factor`
 
@@ -523,9 +531,10 @@ def _rescale_kernel(image, scale_factor, method="linear",
         ratio of input pixel size to output pixel size. Values larger than 1
         zoom the image.
 
-    method : str
-        interpolation method to be passed to ``RegularGridInterpolator``.
-        Default is "linear"
+    spline_order : int
+        spline order for interpolation by ``RectBivariateSpline``.
+        Default is 1 (linear interpolation). The method uses the same order in both
+        directions (i.e. `kx = ky = spline_order`).
 
     image_header : astropy.fits.Header
         an optional image header with a WCS. If ``None``, a WCS is constructed
@@ -537,13 +546,10 @@ def _rescale_kernel(image, scale_factor, method="linear",
     The output kernel size is always odd to avoid half-pixel shift when
     convolved with an image with ``mode="same"``.
     """
-    nxin, nyin = image.shape
-    interp_img = RegularGridInterpolator(
-        (np.arange(nyin), np.arange(nxin)),
-        image, method=method,
-        bounds_error=False,
-        fill_value=0,
-    )
+    nyin, nxin = image.shape
+
+    interp_img = RectBivariateSpline(
+        np.arange(nxin), np.arange(nyin), image, kx=spline_order, ky=spline_order)
 
     # Make the output size odd so that the image remains centred (important
     # for PSF convolution).
@@ -570,12 +576,13 @@ def _rescale_kernel(image, scale_factor, method="linear",
         outwcs.wcs.crpix = [(nxout + 1)/2, (nyout + 1)/2]
         outwcs.wcs.crval = [0., 0.]
         outwcs.wcs.cdelt = inwcs.wcs.cdelt / scale_factor
-    xout, yout = np.meshgrid(np.arange(nxout), np.arange(nyout))
+    xout, yout = np.meshgrid(np.arange(nxout), np.arange(nyout), indexing='ij')
     xworld, yworld = outwcs.all_pix2world(xout, yout, 0)
     xin, yin = inwcs.all_world2pix(xworld, yworld, 0)
-    outimage = interp_img( (yin.ravel(), xin.ravel()) ).reshape((nyout, nxout))
-    logger.info("Interpolating PSF onto %s image", outimage.shape)
 
+    logger.info("Interpolating PSF onto %s image", (nyout, nxout))
+    logger.info("Interpolation order %d", spline_order)
+    outimage = interp_img(xin.ravel(), yin.ravel(), grid=False).reshape((nyout, nxout))
     outimage *= image.sum() / outimage.sum()
 
     return outimage

scopesim/tests/tests_effects/test_ConstPSF.py

@@ -17,7 +17,6 @@
 """
 
 import pytest
-from pytest import approx
 
 import numpy as np
 from numpy import testing as npt
@@ -35,15 +34,17 @@
 # pylint: disable=missing-class-docstring,
 # pylint: disable=missing-function-docstring
 
-@pytest.fixture(scope="function")
-def centre_fov():
+@pytest.fixture(name="centre_fov", scope="function")
+def fixture_centre_fov():
     return _centre_fov()
 
 
-@pytest.fixture(scope="function")
-def const_psf(mock_path):
+@pytest.fixture(name="const_psf", scope="function")
+def fixture_const_psf(mock_path):
+    """Instantiate a field constant psf with linear interpolation"""
     constpsf = FieldConstantPSF(
-        filename=str(mock_path / "test_ConstPSF.fits"))
+        filename=str(mock_path / "test_ConstPSF.fits"),
+        interp_order=1)
     return constpsf
 
 
@@ -98,9 +99,9 @@ def test_returns_array_with_single_kernel_from_fov(
     @pytest.mark.parametrize("scale_factor", (
         pytest.param(1/3, marks=pytest.mark.xfail(reason="PSF center off")),
         1,
-        pytest.param(5, marks=pytest.mark.xfail(reason="1x1 array turned into scalar")),
+        5
     ))
-    def test_kernel_is_scale_properly_if_cdelts_differ(
+    def test_kernel_is_scaled_properly_if_cdelts_differ(
         self,
         scale_factor,
         const_psf,