From c04bfc640b5214b1125e7dafc93c1155122030ed Mon Sep 17 00:00:00 2001 From: oczoske Date: Thu, 19 Feb 2026 13:25:56 +0100 Subject: [PATCH 1/7] Fixes off-by-a-half in FOVs; needs cleaning --- scopesim/effects/psfs/discrete.py | 65 +++++++++++++++++++------------ scopesim/effects/psfs/psf_base.py | 4 ++ 2 files changed, 45 insertions(+), 24 deletions(-) diff --git a/scopesim/effects/psfs/discrete.py b/scopesim/effects/psfs/discrete.py index 3da2da59..c5885242 100644 --- a/scopesim/effects/psfs/discrete.py +++ b/scopesim/effects/psfs/discrete.py @@ -205,7 +205,7 @@ def get_kernel(self, fov): self.kernel = self._file[ext].data self.current_layer_id = ext hdr = self._file[ext].header - + print(hdr['EXTNAME']) self.kernel /= np.sum(self.kernel) # compare kernel and fov pixel scales, rescale if needed @@ -222,7 +222,7 @@ def get_kernel(self, fov): 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) + self.kernel = _rescale_kernel(self.kernel, pix_ratio, hdr) if ((fov.header["NAXIS1"] < hdr["NAXIS1"]) or (fov.header["NAXIS2"] < hdr["NAXIS2"])): @@ -453,7 +453,7 @@ def get_kernel(self, fov): "!SIM.computing.spline_order", cmds=self.cmds) for ii, kern in enumerate(self.kernel): self.kernel[ii][0] = _rescale_kernel( - kern[0], pix_ratio, spline_order) + kern[0], pix_ratio) for i, kern in enumerate(self.kernel): self.kernel[i][0] /= np.sum(kern[0]) @@ -510,31 +510,48 @@ def _make_strehl_map_from_table(tbl, pixel_scale=1*u.arcsec): return map_hdu -def _rescale_kernel(image, scale_factor, spline_order): +def _rescale_kernel(image, scale_factor, image_header=None): sum_image = np.sum(image) - image = zoom(image, scale_factor, order=spline_order) - image = np.nan_to_num(image, copy=False) # numpy version >=1.13 - - # Re-centre kernel - im_shape = image.shape - # TODO: this might be another off-by-something - dy, dx = np.divmod(np.argmax(image), im_shape[1]) - np.array(im_shape) // 2 - if dy > 0: - image = image[2*dy:, :] - elif dy < 0: - image = image[:2*dy, :] - if dx > 0: - image = image[:, 2*dx:] - elif dx < 0: - image = image[:, :2*dx] - - sum_new_image = np.sum(image) - image *= sum_image / sum_new_image - - return image + nxin, nyin = image.shape + iimg = RegularGridInterpolator((np.arange(nyin), np.arange(nxin)), + image, method='linear', bounds_error=False, + fill_value=0) + nxout = int(nxin * scale_factor) + if nxout % 2 != 0: + nxout += 1 + nyout = int(nyin * scale_factor) + if nyout % 2 != 0: + nyout += 1 + + if image_header is not None: + inwcs = WCS(image_header) + outwcs = WCS(image_header) + outwcs.wcs.crpix = [(nxout + 1)/2, (nyout + 1)/2] + outwcs.wcs.cdelt = inwcs.wcs.cdelt / scale_factor + else: + inwcs = WCS(naxis=2) + inwcs.wcs.ctype = ["LINEAR", "LINEAR"] + inwcs.wcs.crpix = [(nxin + 1)/2, (nyin + 1)/2] + inwcs.wcs.crval = [0., 0.] + inwcs.wcs.cdelt = [1, 1] + outwcs = WCS(naxis=2) + outwcs.wcs.ctype = ["LINEAR", "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)) + xworld, yworld = outwcs.all_pix2world(xout, yout, 0) + xin, yin = inwcs.all_world2pix(xworld, yworld, 0) + outimage = iimg( (yin.ravel(), xin.ravel()) ).reshape((nyout, nxout)) + logger.info("Interpolating PSF onto %s image", outimage.shape) + + outimage *= sum_image / outimage.sum() + + return outimage def _cutout_kernel(image, fov_header, kernel_header=None): + logger.info("Going into _cutout_kernel") wk = WCS(kernel_header) h, w = image.shape xcen, ycen = 0.5 * w, 0.5 * h diff --git a/scopesim/effects/psfs/psf_base.py b/scopesim/effects/psfs/psf_base.py index 613ac012..1796c0ef 100644 --- a/scopesim/effects/psfs/psf_base.py +++ b/scopesim/effects/psfs/psf_base.py @@ -101,6 +101,10 @@ def apply_to(self, obj, **kwargs): logger.debug("PSF convolution start") if image.ndim == 2 and kernel.ndim == 2: new_image = convolve(image - bkg_level, kernel, mode=mode) + from astropy.io import fits + fits.writeto(f"image_{obj.meta['wave_min'].value}.fits", data=image.value, overwrite=True) + fits.writeto(f"kernel_{obj.meta['wave_min'].value}.fits", kernel, overwrite=True) + fits.writeto(f"newimage_{obj.meta['wave_min'].value}.fits", new_image, overwrite=True) elif image.ndim == 3 and kernel.ndim == 2: kernel = kernel[None, :, :] bkg_level = bkg_level[:, None, None] From 09e235d9f35c956661d816503e2793b3447a9cb5 Mon Sep 17 00:00:00 2001 From: oczoske Date: Thu, 19 Feb 2026 13:36:36 +0100 Subject: [PATCH 2/7] Fix for tests --- scopesim/effects/psfs/psf_base.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/scopesim/effects/psfs/psf_base.py b/scopesim/effects/psfs/psf_base.py index 1796c0ef..85930b51 100644 --- a/scopesim/effects/psfs/psf_base.py +++ b/scopesim/effects/psfs/psf_base.py @@ -101,10 +101,10 @@ def apply_to(self, obj, **kwargs): logger.debug("PSF convolution start") if image.ndim == 2 and kernel.ndim == 2: new_image = convolve(image - bkg_level, kernel, mode=mode) - from astropy.io import fits - fits.writeto(f"image_{obj.meta['wave_min'].value}.fits", data=image.value, overwrite=True) - fits.writeto(f"kernel_{obj.meta['wave_min'].value}.fits", kernel, overwrite=True) - fits.writeto(f"newimage_{obj.meta['wave_min'].value}.fits", new_image, overwrite=True) + #from astropy.io import fits + #fits.writeto(f"image_{obj.meta['wave_min'].value}.fits", data=image.value, overwrite=True) + #fits.writeto(f"kernel_{obj.meta['wave_min'].value}.fits", kernel, overwrite=True) + #fits.writeto(f"newimage_{obj.meta['wave_min'].value}.fits", new_image, overwrite=True) elif image.ndim == 3 and kernel.ndim == 2: kernel = kernel[None, :, :] bkg_level = bkg_level[:, None, None] From 74d9baa623c60907376023332343ee0f849a1260 Mon Sep 17 00:00:00 2001 From: oczoske Date: Mon, 23 Mar 2026 16:30:38 +0100 Subject: [PATCH 3/7] add method to _rescale_kernel --- scopesim/effects/psfs/discrete.py | 40 ++++++++++++++++++++++++++----- 1 file changed, 34 insertions(+), 6 deletions(-) diff --git a/scopesim/effects/psfs/discrete.py b/scopesim/effects/psfs/discrete.py index c5885242..2ba5fe61 100644 --- a/scopesim/effects/psfs/discrete.py +++ b/scopesim/effects/psfs/discrete.py @@ -205,7 +205,7 @@ def get_kernel(self, fov): self.kernel = self._file[ext].data self.current_layer_id = ext hdr = self._file[ext].header - print(hdr['EXTNAME']) + self.kernel /= np.sum(self.kernel) # compare kernel and fov pixel scales, rescale if needed @@ -222,7 +222,8 @@ def get_kernel(self, fov): 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, hdr) + self.kernel = _rescale_kernel(self.kernel, pix_ratio, + image_header=hdr) if ((fov.header["NAXIS1"] < hdr["NAXIS1"]) or (fov.header["NAXIS2"] < hdr["NAXIS2"])): @@ -452,8 +453,7 @@ def get_kernel(self, fov): spline_order = from_currsys( "!SIM.computing.spline_order", cmds=self.cmds) for ii, kern in enumerate(self.kernel): - self.kernel[ii][0] = _rescale_kernel( - kern[0], pix_ratio) + self.kernel[ii][0] = _rescale_kernel(kern[0], pix_ratio) for i, kern in enumerate(self.kernel): self.kernel[i][0] /= np.sum(kern[0]) @@ -510,12 +510,40 @@ def _make_strehl_map_from_table(tbl, pixel_scale=1*u.arcsec): return map_hdu -def _rescale_kernel(image, scale_factor, image_header=None): +def _rescale_kernel(image, scale_factor, method="linear", image_header=None): + """Rescale `image` by `scale_factor` + + Parameters + ---------- + image : array_like, 2D + the image to be rescaled + + scale_factor : float + 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" + + image_header : astropy.fits.Header + an optional image header with a WCS. If ``None``, a WCS is constructed + that counts image pixels. + + Notes + ----- + The function uses ``scipy.interpolate.RegularGridInterpolator``. + """ sum_image = np.sum(image) nxin, nyin = image.shape iimg = RegularGridInterpolator((np.arange(nyin), np.arange(nxin)), - image, method='linear', bounds_error=False, + image, method=method, bounds_error=False, fill_value=0) + # Adjust the output size so that the image remains centred (important + # for PSF convolution). + # TODO: This assumes that the PSF is applied to an image with even size + # The adjustment for odd sizes requires knowledge about that image + # that this function does not have. nxout = int(nxin * scale_factor) if nxout % 2 != 0: nxout += 1 From 56c468a2a8657147410ffa8a6ceceeca01df3d85 Mon Sep 17 00:00:00 2001 From: oczoske Date: Wed, 1 Apr 2026 15:45:26 +0200 Subject: [PATCH 4/7] Output kernel size must be odd --- scopesim/effects/psfs/discrete.py | 29 ++++++++++++++++------------- 1 file changed, 16 insertions(+), 13 deletions(-) diff --git a/scopesim/effects/psfs/discrete.py b/scopesim/effects/psfs/discrete.py index 2ba5fe61..bd0ee721 100644 --- a/scopesim/effects/psfs/discrete.py +++ b/scopesim/effects/psfs/discrete.py @@ -510,7 +510,8 @@ def _make_strehl_map_from_table(tbl, pixel_scale=1*u.arcsec): return map_hdu -def _rescale_kernel(image, scale_factor, method="linear", image_header=None): +def _rescale_kernel(image, scale_factor, method="linear", + image_header=None): """Rescale `image` by `scale_factor` Parameters @@ -533,22 +534,24 @@ def _rescale_kernel(image, scale_factor, method="linear", image_header=None): Notes ----- The function uses ``scipy.interpolate.RegularGridInterpolator``. + The output kernel size is always odd to avoid half-pixel shift when + convolved with an image with ``mode="same"``. """ - sum_image = np.sum(image) nxin, nyin = image.shape - iimg = RegularGridInterpolator((np.arange(nyin), np.arange(nxin)), - image, method=method, bounds_error=False, - fill_value=0) - # Adjust the output size so that the image remains centred (important + interp_img = RegularGridInterpolator( + (np.arange(nyin), np.arange(nxin)), + image, method=method, + bounds_error=False, + fill_value=0, + ) + + # Make the output size odd so that the image remains centred (important # for PSF convolution). - # TODO: This assumes that the PSF is applied to an image with even size - # The adjustment for odd sizes requires knowledge about that image - # that this function does not have. nxout = int(nxin * scale_factor) - if nxout % 2 != 0: + if nxout % 2 == 0: nxout += 1 nyout = int(nyin * scale_factor) - if nyout % 2 != 0: + if nyout % 2 == 0: nyout += 1 if image_header is not None: @@ -570,10 +573,10 @@ def _rescale_kernel(image, scale_factor, method="linear", image_header=None): xout, yout = np.meshgrid(np.arange(nxout), np.arange(nyout)) xworld, yworld = outwcs.all_pix2world(xout, yout, 0) xin, yin = inwcs.all_world2pix(xworld, yworld, 0) - outimage = iimg( (yin.ravel(), xin.ravel()) ).reshape((nyout, nxout)) + outimage = interp_img( (yin.ravel(), xin.ravel()) ).reshape((nyout, nxout)) logger.info("Interpolating PSF onto %s image", outimage.shape) - outimage *= sum_image / outimage.sum() + outimage *= image.sum() / outimage.sum() return outimage From 37a07e2ad8fde76fd2a9bc4139f3cee60641a227 Mon Sep 17 00:00:00 2001 From: oczoske Date: Wed, 1 Apr 2026 16:00:34 +0200 Subject: [PATCH 5/7] Test for oddness of rescaled kernel --- scopesim/tests/tests_effects/test_PSF.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/scopesim/tests/tests_effects/test_PSF.py b/scopesim/tests/tests_effects/test_PSF.py index 6b70a679..800cd36c 100644 --- a/scopesim/tests/tests_effects/test_PSF.py +++ b/scopesim/tests/tests_effects/test_PSF.py @@ -6,6 +6,7 @@ from scopesim.effects import PSF from scopesim.effects.psfs.psf_base import get_bkg_level +from scopesim.effects.psfs.discrete import _rescale_kernel from scopesim.optics import ImagePlane from scopesim.tests.mocks.py_objects.header_objects import _implane_header @@ -42,6 +43,17 @@ def test_returns_array_from_kernel(self): assert np.sum(psf.kernel) == approx(np.sum(psf.get_kernel(None))) + @pytest.mark.parametrize("nkern, scale_factor", [(128, 1.3), + (129, 1.3), + (512, 2.41), + (511, 2.41)]) + def test_rescale_produces_odd_size_kernel(self, nkern, scale_factor): + psf = PSF() + psf.kernel = basic_kernel(n=nkern) + outkern = _rescale_kernel(psf.kernel, scale_factor) + n_y, n_x = outkern.shape + assert n_y % 2 == 1 + assert n_x % 2 == 1 class TestRotationBlur: @pytest.mark.parametrize("angle", ([1, 5, 15, 60])) From 369108545acba3e99f44369dcae8502293d74885 Mon Sep 17 00:00:00 2001 From: oczoske Date: Thu, 2 Apr 2026 11:04:45 +0200 Subject: [PATCH 6/7] Remove debugging leftovers --- scopesim/effects/psfs/discrete.py | 2 +- scopesim/effects/psfs/psf_base.py | 4 ---- 2 files changed, 1 insertion(+), 5 deletions(-) diff --git a/scopesim/effects/psfs/discrete.py b/scopesim/effects/psfs/discrete.py index bd0ee721..19f7b0bf 100644 --- a/scopesim/effects/psfs/discrete.py +++ b/scopesim/effects/psfs/discrete.py @@ -582,7 +582,7 @@ def _rescale_kernel(image, scale_factor, method="linear", def _cutout_kernel(image, fov_header, kernel_header=None): - logger.info("Going into _cutout_kernel") + logger.debug("Going into _cutout_kernel") wk = WCS(kernel_header) h, w = image.shape xcen, ycen = 0.5 * w, 0.5 * h diff --git a/scopesim/effects/psfs/psf_base.py b/scopesim/effects/psfs/psf_base.py index 85930b51..613ac012 100644 --- a/scopesim/effects/psfs/psf_base.py +++ b/scopesim/effects/psfs/psf_base.py @@ -101,10 +101,6 @@ def apply_to(self, obj, **kwargs): logger.debug("PSF convolution start") if image.ndim == 2 and kernel.ndim == 2: new_image = convolve(image - bkg_level, kernel, mode=mode) - #from astropy.io import fits - #fits.writeto(f"image_{obj.meta['wave_min'].value}.fits", data=image.value, overwrite=True) - #fits.writeto(f"kernel_{obj.meta['wave_min'].value}.fits", kernel, overwrite=True) - #fits.writeto(f"newimage_{obj.meta['wave_min'].value}.fits", new_image, overwrite=True) elif image.ndim == 3 and kernel.ndim == 2: kernel = kernel[None, :, :] bkg_level = bkg_level[:, None, None] From 256434b709ed95c2aa53f0fcd381513524b56886 Mon Sep 17 00:00:00 2001 From: oczoske Date: Thu, 2 Apr 2026 11:34:46 +0200 Subject: [PATCH 7/7] Typo in file name... --- ...opsim_templates_intro.ipynb => scopesim_templates_intro.ipynb} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename docs/source/5_liners/{scopsim_templates_intro.ipynb => scopesim_templates_intro.ipynb} (100%) diff --git a/docs/source/5_liners/scopsim_templates_intro.ipynb b/docs/source/5_liners/scopesim_templates_intro.ipynb similarity index 100% rename from docs/source/5_liners/scopsim_templates_intro.ipynb rename to docs/source/5_liners/scopesim_templates_intro.ipynb