diff --git a/README.md b/README.md index d52eb08..621e582 100644 --- a/README.md +++ b/README.md @@ -132,8 +132,10 @@ PNGs. The command above omits the optional HDR and lens-map exports. This example infers a position facing the hole from its look direction and radius. Adaptive refinement should be configured for the desired image accuracy; it is disabled by default. Camera controls, movie sequences, lens-map reuse, PSF settings, and HDR output -are described in [usage.md](usage.md). Both binaries provide a complete option -list with `--help`. +are described in [usage.md](usage.md). For fast previews, `--fast-mode` replaces +per-event PSF splats with a supersampled delta deposit plus one global PSF +convolution and downsample. Both binaries provide a complete option list with +`--help`. ### Example: Looking outward just above the Schwarzschild horizon diff --git a/benchmarks/fast_mode_cpu_2026-09-18.md b/benchmarks/fast_mode_cpu_2026-09-18.md new file mode 100644 index 0000000..6f2927b --- /dev/null +++ b/benchmarks/fast_mode_cpu_2026-09-18.md @@ -0,0 +1,108 @@ +# Fast-mode CPU benchmark (2026-09-18) + +Compares the default per-event phase-cache PSF splat with `--fast-mode` on a +7.4-degree 2MASS galactic-center field. Commands, raw terminal output, and +wall-clock totals are preserved below. + +## Environment + +- Git revision: `7ddc585` plus the uncommitted fast-mode work. +- Build: `make -j4 PSF_BACKEND=cpu SPACETIME=minkowski backend` + (`-O2 -DNDEBUG -fopenmp`, libpng). +- Host: 12th Gen Intel Core i7-12700K, 16 hardware threads. + +## Commands + +All runs use the same catalog, geometry, and exposure; only the accumulation +mode changes. `--all-sky-catalog assets/2mass/processed/all_sky`, +`--width 320 --height 240 --fov-deg 10 --look-ra-deg 266 --look-dec-deg -29 +--exposure 1e13 --coarse-cell-pixels 16 --refine-max-level 0`. + +```sh +./build/Release/minkowski_sky --all-sky-catalog assets/2mass/processed/all_sky \ + --width 320 --height 240 --fov-deg 10 --look-ra-deg 266 --look-dec-deg -29 \ + --exposure 1e13 --coarse-cell-pixels 16 --refine-max-level 0 \ + --output /tmp/fastmode/bench/normal.png + +./build/Release/minkowski_sky --all-sky-catalog assets/2mass/processed/all_sky \ + --width 320 --height 240 --fov-deg 10 --look-ra-deg 266 --look-dec-deg -29 \ + --exposure 1e13 --coarse-cell-pixels 16 --refine-max-level 0 \ + --fast-mode --fast-supersample 2 \ + --output /tmp/fastmode/bench/fast_n2.png + +./build/Release/minkowski_sky --all-sky-catalog assets/2mass/processed/all_sky \ + --width 320 --height 240 --fov-deg 10 --look-ra-deg 266 --look-dec-deg -29 \ + --exposure 1e13 --coarse-cell-pixels 16 --refine-max-level 0 \ + --fast-mode --fast-supersample 4 \ + --output /tmp/fastmode/bench/fast_n4.png + +./build/Release/minkowski_sky --all-sky-catalog assets/2mass/processed/all_sky \ + --width 320 --height 240 --fov-deg 10 --look-ra-deg 266 --look-dec-deg -29 \ + --exposure 1e13 --coarse-cell-pixels 16 --refine-max-level 0 \ + --fast-mode --fast-supersample 2 --fast-deposit bilinear \ + --output /tmp/fastmode/bench/fast_bil_n2.png +``` + +## Wall-clock totals + +Measured with `date +%s.%N` around each command. + +| run | wall time | images | +| --- | --- | --- | +| normal (phase cache) | 16.645 s | 5,912,840 | +| fast nearest N=2 | 1.663 s | 5,912,840 | +| fast bilinear N=2 | 1.648 s | 5,912,840 | +| fast nearest N=4 | 12.229 s | 5,912,840 | + +Fast nearest N=2 is about 10x faster than the phase-cache path on this +PSF-dominated field. N=4 increases the supersampled buffer and kernel area by +4x each and becomes convolution-bound. + +## Raw terminal output + +```text +===== normal.log ===== +Blackbody LUT: 1024 CIE 1931 2-deg 1 nm-linear XYZ nodes, T=[670.146556, 101408.88] K, loaded assets/blackbody/cie1931_2deg_xyz_1024.grbblut in 0.000056 s; payload fnv1a64=0e39867b0e70a809 +Blackbody backend: lut +PSF cache ready: 64x64 phases, radius 47 px, relative tail 1e-08, tail abs 1e-06, boundary 1e-07, build 0.725 s +Catalog prefetch: 96 requested, 96 newly loaded (6618807 stars), 0 unavailable in 0.696 s; 4 loader workers +Rendered 5912840 images from 6618807 catalog stars to /tmp/fastmode/bench/normal.png (ok) +PSF splats: cached 5912840, cached wing-clipped 0, direct fallbacks 0, discarded below min-Y 0 +===== fast_n2.log ===== +Fast mode is a preview approximation; --max-cache-psf-flux is ignored and one global kernel is used for every event. +Blackbody LUT: 1024 CIE 1931 2-deg 1 nm-linear XYZ nodes, T=[670.146556, 101408.88] K, loaded assets/blackbody/cie1931_2deg_xyz_1024.grbblut in 0.000045 s; payload fnv1a64=0e39867b0e70a809 +Blackbody backend: lut +Fast PSF: supersample 2x, deposit nearest, kernel radius 93 ss px (46.50 final px), retained flux 1.000000, cached 4-point quadrature +Catalog prefetch: 96 requested, 96 newly loaded (6618807 stars), 0 unavailable in 0.708 s; 4 loader workers +Rendered 5912840 images from 6618807 catalog stars to /tmp/fastmode/bench/fast_n2.png (ok) +Fast PSF splats: deposited 5912840, wing-clipped 0, discarded below min-Y 0 +===== fast_n4.log ===== +Fast mode is a preview approximation; --max-cache-psf-flux is ignored and one global kernel is used for every event. +Blackbody LUT: 1024 CIE 1931 2-deg 1 nm-linear XYZ nodes, T=[670.146556, 101408.88] K, loaded assets/blackbody/cie1931_2deg_xyz_1024.grbblut in 0.000058 s; payload fnv1a64=0e39867b0e70a809 +Blackbody backend: lut +Fast PSF: supersample 4x, deposit nearest, kernel radius 185 ss px (46.25 final px), retained flux 1.000000, cached 4-point quadrature +Catalog prefetch: 96 requested, 96 newly loaded (6618807 stars), 0 unavailable in 0.709 s; 4 loader workers +Rendered 5912840 images from 6618807 catalog stars to /tmp/fastmode/bench/fast_n4.png (ok) +Fast PSF splats: deposited 5912840, wing-clipped 0, discarded below min-Y 0 +===== fast_bil_n2.log ===== +Fast mode is a preview approximation; --max-cache-psf-flux is ignored and one global kernel is used for every event. +Blackbody LUT: 1024 CIE 1931 2-deg 1 nm-linear XYZ nodes, T=[670.146556, 101408.88] K, loaded assets/blackbody/cie1931_2deg_xyz_1024.grbblut in 0.000047 s; payload fnv1a64=0e39867b0e70a809 +Blackbody backend: lut +Fast PSF: supersample 2x, deposit bilinear, kernel radius 93 ss px (46.50 final px), retained flux 1.000000, cached 4-point quadrature +Catalog prefetch: 96 requested, 96 newly loaded (6618807 stars), 0 unavailable in 0.712 s; 4 loader workers +Rendered 5912840 images from 6618807 catalog stars to /tmp/fastmode/bench/fast_bil_n2.png (ok) +Fast PSF splats: deposited 5912840, wing-clipped 0, discarded below min-Y 0 +``` + +## Tone-mapped difference vs the normal render + +```text +fast_n2 mean=1.1213 rms=1.7255 max=14 px>16=0/76800 +fast_n4 mean=0.5516 rms=0.9146 max=8 px>16=0/76800 +fast_bil_n2 mean=0.2091 rms=0.4793 max=16 px>16=0/76800 +``` + +`mean`/`rms` are per-channel 8-bit differences; `px>16` counts pixels with any +channel more than 16 levels from the normal render. No pixel exceeds 16 in any +run. These are display-scale differences from the measured deposition tradeoff; +they are not a bit-exact regression baseline. diff --git a/benchmarks/fast_mode_deposit_2026-09-18.md b/benchmarks/fast_mode_deposit_2026-09-18.md new file mode 100644 index 0000000..f33232f --- /dev/null +++ b/benchmarks/fast_mode_deposit_2026-09-18.md @@ -0,0 +1,128 @@ +# Fast-mode deposition error experiment (2026-09-18) + +Decision this note settles: whether `--fast-mode` should deposit each point +source image as a single nearest supersampled pixel or as 4 adjacent +(bilinear) pixels before the one global PSF convolution and downscale. + +## Command + +```sh +python3 scripts/fast_mode_deposit_error.py --phases 12 --grid 48 --tail 1e-6 --max-radius 64 +``` + +## Model + +Reference = pixel-area integral of the normalized circular Moffat +`M(r; α, β)`, `α = FWHM / (2 sqrt(2^(1/β) - 1))`, on the final grid, using +8-point Gauss-Legendre quadrature. Fast image with supersample factor `N`: + +- deposit the event 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] = ∫_ss-cell m M(z; Nα, β) dz` (4-point quadrature); +- downscale by averaging each `N x N` block, which together with the + unnormalised `N`-scaled kernel reproduces the final pixel-area integral. + +Metrics: `maxErr`/`rmsErr` are whole-image errors vs the direct image at the +true position; `coreErr` is the worst central-3x3 pixel error; `FWHMerr` +compares the fast image with a direct image at the same centroid, so it is a +pure shape error (0 when deposition is exact 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). FWHM 1.0 / β 3.5 rows are sub-Nyquist and are not meaningful +for the shape metrics. + +## Raw output + +```text +Fast-mode deposition error experiment +grid=48 phases=12 tail=1e-06 max_radius=64 +kernel quadrature = 4-point Gauss; reference = 8-point Gauss + + FWHM beta N deposit R_ss retain maxErr rmsErr coreErr centErr FWHMerr beterr +--------------------------------------------------------------------------------------------------- + 1.00 8.00 2 nearest 10 1.00000 6.333e-01 1.714e-02 6.333e-01 3.043e-01 1.733e-02 nan + 1.00 8.00 2 bilinear 10 1.00000 1.805e-01 3.972e-03 1.805e-01 2.722e-03 3.187e-01 nan + 1.00 8.00 3 nearest 14 1.00000 3.907e-01 1.061e-02 3.907e-01 1.847e-01 7.491e-03 nan + 1.00 8.00 3 bilinear 14 1.00000 7.587e-02 1.676e-03 7.587e-02 5.562e-03 1.816e-01 nan + 1.00 8.00 4 nearest 18 1.00000 2.608e-01 7.142e-03 2.608e-01 1.237e-01 3.665e-02 nan + 1.00 8.00 4 bilinear 18 1.00000 4.428e-02 9.775e-04 4.428e-02 3.395e-03 3.646e-02 nan + + 2.70 8.00 2 nearest 24 1.00000 1.635e-01 8.693e-03 1.635e-01 2.946e-01 2.325e-08 5.061e-07 + 2.70 8.00 2 bilinear 24 1.00000 4.339e-02 1.476e-03 4.339e-02 6.018e-09 4.046e-02 8.414e-02 + 2.70 8.00 3 nearest 35 1.00000 9.597e-02 5.229e-03 9.597e-02 1.768e-01 1.465e-08 1.067e-06 + 2.70 8.00 3 bilinear 35 1.00000 1.842e-02 6.387e-04 1.842e-02 2.437e-08 1.747e-02 4.029e-02 + 2.70 8.00 4 nearest 46 1.00000 6.308e-02 3.488e-03 6.308e-02 1.179e-01 1.563e-08 2.905e-05 + 2.70 8.00 4 bilinear 46 1.00000 1.011e-02 3.418e-04 1.011e-02 1.778e-08 9.355e-03 3.091e-02 + + 6.00 8.00 2 nearest 51 1.00000 7.012e-02 7.714e-03 7.012e-02 2.946e-01 1.617e-06 3.132e-06 + 6.00 8.00 2 bilinear 51 1.00000 9.563e-03 6.146e-04 9.563e-03 8.749e-08 6.686e-03 1.039e-02 + 6.00 8.00 3 nearest 64 1.00000 4.199e-02 4.631e-03 4.199e-02 1.768e-01 7.456e-10 3.315e-08 + 6.00 8.00 3 bilinear 64 1.00000 4.090e-03 2.640e-04 4.090e-03 9.147e-06 2.683e-03 4.463e-03 + 6.00 8.00 4 nearest 64 0.99994 2.795e-02 3.088e-03 2.795e-02 1.179e-01 2.631e-07 9.889e-06 + 6.00 8.00 4 bilinear 64 0.99994 2.195e-03 1.409e-04 2.195e-03 9.262e-06 1.528e-03 2.354e-03 + + 1.00 4.50 2 nearest 19 1.00000 6.006e-01 1.643e-02 6.006e-01 3.041e-01 1.516e-02 nan + 1.00 4.50 2 bilinear 19 1.00000 1.745e-01 3.823e-03 1.745e-01 2.673e-03 2.485e-01 nan + 1.00 4.50 3 nearest 28 1.00000 3.712e-01 1.016e-02 3.712e-01 1.846e-01 6.881e-03 nan + 1.00 4.50 3 bilinear 28 1.00000 7.354e-02 1.616e-03 7.354e-02 5.429e-03 1.432e-01 nan + 1.00 4.50 4 nearest 36 1.00000 2.480e-01 6.842e-03 2.480e-01 1.236e-01 3.167e-02 nan + 1.00 4.50 4 bilinear 36 1.00000 4.269e-02 9.383e-04 4.269e-02 3.322e-03 3.640e-02 nan + + 2.70 4.50 2 nearest 49 1.00000 1.625e-01 8.537e-03 1.625e-01 2.946e-01 9.993e-07 1.210e-06 + 2.70 4.50 2 bilinear 49 1.00000 4.413e-02 1.439e-03 4.413e-02 5.868e-07 4.453e-02 4.843e-02 + 2.70 4.50 3 nearest 64 1.00000 9.571e-02 5.135e-03 9.571e-02 1.768e-01 1.482e-07 2.493e-07 + 2.70 4.50 3 bilinear 64 1.00000 1.871e-02 6.228e-04 1.871e-02 3.673e-06 1.867e-02 2.057e-02 + 2.70 4.50 4 nearest 64 0.99999 6.299e-02 3.426e-03 6.299e-02 1.179e-01 1.817e-07 4.228e-07 + 2.70 4.50 4 bilinear 64 0.99999 1.030e-02 3.333e-04 1.030e-02 3.726e-06 9.995e-03 1.353e-02 + + 6.00 4.50 2 nearest 64 0.99998 6.924e-02 7.554e-03 6.924e-02 2.945e-01 1.290e-04 1.735e-04 + 6.00 4.50 2 bilinear 64 0.99998 9.859e-03 6.016e-04 9.859e-03 2.462e-06 7.468e-03 9.748e-03 + 6.00 4.50 3 nearest 64 0.99978 4.153e-02 4.535e-03 4.153e-02 1.775e-01 2.110e-08 3.706e-08 + 6.00 4.50 3 bilinear 64 0.99978 4.214e-03 2.585e-04 4.214e-03 7.436e-04 3.203e-03 4.204e-03 + 6.00 4.50 4 nearest 64 0.99872 2.767e-02 3.024e-03 2.767e-02 1.185e-01 3.310e-06 1.891e-04 + 6.00 4.50 4 bilinear 64 0.99872 2.264e-03 1.485e-04 2.264e-03 7.451e-04 1.709e-03 2.399e-03 + + 1.00 3.50 2 nearest 35 1.00000 5.809e-01 1.600e-02 5.809e-01 3.038e-01 1.370e-02 nan + 1.00 3.50 2 bilinear 35 1.00000 1.708e-01 3.732e-03 1.708e-01 2.601e-03 2.192e-01 nan + 1.00 3.50 3 nearest 52 1.00000 3.594e-01 9.895e-03 3.594e-01 1.843e-01 6.403e-03 nan + 1.00 3.50 3 bilinear 52 1.00000 7.210e-02 1.579e-03 7.210e-02 5.255e-03 1.252e-01 nan + 1.00 3.50 4 nearest 64 1.00000 2.402e-01 6.660e-03 2.402e-01 1.234e-01 2.661e-02 nan + 1.00 3.50 4 bilinear 64 1.00000 4.170e-02 9.143e-04 4.170e-02 3.223e-03 3.644e-02 nan + + 2.70 3.50 2 nearest 64 1.00000 1.619e-01 8.440e-03 1.619e-01 2.946e-01 2.790e-05 2.602e-05 + 2.70 3.50 2 bilinear 64 1.00000 4.461e-02 1.417e-03 4.461e-02 3.036e-07 4.689e-02 3.892e-02 + 2.70 3.50 3 nearest 64 0.99997 9.553e-02 5.077e-03 9.553e-02 1.769e-01 3.421e-07 3.302e-07 + 2.70 3.50 3 bilinear 64 0.99997 1.889e-02 6.133e-04 1.889e-02 8.121e-05 2.016e-02 1.646e-02 + 2.70 3.50 4 nearest 64 0.99989 6.293e-02 3.387e-03 6.293e-02 1.179e-01 8.537e-07 5.529e-06 + 2.70 3.50 4 bilinear 64 0.99989 1.042e-02 3.282e-04 1.042e-02 8.142e-05 1.066e-02 1.038e-02 + +Worst case over all tested configurations, by deposit scheme: + nearest: maxErr=6.333e-01 rmsErr=1.714e-02 coreErr=6.333e-01 centErr=3.043e-01 FWHMerr=3.665e-02 beterr=1.891e-04 + bilinear: maxErr=1.805e-01 rmsErr=3.972e-03 coreErr=1.805e-01 centErr=5.562e-03 FWHMerr=3.187e-01 beterr=8.414e-02 +``` + +## Conclusion + +For the production PSF (FWHM 2.7, β 4.5): + +| N | deposit | core HDR err | image RMS | centroid err | FWHM err | beta change | +|---|---|---|---|---|---|---| +| 2 | nearest | 16.3% | 0.85% | 0.295 px | ~0 | ~1.2e-6 | +| 2 | bilinear | 4.4% | 0.14% | ~0 | +4.45% | +4.8% | +| 3 | nearest | 9.6% | 0.51% | 0.177 px | ~0 | ~2.5e-7 | +| 3 | bilinear | 1.9% | 0.06% | ~0 | +1.87% | +2.1% | +| 4 | nearest | 6.3% | 0.34% | 0.118 px | ~0 | ~4.2e-7 | +| 4 | bilinear | 1.0% | 0.03% | ~0 | +1.0% | +1.4% | + +- Nearest preserves the specified FWHM/β exactly (kernel-only shape and beta + error at the numerical floor); its only error is a bounded position + quantization of `0.5/N` output pixels. +- Bilinear recovers the sub-pixel centroid exactly and roughly quarters the + HDR error, but broadens the profile because it interpolates the kernel + samples rather than the continuous kernel; it therefore changes the + specified FWHM and beta unless separately compensated. +- A single fixed global convolution cannot provide both. `--fast-mode` default + is therefore `nearest`; `--fast-deposit bilinear` is available as an + explicit position-over-shape approximation. diff --git a/nr_spacetime_movie_renderer_design.md b/nr_spacetime_movie_renderer_design.md index 1ca24fa..c9ddb2f 100644 --- a/nr_spacetime_movie_renderer_design.md +++ b/nr_spacetime_movie_renderer_design.md @@ -838,7 +838,20 @@ P(x-x_s,y-y_s) 直接 splat 到 HDR framebuffer。 -不要先把星压进 pixel 再卷积，因为那会丢失亚像素位置。 +不要先把星压进最终 pixel 再卷积，因为那会丢失亚像素位置。 + +`--fast-mode` 是显式预览近似：把每个星像作为 delta 累积进 `N×N` 超采样 +buffer（`--fast-deposit nearest` 只写 1 个超采样像素，`bilinear` 写 4 个 +相邻像素），最后对整张 buffer 用单个全局离散核做一次卷积，再以 `N×N` +box 平均降采样。核取目标 Moffat 在超采样尺度（宽度 `Nα`、同一 `beta`）下 +的 pixel-area 积分且不做额外归一化；`1/N²` 的 box 平均正好重建最终 +pixel-area 积分，所以核本身不改变指定的 FWHM/beta。`nearest` 在吸附到的 +超采样格点上精确保持 PSF 形状，代价是至多 `1/(2N)` 个输出像素的位置量化； +`bilinear` 精确保持亚像素质心，但会略微展宽 FWHM，因此默认是 `nearest`。 +沉积误差测量、`N` 与取舍见 +[`benchmarks/fast_mode_deposit_2026-09-18.md`](benchmarks/fast_mode_deposit_2026-09-18.md)。 +该模式只支持 CPU PSF 后端；`--max-cache-psf-flux` 不适用（所有事件共用同一 +个全局核半径）。 PSF 第一版可用 Gaussian； 以后可换成 Airy 或其他相机模型。 diff --git a/scripts/fast_mode_deposit_error.py b/scripts/fast_mode_deposit_error.py new file mode 100644 index 0000000..3e5f07a --- /dev/null +++ b/scripts/fast_mode_deposit_error.py @@ -0,0 +1,349 @@ +#!/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() diff --git a/src/frame.c b/src/frame.c index 19a9c0e..683057b 100644 --- a/src/frame.c +++ b/src/frame.c @@ -1110,6 +1110,7 @@ typedef struct { double psf_relative_tail; double psf_min_y; PsfEventSink *event_sink; + FastPsfAccumulator *fast; size_t images; size_t direct_fallbacks; size_t cached_wing_clipped; @@ -1187,6 +1188,14 @@ static int splat_catalog_tile(const Star *stars, size_t count, weights[2] * context->vertex[2]->log_frequency_ratio; const LinearRgb color = blackbody_to_linear_rgb(star->temperature_K * exp(log_g)); const double flux = context->exposure * star->amplitude * context->magnification; + if (context->fast != NULL) { + const int fast_status = fast_psf_accumulator_deposit( + context->fast, image_x, image_y, color, flux); + context->cached_wing_clipped += fast_status == 2; + context->discarded_below_min_y += fast_status == 3; + ++context->images; + continue; + } PsfCachedEvent event; const int direct_fallback = psf_prepare_cached_event( &event, image_x, image_y, color, flux, context->psf, context->psf_cache, @@ -1272,7 +1281,8 @@ static CatalogSplatStats splat_catalog_triangles( const PsfKernelCache *psf_cache, double max_magnification, double max_cache_psf_flux, double psf_relative_tail, double psf_min_y, - size_t first_triangle, size_t last_triangle, PsfEventSink *event_sink) { + size_t first_triangle, size_t last_triangle, PsfEventSink *event_sink, + FastPsfAccumulator *fast) { CatalogSplatStats stats = {0}; PsfEventSink owned_sink; const int owns_sink = event_sink == NULL; @@ -1321,7 +1331,8 @@ static CatalogSplatStats splat_catalog_triangles( .max_cache_psf_flux = max_cache_psf_flux, .psf_relative_tail = psf_relative_tail, .psf_min_y = psf_min_y, - .event_sink = event_sink}; + .event_sink = event_sink, + .fast = fast}; const int visit_result = catalog_visit_source_triangle( catalog, direction, 0, splat_catalog_tile, &context); if (visit_result == 0) @@ -1383,7 +1394,8 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, int catalog_load_workers, CatalogPrefetchStats *prefetch_stats, PsfSplatStats *psf_stats, - const FrameSplatProgress *progress) { + const FrameSplatProgress *progress, + FastPsfAccumulator *fast) { if (mesh == NULL || catalog == NULL || hdr == NULL || exposure <= 0.0 || psf == NULL || width <= 0 || height <= 0 || catalog_load_workers <= 0 || isnan(max_magnification) || max_magnification <= 0.0 || @@ -1408,6 +1420,63 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, progress->callback(progress->context, FRAME_SPLAT_PROGRESS_BEGIN, 0, mesh->triangle_count); + /* Fast mode replaces the per-event PSF splat with cheap delta deposits into + * one shared supersampled buffer. A single global convolution plus an N x N + * box average then resolves the frame. Deposits are per-cell atomic adds; + * the immutable kernel and input buffers stay shared read-only. */ + if (fast != NULL) { + size_t images = 0, direct_fallbacks = 0, cached_wing_clipped = 0, + discarded_below_min_y = 0; + int failed = 0; + int worker_count = omp_get_max_threads(); + if (worker_count < 1) + worker_count = 1; + fast_psf_accumulator_clear(fast); +#ifdef GR_DEBUG + double max_raw_magnification = 0.0; + size_t magnification_clamped_triangles = 0; +#pragma omp parallel num_threads(worker_count) reduction(+ : images, direct_fallbacks, cached_wing_clipped, discarded_below_min_y, magnification_clamped_triangles) reduction(max : max_raw_magnification, failed) +#else +#pragma omp parallel num_threads(worker_count) reduction(+ : images, direct_fallbacks, cached_wing_clipped, discarded_below_min_y) reduction(max : failed) +#endif + { +#pragma omp for schedule(dynamic, 1) + for (size_t t = 0; t < mesh->triangle_count; ++t) { + const CatalogSplatStats stats = splat_catalog_triangles( + mesh, catalog, hdr, width, height, exposure, psf, NULL, + max_magnification, max_cache_psf_flux, psf_relative_tail, psf_min_y, + t, t + 1, NULL, fast); + images += stats.images; + direct_fallbacks += stats.direct_fallbacks; + cached_wing_clipped += stats.cached_wing_clipped; + discarded_below_min_y += stats.discarded_below_min_y; + failed |= stats.failed; +#ifdef GR_DEBUG + max_raw_magnification = + fmax(max_raw_magnification, stats.max_raw_magnification); + magnification_clamped_triangles += stats.magnification_clamped_triangles; +#endif + } + } + if (!failed && fast_psf_accumulator_resolve(fast, hdr, worker_count)) + failed = 1; + copy_psf_splat_stats(psf_stats, (CatalogSplatStats){ + .images = images, + .direct_fallbacks = direct_fallbacks, + .cached_wing_clipped = cached_wing_clipped, + .discarded_below_min_y = discarded_below_min_y, +#ifdef GR_DEBUG + .max_raw_magnification = max_raw_magnification, + .magnification_clamped_triangles = + magnification_clamped_triangles, +#endif + }); + if (progress != NULL && progress->callback != NULL) + progress->callback(progress->context, FRAME_SPLAT_PROGRESS_END, + mesh->triangle_count, mesh->triangle_count); + return failed ? SIZE_MAX : images; + } + #ifdef PSF_BACKEND_DUMMY PsfEventSink owner; if (psf_event_sink_init(&owner, hdr, width, height, psf_cache)) return SIZE_MAX; @@ -1442,7 +1511,7 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, const CatalogSplatStats stats = splat_catalog_triangles( mesh, catalog, hdr, width, height, exposure, psf, psf_cache, max_magnification, max_cache_psf_flux, psf_relative_tail, psf_min_y, - t, t + 1, &local); + t, t + 1, &local, NULL); dummy_images += stats.images; dummy_direct += stats.direct_fallbacks; dummy_clipped += stats.cached_wing_clipped; @@ -1531,7 +1600,7 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, const CatalogSplatStats stats = splat_catalog_triangles( mesh, catalog, hdr, width, height, exposure, psf, psf_cache, max_magnification, max_cache_psf_flux, psf_relative_tail, psf_min_y, - t, t + 1, &local); + t, t + 1, &local, NULL); hip_images += stats.images; hip_direct += stats.direct_fallbacks; hip_clipped += stats.cached_wing_clipped; @@ -1594,7 +1663,7 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, const CatalogSplatStats stats = splat_catalog_triangles( mesh, catalog, hdr, width, height, exposure, psf, psf_cache, max_magnification, max_cache_psf_flux, psf_relative_tail, psf_min_y, - 0, mesh->triangle_count, NULL); + 0, mesh->triangle_count, NULL, NULL); copy_psf_splat_stats(psf_stats, stats); return stats.images; } @@ -1611,7 +1680,7 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, const CatalogSplatStats stats = splat_catalog_triangles( mesh, catalog, hdr, width, height, exposure, psf, psf_cache, max_magnification, max_cache_psf_flux, psf_relative_tail, psf_min_y, - 0, mesh->triangle_count, NULL); + 0, mesh->triangle_count, NULL, NULL); copy_psf_splat_stats(psf_stats, stats); return stats.images; } @@ -1622,7 +1691,7 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, const CatalogSplatStats stats = splat_catalog_triangles( mesh, catalog, hdr, width, height, exposure, psf, psf_cache, max_magnification, max_cache_psf_flux, psf_relative_tail, psf_min_y, - 0, mesh->triangle_count, NULL); + 0, mesh->triangle_count, NULL, NULL); copy_psf_splat_stats(psf_stats, stats); return stats.images; } @@ -1639,7 +1708,7 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, const CatalogSplatStats stats = splat_catalog_triangles( mesh, catalog, hdr, width, height, exposure, psf, psf_cache, max_magnification, max_cache_psf_flux, psf_relative_tail, psf_min_y, - 0, mesh->triangle_count, NULL); + 0, mesh->triangle_count, NULL, NULL); copy_psf_splat_stats(psf_stats, stats); return stats.images; } @@ -1673,7 +1742,7 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, mesh, catalog, private_hdr[worker], width, height, exposure, psf, psf_cache, max_magnification, max_cache_psf_flux, psf_relative_tail, psf_min_y, - triangle, triangle + 1, &event_sink); + triangle, triangle + 1, &event_sink, NULL); images += stats.images; direct_fallbacks += stats.direct_fallbacks; cached_wing_clipped += stats.cached_wing_clipped; @@ -1710,7 +1779,7 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, mesh, catalog, private_hdr[worker], width, height, exposure, psf, psf_cache, max_magnification, max_cache_psf_flux, psf_relative_tail, psf_min_y, - triangle, triangle + 1, &event_sink); + triangle, triangle + 1, &event_sink, NULL); images += stats.images; direct_fallbacks += stats.direct_fallbacks; cached_wing_clipped += stats.cached_wing_clipped; diff --git a/src/frame.h b/src/frame.h index dec7cfa..de9405d 100644 --- a/src/frame.h +++ b/src/frame.h @@ -122,7 +122,9 @@ int frame_lens_mesh_refine_with_progress( void *context); /* Each locally invertible escaped triangle contributes one image per contained - * star. */ + * star. When `fast` is non-NULL, each image is deposited into its shared + * supersampled buffer and the whole frame is resolved once at the end instead + * of splatting a per-event PSF. */ size_t frame_splat_catalog(const FrameLensMesh *mesh, StarCatalog *catalog, double *hdr, int width, int height, double exposure, @@ -136,7 +138,8 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, int catalog_load_workers, CatalogPrefetchStats *prefetch_stats, PsfSplatStats *psf_stats, - const FrameSplatProgress *progress); + const FrameSplatProgress *progress, + FastPsfAccumulator *fast); void frame_draw_mesh(const FrameLensMesh *mesh, double *hdr, int width, int height, double gray, double opacity); void frame_lens_mesh_destroy(FrameLensMesh *mesh); diff --git a/src/main.c b/src/main.c index 13a573d..8b5be48 100644 --- a/src/main.c +++ b/src/main.c @@ -34,6 +34,10 @@ typedef struct { int radius_specified, roll_specified; PointSpreadFunction psf; PsfKernelCache psf_cache; + int fast_mode; + int fast_supersample; + FastPsfDeposit fast_deposit; + FastPsfAccumulator *fast_psf; const char *catalog_path; const char *all_sky_catalog_path; const char *output_path; @@ -203,6 +207,8 @@ static int parse_args(int argc, char **argv, Settings *s, .psf_min_y = 0.0, .observer_radius = 30.0, .psf = {2.7, 4.5}, + .fast_supersample = 2, + .fast_deposit = FAST_PSF_DEPOSIT_NEAREST, .catalog_path = "assets/sky_grid_5deg.csv", .output_path = default_output_path, .frames_prefix = "frame", @@ -296,6 +302,18 @@ static int parse_args(int argc, char **argv, Settings *s, !parse_nonnegative(argv[++i], &s->psf_min_y)) { } else if (!strcmp(argv[i], "--psf-direct")) { s->psf_direct = 1; + } else if (!strcmp(argv[i], "--fast-mode")) { + s->fast_mode = 1; + } else if (!strcmp(argv[i], "--fast-supersample") && i + 1 < argc && + !parse_int(argv[++i], &s->fast_supersample)) { + } else if (!strcmp(argv[i], "--fast-deposit") && i + 1 < argc) { + const char *mode = argv[++i]; + if (!strcmp(mode, "nearest")) + s->fast_deposit = FAST_PSF_DEPOSIT_NEAREST; + else if (!strcmp(mode, "bilinear")) + s->fast_deposit = FAST_PSF_DEPOSIT_BILINEAR; + else + return -1; } else if (!strcmp(argv[i], "--verbose")) { s->verbose = 1; } else if (!strcmp(argv[i], "--write-catalog") && i + 1 < argc) @@ -380,6 +398,11 @@ static void print_help(const char *program) { " --psf-relative-tail R Maximum omitted relative PSF tail fraction (default: 1e-8)\n" " --psf-min-y Y Skip events below this linear HDR luminance (default: 0, disabled)\n" " --psf-direct Disable the PSF lookup cache (default: disabled)\n" + " --fast-mode Deposit each image as a supersampled delta and run one\n" + " global PSF convolution + downsample (CPU only; preview)\n" + " --fast-supersample N Fast-mode supersample factor, 1..8 (default: 2)\n" + " --fast-deposit MODE nearest (exact FWHM/beta, 0.5/N px quantization) or\n" + " bilinear (exact centroid, broadens FWHM); default: nearest\n" " --catalog-load-workers N All-sky catalog loader workers (default: 4)\n" " --blackbody-table FILE Explicit GRBBLUT3 table\n" " (default: assets/blackbody/cie1931_2deg_xyz_1024.grbblut)\n", @@ -603,6 +626,18 @@ static int build_observer(const Settings *s, const SpacetimeSource *spacetime, return 0; } +static void report_psf_splat(const Settings *s, const PsfSplatStats *stats) { + if (s->fast_mode) { + fprintf(stderr, + "Fast PSF splats: deposited %zu, wing-clipped %zu, discarded " + "below min-Y %zu\n", + stats->cached_splats, stats->cached_wing_clipped, + stats->discarded_below_min_y); + return; + } + psf_kernel_cache_report(&s->psf_cache, stats, stderr); +} + static int render_observer_frame(const Settings *s, StarCatalog *catalog, const SpacetimeSource *spacetime, const ObserverState *observer, @@ -680,7 +715,8 @@ static int render_observer_frame(const Settings *s, StarCatalog *catalog, s->catalog_load_workers, &prefetch, &psf_stats, &(FrameSplatProgress){report_splat_progress, s->verbose ? report_splat_worker_progress : NULL, - &progress}); + &progress}, + s->fast_psf); if (images == SIZE_MAX) { frame_lens_mesh_destroy(&mesh); free(hdr); @@ -702,7 +738,7 @@ static int render_observer_frame(const Settings *s, StarCatalog *catalog, fprintf(stderr, "Rendered %zu images from %zu catalog stars to %s (%s)\n", images, catalog->count, output_path, result == 0 ? "ok" : "write failed"); - psf_kernel_cache_report(&s->psf_cache, &psf_stats, stderr); + report_psf_splat(s, &psf_stats); frame_lens_mesh_destroy(&mesh); free(hdr); return result; @@ -883,7 +919,8 @@ static int render_movie(const Settings *s, StarCatalog *catalog, s->catalog_load_workers, &prefetch, &psf_stats, &(FrameSplatProgress){report_splat_progress, s->verbose ? report_splat_worker_progress : NULL, - &progress}); + &progress}, + s->fast_psf); if (images == SIZE_MAX) { free(hdr); goto done; @@ -894,7 +931,7 @@ static int render_movie(const Settings *s, StarCatalog *catalog, free(hdr); fprintf(stderr, "Rendered %zu images from %zu catalog stars to %s (%s)\n", images, catalog->count, output_path, write_result == 0 ? "ok" : "write failed"); - psf_kernel_cache_report(&s->psf_cache, &psf_stats, stderr); + report_psf_splat(s, &psf_stats); if (write_result) goto done; } @@ -925,6 +962,22 @@ static int render_lens_map(const Settings *s, StarCatalog *catalog) { return -1; } int result = 0; + FastPsfAccumulator local_fast = {0}; + FastPsfAccumulator *fast = NULL; + if (s->fast_mode && + fast_psf_accumulator_init(&local_fast, map.width, map.height, + s->fast_supersample, s->fast_deposit, &s->psf, + s->psf_relative_tail, s->psf_min_y, + s->psf_direct)) { + fputs("Fast-mode accumulator construction failed for the imported map.\n", + stderr); + lens_map_destroy(&map); + return -1; + } + if (s->fast_mode) { + fast = &local_fast; + fast_psf_accumulator_report(&local_fast, stderr); + } for (size_t i = 0; i < map.frame_count; ++i) { const char *output_path = s->output_path; char movie_path[PATH_MAX]; @@ -959,7 +1012,8 @@ static int render_lens_map(const Settings *s, StarCatalog *catalog) { &prefetch, &psf_stats, &(FrameSplatProgress){report_splat_progress, s->verbose ? report_splat_worker_progress : NULL, - &progress}); + &progress}, + fast); if (images == SIZE_MAX) { free(hdr); result = -1; @@ -984,16 +1038,17 @@ static int render_lens_map(const Settings *s, StarCatalog *catalog) { free(hdr); fprintf(stderr, "Rendered %zu images from %zu catalog stars to %s (%s; imported lens map)\n", images, catalog->count, output_path, write_result == 0 ? "ok" : "write failed"); - psf_kernel_cache_report(&s->psf_cache, &psf_stats, stderr); + report_psf_splat(s, &psf_stats); if (write_result) { result = -1; break; } #else free(hdr); fprintf(stderr, "Dummy PSF classified %zu images from %zu catalog stars; no HDR, PNG, or PPM was written.\n", images, catalog->count); - psf_kernel_cache_report(&s->psf_cache, &psf_stats, stderr); + report_psf_splat(s, &psf_stats); #endif } + fast_psf_accumulator_destroy(&local_fast); lens_map_destroy(&map); return result; } @@ -1028,7 +1083,8 @@ int main(int argc, char **argv) { "[--psf-fwhm-pixels N] [--psf-moffat-beta N] " "[--max-magnification M] [--max-cache-psf-flux F] " "[--psf-relative-tail R] [--psf-min-y Y] " - "[--psf-direct] [--verbose] " + "[--psf-direct] [--fast-mode --fast-supersample N " + "--fast-deposit nearest|bilinear] [--verbose] " #ifdef ENABLE_HDR_OUTPUT "[--hdr-output] " #endif @@ -1073,6 +1129,22 @@ int main(int argc, char **argv) { fputs("Dummy PSF backend: --output and --hdr-output are accepted for command parity but no image files will be written.\n", stderr); #endif +#if defined(PSF_BACKEND_HIP) || defined(PSF_BACKEND_DUMMY) + if (settings.fast_mode) { + fputs("--fast-mode is available only in the CPU PSF backend build.\n", + stderr); + return 2; + } +#endif + if (settings.fast_mode && + (settings.fast_supersample < 1 || settings.fast_supersample > 8)) { + fputs("--fast-supersample must be between 1 and 8.\n", stderr); + return 2; + } + if (settings.fast_mode) + fputs("Fast mode is a preview approximation; --max-cache-psf-flux is " + "ignored and one global kernel is used for every event.\n", + stderr); #ifdef ENABLE_HDR_OUTPUT if (settings.frames_dir != NULL && settings.write_hdr_output) { fputs("--hdr-output is available only for a single-frame render.\n", stderr); @@ -1131,8 +1203,28 @@ int main(int argc, char **argv) { return 1; } fprintf(stderr, "Blackbody backend: %s\n", blackbody_backend_name()); - if (!settings.psf_direct && psf_kernel_cache_init(&settings.psf_cache, &settings.psf, - settings.psf_relative_tail)) { + FastPsfAccumulator fast_accumulator = {0}; + if (settings.fast_mode && settings.lens_map_input_path == NULL) { + if (fast_psf_accumulator_init(&fast_accumulator, settings.width, + settings.height, settings.fast_supersample, + settings.fast_deposit, &settings.psf, + settings.psf_relative_tail, settings.psf_min_y, + settings.psf_direct)) { + fputs("Fast-mode accumulator construction failed.\n", stderr); + catalog_destroy(&catalog); + spacetime_destroy(&spacetime); + blackbody_backend_destroy(); + return 1; + } + settings.fast_psf = &fast_accumulator; + fast_psf_accumulator_report(&fast_accumulator, stderr); + } else if (settings.fast_mode) { + fputs("Fast mode on an imported lens map uses the map's own dimensions.\n", + stderr); + } + if (!settings.fast_mode && !settings.psf_direct && + psf_kernel_cache_init(&settings.psf_cache, &settings.psf, + settings.psf_relative_tail)) { #ifdef PSF_BACKEND_DUMMY fputs("PSF cache construction failed; dummy chunk statistics are unavailable.\n", stderr); @@ -1144,11 +1236,13 @@ int main(int argc, char **argv) { fputs("PSF cache construction failed; using direct evaluator.\n", stderr); #endif } - psf_kernel_cache_report_ready(&settings.psf_cache, stderr); + if (!settings.fast_mode) + psf_kernel_cache_report_ready(&settings.psf_cache, stderr); if (settings.lens_map_input_path != NULL) { const int result = render_lens_map(&settings, &catalog); catalog_destroy(&catalog); psf_kernel_cache_destroy(&settings.psf_cache); + fast_psf_accumulator_destroy(&fast_accumulator); blackbody_backend_destroy(); return result == 0 ? 0 : 1; } @@ -1159,6 +1253,7 @@ int main(int argc, char **argv) { spacetime_destroy(&spacetime); catalog_destroy(&catalog); psf_kernel_cache_destroy(&settings.psf_cache); + fast_psf_accumulator_destroy(&fast_accumulator); blackbody_backend_destroy(); return result == 0 ? 0 : 1; } diff --git a/src/optics.c b/src/optics.c index 08a6970..89cdb42 100644 --- a/src/optics.c +++ b/src/optics.c @@ -450,6 +450,246 @@ void splat_moffat(double *hdr, int width, int height, double x, double y, relative_tail_fraction, min_y); } +static void fast_psf_add_cell(FastPsfAccumulator *accumulator, int sx, int sy, + double weight, LinearRgb color, double flux) +{ + const double scale = flux * weight; + double *pixel; + if (sx < 0 || sy < 0 || sx >= accumulator->supersampled_width || + sy >= accumulator->supersampled_height) + return; + pixel = &accumulator->buffer + [3 * ((size_t)sy * accumulator->supersampled_width + sx)]; +#pragma omp atomic + pixel[0] += color.r * scale; +#pragma omp atomic + pixel[1] += color.g * scale; +#pragma omp atomic + pixel[2] += color.b * scale; +} + +int fast_psf_accumulator_init(FastPsfAccumulator *accumulator, int width, + int height, int supersample, + FastPsfDeposit deposit, + const PointSpreadFunction *psf, + double relative_tail_fraction, double min_y, + int use_reference) +{ + if (accumulator == NULL || !valid_psf(psf) || width <= 0 || height <= 0 || + supersample < 1 || supersample > 8 || + (deposit != FAST_PSF_DEPOSIT_NEAREST && + deposit != FAST_PSF_DEPOSIT_BILINEAR) || + !isfinite(relative_tail_fraction) || relative_tail_fraction <= 0.0 || + relative_tail_fraction >= 1.0 || !isfinite(min_y) || min_y < 0.0) + return -1; + *accumulator = (FastPsfAccumulator){0}; + const double alpha = moffat_alpha(psf); + const double alpha_supersampled = alpha * supersample; + const double requested = moffat_support_radius( + alpha_supersampled, psf->moffat_beta, 1.0, relative_tail_fraction); + if (!isfinite(requested) || requested <= 0.0 || + requested > (double)INT_MAX - 1.0) + return -1; + const int radius = (int)ceil(requested) + 1; + const size_t side = (size_t)2 * radius + 1; + const size_t ss_width = (size_t)width * supersample; + const size_t ss_height = (size_t)height * supersample; + if (side > SIZE_MAX / side || ss_width > SIZE_MAX / 3 || + ss_width > SIZE_MAX / ss_height || + ss_width * ss_height > SIZE_MAX / (3 * sizeof(double)) || + side * side > SIZE_MAX / sizeof(float)) + return -1; + float *weights = malloc(side * side * sizeof *weights); + double *buffer = calloc(ss_width * ss_height * 3, sizeof *buffer); + if (weights == NULL || buffer == NULL) { + free(weights); + free(buffer); + return -1; + } + *accumulator = (FastPsfAccumulator){ + .fwhm_pixels = psf->fwhm_pixels, + .moffat_beta = psf->moffat_beta, + .alpha_pixels = alpha, + .alpha_supersampled = alpha_supersampled, + .relative_tail_fraction = relative_tail_fraction, + .min_y = min_y, + .supersample = supersample, + .deposit = deposit, + .width = width, + .height = height, + .supersampled_width = (int)ss_width, + .supersampled_height = (int)ss_height, + .radius_pixels = radius, + .use_reference = use_reference, + .weights = weights, + .buffer = buffer}; + /* k[m] is the pixel-area integral of I(z/N; alpha, beta) over ss cell m, + * i.e. the final-normalized Moffat evaluated at the supersampled scale. + * Its total weight is ~N^2, so the N x N box average restores unit flux + * and the requested FWHM/beta. The star sits at cell 0's centre (0.5). */ + const double *nodes = use_reference ? gauss8_x : gauss4_x; + const double *node_weights = use_reference ? gauss8_w : gauss4_w; + const int order = use_reference ? 8 : PSF_QUADRATURE_ORDER; + const double normalization_scale = (double)supersample * supersample; + for (int dy = -radius; dy <= radius; ++dy) + for (int dx = -radius; dx <= radius; ++dx) { + const double weight = moffat_pixel_integral_quadrature( + alpha_supersampled, psf->moffat_beta, dx, dy, 0.5, 0.5, + nodes, node_weights, order); + weights[(size_t)(dy + radius) * side + (dx + radius)] = + (float)(weight * normalization_scale); + } + return 0; +} + +void fast_psf_accumulator_clear(FastPsfAccumulator *accumulator) +{ + if (accumulator == NULL || accumulator->buffer == NULL) + return; + memset(accumulator->buffer, 0, + (size_t)accumulator->supersampled_width * + accumulator->supersampled_height * 3 * sizeof(double)); +} + +int fast_psf_accumulator_deposit(FastPsfAccumulator *accumulator, double x, + double y, LinearRgb color, double flux) +{ + if (accumulator == NULL || accumulator->buffer == NULL) + return 1; + if (!isfinite(x) || !isfinite(y) || !isfinite(flux) || flux <= 0.0 || + !isfinite(color.r) || !isfinite(color.g) || !isfinite(color.b)) + return 1; + const double support_radius = moffat_effective_support_radius( + accumulator->alpha_pixels, accumulator->moffat_beta, color, flux, + accumulator->relative_tail_fraction, accumulator->min_y); + if (support_radius == 0.0) + return 3; + const int wing_clipped = + support_radius > + (double)accumulator->radius_pixels / accumulator->supersample + ? 2 : 0; + const double zx = accumulator->supersample * x; + const double zy = accumulator->supersample * y; + if (accumulator->deposit == FAST_PSF_DEPOSIT_NEAREST) { + fast_psf_add_cell(accumulator, (int)floor(zx), (int)floor(zy), 1.0, + color, flux); + } else { + const int base_x = (int)floor(zx - 0.5); + const int base_y = (int)floor(zy - 0.5); + const double fx = zx - 0.5 - base_x; + const double fy = zy - 0.5 - base_y; + fast_psf_add_cell(accumulator, base_x, base_y, + (1.0 - fx) * (1.0 - fy), color, flux); + fast_psf_add_cell(accumulator, base_x + 1, base_y, + fx * (1.0 - fy), color, flux); + fast_psf_add_cell(accumulator, base_x, base_y + 1, + (1.0 - fx) * fy, color, flux); + fast_psf_add_cell(accumulator, base_x + 1, base_y + 1, fx * fy, + color, flux); + } + return wing_clipped; +} + +int fast_psf_accumulator_resolve(const FastPsfAccumulator *accumulator, + double *hdr, int worker_count) +{ + const int supersample = accumulator == NULL ? 0 : accumulator->supersample; + const int radius = accumulator == NULL ? 0 : accumulator->radius_pixels; + const size_t side = (size_t)2 * radius + 1; + if (accumulator == NULL || accumulator->weights == NULL || + accumulator->buffer == NULL || hdr == NULL || supersample <= 0) + return -1; + int *row_span = malloc(side * sizeof *row_span); + if (row_span == NULL) + return -1; + for (int dy = -radius; dy <= radius; ++dy) { + const double remaining = + (double)radius * radius - (double)dy * dy; + row_span[dy + radius] = + remaining > 0.0 ? (int)sqrt(remaining) : 0; + } + const double inverse_block = 1.0 / ((double)supersample * supersample); + const int threads = worker_count > 0 ? worker_count : 1; +#pragma omp parallel for schedule(static) num_threads(threads) if (threads > 1) + for (int row = 0; row < accumulator->height; ++row) { + for (int column = 0; column < accumulator->width; ++column) { + double sum[3] = {0.0, 0.0, 0.0}; + for (int sy = supersample * row; + sy < supersample * row + supersample; ++sy) { + if (sy < 0 || sy >= accumulator->supersampled_height) + continue; + for (int sx = supersample * column; + sx < supersample * column + supersample; ++sx) { + if (sx < 0 || sx >= accumulator->supersampled_width) + continue; + for (int dy = -radius; dy <= radius; ++dy) { + const int yy = sy - dy; + if (yy < 0 || yy >= accumulator->supersampled_height) + continue; + const int span = row_span[dy + radius]; + for (int dx = -span; dx <= span; ++dx) { + const int xx = sx - dx; + size_t offset; + const double weight = + accumulator->weights + [(size_t)(dy + radius) * side + + (size_t)(dx + radius)]; + if (xx < 0 || xx >= accumulator->supersampled_width) + continue; + offset = 3 * ((size_t)yy * + accumulator->supersampled_width + + xx); + sum[0] += weight * accumulator->buffer[offset]; + sum[1] += weight * accumulator->buffer[offset + 1]; + sum[2] += weight * accumulator->buffer[offset + 2]; + } + } + } + } + const size_t out = 3 * ((size_t)row * accumulator->width + column); + hdr[out] += sum[0] * inverse_block; + hdr[out + 1] += sum[1] * inverse_block; + hdr[out + 2] += sum[2] * inverse_block; + } + } + free(row_span); + return 0; +} + +void fast_psf_accumulator_destroy(FastPsfAccumulator *accumulator) +{ + if (accumulator == NULL) + return; + free(accumulator->weights); + free(accumulator->buffer); + *accumulator = (FastPsfAccumulator){0}; +} + +void fast_psf_accumulator_report(const FastPsfAccumulator *accumulator, + FILE *stream) +{ + if (stream == NULL || accumulator == NULL || accumulator->weights == NULL) + return; + double sum = 0.0; + const size_t count = + (size_t)(2 * accumulator->radius_pixels + 1) * + (size_t)(2 * accumulator->radius_pixels + 1); + for (size_t i = 0; i < count; ++i) + sum += accumulator->weights[i]; + fprintf(stream, + "Fast PSF: supersample %dx, deposit %s, kernel radius %d ss px " + "(%.2f final px), retained flux %.6f, %s quadrature\n", + accumulator->supersample, + accumulator->deposit == FAST_PSF_DEPOSIT_NEAREST ? "nearest" + : "bilinear", + accumulator->radius_pixels, + (double)accumulator->radius_pixels / accumulator->supersample, + sum / ((double)accumulator->supersample * + accumulator->supersample), + accumulator->use_reference ? "reference 8-point" : "cached 4-point"); +} + + static unsigned char tonemap_channel(double hdr_value) { /* Reinhard tone mapping followed by the sRGB display transfer curve. */ diff --git a/src/optics.h b/src/optics.h index 7f1a6bf..2306a7d 100644 --- a/src/optics.h +++ b/src/optics.h @@ -21,6 +21,33 @@ typedef struct { int ready; } PsfKernelCache; +typedef enum { + FAST_PSF_DEPOSIT_NEAREST = 0, + FAST_PSF_DEPOSIT_BILINEAR = 1, +} FastPsfDeposit; + +/* Fast point-source accumulation: every image event is deposited as a delta + * (one nearest supersampled pixel, or 4 bilinear pixels) into one shared + * supersampled HDR buffer. A single immutable global kernel is convolved once + * over the whole buffer, then an N x N box average downsamples to the final + * image. `nearest` reproduces the current pixel-integrated Moffat exactly at + * the snapped supersampled position; `bilinear` preserves the sub-pixel + * centroid but broadens the profile (see + * benchmarks/fast_mode_deposit_2026-09-18.md). The kernel is the pixel-area + * integral of the Moffat with alpha scaled by `supersample` and the same beta, + * so the box average keeps the requested FWHM/beta semantics. */ +typedef struct { + double fwhm_pixels, moffat_beta, alpha_pixels, alpha_supersampled; + double relative_tail_fraction, min_y; + int supersample; + FastPsfDeposit deposit; + int width, height, supersampled_width, supersampled_height; + int radius_pixels; /* kernel radius in supersampled pixels */ + int use_reference; /* 8-point instead of 4-point kernel quadrature */ + float *weights; /* (2R+1)^2 pixel-area kernel, row-major */ + double *buffer; /* supersampled HDR, 3 channels per pixel */ +} FastPsfAccumulator; + typedef struct { size_t cached_splats; size_t cached_wing_clipped; @@ -67,6 +94,27 @@ void psf_kernel_cache_destroy(PsfKernelCache *cache); void psf_kernel_cache_report_ready(const PsfKernelCache *cache, FILE *stream); void psf_kernel_cache_report(const PsfKernelCache *cache, const PsfSplatStats *stats, FILE *stream); +/* Fast-mode accumulator. init builds the immutable global kernel and the + * shared supersampled buffer; deposit is thread-safe (per-cell atomic add); + * resolve convolves and downsamples, accumulating into the final HDR buffer + * (it does not overwrite it). Returns 0 on success and -1 on invalid input or + * allocation failure. */ +int fast_psf_accumulator_init(FastPsfAccumulator *accumulator, int width, + int height, int supersample, + FastPsfDeposit deposit, + const PointSpreadFunction *psf, + double relative_tail_fraction, double min_y, + int use_reference); +void fast_psf_accumulator_clear(FastPsfAccumulator *accumulator); +/* Returns 0 for a full deposit, 2 when the event's support exceeds the global + * kernel radius (wing clipped), and 3 when discarded by --psf-min-y. */ +int fast_psf_accumulator_deposit(FastPsfAccumulator *accumulator, double x, + double y, LinearRgb color, double flux); +int fast_psf_accumulator_resolve(const FastPsfAccumulator *accumulator, + double *hdr, int worker_count); +void fast_psf_accumulator_destroy(FastPsfAccumulator *accumulator); +void fast_psf_accumulator_report(const FastPsfAccumulator *accumulator, + FILE *stream); /* Reference implementation: pixel-area-integrated Moffat with the same tail * budgets as the cache-aware renderer. */ void splat_moffat_direct(double *hdr, int width, int height, double x, double y, diff --git a/tests/capture_psf.c b/tests/capture_psf.c index 48e6edd..498e752 100644 --- a/tests/capture_psf.c +++ b/tests/capture_psf.c @@ -119,7 +119,7 @@ int main(int argc, char **argv) { #pragma omp for schedule(dynamic,1) for(size_t t=first;t @@ -304,6 +304,140 @@ int main(void) { free(cached_hdr); free(reference_hdr); psf_kernel_cache_destroy(&cache); + /* Fast mode: a nearest deposit must reproduce the current pixel-integrated + * Moffat at the snapped supersampled centre, preserve total flux, and honour + * the min-Y discard rule; bilinear deposition must preserve the centroid. */ + { + int fast_ok = 1; + const int supersample = 2; + FastPsfAccumulator fast = {0}; + FastPsfAccumulator bilinear = {0}; + FastPsfAccumulator min_y_fast = {0}; + FastPsfAccumulator accumulation = {0}; + double *fast_hdr = calloc((size_t)width * height * 3, sizeof *fast_hdr); + double *direct_hdr = calloc((size_t)width * height * 3, sizeof *direct_hdr); + double *background_hdr = + calloc((size_t)width * height * 3, sizeof *background_hdr); + if (fast_hdr == NULL || direct_hdr == NULL || background_hdr == NULL || + fast_psf_accumulator_init(&fast, width, height, supersample, + FAST_PSF_DEPOSIT_NEAREST, &psf, + psf_relative_tail, 0.0, 1) || + fast_psf_accumulator_init(&bilinear, width, height, supersample, + FAST_PSF_DEPOSIT_BILINEAR, &psf, + psf_relative_tail, 0.0, 1) || + fast_psf_accumulator_init(&min_y_fast, width, height, supersample, + FAST_PSF_DEPOSIT_NEAREST, &psf, + psf_relative_tail, 0.5, 1) || + fast_psf_accumulator_init(&accumulation, width, height, supersample, + FAST_PSF_DEPOSIT_NEAREST, &psf, + psf_relative_tail, 0.0, 1)) { + fputs("fast-mode accumulator construction regression failed\n", stderr); + fast_ok = 0; + goto fast_done; + } + /* (50.2, 50.2) snaps to supersampled cell 100, centre (100.5, 100.5) in + * ss coordinates, i.e. final position (50.25, 50.25). */ + if (fast_psf_accumulator_deposit(&fast, 50.2, 50.2, + (LinearRgb){1.0, 1.0, 1.0}, 1.0) == 3 || + fast_psf_accumulator_resolve(&fast, fast_hdr, 4)) { + fputs("fast-mode nearest deposit regression failed\n", stderr); + fast_ok = 0; + goto fast_done; + } + splat_moffat_direct(direct_hdr, width, height, 50.25, 50.25, + (LinearRgb){1.0, 1.0, 1.0}, 1.0, &psf, + psf_relative_tail, 0.0); + double peak = 0.0, max_error = 0.0, fast_flux = 0.0, direct_flux = 0.0; + for (int value = 0; value < width * height * 3; ++value) { + peak = fmax(peak, direct_hdr[value]); + max_error = fmax(max_error, fabs(fast_hdr[value] - direct_hdr[value])); + fast_flux += fast_hdr[value]; + direct_flux += direct_hdr[value]; + } + if (!(peak > 0.0) || max_error > 1e-4 * peak || + fabs(fast_flux - direct_flux) > 1e-4) { + fputs("fast-mode nearest semantics regression failed\n", stderr); + fast_ok = 0; + goto fast_done; + } + /* Bilinear keeps the exact continuous centroid. */ + memset(fast_hdr, 0, (size_t)width * height * 3 * sizeof *fast_hdr); + if (fast_psf_accumulator_deposit(&bilinear, 50.37, 50.62, + (LinearRgb){1.0, 1.0, 1.0}, 1.0) == 3 || + fast_psf_accumulator_resolve(&bilinear, fast_hdr, 4)) { + fputs("fast-mode bilinear deposit regression failed\n", stderr); + fast_ok = 0; + goto fast_done; + } + double weight_sum = 0.0, cx = 0.0, cy = 0.0; + for (int row = 0; row < height; ++row) + for (int column = 0; column < width; ++column) { + const double weight = fast_hdr[3 * (row * width + column)]; + weight_sum += weight; + cx += weight * (column + 0.5); + cy += weight * (row + 0.5); + } + if (!(weight_sum > 0.0) || + hypot(cx / weight_sum - 50.37, cy / weight_sum - 50.62) > 1e-6) { + fputs("fast-mode bilinear centroid regression failed\n", stderr); + fast_ok = 0; + goto fast_done; + } + /* The min-Y cutoff discards an event whose peak luminance is below it. */ + if (fast_psf_accumulator_deposit(&fast, 10.5, 10.5, + (LinearRgb){1.0, 1.0, 1.0}, 1.0) != 0 || + fast_psf_accumulator_deposit(&min_y_fast, 10.5, 10.5, + (LinearRgb){1.0, 1.0, 1.0}, 1.0) != 3) { + fputs("fast-mode min-Y discard regression failed\n", stderr); + fast_ok = 0; + goto fast_done; + } + /* End-to-end plumbing through the frame splat path. */ + memset(fast_hdr, 0, (size_t)width * height * 3 * sizeof *fast_hdr); + PsfSplatStats fast_stats = {0}; + const size_t fast_images = frame_splat_catalog( + &mesh, &catalog, fast_hdr, width, height, test_exposure, &psf, NULL, + INFINITY, 1.0, psf_relative_tail, 0.0, 0, 1, NULL, &fast_stats, NULL, + &fast); + if (fast_images != 1 || fast_stats.discarded_below_min_y != 0) { + fputs("fast-mode frame splat regression failed\n", stderr); + fast_ok = 0; + goto fast_done; + } + /* HDR accumulation semantics: resolve must add onto an existing + * background, not overwrite it. A prefilled buffer plus one deposit must + * preserve the far-field background exactly and add the PSF core. */ + for (int value = 0; value < width * height * 3; ++value) + background_hdr[value] = 0.25; + if (fast_psf_accumulator_deposit(&accumulation, 10.5, 10.5, + (LinearRgb){1.0, 1.0, 1.0}, 1.0) == 3 || + fast_psf_accumulator_resolve(&accumulation, background_hdr, 4)) { + fputs("fast-mode HDR accumulation regression failed\n", stderr); + fast_ok = 0; + goto fast_done; + } + /* (90, 90) is far outside the kernel support of a star at (10.5, 10.5). */ + if (background_hdr[3 * (90 * width + 90)] != 0.25) { + fputs("fast-mode HDR accumulation lost the background\n", stderr); + fast_ok = 0; + goto fast_done; + } + if (!(background_hdr[3 * (10 * width + 10)] > 0.25)) { + fputs("fast-mode HDR accumulation did not add the deposit\n", stderr); + fast_ok = 0; + goto fast_done; + } + fast_done: + free(fast_hdr); + free(direct_hdr); + free(background_hdr); + fast_psf_accumulator_destroy(&fast); + fast_psf_accumulator_destroy(&bilinear); + fast_psf_accumulator_destroy(&min_y_fast); + fast_psf_accumulator_destroy(&accumulation); + if (!fast_ok) + goto done; + } frame_draw_mesh(&mesh, hdr, width, height, 0.5, 0.5); if (hdr[3 * (10 * width + 20)] != 0.25) { fputs("mesh diagnostic overlay regression failed\n", stderr); @@ -325,7 +459,7 @@ int main(void) { frame_splat_catalog(&fine_mesh, &fine_catalog, hdr, width, height, test_exposure, &psf, NULL, INFINITY, 1.0, psf_relative_tail, 0.0, 0, 1, NULL, - NULL, NULL) != 1) { + NULL, NULL, NULL) != 1) { fputs("fine source-triangle containment regression failed\n", stderr); frame_lens_mesh_destroy(&fine_mesh); goto done; @@ -369,7 +503,7 @@ int main(void) { memset(hdr, 0, (size_t)width * height * 3 * sizeof *hdr); if (frame_splat_catalog(&thin_mesh, &thin_catalog, hdr, width, height, test_exposure, &psf, NULL, 1.0, 1.0, - psf_relative_tail, 0.0, 0, 1, NULL, NULL, NULL) != 1 || + psf_relative_tail, 0.0, 0, 1, NULL, NULL, NULL, NULL) != 1 || hdr[3 * (43 * width + 43)] <= 0.0) { fputs("thin source-triangle inverse-map regression failed\n", stderr); goto done; diff --git a/usage.md b/usage.md index 49ef529..2710529 100644 --- a/usage.md +++ b/usage.md @@ -267,6 +267,30 @@ at the existing cache radius. Images above `F` retain the direct fallback. `--max-magnification M` (default unlimited) caps the per-triangle rendering magnification before flux is formed; it is an explicit preview approximation. +### Fast preview mode + +`--fast-mode` replaces the per-event PSF splat with a two-stage approximation: +each point-source image is deposited as a delta into an `N×N` supersampled HDR +buffer (`--fast-supersample N`, default 2), then the whole frame is convolved +once with a single global Moffat kernel and averaged down by the `N×N` block. +The kernel is the pixel-area integral of the requested Moffat at the +supersampled scale, so the requested FWHM and beta are preserved by the +downsample itself. + +`--fast-deposit nearest` (default) deposits into the single nearest +supersampled pixel. It keeps the PSF shape exactly but quantizes the image +position to `1/(2N)` output pixels. `--fast-deposit bilinear` deposits into 4 +adjacent pixels; it preserves the continuous centroid but broadens the profile +(measured +4.5%/+1.9%/+1.0% FWHM at N=2/3/4). The deposition experiment and its +raw output are in +[benchmarks/fast_mode_deposit_2026-09-18.md](benchmarks/fast_mode_deposit_2026-09-18.md). + +Fast mode is available only in the CPU PSF backend. It ignores +`--max-cache-psf-flux` and uses one global kernel radius for every event, so a +bright event whose requested support exceeds that radius is wing-clipped and +counted in the report. `--psf-min-y` is still applied per event. Fast mode is a +preview approximation, not the physically exact per-event PSF path. + ## Progress and diagnostics The PSF-cache completion line is printed before tracing and catalog splatting