diff --git a/README.md b/README.md index 621e582..234c66d 100644 --- a/README.md +++ b/README.md @@ -132,10 +132,12 @@ 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). 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`. +are described in [usage.md](usage.md). For fast single-frame previews, +`--fast-mode` replaces per-event PSF splats with a supersampled delta deposit +plus one global PSF convolution and downsample. Its nearest-deposit path is not +qualified for temporally coherent movie output; see [usage.md](usage.md) for +the spatial and temporal accuracy limits. Both binaries provide a complete +option list with `--help`. ### Example: Looking outward just above the Schwarzschild horizon diff --git a/benchmarks/fast_mode_fftw_2026-09-25.md b/benchmarks/fast_mode_fftw_2026-09-25.md index 2464b7c..122b5de 100644 --- a/benchmarks/fast_mode_fftw_2026-09-25.md +++ b/benchmarks/fast_mode_fftw_2026-09-25.md @@ -5,6 +5,10 @@ convolution. This record separates one-time plan/setup from steady-state execution and repeats the bounded Galactic-center comparison. Raw resolver and end-to-end output is pasted verbatim from the runs below. +The later production-density 4K comparison is recorded separately against the +first committed FFTW implementation, revision +`229f50cd864d9768324042c10edbb56bfc120f08`. + ## Environment - Git revision: `3deebfb` plus the uncommitted FFTW work. @@ -269,6 +273,222 @@ fast_n2.png 1f144aefb4a5520d9d76a13ab69026a72580663e5e45f53128a87a0801176ccf fast_n4.png 679f14a67a085f76c144ecbe36af51b3b2534c152e8ef59708a7884a3f7c8464 ``` +## Production-density 4K Galactic-center comparison + +This is the full 3840x2160, 50-degree-field run over the local processed 2MASS +all-sky catalog, not the bounded 320x240 comparison above. + +### Provenance + +- Git revision: `229f50cd864d9768324042c10edbb56bfc120f08` + (`Feat: Use FFTW linear convolution for CPU fast-mode PSF resolve`). +- The tracked worktree matched that revision when this record was written; the + unrelated untracked `nmesh_spacetime_output_format.md` was not an input. +- Executable: `build/Release/minkowski_sky`, built at + `2026-09-25 23:16:45 -0400`, before the three runs below; size 146,856 bytes; + SHA-256 + `7f6486285b2ff498adb2d4f1acaf62dc96b11125c9cec114f572b9379648962c`. + A Make dry-run after the commit requested only the Makefile's unconditional + `FORCE` relink and no object recompilation, so the recorded binary used the + object set corresponding to the committed tracked sources. +- Host: 12th Gen Intel Core i7-12700K, 16 hardware threads, 128 GB. The commands + did not explicitly set `OMP_NUM_THREADS`; fast-mode initialization reported + `workers=16`. +- Input: `assets/2mass/processed/all_sky`, observed when recording as 64,801 + regular files and 16,207,861,727 bytes. The renderer selected 1,742 tiles and + loaded 54,509,992 stars in each run. No content manifest hash was recorded for + this external processed catalog. +- Timing: zsh built-in `time`; its final timing line is included verbatim for + each command. + +### Standard cached spatial PSF + +```sh +time ./build/Release/minkowski_sky \ + --psf-relative-tail 1e-6 --psf-min-y 1e-8 --max-cache-psf-flux 1e8 \ + --all-sky-catalog assets/2mass/processed/all_sky \ + --width 3840 --height 2160 \ + --look-ra-deg 262.5 --look-dec-deg -30 \ + --fov-deg 50 \ + --coarse-cell-pixels 16 \ + --exposure 1e13 \ + --psf-fwhm-pixels 2.7 --psf-moffat-beta 3 \ + --catalog-load-workers 4 \ + --hdr-output \ + --output output/imgs/2mass_galactic_center_flat_4k_50deg_1e13_beta3-base.png +``` + +```text +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.000130 s; payload fnv1a64=0e39867b0e70a809 +Blackbody backend: lut +PSF cache ready: 64x64 phases, radius 85 px, relative tail 1e-06, tail abs 1e-06, boundary 1e-07, build 2.176 s +Catalog prefetch: 1742 requested, 1742 newly loaded (54509992 stars), 0 unavailable in 5.886 s; 4 loader workers +Rendered 52544535 images from 54509992 catalog stars to output/imgs/2mass_galactic_center_flat_4k_50deg_1e13_beta3-base.png (ok) +PSF splats: cached 52477581, cached wing-clipped 0, direct fallbacks 0, discarded below min-Y 66954 +Warning: --psf-min-y discarded one or more PSF events. +./build/Release/minkowski_sky --psf-relative-tail 1e-6 --psf-min-y 1e-8 1e8 938.60s user 15.80s system 1391% cpu 1:08.61 total +``` + +### Fast mode, N=2 + +```sh +time ./build/Release/minkowski_sky --fast-mode --fast-supersample 2 \ + --psf-relative-tail 1e-6 --psf-min-y 1e-8 --max-cache-psf-flux 1e8 \ + --all-sky-catalog assets/2mass/processed/all_sky \ + --width 3840 --height 2160 \ + --look-ra-deg 262.5 --look-dec-deg -30 \ + --fov-deg 50 \ + --coarse-cell-pixels 16 \ + --exposure 1e13 \ + --psf-fwhm-pixels 2.7 --psf-moffat-beta 3 \ + --catalog-load-workers 4 \ + --hdr-output \ + --output output/imgs/2mass_galactic_center_flat_4k_50deg_1e13_beta3-fast-2.png +``` + +```text +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.000101 s; payload fnv1a64=0e39867b0e70a809 +Blackbody backend: lut +Fast PSF: supersample 2x, deposit nearest, kernel radius 169 ss px (84.50 final px), retained flux 0.999999, cached 4-point quadrature +Fast FFTW: linear min=4658x8018, fft=4704x8064, workers=16, plan=estimate, plan=0.005275 s, kernel_fft=0.141445 s, scratch=2026.1 MiB +Catalog prefetch: 1742 requested, 1742 newly loaded (54509992 stars), 0 unavailable in 5.873 s; 4 loader workers +Rendered 52544535 images from 54509992 catalog stars to output/imgs/2mass_galactic_center_flat_4k_50deg_1e13_beta3-fast-2.png (ok) +Fast PSF splats: deposited 52477581, wing-clipped 0, discarded below min-Y 66954 +./build/Release/minkowski_sky --fast-mode --fast-supersample 2 1e-6 1e-8 115.74s user 2.08s system 763% cpu 15.439 total +``` + +### Fast mode, N=4 + +```sh +time ./build/Release/minkowski_sky --fast-mode --fast-supersample 4 \ + --psf-relative-tail 1e-6 --psf-min-y 1e-8 --max-cache-psf-flux 1e8 \ + --all-sky-catalog assets/2mass/processed/all_sky \ + --width 3840 --height 2160 \ + --look-ra-deg 262.5 --look-dec-deg -30 \ + --fov-deg 50 \ + --coarse-cell-pixels 16 \ + --exposure 1e13 \ + --psf-fwhm-pixels 2.7 --psf-moffat-beta 3 \ + --catalog-load-workers 4 \ + --hdr-output \ + --output output/imgs/2mass_galactic_center_flat_4k_50deg_1e13_beta3-fast-4.png +``` + +```text +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.000043 s; payload fnv1a64=0e39867b0e70a809 +Blackbody backend: lut +Fast PSF: supersample 4x, deposit nearest, kernel radius 336 ss px (84.00 final px), retained flux 0.999999, cached 4-point quadrature +Fast FFTW: linear min=9312x16032, fft=9375x16128, workers=16, plan=estimate, plan=0.005384 s, kernel_fft=0.743673 s, scratch=8075.5 MiB +Catalog prefetch: 1742 requested, 1742 newly loaded (54509992 stars), 0 unavailable in 5.796 s; 4 loader workers +Rendered 52544535 images from 54509992 catalog stars to output/imgs/2mass_galactic_center_flat_4k_50deg_1e13_beta3-fast-4.png (ok) +Fast PSF splats: deposited 52477581, wing-clipped 0, discarded below min-Y 66954 +./build/Release/minkowski_sky --fast-mode --fast-supersample 4 1e-6 1e-8 135.62s user 5.04s system 722% cpu 19.469 total +``` + +### Artifacts and comparison + +```text +05517bde1383327043691ac4b713f638b03ebf6dc72207fe8710a8dce747a03e output/imgs/2mass_galactic_center_flat_4k_50deg_1e13_beta3-base.png +86e9220feac277d77eba87aa8d0febadca5f4eea35cb56164a4eb1578b20fccd output/imgs/2mass_galactic_center_flat_4k_50deg_1e13_beta3-fast-2.png +5dfef3b4bd39d9809463594a1065c32c65ad397396514664c8aa22a5da6d780b output/imgs/2mass_galactic_center_flat_4k_50deg_1e13_beta3-fast-4.png +1ea508d1ee787c2e5f35ccdf37ca5f4a58e9ff5e9984b50928d66ecf88c44c3d output/imgs/2mass_galactic_center_flat_4k_50deg_1e13_beta3-base_HDR.fits +71de3f1e11cb9f811f5d08ea45b51baf5968e8523873873d6f7c414c1849ba53 output/imgs/2mass_galactic_center_flat_4k_50deg_1e13_beta3-fast-2_HDR.fits +554848e39fdfb3b4815fb8f2a591b8d356cbdeaaea235a73f19b46017a18e42a output/imgs/2mass_galactic_center_flat_4k_50deg_1e13_beta3-fast-4_HDR.fits +``` + +All three runs produced 52,544,535 images and applied the same event filtering: +52,477,581 cached/deposited events, no wing clipping or direct fallback, and +66,954 events (0.12742334% of rendered images) discarded by `--psf-min-y`. + +Using the standard HDR buffer as the reference, with +`relative_L2 = ||fast - base||_2 / ||base||_2`, the existing FITS artifacts give: + +| Mode | Wall time | Wall speedup | Approx. speedup after subtracting catalog prefetch | User-CPU speedup | Scratch | HDR relative L2 | HDR total-flux ratio | Tone-mapped PNG normalized RMSE | +|---|---:|---:|---:|---:|---:|---:|---:|---:| +| Standard | 68.610 s | 1.000x | 1.000x | 1.000x | n/a | 0 | 1 | 0 | +| Fast N=2 | 15.439 s | 4.444x | 6.557x | 8.110x | 2026.1 MiB | 10.4384% | 1.000176588 | 0.00678916 | +| Fast N=4 | 19.469 s | 3.524x | 4.587x | 6.921x | 8075.5 MiB | 4.37249% | 1.000178510 | 0.00365798 | + +The post-prefetch figures subtract only the renderer-reported catalog prefetch +duration from wall time; they are useful context, not isolated PSF timings. N=4 +costs 26.1% more wall time and 3.986x the scratch of N=2, while reducing the HDR +relative-L2 difference by about 58% and the tone-mapped PNG RMSE by about 46%. +The nearly identical `+0.0177%` total-flux offset for N=2 and N=4 shows that the +larger N primarily improves spatial placement rather than global photometry. +For the intended rapid single-frame exposure/PSF preview, N=2 is therefore the +default throughput point; N=4 is the higher-fidelity preview option. + +### Error distribution across the image + +The unweighted relative L2 above is dominated by a small number of very bright +PSF cores, so it understates how well the dense faint-star texture is +reproduced in this still frame. `scripts/fast_mode_error_decomposition.py` +reports the error over the flattened scalar RGB channel samples: + +```sh +python3 scripts/fast_mode_error_decomposition.py \ + output/imgs/2mass_galactic_center_flat_4k_50deg_1e13_beta3-base_HDR.fits \ + output/imgs/2mass_galactic_center_flat_4k_50deg_1e13_beta3-fast-2_HDR.fits \ + output/imgs/2mass_galactic_center_flat_4k_50deg_1e13_beta3-fast-4_HDR.fits +``` + +```text +===== 2mass_galactic_center_flat_4k_50deg_1e13_beta3-fast-2_HDR.fits ===== relative_L2=10.4384% +top |d| RGB samples: sq-error energy base energy linear-RGB sum sample share + top 0.0010%: 69.136% 48.948% 2.822% 0.0010% + top 0.0100%: 95.732% 87.994% 9.141% 0.0100% + top 0.1000%: 99.365% 97.633% 16.396% 0.1000% + top 1.0000%: 99.893% 99.224% 25.984% 1.0000% + top 10.0000%: 99.988% 99.679% 46.006% 10.0000% +base-value range (disjoint): rel_RMS signed bias base energy % linear-RGB sum % samples + [1e-06,0.001): 4.7788% +0.7665% 0.000% 0.001% 67026 + [0.001,0.0599): 4.5658% +0.1902% 0.013% 11.063% 12374574 + [0.0599,0.1): 4.6294% +0.1297% 0.027% 10.881% 4486724 + [0.1,0.229): 4.9079% +0.0866% 0.124% 25.365% 5466556 + [0.229,0.769): 6.9340% +0.0100% 0.301% 24.638% 2239488 + [0.769,4.62): 11.4421% -0.0990% 0.615% 10.315% 223948 + [4.62,inf): 10.4475% -0.1789% 98.921% 17.738% 24884 + +===== 2mass_galactic_center_flat_4k_50deg_1e13_beta3-fast-4_HDR.fits ===== relative_L2=4.3725% +top |d| RGB samples: sq-error energy base energy linear-RGB sum sample share + top 0.0010%: 65.911% 35.557% 2.418% 0.0010% + top 0.0100%: 94.540% 86.397% 8.879% 0.0100% + top 0.1000%: 99.108% 96.708% 16.179% 0.1000% + top 1.0000%: 99.848% 99.208% 25.909% 1.0000% + top 10.0000%: 99.982% 99.694% 45.966% 10.0000% +base-value range (disjoint): rel_RMS signed bias base energy % linear-RGB sum % samples + [1e-06,0.001): 2.4187% +0.5689% 0.000% 0.001% 67026 + [0.001,0.0599): 2.2772% +0.0729% 0.013% 11.063% 12374574 + [0.0599,0.1): 2.3129% +0.0550% 0.027% 10.881% 4486724 + [0.1,0.229): 2.4542% +0.0393% 0.124% 25.365% 5466556 + [0.229,0.769): 3.4700% +0.0166% 0.301% 24.638% 2239488 + [0.769,4.62): 5.6752% -0.0189% 0.615% 10.315% 223948 + [4.62,inf): 4.3682% -0.0469% 98.921% 17.738% 24884 +``` + +The FITS RGB channels are flattened into scalar samples: `sq-error energy` is +the fraction of `||fast-base||^2`, `linear-RGB sum` is the sum of baseline +linear-RGB sample values (proportional to total flux, not a per-source +photometry), and the base-value ranges partition all samples and sum to 100%. + +The top 0.01% of RGB channel samples account for 95.7% (`N=2`) or 94.5% +(`N=4`) of squared-error energy and about 9% of the baseline linear-RGB sum. +The flux-bearing faint ranges (`[1e-3,0.0599)`, `[0.0599,0.1)`, `[0.1,0.229)`, +together ~47% of the linear-RGB sum) have 2.3-4.9% local relative RMS and +signed bias of `+0.04%` to `+0.19%`; the many faint PSFs partially cancel. This +is consistent with the concentrated bright-core difference being predominantly a +bounded sub-pixel redistribution rather than a large net linear-RGB loss, and +explains why this tone-mapped still PNG remains smooth despite the global +relative L2. + +This comparison contains one independent frame only. It does not measure +temporal coherence: nearest deposition jumps by `1/N` output pixel per axis +when a moving image crosses a supersampled-cell boundary (`0.5 px` at `N=2`, +`0.25 px` at `N=4`). Both can produce visible movie flicker and are not +qualified here as visually faithful movie-output paths. + ## Reading - Steady-state FFTW is ~109x (N=2) and ~279x (N=4) faster than the spatial @@ -278,3 +498,12 @@ fast_n4.png 679f14a67a085f76c144ecbe36af51b3b2534c152e8ef59708a7884a3f7c8464 - FFTW-versus-spatial differences are floating-point roundoff (`max_abs` at the `1e-15` level against peaks near `0.5`), not crop, wrap, channel, or normalization errors. +- At the committed revision, the production-density 4K Galactic-center run + completes in 15.439 s at N=2 versus 68.610 s for the standard path (4.444x + end-to-end), with identical catalog image/deposit/filter counts. N=4 trades + 26.1% more wall time for materially lower spatial error. +- The 4K still-frame error is concentrated: the top 0.01% of RGB channel + samples hold ~95% of the squared-error energy but only ~9% of the baseline + linear-RGB sum, while the flux-bearing faint-value ranges stay at ~2-5% local + relative RMS with signed bias at or below ~0.2%. This does not establish + temporal fidelity. diff --git a/nr_spacetime_movie_renderer_design.md b/nr_spacetime_movie_renderer_design.md index 46e76fd..9aea84d 100644 --- a/nr_spacetime_movie_renderer_design.md +++ b/nr_spacetime_movie_renderer_design.md @@ -859,6 +859,32 @@ C2R，以及在 `(R,R)` 处的裁剪、`1/(Pwidth·Pheight)` 与 `1/N²` 归一 依赖与 plan 模式取舍见 `build.md` 与 [`benchmarks/fast_mode_fftw_2026-09-25.md`](benchmarks/fast_mode_fftw_2026-09-25.md)。 +fast mode 的单星精度由 deposit 模式与 `N` 决定：`nearest` 的格点间距是 +每轴 `1/N` 个输出像素，单帧瞬时舍入误差至多是 `1/(2N)`；在 +`--psf-fwhm-pixels` 不变时它与图像分辨率无关，只有增大 `N`（或使用保持 +质心、但会展宽 profile 的 `bilinear`）才会降低。 + +但单星偏移不等于一张稠密静态图的整图误差。生产密度 4K 银心单帧对照实测 +fast 对标准路径的 HDR relative L2 为 `N=2` 时 `10.44%`、`N=4` 时 `4.37%`， +线性 RGB 总和偏差 `+0.018%`。把 FITS 三通道展平后，按 `|fast-base|` 排序 +的前 `0.01%` RGB channel samples 占平方误差能量的 `95.7%`（`N=2`）或 +`94.5%`（`N=4`），其 baseline linear-RGB 总和约占 `9%`。不重叠的 base-value +分区中，承载 flux 的暗值区间（`[1e-3,0.06)`、`[0.06,0.1)`、`[0.1,0.23)`， +合计约 `47%` 的 linear-RGB 总和）给出约 `2-5%` local relative RMS 与不高于 +`0.19%` 的 signed bias。这一分布与大量暗星 PSF 的位置误差部分相消、已分辨 +亮星核心保留亚像素位移的解释一致。因此只能说这次高密度 4K 单帧的银河纹理 +明显好于全局 L2 数字所暗示，不能据此外推时间连续性或逐像素精度。 + +对 movie，平滑移动的星像跨越 nearest 格点边界时会跳变完整的 `1/N` 输出 +像素：`N=2` 每轴 `0.5 px`、`N=4` 每轴 `0.25 px`。即使在 4K，这种不连续 +跳变也可能形成可见闪烁，因此 nearest fast mode 尚未获得保持视频视觉特征的 +资格。`bilinear` 可消除这一质心跳变，但存在上述 profile 展宽，也尚未作为 +最终 movie 路径验收。定位上仍是近似路径：对单帧稠密星场，`N=2` 是吞吐 +默认、`N=4` 是更高保真预览；per-event 参考 PSF 路径仍是定量与 movie 基准。 +单帧误差分解可用 +[`scripts/fast_mode_error_decomposition.py`](scripts/fast_mode_error_decomposition.py) +复现。 + PSF 第一版可用 Gaussian； 以后可换成 Airy 或其他相机模型。 diff --git a/scripts/fast_mode_error_decomposition.py b/scripts/fast_mode_error_decomposition.py new file mode 100644 index 0000000..8ae448a --- /dev/null +++ b/scripts/fast_mode_error_decomposition.py @@ -0,0 +1,98 @@ +#!/usr/bin/env python3 +"""Decompose fast-mode-versus-standard HDR error by scalar RGB sample value. + +The unweighted whole-image relative L2 is dominated by a small number of very +bright PSF cores, so it is not a good proxy for how the dense faint-star +texture is reproduced. FITS RGB channels are flattened into scalar samples. +This tool reports, for each fast-mode FITS image: + + * the cumulative fraction of ||fast-base||^2 contributed by the largest |d| + RGB channel samples, alongside their share of the baseline linear-RGB sum; + * local relative RMS, signed bias, and share of base energy and linear-RGB + sum for a disjoint partition of base-value ranges. + +Usage: + python3 scripts/fast_mode_error_decomposition.py BASE.fits FAST.fits [...] +""" +import sys + +import numpy as np + + +def read_primary_fits(path): + with open(path, 'rb') as stream: + header = b'' + while True: + block = stream.read(2880) + if not block: + raise ValueError(f'{path}: truncated FITS header') + header += block + if any(block[i:i + 8] == b'END ' for i in range(0, 2880, 80)): + break + cards = {} + for i in range(0, len(header), 80): + card = header[i:i + 80].decode('ascii', 'replace') + key = card[:8].strip() + if key in ('NAXIS1', 'NAXIS2', 'NAXIS3', 'BITPIX', 'NAXIS'): + cards[key] = int(card[10:30]) + if key == 'END': + break + stream.seek(len(header)) + count = cards['NAXIS1'] * cards.get('NAXIS2', 1) * cards.get('NAXIS3', 1) + dtype = '>f4' if cards['BITPIX'] == -32 else '>f8' + return np.frombuffer(stream.read(count * np.dtype(dtype).itemsize), + dtype=dtype).astype(np.float64) + + +def decompose(name, base, fast, ranks=(1e-5, 1e-4, 1e-3, 1e-2, 1e-1)): + diff = fast - base + abs_diff = np.abs(diff) + base_l2_sq = float(np.dot(base, base)) + diff_l2_sq = float(np.dot(abs_diff, abs_diff)) + base_rgb_sum = float(base.sum()) + print(f'\n===== {name} ===== relative_L2={np.sqrt(diff_l2_sq / base_l2_sq) * 100:.4f}%') + order = np.argsort(abs_diff)[::-1] + print('top |d| RGB samples: sq-error energy base energy linear-RGB sum sample share') + for frac in ranks: + k = max(1, int(frac * base.size)) + idx = order[:k] + print(f' top {frac * 100:8.4f}%: {np.dot(abs_diff[idx], abs_diff[idx]) / diff_l2_sq * 100:7.3f}%' + f' {np.dot(base[idx], base[idx]) / base_l2_sq * 100:9.3f}%' + f' {base[idx].sum() / base_rgb_sum * 100:8.3f}%' + f' {k / base.size * 100:9.4f}%') + nonzero = base[base > 0] + quantiles = np.quantile(nonzero, [0.5, 0.9, 0.99, 0.999]) + edges = np.unique(np.concatenate( + ([-np.inf, 1e-6, 1e-3, 1e-1], quantiles, [np.inf]))) + print('base-value range (disjoint):' + ' rel_RMS signed bias base energy % linear-RGB sum % samples') + for low, high in zip(edges[:-1], edges[1:]): + mask = (base > low) & (base <= high) + if not mask.any(): + continue + bin_diff = diff[mask] + bin_base = base[mask] + rel = np.sqrt(np.dot(bin_diff, bin_diff) / np.dot(bin_base, bin_base)) + bias = bin_diff.sum() / bin_base.sum() + lo = '-inf' if low == -np.inf else f'{low:.3g}' + hi = 'inf' if high == np.inf else f'{high:.3g}' + print(f' [{lo},{hi}): {rel * 100:9.4f}% {bias * 100:+10.4f}%' + f' {np.dot(bin_base, bin_base) / base_l2_sq * 100:8.3f}%' + f' {bin_base.sum() / base_rgb_sum * 100:7.3f}% {int(mask.sum())}') + + +def main(argv): + if len(argv) < 3: + print(__doc__) + return 2 + base = read_primary_fits(argv[1]) + for path in argv[2:]: + fast = read_primary_fits(path) + if fast.shape != base.shape: + raise ValueError(f'{path}: shape {fast.shape} != base {base.shape}') + decompose(path.rsplit('/', 1)[-1], base, fast) + return 0 + + +if __name__ == '__main__': + sys.exit(main(sys.argv)) diff --git a/usage.md b/usage.md index 70f54d1..ef6c597 100644 --- a/usage.md +++ b/usage.md @@ -278,11 +278,12 @@ 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 +supersampled pixel. It keeps the PSF shape exactly on that grid. The grid +spacing is `1/N` output pixel per axis, while the instantaneous rounding error +is bounded by `1/(2N)` per axis. `--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 +(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 @@ -291,6 +292,42 @@ 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. +Its accuracy is governed by the deposit mode and `N`, not by the output +resolution. A nearest deposit shifts each image by at most `1/(2N)` output +pixels per axis. The separate two-dimensional centroid-distance experiment +measured worst cases of ~`0.30 px` at `N=2` and ~`0.12 px` at `N=4` for a FWHM +`2.7`, beta `4.5` PSF. These offsets remain fixed fractions of one output +pixel: with `--psf-fwhm-pixels` held fixed, a larger image does not reduce +them, and only a larger `N` (or the centroid-preserving bilinear deposit, at +the cost of a broadened profile) does. + +That per-star number is not how one dense still image behaves. In the measured +production-density 4K Galactic-center frame, fast-versus-standard HDR relative +L2 is `10.44%` at `N=2` and `4.37%` at `N=4`, with a `+0.018%` total +linear-RGB-sum offset. The difference is concentrated in a small set of scalar +RGB channel samples: the top `0.01%` ranked by `|fast-base|` account for +`95.7%` (`N=2`) or `94.5%` (`N=4`) of squared-error energy, while containing +about `9%` of the summed baseline linear-RGB values. The disjoint base-value +partition puts the flux-bearing faint ranges (`[1e-3,0.06)`, `[0.06,0.1)`, +`[0.1,0.23)`, together ~`47%` of the linear-RGB sum) at ~`2-5%` local relative +RMS with signed bias at or below `0.19%`. This pattern is consistent with +partial cancellation among many faint PSFs, while resolved bright cores retain +bounded sub-pixel displacement. It does not make fast mode pixelwise or +visually lossless. + +These measurements characterize independent still frames only. With nearest +deposition, a smoothly moving image jumps by one complete grid spacing when it +crosses a cell boundary: `1/N` output pixel per axis, or `0.5 px` at `N=2` and +`0.25 px` at `N=4`. Such discontinuities can appear as visible flicker even at +4K. Neither nearest mode is therefore qualified as a temporally coherent movie +output path. Bilinear deposition removes this particular centroid jump but has +the profile-broadening tradeoff above and has not been qualified as a final +movie path either. For single-frame dense-star previews, `N=2` remains the +throughput default and `N=4` the higher-fidelity point; the exact per-event +reference PSF path remains the quantitative and movie reference. The +single-frame error decomposition is reproducible with +[`scripts/fast_mode_error_decomposition.py`](scripts/fast_mode_error_decomposition.py). + In CPU builds the global convolution is a zero-padded FFTW linear convolution (`fftw3`/`fftw3_omp`, double precision). The immutable kernel spectrum, FFTW plans, and scratch buffers are built once when the accumulator is initialized