From c04bfc640b5214b1125e7dafc93c1155122030ed Mon Sep 17 00:00:00 2001 From: oczoske Date: Thu, 19 Feb 2026 13:25:56 +0100 Subject: [PATCH 1/3] 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/3] 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/3] 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