Frame: buffer cached PSF events per worker

src/frame.c

@@ -10,6 +10,10 @@
 #include <stdlib.h>
 #include <string.h>
 
+#ifndef FRAME_PSF_EVENT_SINK
+#define FRAME_PSF_EVENT_SINK 1
+#endif
+
 /* Numerical metric backends may reserve substantial memory for slabs and
  * thread-local evaluators, so they retain this private-HDR allocation budget.
  * Analytic backends deliberately use all OpenMP render workers instead. */
@@ -115,6 +119,48 @@ int frame_lens_mesh_trace(FrameLensMesh *mesh, const SpacetimeSource *spacetime,
   return 0;
 }
 
+enum { PSF_EVENT_SINK_CAPACITY = 16384 };
+
+typedef struct {
+  PsfCachedEvent *events;
+  size_t count;
+  double *hdr;
+  int width, height;
+  const PsfKernelCache *cache;
+} PsfEventSink;
+
+static void psf_event_sink_init(PsfEventSink *sink, double *hdr, int width,
+                                int height, const PsfKernelCache *cache) {
+  *sink = (PsfEventSink){.hdr = hdr, .width = width, .height = height,
+                         .cache = cache};
+  if (cache != NULL)
+    sink->events = malloc(PSF_EVENT_SINK_CAPACITY * sizeof *sink->events);
+}
+
+static void psf_event_sink_flush(PsfEventSink *sink) {
+  for (size_t i = 0; i < sink->count; ++i)
+    splat_prepared_cached_event(sink->hdr, sink->width, sink->height,
+                                &sink->events[i], sink->cache);
+  sink->count = 0;
+}
+
+static void psf_event_sink_destroy(PsfEventSink *sink) {
+  psf_event_sink_flush(sink);
+  free(sink->events);
+}
+
+#if FRAME_PSF_EVENT_SINK
+static void psf_event_sink_emit(PsfEventSink *sink, const PsfCachedEvent *event) {
+  if (sink->events == NULL) {
+    splat_prepared_cached_event(sink->hdr, sink->width, sink->height, event, sink->cache);
+    return;
+  }
+  if (sink->count == PSF_EVENT_SINK_CAPACITY)
+    psf_event_sink_flush(sink);
+  sink->events[sink->count++] = *event;
+}
+#endif
+
 typedef struct {
   size_t a, b, triangle;
   unsigned int side;
@@ -852,6 +898,7 @@ typedef struct {
   double max_cache_psf_flux;
   double psf_relative_tail;
   double psf_min_y;
+  PsfEventSink *event_sink;
   size_t images;
   size_t direct_fallbacks;
   size_t cached_wing_clipped;
@@ -920,11 +967,26 @@ 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;
-    const int direct_fallback = splat_moffat_cached(
-        context->hdr, context->width, context->height, image_x, image_y, color,
-        flux,
-        context->psf, context->psf_cache, context->max_cache_psf_flux,
-        context->psf_relative_tail, context->psf_min_y);
+    PsfCachedEvent event;
+    const int direct_fallback = psf_prepare_cached_event(
+        &event, image_x, image_y, color, flux, context->psf, context->psf_cache,
+        context->max_cache_psf_flux, context->psf_relative_tail, context->psf_min_y);
+    if (direct_fallback == 1) {
+      /* Preserve the old per-event accumulation order exactly while this is
+       * the CPU reference sink.  A future asynchronous GPU sink must instead
+       * finish the submitted chunk before this CPU reference contribution. */
+      psf_event_sink_flush(context->event_sink);
+      splat_moffat_direct(context->hdr, context->width, context->height, image_x,
+                          image_y, color, flux, context->psf,
+                          context->psf_relative_tail, context->psf_min_y);
+    } else if (direct_fallback != 3) {
+#if FRAME_PSF_EVENT_SINK
+      psf_event_sink_emit(context->event_sink, &event);
+#else
+      splat_prepared_cached_event(context->hdr, context->width, context->height,
+                                  &event, context->psf_cache);
+#endif
+    }
     context->direct_fallbacks += direct_fallback == 1;
     context->cached_wing_clipped += direct_fallback == 2;
     context->discarded_below_min_y += direct_fallback == 3;
@@ -951,8 +1013,14 @@ 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) {
+    size_t first_triangle, size_t last_triangle, PsfEventSink *event_sink) {
   CatalogSplatStats stats = {0};
+  PsfEventSink owned_sink;
+  const int owns_sink = event_sink == NULL;
+  if (owns_sink) {
+    psf_event_sink_init(&owned_sink, hdr, width, height, psf_cache);
+    event_sink = &owned_sink;
+  }
   for (size_t t = first_triangle; t < last_triangle; ++t) {
     const LensVertex *vertex[3];
     if (!usable_triangle(mesh, &mesh->triangles[t], vertex))
@@ -986,7 +1054,8 @@ static CatalogSplatStats splat_catalog_triangles(
                                     .psf = psf, .psf_cache = psf_cache,
                                     .max_cache_psf_flux = max_cache_psf_flux,
                                     .psf_relative_tail = psf_relative_tail,
-                                    .psf_min_y = psf_min_y};
+                                    .psf_min_y = psf_min_y,
+                                    .event_sink = event_sink};
     if (catalog_visit_source_triangle(catalog, direction, 0, splat_catalog_tile,
                                       &context) == 0)
       stats.images += context.images;
@@ -994,6 +1063,8 @@ static CatalogSplatStats splat_catalog_triangles(
     stats.cached_wing_clipped += context.cached_wing_clipped;
     stats.discarded_below_min_y += context.discarded_below_min_y;
   }
+  if (owns_sink)
+    psf_event_sink_destroy(&owned_sink);
   return stats;
 }
 
@@ -1061,7 +1132,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);
+          0, mesh->triangle_count, NULL);
       copy_psf_splat_stats(psf_stats, stats);
       return stats.images;
   }
@@ -1078,7 +1149,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);
+          0, mesh->triangle_count, NULL);
       copy_psf_splat_stats(psf_stats, stats);
       return stats.images;
     }
@@ -1089,7 +1160,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);
+          0, mesh->triangle_count, NULL);
       copy_psf_splat_stats(psf_stats, stats);
       return stats.images;
     }
@@ -1106,7 +1177,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);
+          0, mesh->triangle_count, NULL);
     copy_psf_splat_stats(psf_stats, stats);
     return stats.images;
   }
@@ -1126,6 +1197,8 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
 #endif
     {
       const size_t worker = (size_t)omp_get_thread_num();
+      PsfEventSink event_sink;
+      psf_event_sink_init(&event_sink, private_hdr[worker], width, height, psf_cache);
       size_t local_triangles = 0, next_report = 8;
       progress->worker_callback(progress->context, worker, worker_count, 0, 0);
       /* Source density and lens magnification can vary by orders of magnitude
@@ -1138,7 +1211,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);
+            triangle, triangle + 1, &event_sink);
         images += stats.images;
         direct_fallbacks += stats.direct_fallbacks;
         cached_wing_clipped += stats.cached_wing_clipped;
@@ -1155,6 +1228,7 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
             next_report *= 2;
         }
       }
+      psf_event_sink_destroy(&event_sink);
       progress->worker_callback(progress->context, worker, worker_count,
                                 local_triangles, 1);
     }
@@ -1166,13 +1240,15 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
 #endif
     {
       const size_t worker = (size_t)omp_get_thread_num();
+      PsfEventSink event_sink;
+      psf_event_sink_init(&event_sink, private_hdr[worker], width, height, psf_cache);
 #pragma omp for schedule(dynamic, 1)
       for (size_t triangle = 0; triangle < mesh->triangle_count; ++triangle) {
         const CatalogSplatStats stats = splat_catalog_triangles(
             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);
+            triangle, triangle + 1, &event_sink);
         images += stats.images;
         direct_fallbacks += stats.direct_fallbacks;
         cached_wing_clipped += stats.cached_wing_clipped;
@@ -1182,6 +1258,7 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
         magnification_clamped_triangles += stats.magnification_clamped_triangles;
 #endif
       }
+      psf_event_sink_destroy(&event_sink);
     }
   }
 #pragma omp parallel for schedule(static)

src/optics.c

@@ -401,14 +401,14 @@ void splat_moffat_direct(double *hdr, int width, int height, double x, double y,
     }
 }
 
-int splat_moffat_cached(double *hdr, int width, int height, double x, double y,
-                        LinearRgb color, double flux,
-                        const PointSpreadFunction *psf,
-                        const PsfKernelCache *cache,
-                        double max_cache_psf_flux,
-                        double relative_tail_fraction, double min_y)
+int psf_prepare_cached_event(PsfCachedEvent *event, double x, double y,
+                             LinearRgb color, double flux,
+                             const PointSpreadFunction *psf,
+                             const PsfKernelCache *cache,
+                             double max_cache_psf_flux,
+                             double relative_tail_fraction, double min_y)
 {
-    if (hdr == NULL || width <= 0 || height <= 0 || flux <= 0.0 || !valid_psf(psf))
+    if (event == NULL || flux <= 0.0 || !valid_psf(psf))
         return 1;
     const double alpha = moffat_alpha(psf);
     if (!isfinite(min_y) || min_y < 0.0)
@@ -422,24 +422,58 @@ int splat_moffat_cached(double *hdr, int width, int height, double x, double y,
         cache->relative_tail_fraction != relative_tail_fraction ||
         !isfinite(support_radius) ||
         (support_radius > cache->max_radius_pixels && flux > max_cache_psf_flux)) {
-        splat_moffat_direct(hdr, width, height, x, y, color, flux, psf,
-                            relative_tail_fraction, min_y);
         return 1;
     }
-    const double base_x = floor(x), base_y = floor(y);
-    const double fx = x - base_x, fy = y - base_y;
+    *event = (PsfCachedEvent){.x = x, .y = y, .color = color, .flux = flux,
+                              .support_radius = support_radius};
+    return support_radius > cache->max_radius_pixels ? 2 : 0;
+}
+
+int splat_moffat_cached(double *hdr, int width, int height, double x, double y,
+                        LinearRgb color, double flux,
+                        const PointSpreadFunction *psf,
+                        const PsfKernelCache *cache,
+                        double max_cache_psf_flux,
+                        double relative_tail_fraction, double min_y)
+{
+    if (hdr == NULL || width <= 0 || height <= 0)
+        return 1;
+    PsfCachedEvent event;
+    const int status = psf_prepare_cached_event(
+        &event, x, y, color, flux, psf, cache, max_cache_psf_flux,
+        relative_tail_fraction, min_y);
+    if (status == 1) {
+        splat_moffat_direct(hdr, width, height, x, y, color, flux, psf,
+                            relative_tail_fraction, min_y);
+        return status;
+    }
+    if (status == 3)
+        return status;
+    splat_prepared_cached_event(hdr, width, height, &event, cache);
+    return status;
+}
+
+void splat_prepared_cached_event(double *hdr, int width, int height,
+                                 const PsfCachedEvent *event,
+                                 const PsfKernelCache *cache)
+{
+    if (hdr == NULL || width <= 0 || height <= 0 || event == NULL || cache == NULL ||
+        !cache->ready)
+        return;
+    const double base_x = floor(event->x), base_y = floor(event->y);
+    const double fx = event->x - base_x, fy = event->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)fmin(ceil(support_radius), cache->max_radius_pixels);
+    const int support = (int)fmin(ceil(event->support_radius), cache->max_radius_pixels);
     for (int offset_y = -support; offset_y <= support; ++offset_y) {
         const int py = (int)base_y + offset_y;
         if (py < 0 || py >= height)
             continue;
         int first_offset_x, last_offset_x;
-        if (!moffat_row_offset_range(support_radius, fx, fy, support, offset_y,
+        if (!moffat_row_offset_range(event->support_radius, fx, fy, support, offset_y,
                                      &first_offset_x, &last_offset_x))
             continue;
         if (first_offset_x < -(int)base_x)
@@ -455,12 +489,11 @@ int splat_moffat_cached(double *hdr, int width, int height, double x, double 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;
+            pixel[0] += event->color.r * event->flux * weight;
+            pixel[1] += event->color.g * event->flux * weight;
+            pixel[2] += event->color.b * event->flux * weight;
         }
     }
-    return support_radius > cache->max_radius_pixels ? 2 : 0;
 }
 
 void splat_moffat(double *hdr, int width, int height, double x, double y,

src/optics.h

@@ -32,6 +32,15 @@ typedef struct {
 #endif
 } PsfSplatStats;
 
+/* A cache-eligible star image after all catalog/lens/colour decisions.  This
+ * is the lossless CPU -> PSF-backend boundary; values deliberately remain
+ * double because the production HDR path is double. */
+typedef struct {
+  double x, y;
+  LinearRgb color;
+  double flux, support_radius;
+} PsfCachedEvent;
+
 /* 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);
@@ -61,6 +70,18 @@ int splat_moffat_cached(double *hdr, int width, int height, double x, double y,
                         const PsfKernelCache *cache,
                         double max_cache_psf_flux,
                         double relative_tail_fraction, double min_y);
+/* Applies precisely the eligibility rules used by splat_moffat_cached().
+ * Returns its status code; status 0 and 2 populate event for cache backends. */
+int psf_prepare_cached_event(PsfCachedEvent *event, double x, double y,
+                             LinearRgb color, double flux,
+                             const PointSpreadFunction *psf,
+                             const PsfKernelCache *cache,
+                             double max_cache_psf_flux,
+                             double relative_tail_fraction, double min_y);
+/* Consumes an event already accepted by psf_prepare_cached_event(). */
+void splat_prepared_cached_event(double *hdr, int width, int height,
+                                 const PsfCachedEvent *event,
+                                 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);