GR-raytracing/fast_mode_fftw_optimization_plan.md

Fast-mode FFTW convolution optimization plan

1. Status and objective

This plan starts from commit 3deebfb (Feat: Add --fast-mode supersampled point-source accumulation). The current CPU fast mode already:

• maps catalog images to continuous image positions;
• deposits each image into one nearest supersampled cell by default, or four bilinear cells when explicitly requested;
• stores a supersampled double-RGB impulse buffer;
• applies one global pixel-integrated Moffat kernel;
• averages each N x N supersampled block into the final HDR framebuffer;
• preserves additive HDR semantics, --psf-min-y, wing-clipping statistics, and the documented preview-only boundary.

The remaining bottleneck is fast_psf_accumulator_resolve(), which evaluates the global convolution by nested spatial loops. The earlier HIP work already established PSF accumulation as the dominant stage, and the target machine is the local 128 GB host. FFTW is an accepted dependency. This work therefore does not repeat a feasibility study or reconsider the selected nearest-cell visual tradeoff. Its objective is:

Replace the spatial global convolution with a reusable, multithreaded FFTW linear-convolution path while reproducing the current discrete fast-mode result to floating-point tolerance and retaining the spatial implementation as a test/benchmark reference.

The target is the CPU PSF backend. HIP and dummy PSF builds continue to reject --fast-mode and must not acquire an unnecessary FFTW link dependency.

2. Fixed requirements and non-goals

2.1 Required behavior

• Use double-precision FFTW (fftw3), not the float API.
• Use the locally available FFTW OpenMP backend (fftw3_omp); the current host reports FFTW 3.3.10 and pkg-config --libs fftw3_omp returns -lfftw3_omp -lfftw3.
• On Gentoo this requires sci-libs/fftw[openmp], which is already enabled on the target host. It does not require the separate threads USE flag: threads selects FFTW's pthread backend, while this plan deliberately uses the OpenMP backend so it shares the renderer's OpenMP runtime.
• Compute zero-padded linear convolution. Circular wraparound is a failure.
• Match the current circular kernel traversal, not the unused square corners of the allocated weights array.
• Preserve the current impulse-buffer contents so the same buffer can be sent through the spatial reference and FFTW implementations in one process.
• Normalize the unscaled FFTW inverse exactly once.
• Apply the existing 1/N^2 box-average normalization exactly once.
• Crop the full linear convolution at the mathematically correct offset before downsampling.
• Add the result to the caller's HDR buffer; never overwrite existing HDR.
• Cache plans, the kernel spectrum, FFT buffers, and dimension metadata across movie frames.
• Keep FFTW planning outside OpenMP producer regions and avoid nested OpenMP/FFTW oversubscription.
• Report one-time planning/kernel-transform time separately from per-frame FFT execution time.
• Propagate allocation or plan failures as explicit fast-mode initialization failures. Do not silently fall back to a potentially multi-minute spatial convolution.

2.2 Preserved preview semantics

This optimization must not change:

• nearest as the default deposit mode;
• bilinear as the explicit position-over-shape alternative;
• supersample factor or coordinate convention;
• the Moffat FWHM, beta, quadrature, normalization, or support radius;
• flux, color, magnification, frequency-shift, exposure, or min-Y decisions;
• fixed global-kernel wing clipping and its statistics;
• zero-outside-image behavior, including the existing narrow bilinear edge approximation;
• tone mapping, PNG/FITS output, lens-map reuse, or catalog traversal.

2.3 Non-goals

• Do not redesign the producer/catalog path.
• Do not introduce a GPU FFT backend in this change.
• Do not add polyphase convolution as a competing production algorithm.
• Do not change nearest/bilinear acceptance criteria.
• Do not add an artificial memory cap for the 128 GB target host.
• Do not make bitwise identity an acceptance criterion; FFT summation order is different from direct spatial accumulation.
• Do not combine this work with impulse-buffer persistence or PSF/exposure sweep orchestration. The internal state should permit those later, but they are a separate feature.

3. Exact discrete operation to preserve

Let the supersampled impulse buffer be

B[c][y][x], c in {R,G,B}, 0 <= x < Wss, 0 <= y < Hss

with Wss = N * width and Hss = N * height. Let the stored kernel radius be R, its side be K = 2R + 1, and its circularly retained weights be

A[dy + R][dx + R] = weights[(dy + R) * K + (dx + R)]

only where

-R <= dy <= R
|dx| <= floor(sqrt(R^2 - dy^2)).

All other elements of A are zero. The current fine-grid result is

C[c][sy][sx] = sum_dy sum_dx A[dy + R][dx + R]
                                  * B[c][sy - dy][sx - dx],

where samples of B outside its image are zero. The final operation is

hdr[c][row][column] +=
    (1 / N^2) * sum_j=0..N-1 sum_i=0..N-1
        C[c][N*row + j][N*column + i].

The FFTW implementation must preserve this discrete operation, including the current float kernel weights promoted to double during multiplication.

4. Linear-convolution layout and crop convention

This is the most error-prone part of the implementation and must be fixed by tests before performance work.

4.1 Padding bounds

Choose FFT dimensions satisfying

Pwidth  >= Wss + K - 1 = Wss + 2R
Pheight >= Hss + K - 1 = Hss + 2R.

B is copied at padded origin (0,0). The kernel is also stored at padded origin with its array indices unchanged:

padded_kernel[dy + R][dx + R] = A[dy + R][dx + R].

Do not apply an additional fftshift/ifftshift in this convention. The ordinary full convolution produced by FFT multiplication has dimensions Wss + K - 1 by Hss + K - 1, and the desired same-size fine-grid result is

C[sy][sx] = full_convolution[sy + R][sx + R].

Thus the fine-grid crop begins exactly at (R,R) and has extent Wss x Hss.

An alternative center-at-frequency-origin layout is mathematically possible, but it changes the crop/wrap convention and is deliberately excluded from the first implementation.

4.2 FFT-friendly dimensions

Add a checked helper such as next_smooth_size() that returns the first value at least the required extent whose prime factors are limited to small FFTW- friendly factors, initially {2,3,5,7}. The helper must:

• accept and return size_t while checking overflow;
• reject a selected dimension above INT_MAX before passing it to FFTW's int n[2] interface;
• have unit tests for exact, next-size, and near-overflow inputs;
• record both the minimum linear-convolution extent and selected FFT extent.

Do not hard-code powers of two: they can waste substantially more memory and work than a nearby smooth composite size.

5. Source and module boundaries

5.1 New private FFTW module

Keep FFTW details out of the public optics API and out of unrelated backends. Add:

src/fast_psf_fftw.c
src/fast_psf_fftw.h

fast_psf_fftw.h is a private renderer header. It declares an opaque state and functions conceptually equivalent to:

typedef struct FastPsfFftwState FastPsfFftwState;

FastPsfFftwState *fast_psf_fftw_create(
    int ss_width, int ss_height,
    int final_width, int final_height,
    int supersample, int kernel_radius,
    const float *weights,
    int fft_workers,
    FastPsfFftwTiming *timing);

int fast_psf_fftw_resolve(
    FastPsfFftwState *state,
    const double *interleaved_ss_rgb,
    double *interleaved_hdr_rgb,
    FastPsfFftwTiming *timing);

void fast_psf_fftw_destroy(FastPsfFftwState *state);

Exact names may follow project style, but ownership must stay explicit:

• FastPsfAccumulator owns one FFTW state pointer;
• the FFTW state owns every FFTW allocation and plan;
• the input impulse buffer remains owned by FastPsfAccumulator;
• the output HDR remains caller-owned;
• the kernel spectrum is immutable after initialization;
• destruction is safe for a zero/partially initialized state.

Avoid exposing fftw_plan or fftw_complex through optics.h. If the public FastPsfAccumulator must carry the private pointer, use a forward-declared opaque struct.

5.2 Spatial reference path

Extract the current nested-loop resolver into a clearly named reference helper, for example:

fast_psf_accumulator_resolve_spatial_reference(...)

It must remain callable by a dedicated FFTW test/benchmark with the identical impulse buffer. Production fast_psf_accumulator_resolve() dispatches to FFTW in CPU builds after validation. The spatial path is not a silent runtime fallback.

The reference helper should reuse a cached circular row_span array rather than allocating it on every call. This small cleanup also guarantees that the spatial and FFTW kernel masks are built from one definition.

6. FFTW state and memory layout

6.1 Planar FFT scratch

Keep the existing interleaved impulse buffer because the deposit hot path and its atomic RGB updates already use that layout. At resolve time, pack it into FFTW-aligned planar scratch:

real_rgb[channel][padded_y][padded_x]
frequency_rgb[channel][ky][kx], kx extent = Pwidth/2 + 1
kernel_frequency[ky][kx]

Use:

• fftw_alloc_real() for real scratch;
• fftw_alloc_complex() for channel spectra and the kernel spectrum;
• out-of-place R2C and C2R transforms in the first implementation;
• fftw_plan_many_dft_r2c() and fftw_plan_many_dft_c2r() for three planar RGB transforms with unit stride and per-plane distance;
• one shared kernel spectrum multiplied into each of the three channel spectra.

Out-of-place planar buffers are intentionally preferred over an in-place or strided-interleaved first version. The target host has ample memory, while the simpler layout reduces alignment, padding, and plan-many mistakes.

6.2 Checked allocation sizes

Before allocation, check every product used for:

3 * Pheight * Pwidth * sizeof(double)
3 * Pheight * (Pwidth/2 + 1) * sizeof(fftw_complex)
Pheight * (Pwidth/2 + 1) * sizeof(fftw_complex)

Allocation failure is a hard fast-mode initialization error with a diagnostic that includes requested dimensions and byte counts. These numbers are logged for provenance, not used as a policy cap.

6.3 Packing

For each frame:

1. zero the full padded planar real scratch;
2. copy the active Hss x Wss region from interleaved RGB into three planes;
3. leave all right and bottom padding zero;
4. keep the original impulse buffer unchanged.

The zero and pack loops may use one ordinary OpenMP parallel-for region. They must finish before FFTW execution starts.

7. Kernel transform

Build the padded real kernel once after plans and scratch buffers exist:

1. clear a single padded real plane;
2. for dy = -R..R, compute the same circular row_span as the spatial reference;
3. copy only dx = -row_span..row_span to index (dy+R, dx+R);
4. leave square corners and all FFT padding zero;
5. execute one R2C kernel transform;
6. store its spectrum in immutable kernel_frequency;
7. compute the reported retained-flux sum from this same circular mask.

Do not assume the kernel spectrum is purely real. The unshifted full- convolution layout places the kernel center at (R,R), so the spectrum generally has a phase. Complex multiplication must use both real and imaginary components:

(ar + i ai) * (br + i bi)
  = (ar*br - ai*bi) + i(ar*bi + ai*br).

The one-time kernel transform and plan creation must not mutate the production impulse buffer.

8. Per-frame FFTW resolve

The production resolver executes these stages in order:

1. Zero and pack interleaved impulses into padded planar real scratch.
2. Forward FFT all three planes with the cached R2C plan.
3. Frequency multiply each RGB spectrum by the shared kernel spectrum.
4. Inverse FFT all three planes with the cached C2R plan.
5. Crop, normalize, and box-downsample into the caller HDR.

FFTW's transforms are unnormalized. Let

fft_scale = 1 / (Pwidth * Pheight)
box_scale = 1 / (N * N).

For output pixel (row,column) and fine offsets (j,i), read inverse scratch at

y = R + N*row    + j
x = R + N*column + i.

Then accumulate

hdr += inverse[y][x] * fft_scale * box_scale.

Fuse crop, inverse normalization, box summation, RGB interleaving, and HDR addition in one OpenMP parallel-for over final output rows. Do not materialize a separate cropped fine-grid image.

The frequency multiply is embarrassingly parallel over frequency bins and may use OpenMP if measurement shows it is not already hidden by FFT cost. It must not run inside an active FFTW OpenMP region.

9. FFTW planning, threads, and wisdom

9.1 Process-level threading lifecycle

The local host provides libfftw3_omp. Use the FFTW threading API exported by that OpenMP library:

fftw_init_threads();
fftw_plan_with_nthreads(fft_workers);

Requirements:

• initialize FFTW threading before plan creation;
• create and destroy plans only from the serial control thread;
• never create a plan inside the catalog OpenMP region;
• do not call fftw_cleanup_threads() while any accumulator or plan exists;
• keep process-global FFTW runtime ownership/reference counting in one private module rather than letting individual accumulators independently clean up global state;
• choose fft_workers from the same OpenMP worker policy already used by the renderer and report it;
• do not set a nested OpenMP region around fftw_execute();
• after FFTW returns, use a separate OpenMP region for crop/downsample.

OMP_DYNAMIC and the user's OpenMP environment remain authoritative. Do not silently force a different global thread count.

9.2 Planning flags

Implement in two steps:

1. correctness and unit tests use FFTW_ESTIMATE so tests are fast and do not depend on host wisdom;
2. the release benchmark compares FFTW_ESTIMATE with FFTW_MEASURE, recording plan time separately from steady-state execution.

Select the release default only from that measured comparison. If FFTW_MEASURE materially improves repeated movie-frame execution, add optional wisdom import/export in a follow-up or in the same implementation phase:

• wisdom is host/FFTW-version/dimension specific;
• a missing or incompatible wisdom file must be explicit in verbose output;
• wisdom files are local artifacts and are not committed as portable project data;
• planning may destroy scratch contents, so plans are created before the kernel and frame inputs are populated.

Do not mix plan time into the per-frame convolution metric.

10. Build-system integration

10.1 Dependency detection

For PSF_BACKEND=cpu, require:

pkg-config --exists fftw3_omp
pkg-config --modversion fftw3_omp
pkg-config --cflags --libs fftw3_omp

The current host returns version 3.3.10 and libraries -lfftw3_omp -lfftw3.

On Gentoo the package requirement is sci-libs/fftw[openmp]; do not require sci-libs/fftw[threads] and do not link libfftw3_threads. The similarly named FFTW threading API is implemented by both backends, but exactly one backend library should provide it in this executable.

Add CPU-only build variables such as:

FFTW_CPPFLAGS := $(shell pkg-config --cflags fftw3_omp)
FFTW_LDLIBS := $(shell pkg-config --libs fftw3_omp)

and a CPU-only definition such as -DFAST_PSF_FFTW. Fail early with a clear Make error if fftw3_omp is unavailable in a CPU build.

Do not append FFTW flags to HIP or dummy backend links. CPU tests that compile src/optics.c and the new FFTW module must receive the same FFTW compile/link flags.

The current Makefile obtains common sources through a src/*.c wildcard. Do not let that wildcard pull the FFTW module into HIP and dummy binaries. Exclude src/fast_psf_fftw.c from COMMON_SOURCES and append it only to the CPU renderer and CPU-test source lists. This makes the dependency boundary visible. Guard calls from src/optics.c consistently so HIP and dummy links have no unresolved FFTW-module symbols.

Audit at least:

• renderer backend link;
• test_frame and the other CPU tests using COMMON_SOURCES;
• the new FFTW unit test;
• make_psf_fixture if it links the FFTW-enabled optics object;
• capture/replay tools, ensuring HIP-only tools do not accidentally gain a CPU FFTW requirement.

Keep dependency tracking (-MMD -MP) for the new module.

10.2 Build identities

The CPU build always includes FFTW after this change, so a separate executable suffix is not required. If a temporary compile-time A/B switch is introduced, include it in the object-directory tag to prevent stale-object reuse. Prefer a single CPU binary with FFTW production resolution and a test-only spatial reference over persistent compile-time variants.

11. Diagnostics and timing

Extend fast-mode reporting with one-time state:

Fast FFTW: linear min=<Hmin>x<Wmin>, fft=<Pheight>x<Pwidth>,
           workers=<n>, plan=<estimate|measure>, plan=<seconds>,
           kernel_fft=<seconds>, scratch=<bytes>

Under --verbose, report per-frame stages:

Fast FFTW frame: zero_pack=<s>, forward=<s>, multiply=<s>,
                 inverse=<s>, crop_downsample=<s>, total=<s>

Timing rules:

• use wall time, not sums of worker-local time;
• plan/kernel setup is one-time and separate;
• total spans zero/pack through completed HDR addition;
• keep catalog prefetch/deposit timing outside the FFTW resolve metric;
• include selected dimensions, supersample factor, FWHM/beta, kernel radius, FFTW version, planning flag, worker environment, and binary revision in benchmark records.

These timings validate and attribute the implementation; they are not another feasibility gate for fast mode.

12. Correctness test plan

12.1 Dedicated FFTW-versus-spatial unit test

Add tests/test_fast_psf_fftw.c. It constructs one accumulator/impulse buffer, runs both resolvers without re-depositing, and compares their final HDR arrays. Cover:

• non-power-of-two image sizes;
• N = 1, 2, 3, 4;
• nearest and bilinear deposit modes;
• one white center impulse;
• distinct R/G/B values to catch channel-layout mistakes;
• impulses near all four edges and four corners;
• an impulse exactly on a supersampled cell center;
• phases immediately on both sides of a nearest-cell boundary;
• multiple separated impulses;
• many impulses in one supersampled cell;
• a deterministic dense pseudo-random impulse field;
• an empty impulse buffer;
• a nonzero prefilled HDR background;
• a kernel radius that makes the selected FFT padding larger in both axes;
• dimensions where only one axis needs the next smooth FFT size.

For each case compute:

• maximum absolute RGB error;
• maximum relative error above a reference-magnitude floor;
• RMS error;
• total per-channel flux difference;
• worst sample coordinate and channel;
• NaN/Inf count.

Use provisional double-precision acceptance limits of:

max_abs <= 1e-10 * max(1, reference_peak)
max_rel <= 1e-9 for |reference| above 1e-12 * reference_peak
relative total-flux error <= 1e-10

Tighten these after observing stable results on the target host; do not loosen them merely to hide crop, wrap, channel, or normalization errors. Any coherent edge band, shifted peak, factor-of-area scaling, or opposite-edge ghost is an algorithmic failure regardless of aggregate tolerance.

12.2 Analytic impulse checks

Before relying only on spatial comparison, add exact index tests with a single unit impulse:

• verify the full-convolution peak appears at impulse index plus R;
• verify the cropped fine-grid center returns to the original impulse index;
• verify no signal appears at the opposite edge;
• verify a unit impulse's cropped result equals the circular kernel samples;
• verify box-downsample indices for every phase 0..N-1 in both axes;
• verify RGB channels do not cross-contaminate.

These tests diagnose layout errors that a broad image statistic can obscure.

12.3 Existing regression suite

After FFTW becomes the production resolver:

make -j4 BUILD_TYPE=Debug PSF_BACKEND=cpu test

must pass. In particular:

• the existing fast-mode snapped-kernel regression now exercises FFTW;
• the additive-HDR regression must still preserve its 0.25 background;
• ordinary non-fast FITS references remain unchanged;
• Minkowski and Schwarzschild camera/lens-map regressions remain unchanged.

Add CLI smoke coverage for --fast-mode in a CPU build. HIP and dummy builds must continue to reject it with the existing message.

12.4 Scientific/visual metrics

Using identical impulses, compare FFTW and spatial fast-mode outputs for:

• centroid;
• fitted FWHM;
• fitted Moffat beta relative to the same direct-reference fitting procedure;
• total linear-HDR flux;
• center, edge, and corner crops.

The FFTW change should add only floating-point roundoff to the already measured nearest/bilinear approximation. It must not change the deposition-error table or the rationale for nearest as default.

13. Performance benchmark plan

13.1 Resolver-only benchmark

Add tests/benchmark_fast_psf_fftw.c or an equivalent bounded benchmark target that:

• initializes one accumulator;
• deposits or directly fills one deterministic impulse buffer once;
• runs spatial and FFTW resolution on the identical buffer;
• separates first-use plan/setup from repeated execution;
• performs warm-up before measured FFT executions;
• records at least five steady-state executions when runtime permits;
• reports median, minimum, and all raw samples;
• hashes the impulse buffer and both HDR outputs;
• records RSS and FFTW dimensions for provenance.

Required cases:

1. current default PSF, N=2, 320x240;
2. current default PSF, N=4, 320x240;
3. a medium frame that is large enough for FFT execution to dominate startup;
4. the intended 4K N=2 target;
5. 4K N=4 if the user requests that mode to be production-supported.

The 4K cases are bounded synthetic/identical-buffer resolver tests, not another catalog or geodesic run.

13.2 End-to-end comparison

Repeat the recorded 320x240 Galactic-center command from benchmarks/fast_mode_cpu_2026-09-18.md with the same catalog, geometry, exposure, supersample factor, deposit mode, worker environment, and revision metadata. Preserve raw terminal output rather than only a timing table.

Do not claim an end-to-end speedup unless inputs, event/image counts, PSF parameters, and output mode match. Do not add worker-summed times to process wall time.

An expensive full-sky or production 4K catalog render is not automatically run as part of implementation. It requires separate explicit authorization and must use a frozen lens map/input provenance when used for final evidence.

13.3 Performance acceptance

The FFTW implementation is ready to become the default resolver when:

• all correctness tests pass;
• the resolver-only benchmark shows a clear steady-state reduction for default N=2 and the N=4 case that exposed spatial scaling;
• plan/setup time is reported separately and is acceptable for the intended single-frame/movie usage, or wisdom support addresses it;
• the recorded end-to-end bounded command does not regress catalog/deposit behavior or image count;
• FFTW output differences remain within the linear-HDR tolerances and do not change measured centroid/FWHM/beta beyond roundoff.

No new arbitrary global wall-time target is introduced by this plan.

14. Failure handling and cleanup

• Reject invalid/overflowed FFT dimensions before allocation or plan creation.
• Treat fftw_alloc_* == NULL or a null plan as initialization failure.
• Print which allocation/plan failed, dimensions, plan mode, and worker count.
• Never continue with a partially initialized kernel spectrum.
• Never silently switch to spatial convolution after an FFTW failure.
• Destroy forward/inverse/kernel plans before freeing their buffers.
• Release every FFTW allocation in the partial-initialization error path.
• Ensure movie failure exits and imported-lens-map early returns destroy the FFTW state through the existing accumulator cleanup.
• Do not call global FFTW cleanup while another accumulator remains alive.
• Keep fast_psf_accumulator_destroy() idempotent on a zero state.

15. Implementation sequence

Phase 1: Freeze the spatial reference

• Extract the existing resolver without changing its arithmetic.
• Cache/reuse its circular row spans.
• Add a small digest/reference test for the current spatial output.
• Add checked smooth-size and product helpers with unit tests.

Exit criterion: current fast-mode tests and the spatial output digest pass.

Phase 2: Add FFTW build and private state

• Add CPU-only fftw3_omp detection and link flags.
• Add the private FFTW module and opaque ownership.
• Allocate planar real/frequency scratch.
• Initialize process-level FFTW threading.
• Create FFTW_ESTIMATE plans and build the circular kernel spectrum.
• Report plan dimensions and setup time.

Exit criterion: all CPU/HIP/dummy build variants still compile with the correct dependency boundary; state construction/destruction passes allocation tests.

Phase 3: Implement exact FFT resolution

• Implement zero/pack, batched forward FFT, shared-kernel multiplication, batched inverse FFT, (R,R) crop, normalization, box downsample, and additive HDR write.
• Add analytic impulse/crop tests first.
• Add dense FFTW-versus-spatial comparisons.

Exit criterion: dedicated correctness tests pass for all required dimensions, phases, channels, and edge cases.

Phase 4: Integrate the production path

• Make fast_psf_accumulator_resolve() use FFTW in CPU fast mode.
• Retain the spatial helper only for tests/benchmarks.
• Add CLI smoke coverage and verbose timing.
• Run the full CPU Debug suite.

Exit criterion: existing behavior and references pass, and additive HDR semantics remain covered.

Phase 5: Benchmark planning and execution

• Run the resolver-only matrix.
• Compare FFTW_ESTIMATE and FFTW_MEASURE.
• Decide whether wisdom support is required.
• Repeat the bounded Galactic-center command with full provenance.
• Save raw logs and output hashes.

Exit criterion: correctness remains within tolerance and FFTW is demonstrably faster than the spatial resolver for the default and N=4 cases.

Phase 6: Documentation finalization

• Update usage.md with FFTW dependency and planning/wisdom behavior.
• Update nr_spacetime_movie_renderer_design.md to state that the global fast convolution is implemented as zero-padded FFTW linear convolution.
• Add a dated benchmark record with commands, revision, environment, raw logs, hashes, plan/setup time, and steady-state execution samples.
• Keep the existing deposition-error record unchanged except for links to the FFTW validation record; FFTW does not change nearest/bilinear errors.

16. Review checklist

Before considering the implementation complete, review the following explicitly:

• FFTW is linked only for the CPU PSF backend and CPU tests.
• Local dependency is fftw3_omp, version and libraries are reported.
• Padding satisfies Wss + 2R and Hss + 2R before smooth rounding.
• Kernel square corners outside the circular support are zero.
• Kernel is stored unshifted at indices (dy+R, dx+R).
• Fine-grid crop begins at (R,R).
• Inverse normalization uses 1/(Pwidth*Pheight) exactly once.
• Box normalization uses 1/N^2 exactly once.
• Output uses += into HDR.
• RGB batch layout and distances are covered by distinct-channel tests.
• Empty input produces no HDR change.
• No opposite-edge wraparound appears.
• Original impulse buffer is unchanged by resolve.
• Plans and kernel spectrum are reused across frames.
• Planning occurs outside OpenMP producer regions.
• FFTW and renderer OpenMP stages do not nest/oversubscribe.
• Partial initialization and all early exits free FFTW state.
• Spatial reference and FFTW consume the identical impulse buffer.
• Linear-HDR errors, flux, centroid, FWHM, and beta are recorded.
• Plan time and steady-state execution time are reported separately.
• Benchmark commands, revisions, inputs, event counts, and raw output are preserved.
• No expensive full-sky/4K catalog render is run without explicit authorization.

17. Expected final architecture

catalog / inverse lens map / g / blackbody / flux
                       |
                       v
        atomic nearest or bilinear SS delta deposit
                       |
                       v
          interleaved double-RGB impulse buffer
                       |
                       v
       zero + pack to FFTW planar padded real arrays
                       |
                       v
             batched RGB R2C forward FFT
                       |
                       v
        multiply by one cached Moffat kernel spectrum
                       |
                       v
             batched RGB C2R inverse FFT
                       |
                       v
       crop at (R,R) + normalize + N x N downsample
                       |
                       v
             additive double-RGB HDR framebuffer

The physical/catalog decisions remain on the CPU producer side; FFTW replaces only the global convolution implementation inside the explicit preview path.