Lose the PSF light that falls off the frame

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

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
 

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]:

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(