Feat: Add --fast-mode supersampled point-source accumulation

Add an optional CPU preview path that deposits each point-source image as a
supersampled delta and resolves the whole frame with one global Moffat
convolution plus an N x N box average, instead of splatting a per-event PSF.

- optics: FastPsfAccumulator builds the pixel-area-integral kernel of the
  target Moffat at the supersampled scale (width N*alpha, same beta).  The
  1/N^2 box average then reproduces the final pixel-area integral, so the
  requested FWHM and beta are preserved without renormalisation.  Deposits are
  per-cell atomic adds; resolve accumulates into the caller's HDR buffer.
- frame: fast branch in frame_splat_catalog with one shared supersampled
  buffer and a single resolve per frame; the accumulator is reused across
  movie frames and built from the map dimensions on lens-map import.
- main: --fast-mode, --fast-supersample N (1..8, default 2) and
  --fast-deposit nearest|bilinear (default nearest).  CPU-only and rejected in
  the HIP/dummy backends; --psf-min-y still applies per event while
  --max-cache-psf-flux does not.
- The deposition scheme was chosen by scripts/fast_mode_deposit_error.py:
  nearest keeps the PSF shape exactly with <= 0.5/N px position quantization;
  bilinear keeps the exact centroid but broadens FWHM and beta.  Recorded in
  benchmarks/fast_mode_deposit_2026-09-18.md.
- tests/test_frame.c covers fast nearest vs the direct evaluator at the snapped
  centre, flux conservation, bilinear centroid, min-Y discard, frame plumbing,
  and HDR accumulation onto a non-zero background.
- benchmarks/fast_mode_cpu_2026-09-18.md records a ~10x speedup on the 2MASS
  galactic-centre field with small tone-mapped differences.

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;

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);

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;
 }

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. */

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,