Skip to content

Use RegularGridInterpolator for proper extrapolation of small PSF images #798

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
oczoske merged 2 commits into from
Sep 20, 2025
Merged

Use RegularGridInterpolator for proper extrapolation of small PSF images #798

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
0bbaa9e
Use RegularGridInterpolator for proper extrapolation of small psf images
oczoske Sep 17, 2025
f61afe8
Introduce psf_maxsize parameter
oczoske Sep 19, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
26 changes: 15 additions & 11 deletions scopesim/effects/psfs/discrete.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,8 @@
from tqdm.auto import tqdm
from scipy.signal import convolve
from scipy.ndimage import zoom
from scipy.interpolate import RectBivariateSpline, griddata
from scipy.interpolate import (RectBivariateSpline, griddata,
RegularGridInterpolator)

from astropy import units as u
from astropy.io import fits
Expand Down Expand Up @@ -150,10 +151,6 @@ def make_psf_cube(self, fov):
if not self.cmds.get("!OBS.interp_psf", True):
lam = np.array([lam[len(lam)//2]])

# adapt the size of the output cube to the FOV's spatial shape
nxpsf = min(512, 2 * nxfov + 1)
nypsf = min(512, 2 * nyfov + 1)

# Some data from the psf file
ext = self.current_layer_id
hdr = self._file[ext].header
Expand All @@ -168,15 +165,22 @@ def make_psf_cube(self, fov):

psfwcs = WCS(hdr)
psf = self._file[ext].data
psf = psf/psf.sum() # normalisation of the input psf
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.
# Note: According to some basic profiling, this line is one of the
# single largest hits on performance.
ipsf = RectBivariateSpline(np.arange(nyin), np.arange(nxin), psf,
kx=1, ky=1)
# 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)

# adapt the size of the output cube to the FOV's spatial shape
psf_maxsize = self.cmds.get("!SIM.computing.psf_maxsize", np.nan)

nxpsf = min(nxin, 2 * nxfov + 1, psf_maxsize)
nypsf = min(nyin, 2 * nyfov + 1, psf_maxsize)

xcube, ycube = np.meshgrid(np.arange(nxpsf), np.arange(nypsf))
cubewcs = WCS(naxis=2)
Expand All @@ -197,7 +201,7 @@ 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, xpsf, grid=False)
outcube[i,] = (ipsf( (ypsf.ravel(), xpsf.ravel()) ).reshape(outcube[i,].shape)
* fov_pixel_scale**2 / psf_wave_pixscale**2)

# .squeeze() gets rid of any axes with length one
Expand Down
Loading