Optics: cache pixel-area Moffat kernels

src/frame.c

@@ -184,7 +184,9 @@ typedef struct {
   int width, height;
   double exposure, magnification;
   const PointSpreadFunction *psf;
+  const PsfKernelCache *psf_cache;
   size_t images;
+  size_t direct_fallbacks;
 } TriangleSplatContext;
 
 static int splat_catalog_tile(const Star *stars, size_t count,
@@ -218,9 +220,10 @@ static int splat_catalog_tile(const Star *stars, size_t count,
                          weights[1] * context->vertex[1]->log_frequency_ratio +
                          weights[2] * context->vertex[2]->log_frequency_ratio;
     const LinearRgb color = blackbody_to_linear_rgb(star->temperature_K * exp(log_g));
-    splat_moffat(context->hdr, context->width, context->height, image_x, image_y,
-                 color, context->exposure * star->amplitude * context->magnification,
-                 context->psf);
+    context->direct_fallbacks += splat_moffat_cached(
+        context->hdr, context->width, context->height, image_x, image_y, color,
+        context->exposure * star->amplitude * context->magnification,
+        context->psf, context->psf_cache);
     ++context->images;
   }
   return 0;
@@ -230,8 +233,10 @@ static size_t splat_catalog_triangles(const FrameLensMesh *mesh,
                                       StarCatalog *catalog, double *hdr,
                                       int width, int height, double exposure,
                                       const PointSpreadFunction *psf,
+                                      const PsfKernelCache *psf_cache,
                                       size_t first_triangle,
-                                      size_t last_triangle) {
+                                      size_t last_triangle,
+                                      size_t *direct_fallbacks) {
   size_t images = 0;
   for (size_t t = first_triangle; t < last_triangle; ++t) {
     const LensVertex *vertex[3];
@@ -251,10 +256,11 @@ static size_t splat_catalog_triangles(const FrameLensMesh *mesh,
                                     .triangle = &mesh->triangles[t], .hdr = hdr,
                                     .width = width, .height = height,
                                     .exposure = exposure, .magnification = magnification,
-                                    .psf = psf};
+                                    .psf = psf, .psf_cache = psf_cache};
     if (catalog_visit_source_triangle(catalog, direction, 0, splat_catalog_tile,
                                       &context) == 0)
       images += context.images;
+    *direct_fallbacks += context.direct_fallbacks;
   }
   return images;
 }
@@ -283,8 +289,10 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
                            StarCatalog *catalog, double *hdr, int width,
                            int height, double exposure,
                            const PointSpreadFunction *psf,
+                           const PsfKernelCache *psf_cache,
                            int catalog_load_workers,
-                           CatalogPrefetchStats *prefetch_stats) {
+                           CatalogPrefetchStats *prefetch_stats,
+                           PsfSplatStats *psf_stats) {
   if (mesh == NULL || catalog == NULL || hdr == NULL || exposure <= 0.0 ||
       psf == NULL || width <= 0 || height <= 0 || catalog_load_workers <= 0)
     return 0;
@@ -292,25 +300,45 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
   /* A bounded parallel read phase completes before splatting. Its serial cache
    * commit leaves immutable tile data for the OpenMP splat workers. */
   prefetch_catalog_for_mesh(mesh, catalog, catalog_load_workers, prefetch_stats);
+  if (psf_stats != NULL)
+    *psf_stats = (PsfSplatStats){0};
 
   const size_t pixel_count = (size_t)width * height * 3;
   if (pixel_count > SIZE_MAX / sizeof(double) ||
       pixel_count * sizeof(double) > FRAME_SPLAT_MAX_PRIVATE_HDR_BYTES / 2)
-    return splat_catalog_triangles(mesh, catalog, hdr, width, height, exposure,
-                                   psf, 0, mesh->triangle_count);
+    {
+      size_t fallbacks = 0;
+      const size_t images = splat_catalog_triangles(mesh, catalog, hdr, width, height,
+                                                     exposure, psf, psf_cache, 0,
+                                                     mesh->triangle_count, &fallbacks);
+      if (psf_stats != NULL) *psf_stats = (PsfSplatStats){images - fallbacks, fallbacks};
+      return images;
+    }
   const size_t buffer_bytes = pixel_count * sizeof(double);
   size_t worker_count = FRAME_SPLAT_MAX_PRIVATE_HDR_BYTES / buffer_bytes;
   const int max_threads = omp_get_max_threads();
   if (worker_count > (size_t)max_threads)
     worker_count = (size_t)max_threads;
   if (worker_count < 2 || worker_count > INT_MAX)
-    return splat_catalog_triangles(mesh, catalog, hdr, width, height, exposure,
-                                   psf, 0, mesh->triangle_count);
+    {
+      size_t fallbacks = 0;
+      const size_t images = splat_catalog_triangles(mesh, catalog, hdr, width, height,
+                                                     exposure, psf, psf_cache, 0,
+                                                     mesh->triangle_count, &fallbacks);
+      if (psf_stats != NULL) *psf_stats = (PsfSplatStats){images - fallbacks, fallbacks};
+      return images;
+    }
 
   double **private_hdr = calloc(worker_count, sizeof *private_hdr);
   if (private_hdr == NULL)
-    return splat_catalog_triangles(mesh, catalog, hdr, width, height, exposure,
-                                   psf, 0, mesh->triangle_count);
+    {
+      size_t fallbacks = 0;
+      const size_t images = splat_catalog_triangles(mesh, catalog, hdr, width, height,
+                                                     exposure, psf, psf_cache, 0,
+                                                     mesh->triangle_count, &fallbacks);
+      if (psf_stats != NULL) *psf_stats = (PsfSplatStats){images - fallbacks, fallbacks};
+      return images;
+    }
   size_t allocated = 0;
   for (; allocated < worker_count; ++allocated) {
     private_hdr[allocated] = calloc(pixel_count, sizeof **private_hdr);
@@ -321,12 +349,17 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
     while (allocated > 0)
       free(private_hdr[--allocated]);
     free(private_hdr);
-    return splat_catalog_triangles(mesh, catalog, hdr, width, height, exposure,
-                                   psf, 0, mesh->triangle_count);
+    size_t fallbacks = 0;
+    const size_t images = splat_catalog_triangles(mesh, catalog, hdr, width, height,
+                                                   exposure, psf, psf_cache, 0,
+                                                   mesh->triangle_count, &fallbacks);
+    if (psf_stats != NULL)
+      *psf_stats = (PsfSplatStats){images - fallbacks, fallbacks};
+    return images;
   }
 
-  size_t images = 0;
-#pragma omp parallel num_threads((int)worker_count) reduction(+ : images)
+  size_t images = 0, direct_fallbacks = 0;
+#pragma omp parallel num_threads((int)worker_count) reduction(+ : images, direct_fallbacks)
   {
     const size_t worker = (size_t)omp_get_thread_num();
     /* Source density and lens magnification can vary by orders of magnitude
@@ -335,9 +368,9 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
      * idle.  Each worker still owns its HDR buffer exclusively. */
 #pragma omp for schedule(dynamic, 1)
     for (size_t triangle = 0; triangle < mesh->triangle_count; ++triangle)
-      images += splat_catalog_triangles(mesh, catalog, private_hdr[worker],
-                                        width, height, exposure, psf, triangle,
-                                        triangle + 1);
+      images += splat_catalog_triangles(mesh, catalog, private_hdr[worker], width,
+                                        height, exposure, psf, psf_cache, triangle,
+                                        triangle + 1, &direct_fallbacks);
   }
 #pragma omp parallel for schedule(static)
   for (size_t pixel = 0; pixel < pixel_count; ++pixel)
@@ -346,6 +379,8 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
   for (size_t worker = 0; worker < worker_count; ++worker)
     free(private_hdr[worker]);
   free(private_hdr);
+  if (psf_stats != NULL)
+    *psf_stats = (PsfSplatStats){images - direct_fallbacks, direct_fallbacks};
   return images;
 }
 

src/frame.h

@@ -39,8 +39,10 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
                            StarCatalog *catalog, double *hdr, int width,
                            int height, double exposure,
                            const PointSpreadFunction *psf,
+                           const PsfKernelCache *psf_cache,
                            int catalog_load_workers,
-                           CatalogPrefetchStats *prefetch_stats);
+                           CatalogPrefetchStats *prefetch_stats,
+                           PsfSplatStats *psf_stats);
 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

@@ -18,13 +18,18 @@ typedef struct {
   int width, height;
   int coarse_cell_pixels;
   int draw_mesh;
+  int psf_direct;
   double horizontal_fov_deg, look_ra_deg, look_dec_deg, exposure;
   double observer_radius;
   double observer_inward_speed;
   PointSpreadFunction psf;
+  PsfKernelCache psf_cache;
   const char *catalog_path;
   const char *all_sky_catalog_path;
   const char *output_path;
+#ifdef ENABLE_HDR_DEBUG
+  const char *hdr_output_path;
+#endif
   const char *observer_track_path;
   const char *frames_dir;
   const char *frames_prefix;
@@ -126,6 +131,10 @@ static int parse_args(int argc, char **argv, Settings *s,
       s->all_sky_catalog_path = argv[++i];
     else if (!strcmp(argv[i], "--output") && i + 1 < argc)
       s->output_path = argv[++i];
+#ifdef ENABLE_HDR_DEBUG
+    else if (!strcmp(argv[i], "--hdr-output") && i + 1 < argc)
+      s->hdr_output_path = argv[++i];
+#endif
     else if (!strcmp(argv[i], "--width") && i + 1 < argc &&
              !parse_int(argv[++i], &s->width)) {
     } else if (!strcmp(argv[i], "--height") && i + 1 < argc &&
@@ -150,6 +159,8 @@ static int parse_args(int argc, char **argv, Settings *s,
                !parse_positive(argv[++i], &s->psf.fwhm_pixels)) {
     } else if (!strcmp(argv[i], "--psf-moffat-beta") && i + 1 < argc &&
                !parse_moffat_beta(argv[++i], &s->psf.moffat_beta)) {
+    } else if (!strcmp(argv[i], "--psf-direct")) {
+      s->psf_direct = 1;
     } else if (!strcmp(argv[i], "--write-catalog") && i + 1 < argc)
       *write_path = argv[++i];
     else if (!strcmp(argv[i], "--observer-track") && i + 1 < argc)
@@ -219,15 +230,26 @@ static int render_observer_frame(const Settings *s, StarCatalog *catalog,
     return -1;
   }
   CatalogPrefetchStats prefetch = {0};
+  PsfSplatStats psf_stats = {0};
   size_t images = frame_splat_catalog(
       &mesh, catalog, hdr, s->width, s->height, s->exposure, &s->psf,
-      s->catalog_load_workers, &prefetch);
+      &s->psf_cache, s->catalog_load_workers, &prefetch, &psf_stats);
   if (s->draw_mesh)
     frame_draw_mesh(&mesh, hdr, s->width, s->height, 0.5, 0.5);
+#ifdef ENABLE_HDR_DEBUG
+  if (s->hdr_output_path != NULL &&
+      write_hdr_pfm(s->hdr_output_path, hdr, s->width, s->height)) {
+    perror(s->hdr_output_path);
+    frame_lens_mesh_destroy(&mesh);
+    free(hdr);
+    return -1;
+  }
+#endif
   int result = write_tonemapped_image(output_path, hdr, s->width, s->height);
   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);
   if (catalog->kind == STAR_CATALOG_ALL_SKY)
     fprintf(stderr,
             "Catalog prefetch: %zu requested, %zu newly loaded (%zu stars), "
@@ -319,15 +341,17 @@ static int render_movie(const Settings *s, StarCatalog *catalog,
       goto done;
     }
     CatalogPrefetchStats prefetch = {0};
+    PsfSplatStats psf_stats = {0};
     const size_t images = frame_splat_catalog(
         &movie.frames[i].mesh, catalog, hdr, s->width, s->height, s->exposure,
-        &s->psf, s->catalog_load_workers, &prefetch);
+        &s->psf, &s->psf_cache, s->catalog_load_workers, &prefetch, &psf_stats);
     if (s->draw_mesh)
       frame_draw_mesh(&movie.frames[i].mesh, hdr, s->width, s->height, 0.5, 0.5);
     const int write_result = write_tonemapped_image(output_path, hdr, s->width, s->height);
     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);
     if (catalog->kind == STAR_CATALOG_ALL_SKY)
       fprintf(stderr,
               "Catalog prefetch: %zu requested, %zu newly loaded (%zu stars), "
@@ -369,6 +393,10 @@ int main(int argc, char **argv) {
             "N] [--fov-deg D] [--look-ra-deg D] [--look-dec-deg D] "
             "[--exposure E] [--observer-radius R] [--observer-inward-speed V] "
             "[--psf-fwhm-pixels N] [--psf-moffat-beta N] "
+            "[--psf-direct] "
+#ifdef ENABLE_HDR_DEBUG
+            "[--hdr-output PATH] "
+#endif
             "[--coarse-cell-pixels N] [--draw-mesh] [--write-catalog PATH] "
             "[--catalog-load-workers N] "
             "[--observer-track PATH --frames-dir DIR --frames-prefix NAME "
@@ -384,6 +412,12 @@ int main(int argc, char **argv) {
     return write_minkowski_accel_track(&settings) == 0
                ? 0
                : (perror(settings.write_minkowski_accel_track_path), 1);
+#ifdef ENABLE_HDR_DEBUG
+  if (settings.frames_dir != NULL && settings.hdr_output_path != NULL) {
+    fputs("--hdr-output is available only for a single-frame render.\n", stderr);
+    return 2;
+  }
+#endif
   StarCatalog catalog = {0};
   if (settings.all_sky_catalog_path != NULL) {
     if (catalog_load_all_sky(&catalog, settings.all_sky_catalog_path)) {
@@ -404,10 +438,13 @@ int main(int argc, char **argv) {
     catalog_destroy(&catalog);
     return 1;
   }
+  if (!settings.psf_direct && psf_kernel_cache_init(&settings.psf_cache, &settings.psf))
+    fputs("PSF cache construction failed; using direct evaluator.\n", stderr);
   int result = settings.frames_dir != NULL
                    ? render_movie(&settings, &catalog, &spacetime)
                    : render_frame(&settings, &catalog, &spacetime);
   spacetime_destroy(&spacetime);
   catalog_destroy(&catalog);
+  psf_kernel_cache_destroy(&settings.psf_cache);
   return result == 0 ? 0 : 1;
 }

src/optics.c

@@ -1,9 +1,11 @@
 #include "optics.h"
 
 #include <math.h>
+#include <stdint.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
+#include <time.h>
 
 #ifdef ENABLE_PNG
 #include <png.h>
@@ -76,39 +78,272 @@ LinearRgb blackbody_to_linear_rgb(double temperature_K)
               0.0, INFINITY)};
 }
 
+enum {
+    PSF_PHASE_RESOLUTION = 64,
+    PSF_MAX_CACHE_RADIUS_PIXELS = 48,
+    PSF_QUADRATURE_ORDER = 4,
+};
+
+/* These are display-space HDR error budgets, deliberately independent of a
+ * relative flux cut.  The historical 1e-8 relative tail remains a floor for
+ * ordinary images; brighter images grow their support or use the reference
+ * fallback rather than acquiring a visible clipped wing. */
+static const double psf_relative_tail_fraction = 1e-8;
+static const double psf_tail_absolute_hdr_budget = 1e-6;
+static const double psf_boundary_hdr_budget = 1e-7;
+static const double pi = 3.14159265358979323846;
+static const double gauss4_x[PSF_QUADRATURE_ORDER] = {
+    -0.8611363115940526, -0.3399810435848563,
+     0.3399810435848563,  0.8611363115940526};
+static const double gauss4_w[PSF_QUADRATURE_ORDER] = {
+    0.3478548451374539, 0.6521451548625461,
+    0.6521451548625461, 0.3478548451374539};
+static const double gauss8_x[8] = {
+    -0.9602898564975363, -0.7966664774136267, -0.5255324099163290,
+    -0.1834346424956498,  0.1834346424956498,  0.5255324099163290,
+     0.7966664774136267,  0.9602898564975363};
+static const double gauss8_w[8] = {
+    0.1012285362903763, 0.2223810344533745, 0.3137066458778873,
+    0.3626837833783620, 0.3626837833783620, 0.3137066458778873,
+    0.2223810344533745, 0.1012285362903763};
+
+static int valid_psf(const PointSpreadFunction *psf)
+{
+    return psf != NULL && isfinite(psf->fwhm_pixels) &&
+           isfinite(psf->moffat_beta) && psf->fwhm_pixels > 0.0 &&
+           psf->moffat_beta > 1.0;
+}
+
+static double moffat_alpha(const PointSpreadFunction *psf)
+{
+    return psf->fwhm_pixels /
+           (2.0 * sqrt(pow(2.0, 1.0 / psf->moffat_beta) - 1.0));
+}
+
+static double moffat_support_radius(double alpha, double beta,
+                                    double scaled_flux)
+{
+    const double normalization = scaled_flux * (beta - 1.0) / (pi * alpha * alpha);
+    const double relative_tail = fmin(psf_relative_tail_fraction,
+                                      psf_tail_absolute_hdr_budget / scaled_flux);
+    const double tail_radius = alpha * sqrt(pow(relative_tail,
+                                                 1.0 / (1.0 - beta)) - 1.0);
+    const double boundary_ratio = psf_boundary_hdr_budget / normalization;
+    const double boundary_radius = boundary_ratio >= 1.0 ? 0.0 : alpha * sqrt(
+        pow(boundary_ratio, -1.0 / beta) - 1.0);
+    return fmax(tail_radius, boundary_radius);
+}
+
+static double moffat_pixel_integral_quadrature(double alpha, double beta,
+                                               double pixel_x, double pixel_y,
+                                               double star_x, double star_y,
+                                               const double *nodes,
+                                               const double *weights, int order)
+{
+    const double normalization = (beta - 1.0) / (pi * alpha * alpha);
+    double sum = 0.0;
+    for (int iy = 0; iy < order; ++iy)
+        for (int ix = 0; ix < order; ++ix) {
+            const double sx = pixel_x + 0.5 + 0.5 * nodes[ix] - star_x;
+            const double sy = pixel_y + 0.5 + 0.5 * nodes[iy] - star_y;
+            sum += 0.25 * weights[ix] * weights[iy] * normalization *
+                   pow(1.0 + (sx * sx + sy * sy) / (alpha * alpha), -beta);
+        }
+    return sum;
+}
+
+static double moffat_pixel_integral_cached(double alpha, double beta,
+                                           double pixel_x, double pixel_y,
+                                           double star_x, double star_y)
+{
+    return moffat_pixel_integral_quadrature(alpha, beta, pixel_x, pixel_y,
+                                            star_x, star_y, gauss4_x, gauss4_w, 4);
+}
+
+static double moffat_pixel_integral_reference(double alpha, double beta,
+                                              double pixel_x, double pixel_y,
+                                              double star_x, double star_y)
+{
+    return moffat_pixel_integral_quadrature(alpha, beta, pixel_x, pixel_y,
+                                            star_x, star_y, gauss8_x, gauss8_w, 8);
+}
+
+static size_t kernel_index(const PsfKernelCache *cache, int phase_x,
+                           int phase_y, int offset_x, int offset_y)
+{
+    const size_t nodes = (size_t)cache->phase_resolution + 1;
+    const size_t side = (size_t)cache->radius_pixels * 2 + 1;
+    return (((size_t)phase_y * nodes + phase_x) * side +
+            (size_t)(offset_y + cache->radius_pixels)) * side +
+           (size_t)(offset_x + cache->radius_pixels);
+}
+
+void psf_kernel_cache_destroy(PsfKernelCache *cache)
+{
+    if (cache == NULL)
+        return;
+    free(cache->weights);
+    *cache = (PsfKernelCache){0};
+}
+
+int psf_kernel_cache_init(PsfKernelCache *cache, const PointSpreadFunction *psf)
+{
+    if (cache == NULL || !valid_psf(psf))
+        return -1;
+    psf_kernel_cache_destroy(cache);
+    const double alpha = moffat_alpha(psf);
+    const double requested_radius = moffat_support_radius(alpha, psf->moffat_beta, 1.0);
+    if (!isfinite(requested_radius) || requested_radius <= 0.0)
+        return -1;
+    const int radius = (int)fmin((double)PSF_MAX_CACHE_RADIUS_PIXELS,
+                                 ceil(requested_radius) + 1.0);
+    const size_t nodes = PSF_PHASE_RESOLUTION + 1u;
+    const size_t side = (size_t)radius * 2 + 1u;
+    if (nodes > SIZE_MAX / nodes || nodes * nodes > SIZE_MAX / side ||
+        nodes * nodes * side > SIZE_MAX / side ||
+        nodes * nodes * side * side > SIZE_MAX / sizeof(float))
+        return -1;
+    const size_t count = nodes * nodes * side * side;
+    float *weights = malloc(count * sizeof *weights);
+    if (weights == NULL)
+        return -1;
+    const clock_t start = clock();
+    PsfKernelCache building = {.fwhm_pixels = psf->fwhm_pixels,
+                               .moffat_beta = psf->moffat_beta,
+                               .alpha_pixels = alpha,
+                               .max_radius_pixels = radius,
+                               .weights = weights,
+                               .phase_resolution = PSF_PHASE_RESOLUTION,
+                               .radius_pixels = radius};
+    for (int phase_y = 0; phase_y <= PSF_PHASE_RESOLUTION; ++phase_y)
+        for (int phase_x = 0; phase_x <= PSF_PHASE_RESOLUTION; ++phase_x) {
+            const double star_x = (double)phase_x / PSF_PHASE_RESOLUTION;
+            const double star_y = (double)phase_y / PSF_PHASE_RESOLUTION;
+            double sum = 0.0;
+            for (int offset_y = -radius; offset_y <= radius; ++offset_y)
+                for (int offset_x = -radius; offset_x <= radius; ++offset_x) {
+                    const size_t index = kernel_index(&building, phase_x, phase_y,
+                                                      offset_x, offset_y);
+                    const double weight = moffat_pixel_integral_cached(
+                        alpha, psf->moffat_beta, offset_x, offset_y, star_x, star_y);
+                    weights[index] = (float)weight;
+                    sum += weight;
+                }
+            if (!(sum > 0.0) || !isfinite(sum)) {
+                free(weights);
+                return -1;
+            }
+            for (int offset_y = -radius; offset_y <= radius; ++offset_y)
+                for (int offset_x = -radius; offset_x <= radius; ++offset_x) {
+                    const size_t index = kernel_index(&building, phase_x, phase_y,
+                                                      offset_x, offset_y);
+                    weights[index] = (float)(weights[index] / sum);
+                }
+        }
+    building.build_seconds = (double)(clock() - start) / CLOCKS_PER_SEC;
+    building.ready = 1;
+    *cache = building;
+    return 0;
+}
+
+void psf_kernel_cache_report(const PsfKernelCache *cache,
+                             const PsfSplatStats *stats, FILE *stream)
+{
+    if (stream == NULL)
+        return;
+    if (cache == NULL || !cache->ready) {
+        fprintf(stream, "PSF cache: disabled; cached splats 0, direct fallbacks %zu\n",
+                stats == NULL ? 0u : stats->direct_fallbacks);
+        return;
+    }
+    fprintf(stream, "PSF cache: 64x64 phases, radius %.0f px, tail abs %.0e, "
+            "boundary %.0e, build %.3f s; cached splats %zu, direct fallbacks %zu\n",
+            cache->max_radius_pixels, psf_tail_absolute_hdr_budget,
+            psf_boundary_hdr_budget, cache->build_seconds,
+            stats == NULL ? 0u : stats->cached_splats,
+            stats == NULL ? 0u : stats->direct_fallbacks);
+}
+
+void splat_moffat_direct(double *hdr, int width, int height, double x, double y,
+                         LinearRgb color, double flux,
+                         const PointSpreadFunction *psf)
+{
+    if (hdr == NULL || width <= 0 || height <= 0 || flux <= 0.0 || !valid_psf(psf))
+        return;
+    const double alpha = moffat_alpha(psf);
+    if (!isfinite(flux))
+        return;
+    const double support_radius = moffat_support_radius(alpha, psf->moffat_beta, flux);
+    const int min_x = fmax(0.0, floor(x - support_radius));
+    const int max_x = fmin((double)width - 1.0, ceil(x + support_radius));
+    const int min_y = fmax(0.0, floor(y - support_radius));
+    const int max_y = fmin((double)height - 1.0, ceil(y + support_radius));
+    for (int py = min_y; py <= max_y; ++py) for (int px = min_x; px <= max_x; ++px) {
+        if (px < 0 || px >= width || py < 0 || py >= height)
+            continue;
+        const double dx = (px + 0.5) - x, dy = (py + 0.5) - y;
+        if (dx * dx + dy * dy > support_radius * support_radius)
+            continue;
+        const double weight = moffat_pixel_integral_reference(alpha, psf->moffat_beta,
+                                                               px, py, x, y);
+        double *pixel = &hdr[3 * (py * width + px)];
+        pixel[0] += color.r * flux * weight;
+        pixel[1] += color.g * flux * weight;
+        pixel[2] += color.b * flux * weight;
+    }
+}
+
+int splat_moffat_cached(double *hdr, int width, int height, double x, double y,
+                        LinearRgb color, double flux,
+                        const PointSpreadFunction *psf,
+                        const PsfKernelCache *cache)
+{
+    if (hdr == NULL || width <= 0 || height <= 0 || flux <= 0.0 || !valid_psf(psf))
+        return 1;
+    const double alpha = moffat_alpha(psf);
+    const double support_radius = moffat_support_radius(alpha, psf->moffat_beta,
+                                                        flux);
+    if (cache == NULL || !cache->ready || cache->fwhm_pixels != psf->fwhm_pixels ||
+        cache->moffat_beta != psf->moffat_beta || !isfinite(support_radius) ||
+        support_radius > cache->max_radius_pixels) {
+        splat_moffat_direct(hdr, width, height, x, y, color, flux, psf);
+        return 1;
+    }
+    const double base_x = floor(x), base_y = floor(y);
+    const double fx = x - base_x, fy = y - base_y;
+    const double phase_x = fx * cache->phase_resolution;
+    const double phase_y = fy * cache->phase_resolution;
+    const int x0 = (int)floor(phase_x), y0 = (int)floor(phase_y);
+    const int x1 = x0 + 1, y1 = y0 + 1;
+    const double tx = phase_x - x0, ty = phase_y - y0;
+    const int support = (int)ceil(support_radius);
+    for (int offset_y = -support; offset_y <= support; ++offset_y)
+        for (int offset_x = -support; offset_x <= support; ++offset_x) {
+            const int px = (int)base_x + offset_x, py = (int)base_y + offset_y;
+            if (px < 0 || px >= width || py < 0 || py >= height)
+                continue;
+            const double dx = offset_x + 0.5 - fx, dy = offset_y + 0.5 - fy;
+            if (dx * dx + dy * dy > support_radius * support_radius)
+                continue;
+            const double w00 = cache->weights[kernel_index(cache, x0, y0, offset_x, offset_y)];
+            const double w10 = cache->weights[kernel_index(cache, x1, y0, offset_x, offset_y)];
+            const double w01 = cache->weights[kernel_index(cache, x0, y1, offset_x, offset_y)];
+            const double w11 = cache->weights[kernel_index(cache, x1, y1, offset_x, offset_y)];
+            const double weight = (1.0 - ty) * ((1.0 - tx) * w00 + tx * w10) +
+                                  ty * ((1.0 - tx) * w01 + tx * w11);
+            double *pixel = &hdr[3 * (py * width + px)];
+            pixel[0] += color.r * flux * weight;
+            pixel[1] += color.g * flux * weight;
+            pixel[2] += color.b * flux * weight;
+        }
+    return 0;
+}
+
 void splat_moffat(double *hdr, int width, int height, double x, double y,
                   LinearRgb color, double flux,
                   const PointSpreadFunction *psf)
 {
-    /* The tail omitted outside this radius contains 1e-8 of the Moffat's
-     * total flux. Unlike the old fixed 3-sigma box, this is both circular and
-     * far below the displayed HDR precision for the chosen beta. */
-    const double tail_fraction = 1e-8;
-    if (hdr == NULL || psf == NULL || flux <= 0.0 ||
-        psf->fwhm_pixels <= 0.0 || psf->moffat_beta <= 1.0)
-        return;
-    const double beta = psf->moffat_beta;
-    const double alpha = psf->fwhm_pixels /
-                         (2.0 * sqrt(pow(2.0, 1.0 / beta) - 1.0));
-    const double support_radius = alpha * sqrt(
-        pow(tail_fraction, 1.0 / (1.0 - beta)) - 1.0);
-    const double support_radius_squared = support_radius * support_radius;
-    const int min_x = (int)floor(x - support_radius);
-    const int max_x = (int)ceil(x + support_radius);
-    const int min_y = (int)floor(y - support_radius);
-    const int max_y = (int)ceil(y + support_radius);
-    const double normalization = flux * (beta - 1.0) /
-                                 (3.14159265358979323846 * alpha * alpha);
-    for (int py = min_y; py <= max_y; ++py) for (int px = min_x; px <= max_x; ++px) {
-        if (px < 0 || px >= width || py < 0 || py >= height) continue;
-        double dx = (px + 0.5) - x, dy = (py + 0.5) - y;
-        const double radius_squared = dx * dx + dy * dy;
-        if (radius_squared > support_radius_squared) continue;
-        const double w = normalization *
-                         pow(1.0 + radius_squared / (alpha * alpha), -beta);
-        double *pixel = &hdr[3 * (py * width + px)];
-        pixel[0] += color.r * w; pixel[1] += color.g * w; pixel[2] += color.b * w;
-    }
+    splat_moffat_direct(hdr, width, height, x, y, color, flux, psf);
 }
 
 static unsigned char tonemap_channel(double hdr_value)
@@ -186,3 +421,30 @@ int write_tonemapped_image(const char *path, const double *hdr, int width, int h
     return -1;
 #endif
 }
+
+#ifdef ENABLE_HDR_DEBUG
+int write_hdr_pfm(const char *path, const double *hdr, int width, int height)
+{
+    if (path == NULL || hdr == NULL || width <= 0 || height <= 0)
+        return -1;
+    FILE *file = fopen(path, "wb");
+    if (file == NULL)
+        return -1;
+    const uint16_t endian_probe = 1;
+    const char *scale = *(const unsigned char *)&endian_probe == 1 ? "-1.0" : "1.0";
+    int result = fprintf(file, "PF\n%d %d\n%s\n", width, height, scale) < 0 ? -1 : 0;
+    /* PFM rows are stored bottom-to-top.  Its negative scale declares little-endian
+     * float samples, avoiding an unnecessary byte swap on the normal test host. */
+    for (int row = height - 1; result == 0 && row >= 0; --row)
+        for (int column = 0; column < width * 3; ++column) {
+            const float sample = (float)hdr[(size_t)row * width * 3 + column];
+            if (fwrite(&sample, sizeof sample, 1, file) != 1) {
+                result = -1;
+                break;
+            }
+        }
+    if (fclose(file) != 0)
+        result = -1;
+    return result;
+}
+#endif

src/optics.h

@@ -1,20 +1,57 @@
 #ifndef OPTICS_H
 #define OPTICS_H
 
+#include <stddef.h>
+#include <stdio.h>
+
 typedef struct { double r, g, b; } LinearRgb;
 typedef struct {
   double fwhm_pixels;
   double moffat_beta;
 } PointSpreadFunction;
 
+/* One immutable process-wide kernel for the one PSF parameter pair accepted
+ * by the current renderer.  Its storage remains private to optics.c. */
+typedef struct {
+  double fwhm_pixels, moffat_beta, alpha_pixels;
+  double max_radius_pixels, build_seconds;
+  float *weights;
+  int phase_resolution, radius_pixels;
+  int ready;
+} PsfKernelCache;
+
+typedef struct {
+  size_t cached_splats;
+  size_t direct_fallbacks;
+} PsfSplatStats;
+
 /* Integrate a Planck spectrum into absolute linear-sRGB spectral radiance
  * (W m^-2 sr^-1), before catalog amplitude and display exposure. */
 LinearRgb blackbody_to_linear_rgb(double temperature_K);
+int psf_kernel_cache_init(PsfKernelCache *cache,
+                          const PointSpreadFunction *psf);
+void psf_kernel_cache_destroy(PsfKernelCache *cache);
+void psf_kernel_cache_report(const PsfKernelCache *cache,
+                             const PsfSplatStats *stats, 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,
+                         LinearRgb color, double flux,
+                         const PointSpreadFunction *psf);
 void splat_moffat(double *hdr, int width, int height, double x, double y,
                   LinearRgb color, double flux,
                   const PointSpreadFunction *psf);
+/* Returns nonzero when this event used the direct reference fallback. */
+int splat_moffat_cached(double *hdr, int width, int height, double x, double y,
+                        LinearRgb color, double flux,
+                        const PointSpreadFunction *psf,
+                        const PsfKernelCache *cache);
 /* Writes PNG when built with libpng; non-libpng builds use PPM fallback. */
 int write_tonemapped_image(const char *path, const double *hdr, int width,
                            int height);
+#ifdef ENABLE_HDR_DEBUG
+/* Test-build-only: writes the pre-tone-mapping framebuffer as RGB float PFM. */
+int write_hdr_pfm(const char *path, const double *hdr, int width, int height);
+#endif
 
 #endif