From a53515aa02faa20b4e582f16dbe9439f681b832d Mon Sep 17 00:00:00 2001 From: oczoske Date: Thu, 25 Jun 2026 12:53:46 +0200 Subject: [PATCH 1/3] PSF interpolation method and parameters --- scopesim/effects/psfs/discrete.py | 29 ++++++++++++++--------------- 1 file changed, 14 insertions(+), 15 deletions(-) diff --git a/scopesim/effects/psfs/discrete.py b/scopesim/effects/psfs/discrete.py index 78a1a42b..92f7cfab 100644 --- a/scopesim/effects/psfs/discrete.py +++ b/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 @@ -221,8 +221,9 @@ def get_kernel(self, 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.meta["interp_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 @@ -510,7 +511,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 +524,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 +539,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 + print("_rescale_kernel: spline_order = ", spline_order) + 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,10 +569,10 @@ 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)) + outimage = interp_img(xin.ravel(), yin.ravel(), grid=False).reshape((nyout, nxout)) logger.info("Interpolating PSF onto %s image", outimage.shape) outimage *= image.sum() / outimage.sum() From 51fbb884cdf4ee03080d7b0678e9d4145986a9b2 Mon Sep 17 00:00:00 2001 From: oczoske Date: Sun, 5 Jul 2026 11:22:14 +0200 Subject: [PATCH 2/3] Interpolation parameters for PSF cube --- scopesim/effects/psfs/discrete.py | 46 ++++++++++++++++++------------- 1 file changed, 27 insertions(+), 19 deletions(-) diff --git a/scopesim/effects/psfs/discrete.py b/scopesim/effects/psfs/discrete.py index 92f7cfab..a5e84a07 100644 --- a/scopesim/effects/psfs/discrete.py +++ b/scopesim/effects/psfs/discrete.py @@ -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,8 +222,6 @@ 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( - self.meta["interp_order"], cmds=self.cmds) self.kernel = _rescale_kernel(self.kernel, pix_ratio, spline_order=spline_order, image_header=hdr) @@ -233,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"] @@ -246,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 @@ -266,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) @@ -291,7 +297,7 @@ 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)): @@ -299,11 +305,12 @@ def make_psf_cube(self, fov): 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): @@ -540,7 +547,7 @@ def _rescale_kernel(image, scale_factor, spline_order=1, convolved with an image with ``mode="same"``. """ nyin, nxin = image.shape - print("_rescale_kernel: spline_order = ", spline_order) + interp_img = RectBivariateSpline( np.arange(nxin), np.arange(nyin), image, kx=spline_order, ky=spline_order) @@ -572,9 +579,10 @@ def _rescale_kernel(image, scale_factor, spline_order=1, 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(xin.ravel(), yin.ravel(), grid=False).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 From 4e5b2f278d59d9e960e78dfe6fb7c1b95f86685b Mon Sep 17 00:00:00 2001 From: oczoske Date: Sun, 5 Jul 2026 15:15:56 +0200 Subject: [PATCH 3/3] Fix test (unmark one xpass); add default parameters --- scopesim/defaults.yaml | 4 ++++ scopesim/tests/tests_effects/test_ConstPSF.py | 17 +++++++++-------- 2 files changed, 13 insertions(+), 8 deletions(-) diff --git a/scopesim/defaults.yaml b/scopesim/defaults.yaml index 9c9cf8f7..b561bfa1 100644 --- a/scopesim/defaults.yaml +++ b/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/" diff --git a/scopesim/tests/tests_effects/test_ConstPSF.py b/scopesim/tests/tests_effects/test_ConstPSF.py index 7f6ae9a8..d6270dd1 100644 --- a/scopesim/tests/tests_effects/test_ConstPSF.py +++ b/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,