#!/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()