Feat: Use FFTW linear convolution for CPU fast-mode PSF resolve

Replace the nested spatial global convolution in fast_psf_accumulator_resolve with a reusable double-precision FFTW linear convolution on the CPU PSF backend.

- Add the private src/fast_psf_fftw.{c,h} module: zero-padded R2C/C2R plans, cached kernel spectrum and planar scratch, exact 1/(Pwidth*Pheight) and 1/N^2 normalization, (R,R) crop, and additive HDR output.
- Keep the previous nested loops as fast_psf_accumulator_resolve_spatial_reference for tests/benchmarks only; it is not a runtime fallback.
- Cache the circular row spans on FastPsfAccumulator and report one-time plan, kernel transform, scratch, and per-frame stage timings.
- Require fftw3_omp for CPU builds; HIP and dummy builds do not link FFTW.
- Namespace test/helper binaries by spacetime and build tag, and reject make test / psf-capture for non-CPU backends.
- Add tests/test_fast_psf_fftw.c (FFTW versus spatial), tests/benchmark_fast_psf_fftw.c, an FFTW CLI smoke check, and the 2026-09-25 benchmark record.

src/fast_psf_fftw.c

@@ -0,0 +1,393 @@
+#include "fast_psf_fftw.h"
+
+#include <fftw3.h>
+#include <limits.h>
+#include <math.h>
+#include <omp.h>
+#include <stdint.h>
+#include <stdlib.h>
+#include <string.h>
+
+/* Process-global FFTW threading lifecycle.  Plans are created only from the
+ * serial control thread, and global cleanup happens only after the last live
+ * state is gone.  This keeps the individual states from racing on FFTW's
+ * process-wide planner state. */
+static int g_threading_initialized = 0;
+static int g_live_states = 0;
+static int g_plan_measure = 0;
+
+void fast_psf_fftw_set_plan_mode(int measure)
+{
+  g_plan_measure = measure ? 1 : 0;
+}
+
+static unsigned plan_flags(void)
+{
+  return g_plan_measure ? FFTW_MEASURE : FFTW_ESTIMATE;
+}
+
+struct FastPsfFftwState {
+  int fft_width, fft_height;
+  int ss_width, ss_height;
+  int final_width, final_height;
+  int supersample;
+  int kernel_radius;
+  int workers;
+  int registered;
+  int plan_measure;
+  size_t scratch_bytes;
+  double setup_seconds;
+  double kernel_seconds;
+  double fft_scale;
+  double box_scale;
+
+  double *real_rgb;          /* 3 planes of fft_height * fft_width */
+  fftw_complex *freq_rgb;    /* 3 planes of fft_height * (fft_width/2 + 1) */
+  fftw_complex *kernel_freq; /* one plane of fft_height * (fft_width/2 + 1) */
+  fftw_plan forward_plan;
+  fftw_plan inverse_plan;
+};
+
+static int checked_mul_size(size_t a, size_t b, size_t *out)
+{
+  if (a != 0 && b > SIZE_MAX / a)
+    return -1;
+  *out = a * b;
+  return 0;
+}
+
+int fast_psf_fftw_next_smooth_size(size_t min_extent, size_t *out)
+{
+  if (out == NULL || min_extent == 0)
+    return -1;
+  if (min_extent > (size_t)INT_MAX)
+    return -1;
+  static const unsigned factors[] = {2, 3, 5, 7};
+  for (size_t candidate = min_extent; candidate <= (size_t)INT_MAX;
+       ++candidate) {
+    size_t remaining = candidate;
+    for (size_t i = 0; i < sizeof factors / sizeof factors[0]; ++i)
+      while (remaining % factors[i] == 0)
+        remaining /= factors[i];
+    if (remaining == 1) {
+      *out = candidate;
+      return 0;
+    }
+  }
+  return -1;
+}
+
+void fast_psf_fftw_destroy(FastPsfFftwState *state)
+{
+  if (state == NULL)
+    return;
+  if (state->forward_plan != NULL)
+    fftw_destroy_plan(state->forward_plan);
+  if (state->inverse_plan != NULL)
+    fftw_destroy_plan(state->inverse_plan);
+  if (state->real_rgb != NULL)
+    fftw_free(state->real_rgb);
+  if (state->freq_rgb != NULL)
+    fftw_free(state->freq_rgb);
+  if (state->kernel_freq != NULL)
+    fftw_free(state->kernel_freq);
+  const int registered = state->registered;
+  free(state);
+
+  if (!registered)
+    return;
+  g_live_states--;
+  if (g_live_states <= 0) {
+    g_live_states = 0;
+    if (g_threading_initialized) {
+      fftw_cleanup_threads();
+      g_threading_initialized = 0;
+    }
+  }
+}
+
+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,
+                                       const int *row_span, int fft_workers,
+                                       double *setup_seconds,
+                                       double *kernel_seconds)
+{
+  if (setup_seconds != NULL)
+    *setup_seconds = 0.0;
+  if (kernel_seconds != NULL)
+    *kernel_seconds = 0.0;
+  if (ss_width <= 0 || ss_height <= 0 || final_width <= 0 ||
+      final_height <= 0 || supersample < 1 || kernel_radius < 0 ||
+      weights == NULL || row_span == NULL || fft_workers < 1) {
+    fprintf(stderr,
+            "Fast FFTW init failed (arguments): ss=%dx%d final=%dx%d "
+            "supersample=%d kernel_radius=%d workers=%d\n",
+            ss_width, ss_height, final_width, final_height, supersample,
+            kernel_radius, fft_workers);
+    return NULL;
+  }
+  if ((size_t)final_width * (size_t)supersample != (size_t)ss_width ||
+      (size_t)final_height * (size_t)supersample != (size_t)ss_height) {
+    fprintf(stderr,
+            "Fast FFTW init failed (dimensions): ss=%dx%d is not "
+            "final=%dx%d times supersample=%d\n",
+            ss_width, ss_height, final_width, final_height, supersample);
+    return NULL;
+  }
+
+  const char *stage = "state allocation";
+  FastPsfFftwState *state = calloc(1, sizeof *state);
+  if (state == NULL) {
+    fprintf(stderr, "Fast FFTW init failed (state allocation): %zu bytes\n",
+            sizeof *state);
+    return NULL;
+  }
+  state->ss_width = ss_width;
+  state->ss_height = ss_height;
+  state->final_width = final_width;
+  state->final_height = final_height;
+  state->supersample = supersample;
+  state->kernel_radius = kernel_radius;
+  state->workers = fft_workers;
+  state->plan_measure = g_plan_measure;
+
+  size_t real_bytes = 0, freq_bytes = 0, kernel_bytes = 0;
+  size_t fft_width = 0, fft_height = 0;
+  const size_t radius = (size_t)kernel_radius;
+  size_t width_min, height_min;
+  stage = "kernel extent";
+  if (checked_mul_size(radius, 2, &width_min) ||
+      checked_mul_size(radius, 2, &height_min))
+    goto fail;
+  if ((size_t)ss_width > SIZE_MAX - width_min ||
+      (size_t)ss_height > SIZE_MAX - height_min)
+    goto fail;
+  width_min += (size_t)ss_width;
+  height_min += (size_t)ss_height;
+  stage = "smooth FFT size";
+  if (fast_psf_fftw_next_smooth_size(width_min, &fft_width) ||
+      fast_psf_fftw_next_smooth_size(height_min, &fft_height))
+    goto fail;
+  state->fft_width = (int)fft_width;
+  state->fft_height = (int)fft_height;
+
+  const size_t complex_width = fft_width / 2 + 1;
+  size_t real_plane, real_count;
+  size_t freq_plane, freq_count;
+  stage = "allocation size";
+  if (checked_mul_size(fft_height, fft_width, &real_plane) ||
+      checked_mul_size(real_plane, 3, &real_count) ||
+      checked_mul_size(real_count, sizeof(double), &real_bytes) ||
+      checked_mul_size(fft_height, complex_width, &freq_plane) ||
+      checked_mul_size(freq_plane, 3, &freq_count) ||
+      checked_mul_size(freq_count, sizeof(fftw_complex), &freq_bytes) ||
+      checked_mul_size(freq_plane, sizeof(fftw_complex), &kernel_bytes) ||
+      real_bytes > SIZE_MAX - freq_bytes - kernel_bytes)
+    goto fail;
+  /* FFTW's plan_many interface takes idist/odist as int, so the per-plane
+   * element counts must fit before the cast below even when each axis is a
+   * valid int. */
+  stage = "plan_many stride";
+  if (real_plane > (size_t)INT_MAX || freq_plane > (size_t)INT_MAX)
+    goto fail;
+  state->scratch_bytes = real_bytes + freq_bytes + kernel_bytes;
+
+  stage = "FFTW threading initialization";
+  if (g_live_states == 0) {
+    if (!fftw_init_threads())
+      goto fail;
+    g_threading_initialized = 1;
+  }
+  g_live_states++;
+  state->registered = 1;
+  fftw_plan_with_nthreads(fft_workers);
+
+  const double plan_start = omp_get_wtime();
+  stage = "FFTW scratch allocation";
+  state->real_rgb = fftw_alloc_real(real_count);
+  state->freq_rgb = fftw_alloc_complex(freq_count);
+  state->kernel_freq = fftw_alloc_complex(freq_plane);
+  if (state->real_rgb == NULL || state->freq_rgb == NULL ||
+      state->kernel_freq == NULL)
+    goto fail;
+
+  int n[2] = {state->fft_height, state->fft_width};
+  int real_embed[2] = {state->fft_height, state->fft_width};
+  int complex_embed[2] = {state->fft_height, (int)complex_width};
+  const int real_dist = (int)real_plane;
+  const int complex_dist = (int)freq_plane;
+  stage = "FFTW plan creation";
+  state->forward_plan = fftw_plan_many_dft_r2c(
+      2, n, 3, state->real_rgb, real_embed, 1, real_dist, state->freq_rgb,
+      complex_embed, 1, complex_dist, plan_flags());
+  state->inverse_plan = fftw_plan_many_dft_c2r(
+      2, n, 3, state->freq_rgb, complex_embed, 1, complex_dist,
+      state->real_rgb, real_embed, 1, real_dist, plan_flags());
+  if (state->forward_plan == NULL || state->inverse_plan == NULL)
+    goto fail;
+  state->setup_seconds = omp_get_wtime() - plan_start;
+  if (setup_seconds != NULL)
+    *setup_seconds = state->setup_seconds;
+
+  /* One-time kernel transform.  The kernel is stored unshifted at padded origin
+   * indices (dy + R, dx + R), retaining only the circular mask the spatial
+   * reference traverses; square corners and FFT padding stay zero.  The plan is
+   * created before the kernel is written because FFTW_MEASURE overwrites its
+   * input during planning. */
+  stage = "kernel scratch allocation";
+  double *kernel_real = fftw_alloc_real(real_plane);
+  fftw_plan kernel_plan = NULL;
+  if (kernel_real == NULL)
+    goto fail;
+  const double kernel_start = omp_get_wtime();
+  stage = "kernel plan creation";
+  kernel_plan = fftw_plan_dft_r2c_2d(state->fft_height, state->fft_width,
+                                     kernel_real, state->kernel_freq,
+                                     plan_flags());
+  if (kernel_plan == NULL) {
+    fftw_free(kernel_real);
+    goto fail;
+  }
+  memset(kernel_real, 0, real_plane * sizeof *kernel_real);
+  const size_t side = 2 * radius + 1;
+  for (int dy = -kernel_radius; dy <= kernel_radius; ++dy) {
+    const int span = row_span[dy + kernel_radius];
+    for (int dx = -span; dx <= span; ++dx)
+      kernel_real[(size_t)(dy + kernel_radius) * fft_width +
+                  (size_t)(dx + kernel_radius)] =
+          (double)weights[(size_t)(dy + kernel_radius) * side +
+                          (size_t)(dx + kernel_radius)];
+  }
+  fftw_execute(kernel_plan);
+  state->kernel_seconds = omp_get_wtime() - kernel_start;
+  if (kernel_seconds != NULL)
+    *kernel_seconds = state->kernel_seconds;
+  fftw_destroy_plan(kernel_plan);
+  fftw_free(kernel_real);
+
+  state->fft_scale = 1.0 / ((double)fft_width * (double)fft_height);
+  state->box_scale = 1.0 / ((double)supersample * (double)supersample);
+  return state;
+
+fail:
+  fprintf(stderr,
+          "Fast FFTW init failed at '%s': ss=%dx%d, linear min=%zux%zu, "
+          "fft=%zux%zu, plan=%s, workers=%d, bytes real=%zu freq=%zu "
+          "kernel=%zu total=%zu\n",
+          stage, ss_width, ss_height,
+          (size_t)ss_width + 2 * radius,
+          (size_t)ss_height + 2 * radius, fft_width, fft_height,
+          state->plan_measure ? "measure" : "estimate", fft_workers,
+          real_bytes, freq_bytes, kernel_bytes,
+          real_bytes + freq_bytes + kernel_bytes);
+  fast_psf_fftw_destroy(state);
+  return NULL;
+}
+
+int fast_psf_fftw_resolve(FastPsfFftwState *state,
+                          const double *interleaved_ss_rgb,
+                          double *interleaved_hdr_rgb,
+                          FastPsfFftwFrameTiming *timing)
+{
+  if (state == NULL || interleaved_ss_rgb == NULL ||
+      interleaved_hdr_rgb == NULL)
+    return -1;
+  FastPsfFftwFrameTiming local = {0};
+  const double total_start = omp_get_wtime();
+
+  const size_t plane = (size_t)state->fft_height * state->fft_width;
+  const size_t complex_width = (size_t)state->fft_width / 2 + 1;
+  const size_t complex_plane = (size_t)state->fft_height * complex_width;
+  const size_t real_count = 3 * plane;
+  const size_t freq_count = 3 * complex_plane;
+
+  double start = omp_get_wtime();
+#pragma omp parallel for schedule(static) num_threads(state->workers) \
+    if (state->workers > 1)
+  for (size_t i = 0; i < real_count; ++i)
+    state->real_rgb[i] = 0.0;
+#pragma omp parallel for schedule(static) num_threads(state->workers) \
+    if (state->workers > 1)
+  for (int y = 0; y < state->ss_height; ++y) {
+    for (int x = 0; x < state->ss_width; ++x) {
+      const size_t source =
+          3 * ((size_t)y * state->ss_width + (size_t)x);
+      const size_t target = (size_t)y * state->fft_width + (size_t)x;
+      state->real_rgb[target] = interleaved_ss_rgb[source];
+      state->real_rgb[plane + target] = interleaved_ss_rgb[source + 1];
+      state->real_rgb[2 * plane + target] = interleaved_ss_rgb[source + 2];
+    }
+  }
+  local.zero_pack_seconds = omp_get_wtime() - start;
+
+  start = omp_get_wtime();
+  fftw_execute(state->forward_plan);
+  local.forward_seconds = omp_get_wtime() - start;
+
+  start = omp_get_wtime();
+#pragma omp parallel for schedule(static) num_threads(state->workers) \
+    if (state->workers > 1)
+  for (size_t i = 0; i < freq_count; ++i) {
+    const double ar = state->freq_rgb[i][0];
+    const double ai = state->freq_rgb[i][1];
+    const double br = state->kernel_freq[i % complex_plane][0];
+    const double bi = state->kernel_freq[i % complex_plane][1];
+    state->freq_rgb[i][0] = ar * br - ai * bi;
+    state->freq_rgb[i][1] = ar * bi + ai * br;
+  }
+  local.multiply_seconds = omp_get_wtime() - start;
+
+  start = omp_get_wtime();
+  fftw_execute(state->inverse_plan);
+  local.inverse_seconds = omp_get_wtime() - start;
+
+  start = omp_get_wtime();
+  const int radius = state->kernel_radius;
+  const int supersample = state->supersample;
+  const double scale = state->fft_scale * state->box_scale;
+#pragma omp parallel for schedule(static) num_threads(state->workers) \
+    if (state->workers > 1)
+  for (int row = 0; row < state->final_height; ++row) {
+    for (int column = 0; column < state->final_width; ++column) {
+      double sum[3] = {0.0, 0.0, 0.0};
+      const int base_y = radius + supersample * row;
+      const int base_x = radius + supersample * column;
+      for (int j = 0; j < supersample; ++j)
+        for (int i = 0; i < supersample; ++i) {
+          const size_t offset =
+              (size_t)(base_y + j) * state->fft_width +
+              (size_t)(base_x + i);
+          sum[0] += state->real_rgb[offset];
+          sum[1] += state->real_rgb[plane + offset];
+          sum[2] += state->real_rgb[2 * plane + offset];
+        }
+      const size_t out =
+          3 * ((size_t)row * state->final_width + (size_t)column);
+      interleaved_hdr_rgb[out] += sum[0] * scale;
+      interleaved_hdr_rgb[out + 1] += sum[1] * scale;
+      interleaved_hdr_rgb[out + 2] += sum[2] * scale;
+    }
+  }
+  local.crop_downsample_seconds = omp_get_wtime() - start;
+  local.total_seconds = omp_get_wtime() - total_start;
+  if (timing != NULL)
+    *timing = local;
+  return 0;
+}
+
+void fast_psf_fftw_report(const FastPsfFftwState *state, FILE *stream)
+{
+  if (state == NULL || stream == NULL)
+    return;
+  fprintf(stream,
+          "Fast FFTW: linear min=%dx%d, fft=%dx%d, workers=%d, "
+          "plan=%s, plan=%.6f s, kernel_fft=%.6f s, scratch=%.1f MiB\n",
+          state->ss_height + 2 * state->kernel_radius,
+          state->ss_width + 2 * state->kernel_radius, state->fft_height,
+          state->fft_width, state->workers,
+          state->plan_measure ? "measure" : "estimate", state->setup_seconds,
+          state->kernel_seconds,
+          (double)state->scratch_bytes / (1024.0 * 1024.0));
+}

src/fast_psf_fftw.h

@@ -0,0 +1,67 @@
+#ifndef FAST_PSF_FFTW_H
+#define FAST_PSF_FFTW_H
+
+#include <stddef.h>
+#include <stdio.h>
+
+/* Private FFTW-backed global convolution for fast-mode PSF accumulation.
+ *
+ * This module owns every FFTW allocation and plan.  It replaces only the
+ * spatial convolution inside the explicit fast-mode preview path; the impulse
+ * buffer stays owned by the caller, and the HDR framebuffer is caller-owned and
+ * updated with +=.  It is compiled only for the CPU PSF backend (FAST_PSF_FFTW)
+ * so the HIP and dummy backends do not link FFTW. */
+
+typedef struct FastPsfFftwState FastPsfFftwState;
+
+typedef struct {
+  double zero_pack_seconds;
+  double forward_seconds;
+  double multiply_seconds;
+  double inverse_seconds;
+  double crop_downsample_seconds;
+  double total_seconds;
+} FastPsfFftwFrameTiming;
+
+/* Builds one reusable convolution state from the immutable kernel `weights`
+ * (side 2*kernel_radius+1, row-major) and the cached circular `row_span`
+ * (length 2*kernel_radius+1, the number of retained dx at each dy).  The state
+ * computes a zero-padded linear convolution over Wss x Hss and downsamples by
+ * `supersample`.  `setup_seconds` receives plan creation time and
+ * `kernel_seconds` receives the one-time kernel transform time.  Returns NULL
+ * on any allocation or plan failure; nothing is leaked. */
+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,
+                                       const int *row_span, int fft_workers,
+                                       double *setup_seconds,
+                                       double *kernel_seconds);
+
+/* Convolves `interleaved_ss_rgb` (3 doubles per supersampled pixel, row-major)
+ * into `interleaved_hdr_rgb` using +=.  The input buffer is not modified.
+ * Returns 0 on success and -1 on invalid state. */
+int fast_psf_fftw_resolve(FastPsfFftwState *state,
+                          const double *interleaved_ss_rgb,
+                          double *interleaved_hdr_rgb,
+                          FastPsfFftwFrameTiming *timing);
+
+/* Safe for NULL and partially initialized states; releases plans before the
+ * buffers they reference.  Does not touch process-global FFTW state except on
+ * the final live state. */
+void fast_psf_fftw_destroy(FastPsfFftwState *state);
+
+/* Prints one diagnostic line when `state` is non-NULL. */
+void fast_psf_fftw_report(const FastPsfFftwState *state, FILE *stream);
+
+/* First value >= min_extent whose prime factors are limited to {2,3,5,7}.
+ * Returns 0 and stores the result in *out on success; returns -1 when the
+ * result would exceed INT_MAX or on overflow. */
+int fast_psf_fftw_next_smooth_size(size_t min_extent, size_t *out);
+
+/* Selects FFTW_MEASURE (measure != 0) instead of FFTW_ESTIMATE for plans
+ * created after this call.  Intended for the benchmark's plan-mode comparison;
+ * production uses the default FFTW_ESTIMATE. */
+void fast_psf_fftw_set_plan_mode(int measure);
+
+#endif

src/main.c

@@ -977,6 +977,7 @@ static int render_lens_map(const Settings *s, StarCatalog *catalog) {
   if (s->fast_mode) {
     fast = &local_fast;
     fast_psf_accumulator_report(&local_fast, stderr);
+    fast_psf_accumulator_set_verbose(&local_fast, s->verbose);
   }
   for (size_t i = 0; i < map.frame_count; ++i) {
     const char *output_path = s->output_path;
@@ -1218,6 +1219,7 @@ int main(int argc, char **argv) {
     }
     settings.fast_psf = &fast_accumulator;
     fast_psf_accumulator_report(&fast_accumulator, stderr);
+    fast_psf_accumulator_set_verbose(&fast_accumulator, settings.verbose);
   } else if (settings.fast_mode) {
     fputs("Fast mode on an imported lens map uses the map's own dimensions.\n",
           stderr);

src/optics.c

@@ -1,5 +1,9 @@
 #include "optics.h"
 
+#ifdef FAST_PSF_FFTW
+#include "fast_psf_fftw.h"
+#endif
+
 #include <limits.h>
 #include <math.h>
 #include <omp.h>
@@ -500,12 +504,19 @@ int fast_psf_accumulator_init(FastPsfAccumulator *accumulator, int width,
         side * side > SIZE_MAX / sizeof(float))
         return -1;
     float *weights = malloc(side * side * sizeof *weights);
+    int *row_span = malloc(side * sizeof *row_span);
     double *buffer = calloc(ss_width * ss_height * 3, sizeof *buffer);
-    if (weights == NULL || buffer == NULL) {
+    if (weights == NULL || row_span == NULL || buffer == NULL) {
         free(weights);
+        free(row_span);
         free(buffer);
         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;
+    }
     *accumulator = (FastPsfAccumulator){
         .fwhm_pixels = psf->fwhm_pixels,
         .moffat_beta = psf->moffat_beta,
@@ -522,6 +533,7 @@ int fast_psf_accumulator_init(FastPsfAccumulator *accumulator, int width,
         .radius_pixels = radius,
         .use_reference = use_reference,
         .weights = weights,
+        .row_span = row_span,
         .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.
@@ -539,6 +551,23 @@ int fast_psf_accumulator_init(FastPsfAccumulator *accumulator, int width,
             weights[(size_t)(dy + radius) * side + (dx + radius)] =
                 (float)(weight * normalization_scale);
         }
+#ifdef FAST_PSF_FFTW
+    int fft_workers = omp_get_max_threads();
+    if (fft_workers < 1)
+        fft_workers = 1;
+    accumulator->fftw = fast_psf_fftw_create(
+        (int)ss_width, (int)ss_height, width, height, supersample, radius,
+        weights, row_span, fft_workers, &accumulator->fftw_setup_seconds,
+        &accumulator->fftw_kernel_seconds);
+    if (accumulator->fftw == NULL) {
+        free(weights);
+        free(row_span);
+        free(buffer);
+        *accumulator = (FastPsfAccumulator){0};
+        return -1;
+    }
+    accumulator->fftw_enabled = 1;
+#endif
     return 0;
 }
 
@@ -590,24 +619,20 @@ int fast_psf_accumulator_deposit(FastPsfAccumulator *accumulator, double x,
     return wing_clipped;
 }
 
-int fast_psf_accumulator_resolve(const FastPsfAccumulator *accumulator,
-                                 double *hdr, int worker_count)
+/* Spatial reference global convolution.  Kept as a test/benchmark reference
+ * with the identical arithmetic the FFTW path must reproduce; it is never a
+ * silent runtime fallback. */
+int fast_psf_accumulator_resolve_spatial_reference(
+    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)
+        accumulator->row_span == 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 int *row_span = accumulator->row_span;
     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)
@@ -652,15 +677,54 @@ int fast_psf_accumulator_resolve(const FastPsfAccumulator *accumulator,
             hdr[out + 2] += sum[2] * inverse_block;
         }
     }
-    free(row_span);
     return 0;
 }
 
+int fast_psf_accumulator_resolve(FastPsfAccumulator *accumulator,
+                                 double *hdr, int worker_count)
+{
+    if (accumulator == NULL || accumulator->buffer == NULL || hdr == NULL)
+        return -1;
+#ifdef FAST_PSF_FFTW
+    if (accumulator->fftw_enabled && accumulator->fftw != NULL) {
+        FastPsfFftwFrameTiming timing = {0};
+        const double start = omp_get_wtime();
+        if (fast_psf_fftw_resolve(accumulator->fftw, accumulator->buffer, hdr,
+                                  &timing))
+            return -1;
+        accumulator->fftw_frame_seconds = omp_get_wtime() - start;
+        if (accumulator->verbose)
+            fprintf(stderr,
+                    "Fast FFTW frame: zero_pack=%.6f forward=%.6f "
+                    "multiply=%.6f inverse=%.6f crop_downsample=%.6f "
+                    "total=%.6f\n",
+                    timing.zero_pack_seconds, timing.forward_seconds,
+                    timing.multiply_seconds, timing.inverse_seconds,
+                    timing.crop_downsample_seconds, timing.total_seconds);
+        return 0;
+    }
+#endif
+    return fast_psf_accumulator_resolve_spatial_reference(accumulator, hdr,
+                                                          worker_count);
+}
+
+void fast_psf_accumulator_set_verbose(FastPsfAccumulator *accumulator,
+                                      int verbose)
+{
+    if (accumulator != NULL)
+        accumulator->verbose = verbose ? 1 : 0;
+}
+
 void fast_psf_accumulator_destroy(FastPsfAccumulator *accumulator)
 {
     if (accumulator == NULL)
         return;
+#ifdef FAST_PSF_FFTW
+    fast_psf_fftw_destroy(accumulator->fftw);
+    accumulator->fftw = NULL;
+#endif
     free(accumulator->weights);
+    free(accumulator->row_span);
     free(accumulator->buffer);
     *accumulator = (FastPsfAccumulator){0};
 }
@@ -687,6 +751,9 @@ void fast_psf_accumulator_report(const FastPsfAccumulator *accumulator,
             sum / ((double)accumulator->supersample *
                    accumulator->supersample),
             accumulator->use_reference ? "reference 8-point" : "cached 4-point");
+#ifdef FAST_PSF_FFTW
+    fast_psf_fftw_report(accumulator->fftw, stream);
+#endif
 }
 
 

src/optics.h

@@ -26,6 +26,10 @@ typedef enum {
   FAST_PSF_DEPOSIT_BILINEAR = 1,
 } FastPsfDeposit;
 
+/* Private FFTW convolution state.  Defined in fast_psf_fftw.c; only the CPU
+ * PSF backend may reference the implementation. */
+typedef struct FastPsfFftwState FastPsfFftwState;
+
 /* 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
@@ -45,7 +49,16 @@ typedef struct {
   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 */
+  int *row_span;             /* cached circular half-width per dy, length 2R+1 */
   double *buffer;            /* supersampled HDR, 3 channels per pixel */
+#ifdef FAST_PSF_FFTW
+  FastPsfFftwState *fftw;    /* owned private convolution state */
+  int fftw_enabled;
+  double fftw_setup_seconds;
+  double fftw_kernel_seconds;
+  double fftw_frame_seconds; /* most recent resolve */
+#endif
+  int verbose;
 } FastPsfAccumulator;
 
 typedef struct {
@@ -110,8 +123,14 @@ void fast_psf_accumulator_clear(FastPsfAccumulator *accumulator);
  * 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,
+int fast_psf_accumulator_resolve(FastPsfAccumulator *accumulator,
                                  double *hdr, int worker_count);
+/* Test/benchmark reference only: the original nested-loop global convolution.
+ * Production resolve never falls back to it. */
+int fast_psf_accumulator_resolve_spatial_reference(
+    FastPsfAccumulator *accumulator, double *hdr, int worker_count);
+void fast_psf_accumulator_set_verbose(FastPsfAccumulator *accumulator,
+                                      int verbose);
 void fast_psf_accumulator_destroy(FastPsfAccumulator *accumulator);
 void fast_psf_accumulator_report(const FastPsfAccumulator *accumulator,
                                  FILE *stream);