Optics: bound cache-only bright PSF splats

src/frame.c

@@ -787,10 +787,39 @@ typedef struct {
   double exposure, magnification;
   const PointSpreadFunction *psf;
   const PsfKernelCache *psf_cache;
+  double max_cache_psf_flux;
   size_t images;
   size_t direct_fallbacks;
+  size_t cached_wing_clipped;
 } TriangleSplatContext;
 
+typedef struct {
+  size_t images;
+  size_t direct_fallbacks;
+  size_t cached_wing_clipped;
+#ifdef GR_DEBUG
+  double max_raw_magnification;
+  size_t magnification_clamped_triangles;
+#endif
+} CatalogSplatStats;
+
+static void copy_psf_splat_stats(PsfSplatStats *destination,
+                                 CatalogSplatStats source)
+{
+  if (destination == NULL)
+    return;
+  *destination = (PsfSplatStats){.cached_splats =
+                                     source.images - source.direct_fallbacks,
+                                 .cached_wing_clipped = source.cached_wing_clipped,
+                                 .direct_fallbacks = source.direct_fallbacks,
+#ifdef GR_DEBUG
+                                 .max_raw_magnification = source.max_raw_magnification,
+                                 .magnification_clamped_triangles =
+                                     source.magnification_clamped_triangles,
+#endif
+  };
+}
+
 static int splat_catalog_tile(const Star *stars, size_t count,
                               int fully_contained, void *opaque) {
   TriangleSplatContext *context = opaque;
@@ -826,10 +855,11 @@ static int splat_catalog_tile(const Star *stars, size_t count,
     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->direct_fallbacks += direct_fallback;
+        context->psf, context->psf_cache, context->max_cache_psf_flux);
+    context->direct_fallbacks += direct_fallback == 1;
+    context->cached_wing_clipped += direct_fallback == 2;
 #ifdef GR_DEBUG
-    if (direct_fallback)
+    if (direct_fallback == 1)
       /* This is deliberately emitted by the active splat worker: a direct
        * fallback can be the long-running work a Debug render is waiting on.
        * Do not add a critical section here; interleaved Debug lines are more
@@ -845,15 +875,12 @@ static int splat_catalog_tile(const Star *stars, size_t count,
   return 0;
 }
 
-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 *direct_fallbacks) {
-  size_t images = 0;
+static CatalogSplatStats splat_catalog_triangles(
+    const FrameLensMesh *mesh, StarCatalog *catalog, double *hdr, int width,
+    int height, double exposure, const PointSpreadFunction *psf,
+    const PsfKernelCache *psf_cache, double max_magnification,
+    double max_cache_psf_flux, size_t first_triangle, size_t last_triangle) {
+  CatalogSplatStats stats = {0};
   for (size_t t = first_triangle; t < last_triangle; ++t) {
     const LensVertex *vertex[3];
     if (!usable_triangle(mesh, &mesh->triangles[t], vertex))
@@ -863,7 +890,18 @@ static size_t splat_catalog_triangles(const FrameLensMesh *mesh,
     const double image_area =
         spherical_area(vertex[0]->camera_direction, vertex[1]->camera_direction,
                        vertex[2]->camera_direction);
-    const double magnification = image_area / source_area;
+    const double raw_magnification = image_area / source_area;
+    double magnification = raw_magnification;
+#ifdef GR_DEBUG
+    stats.max_raw_magnification = fmax(stats.max_raw_magnification,
+                                       raw_magnification);
+#endif
+    if (raw_magnification > max_magnification) {
+      magnification = max_magnification;
+#ifdef GR_DEBUG
+      ++stats.magnification_clamped_triangles;
+#endif
+    }
     const double direction[3][3] = {
         {vertex[0]->n_infinity[0], vertex[0]->n_infinity[1], vertex[0]->n_infinity[2]},
         {vertex[1]->n_infinity[0], vertex[1]->n_infinity[1], vertex[1]->n_infinity[2]},
@@ -873,13 +911,15 @@ static size_t splat_catalog_triangles(const FrameLensMesh *mesh,
                                     .triangle_index = t, .hdr = hdr,
                                     .width = width, .height = height,
                                     .exposure = exposure, .magnification = magnification,
-                                    .psf = psf, .psf_cache = psf_cache};
+                                    .psf = psf, .psf_cache = psf_cache,
+                                    .max_cache_psf_flux = max_cache_psf_flux};
     if (catalog_visit_source_triangle(catalog, direction, 0, splat_catalog_tile,
                                       &context) == 0)
-      images += context.images;
-    *direct_fallbacks += context.direct_fallbacks;
+      stats.images += context.images;
+    stats.direct_fallbacks += context.direct_fallbacks;
+    stats.cached_wing_clipped += context.cached_wing_clipped;
   }
-  return images;
+  return stats;
 }
 
 static void prefetch_catalog_for_mesh(const FrameLensMesh *mesh,
@@ -907,13 +947,17 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
                            int height, double exposure,
                            const PointSpreadFunction *psf,
                            const PsfKernelCache *psf_cache,
+                           double max_magnification,
+                           double max_cache_psf_flux,
                            int limit_workers_by_memory,
                            int catalog_load_workers,
                            CatalogPrefetchStats *prefetch_stats,
                            PsfSplatStats *psf_stats,
                            const FrameSplatProgress *progress) {
   if (mesh == NULL || catalog == NULL || hdr == NULL || exposure <= 0.0 ||
-      psf == NULL || width <= 0 || height <= 0 || catalog_load_workers <= 0)
+      psf == NULL || width <= 0 || height <= 0 || catalog_load_workers <= 0 ||
+      isnan(max_magnification) || max_magnification <= 0.0 ||
+      !isfinite(max_cache_psf_flux) || max_cache_psf_flux < 1.0)
     return 0;
 
   /* A bounded parallel read phase completes before splatting. Its serial cache
@@ -935,12 +979,11 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
   const size_t pixel_count = (size_t)width * height * 3;
   if (pixel_count > SIZE_MAX / sizeof(double))
     {
-      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 CatalogSplatStats stats = splat_catalog_triangles(
+          mesh, catalog, hdr, width, height, exposure, psf, psf_cache,
+          max_magnification, max_cache_psf_flux, 0, mesh->triangle_count);
+      copy_psf_splat_stats(psf_stats, stats);
+      return stats.images;
   }
   const size_t buffer_bytes = pixel_count * sizeof(double);
   const int max_threads = omp_get_max_threads();
@@ -952,23 +995,21 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
   }
   if (worker_count < 2 || worker_count > INT_MAX)
     {
-      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 CatalogSplatStats stats = splat_catalog_triangles(
+          mesh, catalog, hdr, width, height, exposure, psf, psf_cache,
+          max_magnification, max_cache_psf_flux, 0, mesh->triangle_count);
+      copy_psf_splat_stats(psf_stats, stats);
+      return stats.images;
     }
 
   double **private_hdr = calloc(worker_count, sizeof *private_hdr);
   if (private_hdr == NULL)
     {
-      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 CatalogSplatStats stats = splat_catalog_triangles(
+          mesh, catalog, hdr, width, height, exposure, psf, psf_cache,
+          max_magnification, max_cache_psf_flux, 0, mesh->triangle_count);
+      copy_psf_splat_stats(psf_stats, stats);
+      return stats.images;
     }
   size_t allocated = 0;
   for (; allocated < worker_count; ++allocated) {
@@ -980,19 +1021,25 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
     while (allocated > 0)
       free(private_hdr[--allocated]);
     free(private_hdr);
-    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 CatalogSplatStats stats = splat_catalog_triangles(
+        mesh, catalog, hdr, width, height, exposure, psf, psf_cache,
+        max_magnification, max_cache_psf_flux, 0, mesh->triangle_count);
+    copy_psf_splat_stats(psf_stats, stats);
+    return stats.images;
   }
 
-  size_t images = 0, direct_fallbacks = 0;
+  size_t images = 0, direct_fallbacks = 0, cached_wing_clipped = 0;
+#ifdef GR_DEBUG
+  double max_raw_magnification = 0.0;
+  size_t magnification_clamped_triangles = 0;
+#endif
   /* Keep the ordinary render loop byte-for-byte free of progress checks. */
   if (progress != NULL && progress->worker_callback != NULL) {
-#pragma omp parallel num_threads((int)worker_count) reduction(+ : images, direct_fallbacks)
+#ifdef GR_DEBUG
+#pragma omp parallel num_threads((int)worker_count) reduction(+ : images, direct_fallbacks, cached_wing_clipped, magnification_clamped_triangles) reduction(max : max_raw_magnification)
+#else
+#pragma omp parallel num_threads((int)worker_count) reduction(+ : images, direct_fallbacks, cached_wing_clipped)
+#endif
     {
       const size_t worker = (size_t)omp_get_thread_num();
       size_t local_triangles = 0, next_report = 8;
@@ -1003,9 +1050,16 @@ 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, psf_cache, triangle,
-                                          triangle + 1, &direct_fallbacks);
+        const CatalogSplatStats stats = splat_catalog_triangles(
+            mesh, catalog, private_hdr[worker], width, height, exposure, psf,
+            psf_cache, max_magnification, max_cache_psf_flux, triangle, triangle + 1);
+        images += stats.images;
+        direct_fallbacks += stats.direct_fallbacks;
+        cached_wing_clipped += stats.cached_wing_clipped;
+#ifdef GR_DEBUG
+        max_raw_magnification = fmax(max_raw_magnification, stats.max_raw_magnification);
+        magnification_clamped_triangles += stats.magnification_clamped_triangles;
+#endif
         ++local_triangles;
         if (local_triangles == next_report) {
           progress->worker_callback(progress->context, worker, worker_count,
@@ -1018,14 +1072,26 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
                                 local_triangles, 1);
     }
   } else {
-#pragma omp parallel num_threads((int)worker_count) reduction(+ : images, direct_fallbacks)
+#ifdef GR_DEBUG
+#pragma omp parallel num_threads((int)worker_count) reduction(+ : images, direct_fallbacks, cached_wing_clipped, magnification_clamped_triangles) reduction(max : max_raw_magnification)
+#else
+#pragma omp parallel num_threads((int)worker_count) reduction(+ : images, direct_fallbacks, cached_wing_clipped)
+#endif
     {
       const size_t worker = (size_t)omp_get_thread_num();
 #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, psf_cache, triangle,
-                                          triangle + 1, &direct_fallbacks);
+      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, triangle, triangle + 1);
+        images += stats.images;
+        direct_fallbacks += stats.direct_fallbacks;
+        cached_wing_clipped += stats.cached_wing_clipped;
+#ifdef GR_DEBUG
+        max_raw_magnification = fmax(max_raw_magnification, stats.max_raw_magnification);
+        magnification_clamped_triangles += stats.magnification_clamped_triangles;
+#endif
+      }
     }
   }
 #pragma omp parallel for schedule(static)
@@ -1035,8 +1101,16 @@ 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};
+  copy_psf_splat_stats(psf_stats, (CatalogSplatStats){
+                                        .images = images,
+                                        .direct_fallbacks = direct_fallbacks,
+                                        .cached_wing_clipped = cached_wing_clipped,
+#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);

src/frame.h

@@ -124,6 +124,8 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
                            int height, double exposure,
                            const PointSpreadFunction *psf,
                            const PsfKernelCache *psf_cache,
+                           double max_magnification,
+                           double max_cache_psf_flux,
                            int limit_workers_by_memory,
                            int catalog_load_workers,
                            CatalogPrefetchStats *prefetch_stats,

src/main.c

@@ -22,6 +22,8 @@ typedef struct {
   int psf_direct;
   int verbose;
   double horizontal_fov_deg, look_ra_deg, look_dec_deg, exposure;
+  double max_magnification;
+  double max_cache_psf_flux;
   double observer_radius;
   double observer_inward_speed;
   PointSpreadFunction psf;
@@ -112,6 +114,20 @@ static int parse_moffat_beta(const char *text, double *value) {
   return errno || *end || *value <= 1.0 ? -1 : 0;
 }
 
+static int parse_at_least_one(const char *text, double *value) {
+  char *end;
+  errno = 0;
+  *value = strtod(text, &end);
+  return errno || *end || !isfinite(*value) || *value < 1.0 ? -1 : 0;
+}
+
+static int parse_finite_positive(const char *text, double *value) {
+  char *end;
+  errno = 0;
+  *value = strtod(text, &end);
+  return errno || *end || !isfinite(*value) || *value <= 0.0 ? -1 : 0;
+}
+
 static int validate_tonemapped_output_path(const char *path) {
   const size_t path_length = strlen(path);
 #ifdef ENABLE_PNG
@@ -147,6 +163,8 @@ static int parse_args(int argc, char **argv, Settings *s,
                   .look_ra_deg = 270.0,
                   .look_dec_deg = 0.0,
                   .exposure = 1e-3,
+                  .max_magnification = INFINITY,
+                  .max_cache_psf_flux = 1.0,
                   .observer_radius = 30.0,
                   .psf = {2.7, 4.5},
                   .catalog_path = "assets/sky_grid_5deg.csv",
@@ -212,6 +230,10 @@ 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], "--max-magnification") && i + 1 < argc &&
+               !parse_finite_positive(argv[++i], &s->max_magnification)) {
+    } else if (!strcmp(argv[i], "--max-cache-psf-flux") && i + 1 < argc &&
+               !parse_at_least_one(argv[++i], &s->max_cache_psf_flux)) {
     } else if (!strcmp(argv[i], "--psf-direct")) {
       s->psf_direct = 1;
     } else if (!strcmp(argv[i], "--verbose")) {
@@ -419,7 +441,8 @@ static int render_observer_frame(const Settings *s, StarCatalog *catalog,
                              .all_sky_catalog = catalog->kind == STAR_CATALOG_ALL_SKY};
   size_t images = frame_splat_catalog(
       &mesh, catalog, hdr, s->width, s->height, s->exposure, &s->psf,
-      &s->psf_cache, spacetime_limits_render_workers_by_memory(spacetime),
+      &s->psf_cache, s->max_magnification, s->max_cache_psf_flux,
+      spacetime_limits_render_workers_by_memory(spacetime),
       s->catalog_load_workers, &prefetch, &psf_stats,
       &(FrameSplatProgress){report_splat_progress,
                             s->verbose ? report_splat_worker_progress : NULL,
@@ -596,7 +619,7 @@ static int render_movie(const Settings *s, StarCatalog *catalog,
         .all_sky_catalog = catalog->kind == STAR_CATALOG_ALL_SKY};
     const size_t images = frame_splat_catalog(
         &movie.frames[i].mesh, catalog, hdr, s->width, s->height, s->exposure,
-        &s->psf, &s->psf_cache,
+        &s->psf, &s->psf_cache, s->max_magnification, s->max_cache_psf_flux,
         spacetime_limits_render_workers_by_memory(spacetime),
         s->catalog_load_workers, &prefetch, &psf_stats,
         &(FrameSplatProgress){report_splat_progress,
@@ -641,6 +664,7 @@ 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] "
+            "[--max-magnification M] [--max-cache-psf-flux F] "
             "[--psf-direct] [--verbose] "
 #ifdef ENABLE_HDR_DEBUG
             "[--hdr-output PATH] "

src/optics.c

@@ -286,9 +286,16 @@ void psf_kernel_cache_report(const PsfKernelCache *cache,
   if (stream == NULL)
     return;
   (void)cache;
-  fprintf(stream, "PSF splats: cached %zu, direct fallbacks %zu\n",
+  fprintf(stream, "PSF splats: cached %zu, cached wing-clipped %zu, direct fallbacks %zu\n",
           stats == NULL ? 0u : stats->cached_splats,
+          stats == NULL ? 0u : stats->cached_wing_clipped,
           stats == NULL ? 0u : stats->direct_fallbacks);
+#ifdef GR_DEBUG
+  if (stats != NULL)
+    fprintf(stream, "Debug: max raw magnification %.6g; magnification-clamped "
+                    "triangles %zu\n",
+            stats->max_raw_magnification, stats->magnification_clamped_triangles);
+#endif
 }
 
 void splat_moffat_direct(double *hdr, int width, int height, double x, double y,
@@ -323,7 +330,8 @@ 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)
+                        const PsfKernelCache *cache,
+                        double max_cache_psf_flux)
 {
     if (hdr == NULL || width <= 0 || height <= 0 || flux <= 0.0 || !valid_psf(psf))
         return 1;
@@ -332,7 +340,7 @@ int splat_moffat_cached(double *hdr, int width, int height, double x, double y,
                                                         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) {
+        (support_radius > cache->max_radius_pixels && flux > max_cache_psf_flux)) {
         splat_moffat_direct(hdr, width, height, x, y, color, flux, psf);
         return 1;
     }
@@ -343,7 +351,7 @@ int splat_moffat_cached(double *hdr, int width, int height, double x, double y,
     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);
+    const int support = (int)fmin(ceil(support_radius), cache->max_radius_pixels);
     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;
@@ -363,7 +371,7 @@ int splat_moffat_cached(double *hdr, int width, int height, double x, double y,
             pixel[1] += color.g * flux * weight;
             pixel[2] += color.b * flux * weight;
         }
-    return 0;
+    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

@@ -22,7 +22,12 @@ typedef struct {
 
 typedef struct {
   size_t cached_splats;
+  size_t cached_wing_clipped;
   size_t direct_fallbacks;
+#ifdef GR_DEBUG
+  double max_raw_magnification;
+  size_t magnification_clamped_triangles;
+#endif
 } PsfSplatStats;
 
 /* Integrate a Planck spectrum into absolute linear-sRGB spectral radiance
@@ -42,11 +47,14 @@ void splat_moffat_direct(double *hdr, int width, int height, double x, double y,
 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. */
+/* Returns 0 for a complete cached splat, 1 for the direct reference fallback,
+ * and 2 when the cached core was used with its outer wing intentionally
+ * clipped.  max_cache_psf_flux == 1 preserves the historical behavior. */
 int splat_moffat_cached(double *hdr, int width, int height, double x, double y,
                         LinearRgb color, double flux,
                         const PointSpreadFunction *psf,
-                        const PsfKernelCache *cache);
+                        const PsfKernelCache *cache,
+                        double max_cache_psf_flux);
 /* 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);