diff --git a/CHANGELOG.md b/CHANGELOG.md index 4552fff..91eaaf1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,17 @@ to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ## [Unreleased] +### Fixed + +- **`MoffatPSF`, `EllipticalGaussianPSF` and `AiryPSF` no longer put light that + falls off the frame back onto it.** They normalised their stamp *after* + clipping it to the frame, so a source on the edge column deposited its full + flux instead of the part that lands on the detector (a Moffat star centred on + column 0 kept 100% of its light, where `GaussianPSF` correctly keeps about + half). The stamp is now normalised before clipping, which follows the AO stack + convention that light lost at a detector edge is lost, not renormalised + (aocore CONVENTIONS 3.3). Sources wholly inside the frame are unchanged. + ### Changed - **Generic primitives now come from `aocore`.** getframes is part of an diff --git a/docs/guides/scenes.md b/docs/guides/scenes.md index cbf5d2f..cd4f4a0 100644 --- a/docs/guides/scenes.md +++ b/docs/guides/scenes.md @@ -62,8 +62,9 @@ the camera a `qe_curve`. - [`MoffatPSF(fwhm_arcsec, beta)`][getframes.scene.psf.MoffatPSF] --- a better match to seeing-limited stars; smaller `beta` means broader wings. -Both place a flux-conserving stamp at the source's sub-pixel position. Custom PSFs -can subclass [`PSF`][getframes.scene.psf.PSF] and implement `add_source`. +Both place a flux-conserving stamp at the source's sub-pixel position. Light that +falls beyond the frame edge is lost, as it would be on a real detector; it is +never renormalised back onto the pixels that remain. Custom PSFs can subclass [`PSF`][getframes.scene.psf.PSF] and implement `add_source`. ## Working with the photon map directly diff --git a/src/getframes/scene/psf.py b/src/getframes/scene/psf.py index 98c390c..b7c1cb6 100644 --- a/src/getframes/scene/psf.py +++ b/src/getframes/scene/psf.py @@ -11,6 +11,7 @@ from __future__ import annotations import math +from collections.abc import Callable from dataclasses import dataclass import numpy as np @@ -38,6 +39,33 @@ def _stamp_bounds( return x0, x1, y0, y1 +def _deposit_sampled( + image: NDArray[np.float64], + x: float, + y: float, + radius: int, + flux: float, + profile_at: Callable[[NDArray[np.float64], NDArray[np.float64]], NDArray[np.float64]], +) -> None: + """Add a sampled, normalised stamp, losing whatever falls off the frame. + + ``profile_at(dx, dy)`` evaluates the PSF at column/row offsets from the source. + The profile is normalised over the *whole* stamp before it is clipped to the + frame, so light off the edge is lost rather than renormalised back in. + """ + x0, x1, y0, y1 = _stamp_bounds(x, y, radius, image.shape) + if x0 >= x1 or y0 >= y1: + return # source falls entirely off the frame + ix, iy = round(x), round(y) + offsets = np.arange(-radius, radius + 1, dtype=np.float64) + profile = profile_at(offsets[None, :] + (ix - x), offsets[:, None] + (iy - y)) + total = profile.sum() + if total > 0: + sx, sy = x0 - (ix - radius), y0 - (iy - radius) + stamp = profile[sy : sy + (y1 - y0), sx : sx + (x1 - x0)] + image[y0:y1, x0:x1] += flux * stamp / total + + class PSF: """Base class for point-spread functions.""" @@ -194,17 +222,12 @@ def add_source( alpha = fwhm_pix / (2.0 * np.sqrt(2.0 ** (1.0 / self.beta) - 1.0)) radius = int(np.ceil(6.0 * alpha)) + 1 - x0, x1, y0, y1 = _stamp_bounds(x, y, radius, image.shape) - if x0 >= x1 or y0 >= y1: - return + beta = self.beta - xs = np.arange(x0, x1) - x - ys = np.arange(y0, y1) - y - rr = xs[None, :] ** 2 + ys[:, None] ** 2 - profile = (1.0 + rr / alpha**2) ** (-self.beta) - total = profile.sum() - if total > 0: - image[y0:y1, x0:x1] += flux * profile / total + def profile(dx: NDArray[np.float64], dy: NDArray[np.float64]) -> NDArray[np.float64]: + return np.asarray((1.0 + (dx**2 + dy**2) / alpha**2) ** (-beta)) + + _deposit_sampled(image, x, y, radius, flux, profile) @dataclass(frozen=True) @@ -239,20 +262,15 @@ def add_source( raise ValueError("PSF FWHM and plate scale must be positive.") radius = int(np.ceil(5.0 * sigma_major)) + 1 - x0, x1, y0, y1 = _stamp_bounds(x, y, radius, image.shape) - if x0 >= x1 or y0 >= y1: - return - - xs = np.arange(x0, x1) - x - ys = np.arange(y0, y1) - y theta = math.radians(self.position_angle_deg) cos_t, sin_t = math.cos(theta), math.sin(theta) - u = xs[None, :] * cos_t + ys[:, None] * sin_t - v = -xs[None, :] * sin_t + ys[:, None] * cos_t - profile = np.exp(-0.5 * ((u / sigma_major) ** 2 + (v / sigma_minor) ** 2)) - total = profile.sum() - if total > 0: - image[y0:y1, x0:x1] += flux * profile / total + + def profile(dx: NDArray[np.float64], dy: NDArray[np.float64]) -> NDArray[np.float64]: + u = dx * cos_t + dy * sin_t + v = -dx * sin_t + dy * cos_t + return np.exp(-0.5 * ((u / sigma_major) ** 2 + (v / sigma_minor) ** 2)) + + _deposit_sampled(image, x, y, radius, flux, profile) @dataclass(frozen=True) @@ -299,18 +317,12 @@ def add_source( # First null at 1.22 lambda / D; size the stamp to a few Airy rings. first_null_pix = 1.22 / (arg_per_pixel / math.pi) if arg_per_pixel > 0 else 1.0 radius = int(np.ceil(5.0 * first_null_pix)) + 1 - x0, x1, y0, y1 = _stamp_bounds(x, y, radius, image.shape) - if x0 >= x1 or y0 >= y1: - return + obstruction = self.obstruction + + def profile(dx: NDArray[np.float64], dy: NDArray[np.float64]) -> NDArray[np.float64]: + return _airy_intensity(arg_per_pixel * np.hypot(dx, dy), obstruction) - xs = np.arange(x0, x1) - x - ys = np.arange(y0, y1) - y - rr = np.sqrt(xs[None, :] ** 2 + ys[:, None] ** 2) - arg = arg_per_pixel * rr - profile = _airy_intensity(arg, self.obstruction) - total = profile.sum() - if total > 0: - image[y0:y1, x0:x1] += flux * profile / total + _deposit_sampled(image, x, y, radius, flux, profile) def _airy_intensity(arg: NDArray[np.float64], obstruction: float) -> NDArray[np.float64]: diff --git a/tests/test_scene.py b/tests/test_scene.py index f585a71..8894fd7 100644 --- a/tests/test_scene.py +++ b/tests/test_scene.py @@ -61,6 +61,29 @@ def test_psf_conserves_flux(psf, optics): assert img.sum() == pytest.approx(1000.0, rel=1e-3) +@pytest.mark.parametrize( + "psf", + [ + gf.GaussianPSF(fwhm_arcsec=1.2), + gf.MoffatPSF(fwhm_arcsec=1.2), + gf.EllipticalGaussianPSF(fwhm_major_arcsec=1.6, fwhm_minor_arcsec=0.8), + # lambda / D = 2 arcsec, five pixels. + gf.AiryPSF(aperture_diameter_m=0.1, wavelength_m=1.0e-6), + ], + ids=lambda psf: type(psf).__name__, +) +def test_psf_loses_the_flux_that_falls_off_the_frame(psf): + # Light beyond the detector edge is lost, never renormalised back in + # (aocore CONVENTIONS 3.3). A source on column 0 of a frame must deposit + # exactly the part of a fully captured image that lies at x >= 0. + full = np.zeros((128, 128)) + psf.add_source(full, x=64.3, y=64.0, flux=1000.0, plate_scale_arcsec_per_pixel=0.4) + edge = np.zeros((128, 64)) + psf.add_source(edge, x=0.3, y=64.0, flux=1000.0, plate_scale_arcsec_per_pixel=0.4) + assert edge.sum() == pytest.approx(full[:, 64:].sum(), rel=1e-9) + assert edge.sum() < 0.75 * full.sum() + + def test_psf_peaks_at_source_position(optics): img = np.zeros((64, 64)) gf.GaussianPSF(1.0).add_source(