Optics: parameterize PSF relative tail

src/frame.c

@@ -850,6 +850,7 @@ typedef struct {
   const PointSpreadFunction *psf;
   const PsfKernelCache *psf_cache;
   double max_cache_psf_flux;
+  double psf_relative_tail;
   size_t images;
   size_t direct_fallbacks;
   size_t cached_wing_clipped;
@@ -917,7 +918,8 @@ 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->max_cache_psf_flux);
+        context->psf, context->psf_cache, context->max_cache_psf_flux,
+        context->psf_relative_tail);
     context->direct_fallbacks += direct_fallback == 1;
     context->cached_wing_clipped += direct_fallback == 2;
 #ifdef GR_DEBUG
@@ -941,7 +943,8 @@ 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) {
+    double max_cache_psf_flux, double psf_relative_tail,
+    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];
@@ -974,7 +977,8 @@ static CatalogSplatStats splat_catalog_triangles(
                                     .width = width, .height = height,
                                     .exposure = exposure, .magnification = magnification,
                                     .psf = psf, .psf_cache = psf_cache,
-                                    .max_cache_psf_flux = max_cache_psf_flux};
+                                    .max_cache_psf_flux = max_cache_psf_flux,
+                                    .psf_relative_tail = psf_relative_tail};
     if (catalog_visit_source_triangle(catalog, direction, 0, splat_catalog_tile,
                                       &context) == 0)
       stats.images += context.images;
@@ -1011,6 +1015,7 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
                            const PsfKernelCache *psf_cache,
                            double max_magnification,
                            double max_cache_psf_flux,
+                           double psf_relative_tail,
                            int limit_workers_by_memory,
                            int catalog_load_workers,
                            CatalogPrefetchStats *prefetch_stats,
@@ -1019,7 +1024,9 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
   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 ||
-      !isfinite(max_cache_psf_flux) || max_cache_psf_flux < 1.0)
+      !isfinite(max_cache_psf_flux) || max_cache_psf_flux < 1.0 ||
+      !isfinite(psf_relative_tail) || psf_relative_tail <= 0.0 ||
+         psf_relative_tail >= 1.0)
     return 0;
 
   /* A bounded parallel read phase completes before splatting. Its serial cache
@@ -1043,7 +1050,8 @@ 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, 0, mesh->triangle_count);
+          max_magnification, max_cache_psf_flux, psf_relative_tail,
+          0, mesh->triangle_count);
       copy_psf_splat_stats(psf_stats, stats);
       return stats.images;
   }
@@ -1059,7 +1067,8 @@ 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, 0, mesh->triangle_count);
+          max_magnification, max_cache_psf_flux, psf_relative_tail,
+          0, mesh->triangle_count);
       copy_psf_splat_stats(psf_stats, stats);
       return stats.images;
     }
@@ -1069,7 +1078,8 @@ 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, 0, mesh->triangle_count);
+          max_magnification, max_cache_psf_flux, psf_relative_tail,
+          0, mesh->triangle_count);
       copy_psf_splat_stats(psf_stats, stats);
       return stats.images;
     }
@@ -1084,8 +1094,9 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
       free(private_hdr[--allocated]);
     free(private_hdr);
     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);
+          mesh, catalog, hdr, width, height, exposure, psf, psf_cache,
+          max_magnification, max_cache_psf_flux, psf_relative_tail,
+          0, mesh->triangle_count);
     copy_psf_splat_stats(psf_stats, stats);
     return stats.images;
   }
@@ -1114,7 +1125,8 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
       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);
+            psf_cache, max_magnification, max_cache_psf_flux, psf_relative_tail,
+            triangle, triangle + 1);
         images += stats.images;
         direct_fallbacks += stats.direct_fallbacks;
         cached_wing_clipped += stats.cached_wing_clipped;
@@ -1145,7 +1157,8 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
       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);
+            psf_cache, max_magnification, max_cache_psf_flux, psf_relative_tail,
+            triangle, triangle + 1);
         images += stats.images;
         direct_fallbacks += stats.direct_fallbacks;
         cached_wing_clipped += stats.cached_wing_clipped;

src/frame.h

@@ -126,6 +126,7 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
                            const PsfKernelCache *psf_cache,
                            double max_magnification,
                            double max_cache_psf_flux,
+                           double psf_relative_tail,
                            int limit_workers_by_memory,
                            int catalog_load_workers,
                            CatalogPrefetchStats *prefetch_stats,

src/main.c

@@ -24,6 +24,7 @@ typedef struct {
   double horizontal_fov_deg, look_ra_deg, look_dec_deg, exposure;
   double max_magnification;
   double max_cache_psf_flux;
+  double psf_relative_tail;
   double observer_radius;
   double observer_inward_speed;
   PointSpreadFunction psf;
@@ -186,6 +187,7 @@ static int parse_args(int argc, char **argv, Settings *s,
                   .exposure = 1e-3,
                   .max_magnification = INFINITY,
                   .max_cache_psf_flux = 1.0,
+                  .psf_relative_tail = 1e-8,
                   .observer_radius = 30.0,
                   .psf = {2.7, 4.5},
                   .catalog_path = "assets/sky_grid_5deg.csv",
@@ -255,6 +257,9 @@ static int parse_args(int argc, char **argv, Settings *s,
                !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-relative-tail") && i + 1 < argc &&
+               !parse_finite_positive(argv[++i], &s->psf_relative_tail) &&
+               s->psf_relative_tail < 1.0) {
     } else if (!strcmp(argv[i], "--psf-direct")) {
       s->psf_direct = 1;
     } else if (!strcmp(argv[i], "--verbose")) {
@@ -463,6 +468,7 @@ static int render_observer_frame(const Settings *s, StarCatalog *catalog,
   size_t images = frame_splat_catalog(
       &mesh, catalog, hdr, s->width, s->height, s->exposure, &s->psf,
       &s->psf_cache, s->max_magnification, s->max_cache_psf_flux,
+      s->psf_relative_tail,
       spacetime_limits_render_workers_by_memory(spacetime),
       s->catalog_load_workers, &prefetch, &psf_stats,
       &(FrameSplatProgress){report_splat_progress,
@@ -641,6 +647,7 @@ static int render_movie(const Settings *s, StarCatalog *catalog,
     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->max_magnification, s->max_cache_psf_flux,
+        s->psf_relative_tail,
         spacetime_limits_render_workers_by_memory(spacetime),
         s->catalog_load_workers, &prefetch, &psf_stats,
         &(FrameSplatProgress){report_splat_progress,
@@ -686,6 +693,7 @@ int main(int argc, char **argv) {
             "[--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-relative-tail R] "
             "[--psf-direct] [--verbose] "
 #ifdef ENABLE_HDR_OUTPUT
             "[--hdr-output] "
@@ -749,7 +757,8 @@ int main(int argc, char **argv) {
     catalog_destroy(&catalog);
     return 1;
   }
-  if (!settings.psf_direct && psf_kernel_cache_init(&settings.psf_cache, &settings.psf))
+  if (!settings.psf_direct && psf_kernel_cache_init(&settings.psf_cache, &settings.psf,
+                                                     settings.psf_relative_tail))
     fputs("PSF cache construction failed; using direct evaluator.\n", stderr);
   psf_kernel_cache_report_ready(&settings.psf_cache, stderr);
   int result = settings.frames_dir != NULL

src/optics.c

@@ -92,7 +92,6 @@ enum {
  * 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;
@@ -125,10 +124,14 @@ static double moffat_alpha(const PointSpreadFunction *psf)
 }
 
 static double moffat_support_radius(double alpha, double beta,
-                                    double scaled_flux)
+                                    double scaled_flux,
+                                    double relative_tail_fraction)
 {
+    if (!isfinite(relative_tail_fraction) || relative_tail_fraction <= 0.0 ||
+        relative_tail_fraction >= 1.0)
+        return NAN;
     const double normalization = scaled_flux * (beta - 1.0) / (pi * alpha * alpha);
-    const double relative_tail = fmin(psf_relative_tail_fraction,
+    const double relative_tail = fmin(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);
@@ -190,13 +193,17 @@ void psf_kernel_cache_destroy(PsfKernelCache *cache)
     *cache = (PsfKernelCache){0};
 }
 
-int psf_kernel_cache_init(PsfKernelCache *cache, const PointSpreadFunction *psf)
+int psf_kernel_cache_init(PsfKernelCache *cache, const PointSpreadFunction *psf,
+                          double relative_tail_fraction)
 {
-    if (cache == NULL || !valid_psf(psf))
+    if (cache == NULL || !valid_psf(psf) ||
+        !isfinite(relative_tail_fraction) || relative_tail_fraction <= 0.0 ||
+        relative_tail_fraction >= 1.0)
         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);
+    const double requested_radius = moffat_support_radius(
+        alpha, psf->moffat_beta, 1.0, relative_tail_fraction);
     if (!isfinite(requested_radius) || requested_radius <= 0.0)
         return -1;
     /* The ordinary cache must cover the support implied by this invocation's
@@ -221,6 +228,7 @@ int psf_kernel_cache_init(PsfKernelCache *cache, const PointSpreadFunction *psf)
     PsfKernelCache building = {.fwhm_pixels = psf->fwhm_pixels,
                                .moffat_beta = psf->moffat_beta,
                                .alpha_pixels = alpha,
+                               .relative_tail_fraction = relative_tail_fraction,
                                .max_radius_pixels = radius,
                                .weights = weights,
                                .phase_resolution = PSF_PHASE_RESOLUTION,
@@ -274,9 +282,10 @@ void psf_kernel_cache_report_ready(const PsfKernelCache *cache, FILE *stream) {
         fputs("PSF cache: disabled; using the direct evaluator.\n", stream);
         return;
     }
-    fprintf(stream, "PSF cache ready: 64x64 phases, radius %.0f px, tail abs %.0e, "
-            "boundary %.0e, build %.3f s\n",
-            cache->max_radius_pixels, psf_tail_absolute_hdr_budget,
+    fprintf(stream, "PSF cache ready: 64x64 phases, radius %.0f px, relative tail %.0e, "
+            "tail abs %.0e, boundary %.0e, build %.3f s\n",
+            cache->max_radius_pixels, cache->relative_tail_fraction,
+            psf_tail_absolute_hdr_budget,
             psf_boundary_hdr_budget, cache->build_seconds);
 }
 
@@ -300,14 +309,16 @@ void psf_kernel_cache_report(const PsfKernelCache *cache,
 
 void splat_moffat_direct(double *hdr, int width, int height, double x, double y,
                          LinearRgb color, double flux,
-                         const PointSpreadFunction *psf)
+                         const PointSpreadFunction *psf,
+                         double relative_tail_fraction)
 {
     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 double support_radius = moffat_support_radius(
+        alpha, psf->moffat_beta, flux, relative_tail_fraction);
     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));
@@ -331,17 +342,21 @@ 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 max_cache_psf_flux,
+                        double relative_tail_fraction)
 {
     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);
+    const double support_radius = moffat_support_radius(
+        alpha, psf->moffat_beta, flux, relative_tail_fraction);
     if (cache == NULL || !cache->ready || cache->fwhm_pixels != psf->fwhm_pixels ||
-        cache->moffat_beta != psf->moffat_beta || !isfinite(support_radius) ||
+        cache->moffat_beta != psf->moffat_beta ||
+        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);
+        splat_moffat_direct(hdr, width, height, x, y, color, flux, psf,
+                            relative_tail_fraction);
         return 1;
     }
     const double base_x = floor(x), base_y = floor(y);
@@ -376,9 +391,11 @@ int splat_moffat_cached(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)
+                  const PointSpreadFunction *psf,
+                  double relative_tail_fraction)
 {
-    splat_moffat_direct(hdr, width, height, x, y, color, flux, psf);
+    splat_moffat_direct(hdr, width, height, x, y, color, flux, psf,
+                        relative_tail_fraction);
 }
 
 static unsigned char tonemap_channel(double hdr_value)

src/optics.h

@@ -10,10 +10,11 @@ typedef struct {
   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. */
+/* One immutable process-wide kernel for the one PSF parameter/tail-policy pair
+ * accepted by the current renderer.  Its storage remains private to optics.c. */
 typedef struct {
   double fwhm_pixels, moffat_beta, alpha_pixels;
+  double relative_tail_fraction;
   double max_radius_pixels, build_seconds;
   float *weights;
   int phase_resolution, radius_pixels;
@@ -34,7 +35,8 @@ typedef struct {
  * (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);
+                          const PointSpreadFunction *psf,
+                          double relative_tail_fraction);
 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,
@@ -43,10 +45,12 @@ void psf_kernel_cache_report(const PsfKernelCache *cache,
  * 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);
+                         const PointSpreadFunction *psf,
+                         double relative_tail_fraction);
 void splat_moffat(double *hdr, int width, int height, double x, double y,
                   LinearRgb color, double flux,
-                  const PointSpreadFunction *psf);
+                  const PointSpreadFunction *psf,
+                  double relative_tail_fraction);
 /* 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. */
@@ -54,7 +58,8 @@ 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 max_cache_psf_flux,
+                        double relative_tail_fraction);
 /* 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);