Files
GR-raytracing/psf_lookup_optimization_plan.md

136 lines
6.7 KiB
Markdown

# CPU PSF lookup-table optimization plan
## Goal
Speed up the current CPU point-source rendering path when the physically
correct lens map produces very many star images. This plan changes only how a
known image's Moffat PSF is evaluated and accumulated. It does **not** change
the catalog representation, inverse lens map, number of images, frequency
shift, magnification, adaptive-mesh work, or use of the GPU.
The target is to replace the per-covered-pixel Moffat `pow()` evaluation with
a lookup of the same continuously translated PSF. It must preserve continuous
sub-pixel image positions and have an explicit, testable bound on PSF-tail
error, including for bright sources.
## Interface and ownership
Add a `PsfKernelCache` owned by the rendering settings for the lifetime of one
program invocation. Exactly one cache is built after command-line parsing,
using the process-wide values of:
- `--psf-fwhm-pixels`;
- `--psf-moffat-beta`;
- a named PSF-tail/display-error configuration;
- a fixed sub-pixel phase resolution.
It is immutable after construction, shared read-only by splat workers, and is
destroyed after rendering. It is not serialized to disk and no multi-parameter
cache or LRU is required: the current program has one PSF parameter pair for
the entire image or movie.
`splat_moffat()` gains a cache-aware path (or is replaced by an equivalently
named cache-aware function). A direct-evaluation implementation remains
available as a reference and as a fallback if cache construction fails.
## Kernel construction
For each sub-pixel phase `(fx, fy)`, where
`fx = x - floor(x)` and `fy = y - floor(y)`, build a scalar kernel containing
the **pixel-area integral** of the normalized circular Moffat profile over each
integer pixel footprint. Kernel weights are independent of color and flux.
Use a fixed, documented numerical quadrature whose convergence is verified
against a higher-accuracy reference.
At runtime, retain the image event's full-precision `(x, y)`. Select the four
neighboring phase tables and bilinearly interpolate their weights. Do not
round to a nearest phase table. Normalize every interior phase kernel so its
weight sum is one; an image clipped by the framebuffer boundary retains the
current behavior of losing the off-frame contribution rather than being
renormalized.
Start with a 64 by 64 phase grid. This is an initial implementation value,
not an accuracy claim: increase it if the phase-continuity and reference-image
tests below do not meet their thresholds.
## Tail support and bright sources
A Moffat profile has infinite support, so a finite renderer needs an explicit
finite-support approximation. Do not rely solely on the current fixed
relative omitted-flux fraction for every source. For each image event, choose
the support radius from its final post-lens, post-frequency-shift, post-exposure
flux so that both of these hold:
- omitted integrated flux is below a configured absolute frame error budget;
- the PSF value at the support boundary is below a configured per-pixel HDR
error budget.
The cache contains weights through the radius implied by the selected PSF
parameters at the reference scaled flux of one. A dim image uses a cropped
subset; an unusually bright image uses a larger subset. If an event requires
a radius beyond that cache maximum, use the direct Moffat evaluator for that
event and record it in diagnostics. The radius is not a fixed pixel cap:
low-beta Moffat profiles need broad wings even when their FWHM is small.
This prevents a visible hard cutoff in bright-star wings without silently
discarding their flux.
The selected tail budgets, maximum cache radius, and fallback count must be
reported with the render diagnostics. Their final numerical values are to be
chosen from the validation experiments, not inferred from a performance goal.
## Accumulation and parallel behavior
For every image event, compute color and final flux exactly as in the existing
path. Then multiply the scalar cached PSF weights by that color and flux and
add them to the worker's existing private HDR buffer. Preserve the current
private-HDR ownership, reduction ordering, and serial fallback. Analytic
backends use all OpenMP render workers; a future numerical backend must opt in
to the private-HDR memory cap when budgeting its metric slabs and evaluator
workspaces;
this plan does not introduce atomics, locks, or a different parallel
partitioning.
## Validation and benchmark
Add tests and a reproducible benchmark before making the lookup path default:
1. Compare cached and high-accuracy direct pixel-area integration for a grid
of FWHM, beta, phase, flux, and support-radius cases. Check kernel sum,
per-pixel error, and omitted-flux bound.
2. Translate a star in sub-pixel increments across pixel boundaries. Verify
continuous image values and absence of phase quantization jumps.
3. Test bright stars whose required support is larger than the ordinary cache
crop, including the direct-fallback case. Verify the configured tail-error
bound and no visible-radius discontinuity relative to the reference.
4. Compare serial and OpenMP outputs under the existing private-HDR contract.
Preserve the current expected identity/tolerance policy explicitly.
5. Benchmark the Schwarzschild 1080p, 60-degree, FWHM 1.0 all-sky case and
separately report cache-build time, number of image events, cached-splat
time, direct-fallback count, and total render time. Compare output against
the direct evaluator with the same error budgets.
The lookup path becomes the default only after it meets the agreed visual and
numerical tolerances and demonstrates a material speedup in the PSF-dominated
case. The direct evaluator remains a selectable regression reference.
## Implementation record (2026-08-28)
Implemented as one immutable `PsfKernelCache` per invocation, built after CLI
parsing and passed read-only through the existing frame splat workers. The
cache uses 64 by 64 phase intervals (65 nodes on each axis), 4-point
Gauss-Legendre pixel-area quadrature, bilinear phase interpolation, and a
static OpenMP phase-plane construction loop. The
`--psf-direct` regression mode uses an independent 8-point quadrature
reference; events whose requested support exceeds the cache use that path
automatically.
The selected HDR budgets are omitted flux `1e-6` and boundary contribution
`1e-7`; the former `1e-8` relative tail remains part of the support selection.
The cache reports build time, maximum radius, cached splat count, and direct
fallback count. `tests/test_frame.c` covers cache construction, phases near
pixel boundaries, comparison with the direct reference, and bright-event
fallback. The float-HDR default-PSF acceptance result and the 1080p full-sky
performance result are recorded in
`benchmarks/2mass_all_sky_psf_lookup_2026-08-28.md`.