"""Render the WCC PSF (in-focus Airy or 1/2-wave defocus) on an oversampled grid.
wcc_etc PSFSource.render(ctx) returns a detector-sampled PSF; we get an
oversampled one by describing a finer detector to the context (pixel size and
plate scale divided by `oversample`, npix multiplied by it). Jitter blur is
applied by wcc_etc at the fine plate scale.
"""
import numpy as np
from scipy.ndimage import shift as _nd_shift
from wcc_etc import DEFOCUS_1WAVE_PATH, DEFOCUS_2WAVE_PATH, AiryPSF, DefocusPSF
from wcc_etc.psfsim import DetectorPSFContext
# Detector-pixel stamp sizes (odd). Defocus Huygens data spans ~272 IMX px.
# The in-focus Airy PSF needs ~129 px to hold >=99.5% of the (jitter-blurred)
# energy; a 65 px stamp truncates ~1.3% of the wing energy.
DEFAULT_STAMP = {0: 129, 1: 257, 2: 257}
[docs]
def make_psf_source(focus):
if focus == 0:
return AiryPSF()
if focus == 1:
return DefocusPSF(DEFOCUS_1WAVE_PATH)
if focus == 2:
return DefocusPSF(DEFOCUS_2WAVE_PATH)
raise ValueError(f"focus must be 0, 1, or 2 (waves of defocus), got {focus!r}")
[docs]
def centroid_offset(psf):
"""(dx, dy) of the flux centroid from the grid centre (n-1)/2, in pixels."""
n = psf.shape[0]
i = np.arange(n, dtype=float)
total = psf.sum()
cx = (psf.sum(axis=0) @ i) / total
cy = (psf.sum(axis=1) @ i) / total
c = (n - 1) / 2.0
return float(cx - c), float(cy - c)
[docs]
def recenter_on_centroid(psf):
"""Shift `psf` so its flux centroid sits on the grid centre; renormalize.
wcc_etc renders the Zemax Huygens defocus images about the chief ray
(its fix for issue #18). The flux centroid of a defocused PSF still sits
a little off the chief ray because the ring is asymmetric: +0.10 px for
1 wave and +0.19 px for 2 waves along +Y on zwo. A centroider recovers
the centroid, not the chief ray, so recentring on it makes the catalog
position the quantity a centroider recovers (issue #18).
"""
dx, dy = centroid_offset(psf)
out = _nd_shift(psf, (-dy, -dx), order=1, mode="constant", cval=0.0)
return out / out.sum()
[docs]
def render_oversampled_psf(
sim, focus, oversample=11, stamp_npix=None, jitter_sigma_mas=None,
wavelength_m=None,
):
"""Normalized PSF on a (stamp_npix*oversample)^2 fine grid, centered.
`wavelength_m` overrides the sensor central wavelength (Airy path only;
DefocusPSF is a fixed Huygens image and ignores wavelength). Defocus
PSFs are recentred on their flux centroid (see recenter_on_centroid).
"""
source = make_psf_source(focus)
if stamp_npix is None:
stamp_npix = DEFAULT_STAMP[focus]
if stamp_npix % 2 == 0:
raise ValueError(f"stamp_npix must be odd, got {stamp_npix}")
sensor, telescope = sim.sensor, sim.telescope
if jitter_sigma_mas is None:
jitter_sigma_mas = float(telescope.jitter_sigma.to("mas").value)
plate_mas = (
sensor.get_plate_scale(telescope).to("arcsec/pix").value * 1000.0
)
n_fine = stamp_npix * oversample
ctx = DetectorPSFContext(
npix=n_fine,
pixel_size_um=float(sensor.pixel_size.value) / oversample,
plate_scale_mas=plate_mas / oversample,
wavelength_m=(
float(sensor.wavelength.to("m").value)
if wavelength_m is None
else float(wavelength_m)
),
diameter_m=float(telescope.diameter_primary.to("m").value),
fnum=float(telescope.f_num),
jitter_sigma_mas=jitter_sigma_mas,
oversample=3, # internal supersampling of the fine grid (Airy path)
)
psf = source.render(ctx)
if focus != 0:
psf = recenter_on_centroid(psf)
return psf