Note

This page is generated from the Jupyter notebook notebooks/04_psf_and_image_snr.ipynb in the repository. You can download it and run it interactively.

04 — PSF simulator and PSF-aware SNR

This notebook tours the detector-grid PSF simulator and the PSF-aware SNR features of wcc_etc:

  • PSF sources — AiryPSF (diffraction limited), DefocusPSF (+1 / +2 wave from Zemax Huygens data), and CustomPSF (any image).

  • ``ImageSimulator`` — renders a realistic detector image: source distributed by the PSF, sky + dark added per pixel, Poisson + read noise, and a saturation mask.

  • ``Simulation.get_image_snr`` — an aperture SNR computed from a rendered PSF, so it works for any PSF, not just the analytic Airy disk that get_snr assumes. For the in-focus case it reproduces get_snr; for defocus it quantifies the SNR penalty and the optimal aperture.

Everything is in electrons.

Setup

Imports, a G5V r=15 scene with zodiacal background, and a Simulation on the Sony sensor in the r band.

[ ]:

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.colors import LogNorm

import wcc_etc

wcc_etc.set_wcc_style()

np.random.seed(0)

scene = wcc_etc.get_scene(
    name="G5V",
    mag=15,
    background="zodi",
    bandpass="johnson_r",
    background_prop={"bandpass": "johnson_r", "mag": 22.5},
)

sim = wcc_etc.Simulation.from_sensor_and_scene("sony:r", scene)
print("sensor          :", "sony:r")
print("plate scale     :", round(sim.meta["sensor"]["pixel_size"], 3), "um pixels")
print("default aperture:", sim.meta["r_aper_mas"], "mas")
print("SNR @60s, r=15:", round(float(sim.get_snr(60)["snr"]), 2))

2. Full image simulation

ImageSimulator.simulate(time, psf=...) distributes the source electrons over the PSF, adds the per-pixel sky background and dark current, then Poisson + Gaussian read noise. It returns a SimulatedImage with image_clean (noise-free), image_e (noisy electrons), a saturation_mask, and converters to_adu() / to_fitsimg(). Here: a G5V r=15 star, 30 s, in-focus Airy PSF.

[ ]:

res = imsim.simulate(time=30, psf=wcc_etc.AiryPSF(), add_noise=True, seed=0)

# SimulatedImage.plot_image_row gives the standard 3-panel view:
# PSF + noise, PSF (no noise), and the saturation mask, on a shared color scale.
fig, axes = res.plot_image_row(stretch="log", units="mas")
fig.suptitle("G5V r=15, 30 s, Sony, Airy PSF", y=1.02)
fig.tight_layout()
plt.show()

print(f"clean total electrons : {res.image_clean.sum():.3e}")
print(f"peak pixel (e-)       : {res.image_e.max():.1f}")
print(f"saturated pixels      : {res.saturation_mask.sum()}")

Radial profile and encircled energy

The same SimulatedImage exposes convenience plots for the azimuthally-averaged radial profile and the encircled-energy (EE) curve. Both take units='mas' or 'pix'; plot_radial marks the HWHM and plot_encircled_energy can mark a target EE fraction (here 80%).

[ ]:

fig, (axL, axR) = plt.subplots(1, 2, figsize=(11, 4))
res.plot_radial(ax=axL, units="mas", title="Radial Profile")
res.plot_encircled_energy(ax=axR, units="mas", ee_target=0.8, title="Encircled Energy")
fig.tight_layout()
plt.show()

Comparing PSFs across focus and filter

Azimuthally-averaged normalized PSF (sum = 1) and encircled energy for three configurations: in-focus (Sony r), +1 wave defocus (Sony r), and +2 wave defocus (Sony bb). Defocus spreads the light into a broad doughnut, so the radial peak drops (note the log scale) and the encircled energy needs a much larger radius to reach 80%.

[ ]:

configs = [
    ("In-focus, Sony r", "sony:r", wcc_etc.AiryPSF()),
    (
        "+1 wave defocus, Sony r",
        "sony:r",
        wcc_etc.DefocusPSF(wcc_etc.DEFOCUS_1WAVE_PATH),
    ),
    (
        "+2 wave defocus, Sony bb",
        "sony:bb",
        wcc_etc.DefocusPSF(wcc_etc.DEFOCUS_2WAVE_PATH),
    ),
]

fig, (axR, axE) = plt.subplots(1, 2, figsize=(12, 4.5))
for label, sensorlabel, psf in configs:
    isim = wcc_etc.ImageSimulator.from_sensor_and_scene(sensorlabel, scene, npix=300)
    ctx = isim._context()
    psf_img = psf.render(ctx)  # normalized PSF (sum = 1)
    ps = ctx.plate_scale_mas
    wcc_etc.plot_radial_mpl(
        image_clean=psf_img, pixel_scale_mas=ps, ax=axR, units="mas", show_hwhm=False
    )
    axR.lines[-1].set_label(label)
    wcc_etc.plot_encircled_energy_mpl(
        image_clean=psf_img, pixel_scale_mas=ps, ax=axE, units="mas"
    )
    axE.lines[-1].set_label(label)

axR.set_yscale("log")
axR.set_title("Azimuthally-Averaged PSF")
axE.set_title("Encircled Energy")
axR.legend(fontsize=9)
axE.legend(fontsize=9, loc="lower right")
fig.tight_layout()
plt.show()

Bokeh backend

Every plot has a bokeh variant via backend='bokeh'. In a notebook, call output_notebook() once and show(...) the returned figure. For the web portal you would instead pass return_='components' (script + div) or return_='html' to embed it.

[ ]:

from bokeh.io import output_notebook, show

output_notebook()
show(res.plot_image(backend="bokeh", units="mas", title="Airy PSF (Bokeh)"))

The defocused PSFs render the same way. Here the +1 and +2 wave defocus images are shown side by side as a bokeh row (return_='obj' figures combined with bokeh.layouts.row).

[ ]:

from bokeh.io import show
from bokeh.layouts import row

defocus_psfs = [
    ("+1 wave defocus", wcc_etc.DefocusPSF(wcc_etc.DEFOCUS_1WAVE_PATH)),
    ("+2 wave defocus", wcc_etc.DefocusPSF(wcc_etc.DEFOCUS_2WAVE_PATH)),
]
panels = []
for name, psf in defocus_psfs:
    sim_img = imsim.simulate(time=30, psf=psf, add_noise=True, seed=0)
    panels.append(
        sim_img.plot_image(
            backend="bokeh", units="mas", title=f"{name} (Bokeh)", return_="obj"
        )
    )
show(row(*panels))

3. Effect of pointing jitter

Jitter is applied as a Gaussian blur of the PSF. More jitter -> broader PSF -> lower peak fraction and a more spread-out image.

[ ]:

airy = wcc_etc.AiryPSF()
fig, axes = plt.subplots(1, 3, figsize=(13, 4.2))
for ax, jit in zip(axes, [0.0, 20.0, 50.0]):
    img = render(airy, jitter=jit, npix=96)
    im = ax.imshow(img, origin="lower", cmap="viridis")
    ax.set_title(f"jitter = {jit:.0f} mas\npeak fraction = {img.max():.3f}")
    ax.grid(False)  # 2D image: no gridlines
    fig.colorbar(im, ax=ax, fraction=0.046)
fig.suptitle("Airy PSF vs pointing jitter (Sony)", y=1.03)
fig.tight_layout()
plt.show()

4. PSF-aware SNR cross-check (rendered Airy vs analytic)

get_snr now runs the 2-D image simulation (it delegates to get_image_snr). With the default aperture and the diffraction-limited Airy PSF it should reproduce the legacy analytic formula — now the deprecated get_snr_airy (the small residual is 2-D pixelation vs the 1-D Airy curve). We assert they agree to within ~3%.

[ ]:

import warnings

sim.update(source__mag=15)
with warnings.catch_warnings():
    warnings.simplefilter("ignore", DeprecationWarning)
    etc_snr = float(sim.get_snr_airy(60).value)  # legacy analytic path
img_snr = sim.get_snr(60)["snr"]  # 2-D default (delegates to get_image_snr)
diff_pct = 100.0 * (img_snr - etc_snr) / etc_snr

print(f"get_snr_airy (analytic Airy) : {etc_snr:.3f}")
print(f"get_snr (rendered Airy, 2-D) : {img_snr:.3f}")
print(f"difference                   : {diff_pct:+.2f} %")

assert abs(diff_pct) < 3.0, f"cross-check off by {diff_pct:.2f}%"
print("OK: agree within 3%")

5. Defocus SNR penalty (same aperture)

At a fixed aperture, a defocused PSF spreads the same light over more noise pixels and loses less of it inside the aperture only at large radii — at the default aperture the enclosed fraction drops, so the SNR drops. Comparing the +2 wave defocus to the in-focus Airy at the same default aperture quantifies the penalty.

[ ]:

defocus2 = wcc_etc.DefocusPSF(wcc_etc.DEFOCUS_2WAVE_PATH)

airy_res = sim.get_image_snr(time=60, psf=wcc_etc.AiryPSF())
def2_res = sim.get_image_snr(time=60, psf=defocus2)
drop_pct = 100.0 * (def2_res["snr"] - airy_res["snr"]) / airy_res["snr"]

print(f"aperture radius (mas)       : {airy_res['r_aper_mas']:.1f} (same for both)")
print()
for label, r in [("in-focus Airy", airy_res), ("+2 wave defocus", def2_res)]:
    print(
        f"{label:16s}: SNR = {r['snr']:6.2f}  enclosed = {r['enclosed_fraction']:.3f}  "
        f"n_pix = {r['n_pix']}"
    )
print()
print(f"SNR drop from defocus       : {drop_pct:+.1f} %")

6. Aperture optimization

optimize=True finds the SNR-maximizing aperture radius. For the compact Airy the optimum is small; for the spread-out defocus PSF it moves outward (and recovers some, but not all, of the lost SNR). We compare the optimized result to the default-aperture result for each PSF.

[ ]:

for name, psf in [("in-focus Airy", wcc_etc.AiryPSF()), ("+2 wave defocus", defocus2)]:
    default = sim.get_image_snr(time=60, psf=psf)
    opt = sim.get_image_snr(time=60, psf=psf, optimize=True)
    print(f"{name}:")
    print(
        f"  default aperture : r = {default['r_aper_mas']:6.1f} mas  "
        f"SNR = {default['snr']:6.2f}  n_pix = {default['n_pix']}"
    )
    print(
        f"  optimized        : r = {opt['r_aper_mas']:6.1f} mas  "
        f"SNR = {opt['snr']:6.2f}  n_pix = {opt['n_pix']}  "
        f"enclosed = {opt['enclosed_fraction']:.3f}"
    )
    gain = 100.0 * (opt["snr"] - default["snr"]) / default["snr"]
    print(f"  SNR gain from optimizing aperture: {gain:+.1f} %")
    print()

SNR as a function of aperture radius for each PSF (r=15, 60 s); the marker shows the optimize=True choice — small for the compact Airy, larger for the defocus PSF.

[ ]:

radii = np.linspace(20, 800, 50)
psfs_opt = {
    "in-focus Airy": wcc_etc.AiryPSF(),
    "+1 wave defocus": wcc_etc.DefocusPSF(wcc_etc.DEFOCUS_1WAVE_PATH),
    "+2 wave defocus": defocus2,
}

fig, ax = plt.subplots(figsize=(7.6, 4.6))
for name, psf in psfs_opt.items():
    snr = [sim.get_image_snr(time=60, psf=psf, r_aper_mas=r)["snr"] for r in radii]
    (line,) = ax.plot(radii, snr, label=name)
    opt = sim.get_image_snr(time=60, psf=psf, optimize=True)
    ax.plot(opt["r_aper_mas"], opt["snr"], "o", color=line.get_color(), ms=8)
ax.set_xlabel("aperture radius [mas]")
ax.set_ylabel("SNR (r=15, 60 s)")
ax.set_title("SNR vs aperture radius (markers = optimized aperture)")
ax.legend()
ax.grid(lw=0.3, alpha=0.4)
plt.show()

Summary

  • AiryPSF, DefocusPSF, CustomPSF render normalized PSFs onto the detector grid (pixel-size aware, jitter aware).

  • ImageSimulator.from_sensor_and_scene(...).simulate(time, psf=...) produces a realistic noisy electron image driven by source magnitude, with a saturation mask; to_adu() / to_fitsimg() bridge to ADU and the photometry tools.

  • Simulation.get_image_snr computes an aperture SNR from a rendered PSF. For the in-focus Airy case it reproduces the analytic get_snr (cross-check within a few percent); for defocus it quantifies the SNR penalty and optimize=True finds the SNR-maximizing aperture.

  • get_snr itself (the analytic path) is unchanged.