wyj
3deebfb2fa
Add an optional CPU preview path that deposits each point-source image as a supersampled delta and resolves the whole frame with one global Moffat convolution plus an N x N box average, instead of splatting a per-event PSF. - optics: FastPsfAccumulator builds the pixel-area-integral kernel of the target Moffat at the supersampled scale (width N*alpha, same beta). The 1/N^2 box average then reproduces the final pixel-area integral, so the requested FWHM and beta are preserved without renormalisation. Deposits are per-cell atomic adds; resolve accumulates into the caller's HDR buffer. - frame: fast branch in frame_splat_catalog with one shared supersampled buffer and a single resolve per frame; the accumulator is reused across movie frames and built from the map dimensions on lens-map import. - main: --fast-mode, --fast-supersample N (1..8, default 2) and --fast-deposit nearest|bilinear (default nearest). CPU-only and rejected in the HIP/dummy backends; --psf-min-y still applies per event while --max-cache-psf-flux does not. - The deposition scheme was chosen by scripts/fast_mode_deposit_error.py: nearest keeps the PSF shape exactly with <= 0.5/N px position quantization; bilinear keeps the exact centroid but broadens FWHM and beta. Recorded in benchmarks/fast_mode_deposit_2026-09-18.md. - tests/test_frame.c covers fast nearest vs the direct evaluator at the snapped centre, flux conservation, bilinear centroid, min-Y discard, frame plumbing, and HDR accumulation onto a non-zero background. - benchmarks/fast_mode_cpu_2026-09-18.md records a ~10x speedup on the 2MASS galactic-centre field with small tone-mapped differences.
350 lines
14 KiB
Python
350 lines
14 KiB
Python
#!/usr/bin/env python3
|
|
"""Fast-mode point-source deposition error experiment.
|
|
|
|
Question this script answers: when a single star is deposited into a
|
|
supersampled buffer at its exact continuous position, does a single nearest
|
|
supersampled pixel lose enough HDR information that we must instead deposit
|
|
into the 4 adjacent pixels (bilinear)?
|
|
|
|
Model (must match the intended C implementation):
|
|
* Target PSF is the normalized circular Moffat
|
|
M(r; a, b) = (b-1)/(pi a^2) * (1 + r^2/a^2)^(-b)
|
|
with a = FWHM / (2 sqrt(2^(1/b) - 1)).
|
|
* Reference image = pixel-area integral of M at the true star position
|
|
(8-point Gauss-Legendre quadrature, as used by the direct evaluator).
|
|
* Fast image with supersample factor N:
|
|
- deposit the event delta into the N x N supersampled grid
|
|
(nearest = 1 cell, bilinear = 4 adjacent cells);
|
|
- convolve the whole ss buffer with the single global kernel
|
|
k[m] = integral over ss cell m of M(z; N*a, b) (4-point quadrature);
|
|
- downsample by averaging each N x N block (1/N^2), which together with
|
|
the unnormalised N-scaled kernel reproduces the final pixel integral.
|
|
|
|
Usage: python3 scripts/fast_mode_deposit_error.py [--phases K] [--grid W]
|
|
[--tail R] [--max-radius R]
|
|
"""
|
|
|
|
import argparse
|
|
from math import pi, sqrt, ceil, log
|
|
|
|
import numpy as np
|
|
|
|
GAUSS4_X = np.array([-0.8611363115940526, -0.3399810435848563,
|
|
0.3399810435848563, 0.8611363115940526])
|
|
GAUSS4_W = np.array([0.3478548451374539, 0.6521451548625461,
|
|
0.6521451548625461, 0.3478548451374539])
|
|
GAUSS8_X = np.array([-0.9602898564975363, -0.7966664774136267,
|
|
-0.5255324099163290, -0.1834346424956498,
|
|
0.1834346424956498, 0.5255324099163290,
|
|
0.7966664774136267, 0.9602898564975363])
|
|
GAUSS8_W = np.array([0.1012285362903763, 0.2223810344533745,
|
|
0.3137066458778873, 0.3626837833783620,
|
|
0.3626837833783620, 0.3137066458778873,
|
|
0.2223810344533745, 0.1012285362903763])
|
|
|
|
TAIL_ABS_HDR = 1e-6
|
|
BOUNDARY_HDR = 1e-7
|
|
|
|
|
|
def moffat_alpha(fwhm, beta):
|
|
return fwhm / (2.0 * sqrt(2.0 ** (1.0 / beta) - 1.0))
|
|
|
|
|
|
def support_radius(alpha, beta, scaled_flux=1.0, rel_tail=1e-8):
|
|
norm = scaled_flux * (beta - 1.0) / (pi * alpha * alpha)
|
|
rt = min(rel_tail, TAIL_ABS_HDR / scaled_flux)
|
|
tail = alpha * sqrt(rt ** (1.0 / (1.0 - beta)) - 1.0)
|
|
br = BOUNDARY_HDR / norm
|
|
bound = 0.0 if br >= 1.0 else alpha * sqrt(br ** (-1.0 / beta) - 1.0)
|
|
return max(tail, bound)
|
|
|
|
|
|
def pixel_integral_image(alpha, beta, width, height, xs, ys, nodes, weights):
|
|
"""Reference: Moffat pixel-area integral on a final-resolution grid."""
|
|
img = np.zeros((height, width), dtype=np.float64)
|
|
norm = (beta - 1.0) / (pi * alpha * alpha)
|
|
x = np.arange(width)[None, :]
|
|
y = np.arange(height)[:, None]
|
|
inv = 1.0 / (alpha * alpha)
|
|
for iy in range(len(nodes)):
|
|
sy = y + 0.5 + 0.5 * nodes[iy] - ys
|
|
for ix in range(len(nodes)):
|
|
sx = x + 0.5 + 0.5 * nodes[ix] - xs
|
|
img += (0.25 * weights[ix] * weights[iy] * norm *
|
|
(1.0 + (sx * sx + sy * sy) * inv) ** (-beta))
|
|
return img
|
|
|
|
|
|
def build_kernel(alpha, n, beta, radius, nodes, weights):
|
|
"""Global ss kernel indexed by offset m in [-radius, radius].
|
|
|
|
k[m] = integral over ss cell m of I(z/N; alpha, beta), where I is the
|
|
final-unit normalized Moffat. Its total weight is ~N^2, so averaging each
|
|
N x N downsampled block reproduces the final pixel-area integral exactly.
|
|
"""
|
|
offsets = np.arange(-radius, radius + 1)
|
|
k = np.zeros((2 * radius + 1, 2 * radius + 1), dtype=np.float64)
|
|
alpha_ss = n * alpha
|
|
norm = (beta - 1.0) / (pi * alpha * alpha)
|
|
inv = 1.0 / (alpha_ss * alpha_ss)
|
|
mx = offsets[None, :]
|
|
my = offsets[:, None]
|
|
for iy in range(len(nodes)):
|
|
sy = my + 0.5 * nodes[iy]
|
|
for ix in range(len(nodes)):
|
|
sx = mx + 0.5 * nodes[ix]
|
|
k += (0.25 * weights[ix] * weights[iy] * norm *
|
|
(1.0 + (sx * sx + sy * sy) * inv) ** (-beta))
|
|
return k
|
|
|
|
|
|
def place_kernel(out, kernel, radius, u0x, u0y, scale):
|
|
"""out[p] += scale * kernel[p - u0]; clip to the buffer."""
|
|
h, w = out.shape
|
|
x0, x1 = u0x - radius, u0x + radius
|
|
y0, y1 = u0y - radius, u0y + radius
|
|
kx0 = max(0, -x0)
|
|
ky0 = max(0, -y0)
|
|
kx1 = kernel.shape[1] - max(0, x1 - (w - 1))
|
|
ky1 = kernel.shape[0] - max(0, y1 - (h - 1))
|
|
ix0 = max(0, x0)
|
|
iy0 = max(0, y0)
|
|
out[iy0:iy0 + (ky1 - ky0), ix0:ix0 + (kx1 - kx0)] += \
|
|
scale * kernel[ky0:ky1, kx0:kx1]
|
|
|
|
|
|
def fast_image(width, height, n, xs, ys, kernel, radius, deposit):
|
|
ss_w, ss_h = n * width, n * height
|
|
out = np.zeros((ss_h, ss_w), dtype=np.float64)
|
|
zx, zy = n * xs, n * ys
|
|
ix = int(np.floor(zx))
|
|
iy = int(np.floor(zy))
|
|
if deposit == "nearest":
|
|
place_kernel(out, kernel, radius, ix, iy, 1.0)
|
|
else:
|
|
# Interpolate between the two ss cell centres bracketing the star.
|
|
bx = int(np.floor(zx - 0.5))
|
|
by = int(np.floor(zy - 0.5))
|
|
fx, fy = zx - 0.5 - bx, zy - 0.5 - by
|
|
for dy, wy in ((0, 1 - fy), (1, fy)):
|
|
for dx, wx in ((0, 1 - fx), (1, fx)):
|
|
if wx * wy != 0.0:
|
|
place_kernel(out, kernel, radius, bx + dx, by + dy,
|
|
wx * wy)
|
|
final = out.reshape(height, n, width, n).mean(axis=(1, 3))
|
|
return final
|
|
|
|
|
|
def bilinear_sample(img, x, y):
|
|
h, w = img.shape
|
|
x = min(max(x, 0.0), w - 1.0)
|
|
y = min(max(y, 0.0), h - 1.0)
|
|
x0, y0 = int(np.floor(x)), int(np.floor(y))
|
|
x1, y1 = min(x0 + 1, w - 1), min(y0 + 1, h - 1)
|
|
tx, ty = x - x0, y - y0
|
|
return ((1 - ty) * ((1 - tx) * img[y0, x0] + tx * img[y0, x1]) +
|
|
ty * ((1 - tx) * img[y1, x0] + tx * img[y1, x1]))
|
|
|
|
|
|
def image_centroid(img):
|
|
h, w = img.shape
|
|
xs = np.arange(w)[None, :] + 0.5
|
|
ys = np.arange(h)[:, None] + 0.5
|
|
total = img.sum()
|
|
if total <= 0.0:
|
|
return 0.5 * w, 0.5 * h
|
|
return float((img * xs).sum() / total), float((img * ys).sum() / total)
|
|
|
|
|
|
def radial_fwhm(img, center, step=0.02, r_max=None):
|
|
h, w = img.shape
|
|
cx, cy = center
|
|
if r_max is None:
|
|
r_max = 0.5 * min(w, h) - 1.0
|
|
peak = bilinear_sample(img, cx, cy)
|
|
if peak <= 0.0:
|
|
return float("nan"), float("nan")
|
|
radii = np.arange(0.0, r_max, step)
|
|
prof = np.array([bilinear_sample(img, cx + r, cy) for r in radii])
|
|
half = 0.5 * peak
|
|
idx = np.where(prof <= half)[0]
|
|
if len(idx) == 0:
|
|
return float("nan"), peak
|
|
i = idx[0]
|
|
if i == 0:
|
|
return float("nan"), peak
|
|
r0, r1 = radii[i - 1], radii[i]
|
|
v0, v1 = prof[i - 1], prof[i]
|
|
r_half = r0 + (half - v0) * (r1 - r0) / (v1 - v0)
|
|
return 2.0 * r_half, peak
|
|
|
|
|
|
def moffat_profile(r, amp, alpha, beta):
|
|
return amp * (1.0 + (r / alpha) ** 2) ** (-beta)
|
|
|
|
|
|
def fit_beta(img, center, alpha_hint):
|
|
from scipy.optimize import curve_fit
|
|
h, w = img.shape
|
|
cx, cy = center
|
|
r_max = 0.4 * min(w, h)
|
|
step = 0.05
|
|
radii = np.arange(1.5, r_max, step)
|
|
prof = np.array([bilinear_sample(img, cx + r, cy) for r in radii])
|
|
peak = bilinear_sample(img, cx, cy)
|
|
good = prof > 1e-6 * peak
|
|
if good.sum() < 8:
|
|
return float("nan")
|
|
try:
|
|
popt, _ = curve_fit(moffat_profile, radii[good], prof[good],
|
|
p0=[peak, alpha_hint, 4.0],
|
|
bounds=([0.0, 0.2 * alpha_hint, 1.05],
|
|
[np.inf, 5.0 * alpha_hint, 20.0]),
|
|
maxfev=20000)
|
|
return float(popt[2])
|
|
except Exception:
|
|
return float("nan")
|
|
|
|
|
|
def phase_errors(ref, fast, peak):
|
|
diff = fast - ref
|
|
return (float(np.abs(diff).max() / peak),
|
|
float(np.sqrt(np.mean(diff * diff)) / peak))
|
|
|
|
|
|
def run_config(fwhm, beta, n, phases, grid, nodes_k, weights_k, tail,
|
|
max_radius, deposit):
|
|
alpha = moffat_alpha(fwhm, beta)
|
|
alpha_ss = n * alpha
|
|
r_full = support_radius(alpha_ss, beta, 1.0, tail)
|
|
radius = int(min(ceil(r_full) + 1, max_radius))
|
|
kernel = build_kernel(alpha, n, beta, radius, GAUSS4_X, GAUSS4_W)
|
|
retained = float(kernel.sum() / (n * n))
|
|
|
|
dxs = (np.arange(phases) + 0.5) / phases
|
|
max_err = rms_err = 0.0
|
|
central_err = 0.0
|
|
max_centroid = 0.0
|
|
fwhm_err = 0.0
|
|
beta_bias = None
|
|
ref_fwhm = 0.0
|
|
for dx in dxs:
|
|
for dy in dxs:
|
|
xs = grid / 2 + 0.5 + dx
|
|
ys = grid / 2 + 0.5 + dy
|
|
ref = pixel_integral_image(alpha, beta, grid, grid, xs, ys,
|
|
GAUSS8_X, GAUSS8_W)
|
|
peak = ref.max()
|
|
fast = fast_image(grid, grid, n, xs, ys, kernel, radius, deposit)
|
|
e_max, e_rms = phase_errors(ref, fast, peak)
|
|
max_err = max(max_err, e_max)
|
|
rms_err = max(rms_err, e_rms)
|
|
# HDR error of the star's own core pixels (central 3x3).
|
|
ci = int(np.clip(grid / 2, 0, grid - 3))
|
|
central_err = max(central_err,
|
|
float(np.abs(fast[ci:ci + 3, ci:ci + 3] -
|
|
ref[ci:ci + 3, ci:ci + 3]).max() /
|
|
peak))
|
|
cref = image_centroid(ref)
|
|
cfast = image_centroid(fast)
|
|
max_centroid = max(max_centroid,
|
|
sqrt((cref[0] - cfast[0]) ** 2 +
|
|
(cref[1] - cfast[1]) ** 2))
|
|
# Shape error is isolated by comparing both profiles about the same
|
|
# centre: the fast image's own centroid. Bilinear already keeps
|
|
# that centre at the true position, so its reference is ref.
|
|
if deposit == "bilinear":
|
|
ref_same = ref
|
|
else:
|
|
ref_same = pixel_integral_image(alpha, beta, grid, grid,
|
|
cfast[0], cfast[1],
|
|
GAUSS8_X, GAUSS8_W)
|
|
f_ref, _ = radial_fwhm(ref_same, cfast)
|
|
f_fast, _ = radial_fwhm(fast, cfast)
|
|
if np.isfinite(f_ref) and np.isfinite(f_fast):
|
|
ref_fwhm = f_ref
|
|
fwhm_err = max(fwhm_err, abs(f_fast - f_ref) / f_ref)
|
|
if fwhm >= 2.0:
|
|
# fit_beta() has its own pixel-integration bias; compare the
|
|
# fast fit with the reference fit computed the same way so
|
|
# beterr isolates the deposition-induced shape change.
|
|
b_ref = fit_beta(ref_same, cfast, alpha)
|
|
b_fast = fit_beta(fast, cfast, alpha)
|
|
if np.isfinite(b_ref) and np.isfinite(b_fast) and b_ref > 0.0:
|
|
error = abs(b_fast - b_ref) / b_ref
|
|
beta_bias = error if beta_bias is None else max(beta_bias, error)
|
|
return {
|
|
"fwhm": fwhm, "beta": beta, "n": n, "deposit": deposit,
|
|
"alpha": alpha, "radius": radius, "retained": retained,
|
|
"max_err": max_err, "rms_err": rms_err, "central_err": central_err,
|
|
"centroid": max_centroid, "ref_fwhm": ref_fwhm,
|
|
"fwhm_err": fwhm_err,
|
|
"beta_bias": float("nan") if beta_bias is None else beta_bias,
|
|
}
|
|
|
|
|
|
def main():
|
|
ap = argparse.ArgumentParser()
|
|
ap.add_argument("--phases", type=int, default=24)
|
|
ap.add_argument("--grid", type=int, default=48)
|
|
ap.add_argument("--tail", type=float, default=1e-6)
|
|
ap.add_argument("--max-radius", type=int, default=64)
|
|
ap.add_argument("--deposits", default="nearest,bilinear")
|
|
args = ap.parse_args()
|
|
|
|
configs = [(1.0, 8.0), (2.7, 8.0), (6.0, 8.0),
|
|
(1.0, 4.5), (2.7, 4.5), (6.0, 4.5),
|
|
(1.0, 3.5), (2.7, 3.5)]
|
|
ns = [2, 3, 4]
|
|
deposits = [d for d in args.deposits.split(",") if d]
|
|
|
|
print("Fast-mode deposition error experiment")
|
|
print(f"grid={args.grid} phases={args.phases} tail={args.tail:g} "
|
|
f"max_radius={args.max_radius}")
|
|
print("kernel quadrature = 4-point Gauss; reference = 8-point Gauss")
|
|
print()
|
|
header = (f"{'FWHM':>5} {'beta':>5} {'N':>2} {'deposit':>9} {'R_ss':>5} "
|
|
f"{'retain':>8} {'maxErr':>9} {'rmsErr':>9} {'coreErr':>9} "
|
|
f"{'centErr':>9} {'FWHMerr':>9} {'beterr':>9}")
|
|
print(header)
|
|
print("-" * len(header))
|
|
|
|
rows = []
|
|
for fwhm, beta in configs:
|
|
for n in ns:
|
|
for deposit in deposits:
|
|
r = run_config(fwhm, beta, n, args.phases, args.grid,
|
|
GAUSS4_X, GAUSS4_W, args.tail, args.max_radius,
|
|
deposit)
|
|
rows.append(r)
|
|
print(f"{fwhm:5.2f} {beta:5.2f} {n:2d} {deposit:>9} "
|
|
f"{r['radius']:5d} {r['retained']:8.5f} "
|
|
f"{r['max_err']:9.3e} {r['rms_err']:9.3e} "
|
|
f"{r['central_err']:9.3e} {r['centroid']:9.3e} "
|
|
f"{r['fwhm_err']:9.3e} {r['beta_bias']:9.3e}")
|
|
print()
|
|
|
|
print("Worst case over all tested configurations, by deposit scheme:")
|
|
for deposit in deposits:
|
|
sel = [r for r in rows if r["deposit"] == deposit]
|
|
if not sel:
|
|
continue
|
|
print(f" {deposit:>9}: maxErr={max(r['max_err'] for r in sel):.3e} "
|
|
f"rmsErr={max(r['rms_err'] for r in sel):.3e} "
|
|
f"coreErr={max(r['central_err'] for r in sel):.3e} "
|
|
f"centErr={max(r['centroid'] for r in sel):.3e} "
|
|
f"FWHMerr={max(r['fwhm_err'] for r in sel):.3e} "
|
|
f"beterr={np.nanmax([r['beta_bias'] for r in sel]):.3e}")
|
|
print()
|
|
print("FWHMerr compares the fast image with a direct image at the same "
|
|
"centroid, so it is a pure shape error (0 for exact position "
|
|
"deposition at its snapped centre). centErr is the position "
|
|
"quantization vs the true continuous star position. beterr is the "
|
|
"deposition-induced Moffat beta change, isolated from the radial "
|
|
"fit's own pixel-integration bias by comparing with a reference fit; "
|
|
"reported only for FWHM >= 2 px (nan otherwise).")
|
|
|
|
|
|
if __name__ == "__main__":
|
|
main()
|