From e8deb1ff190f20c662c405e6df303fe9181853e8 Mon Sep 17 00:00:00 2001 From: Yingjie Wang Date: Sun, 30 Aug 2026 15:10:14 -0400 Subject: [PATCH] Optics: add PSF luminance cutoff --- README.md | 5 ++++ src/frame.c | 42 +++++++++++++++++++++---------- src/frame.h | 1 + src/main.c | 8 +++++- src/optics.c | 61 +++++++++++++++++++++++++++++++++++++--------- src/optics.h | 7 +++--- tests/test_frame.c | 37 ++++++++++++++++++++-------- 7 files changed, 123 insertions(+), 38 deletions(-) diff --git a/README.md b/README.md index aee3824..48b6012 100644 --- a/README.md +++ b/README.md @@ -239,6 +239,11 @@ image; it defaults to `1e-8`. Larger values intentionally shorten the Moffat support and rebuild the immutable cache at the corresponding radius, which is useful when a faster, lower-fidelity render is acceptable. +`--psf-min-y Y` defaults to `0` (disabled). A positive value is a linear-HDR +luminance cutoff: a PSF stops where its continuous Moffat Y profile falls below +`Y`, and an image whose central value is already below `Y` is omitted. The +final PSF report counts such omitted events and emits a warning when any occur. + For bounded preview renders, `--max-cache-psf-flux F` (default `1`) allows images with `1 < flux <= F` to use the existing cache instead of the direct evaluator. This does not rebuild or enlarge the cache: the image retains its diff --git a/src/frame.c b/src/frame.c index 393e4f7..0f2eaf1 100644 --- a/src/frame.c +++ b/src/frame.c @@ -851,15 +851,18 @@ typedef struct { const PsfKernelCache *psf_cache; double max_cache_psf_flux; double psf_relative_tail; + double psf_min_y; size_t images; size_t direct_fallbacks; size_t cached_wing_clipped; + size_t discarded_below_min_y; } TriangleSplatContext; typedef struct { size_t images; size_t direct_fallbacks; size_t cached_wing_clipped; + size_t discarded_below_min_y; #ifdef GR_DEBUG double max_raw_magnification; size_t magnification_clamped_triangles; @@ -872,9 +875,11 @@ static void copy_psf_splat_stats(PsfSplatStats *destination, if (destination == NULL) return; *destination = (PsfSplatStats){.cached_splats = - source.images - source.direct_fallbacks, + source.images - source.direct_fallbacks - + source.discarded_below_min_y, .cached_wing_clipped = source.cached_wing_clipped, .direct_fallbacks = source.direct_fallbacks, + .discarded_below_min_y = source.discarded_below_min_y, #ifdef GR_DEBUG .max_raw_magnification = source.max_raw_magnification, .magnification_clamped_triangles = @@ -919,9 +924,10 @@ static int splat_catalog_tile(const Star *stars, size_t count, 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_relative_tail, context->psf_min_y); context->direct_fallbacks += direct_fallback == 1; context->cached_wing_clipped += direct_fallback == 2; + context->discarded_below_min_y += direct_fallback == 3; #ifdef GR_DEBUG if (direct_fallback == 1) /* This is deliberately emitted by the active splat worker: a direct @@ -944,6 +950,7 @@ static CatalogSplatStats splat_catalog_triangles( int height, double exposure, const PointSpreadFunction *psf, 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) { CatalogSplatStats stats = {0}; for (size_t t = first_triangle; t < last_triangle; ++t) { @@ -978,12 +985,14 @@ static CatalogSplatStats splat_catalog_triangles( .exposure = exposure, .magnification = magnification, .psf = psf, .psf_cache = psf_cache, .max_cache_psf_flux = max_cache_psf_flux, - .psf_relative_tail = psf_relative_tail}; + .psf_relative_tail = psf_relative_tail, + .psf_min_y = psf_min_y}; if (catalog_visit_source_triangle(catalog, direction, 0, splat_catalog_tile, &context) == 0) stats.images += context.images; stats.direct_fallbacks += context.direct_fallbacks; stats.cached_wing_clipped += context.cached_wing_clipped; + stats.discarded_below_min_y += context.discarded_below_min_y; } return stats; } @@ -1016,6 +1025,7 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, double max_magnification, double max_cache_psf_flux, double psf_relative_tail, + double psf_min_y, int limit_workers_by_memory, int catalog_load_workers, CatalogPrefetchStats *prefetch_stats, @@ -1026,7 +1036,7 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, isnan(max_magnification) || max_magnification <= 0.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) + psf_relative_tail >= 1.0 || !isfinite(psf_min_y) || psf_min_y < 0.0) return 0; /* A bounded parallel read phase completes before splatting. Its serial cache @@ -1050,7 +1060,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, + max_magnification, max_cache_psf_flux, psf_relative_tail, psf_min_y, 0, mesh->triangle_count); copy_psf_splat_stats(psf_stats, stats); return stats.images; @@ -1067,7 +1077,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, + max_magnification, max_cache_psf_flux, psf_relative_tail, psf_min_y, 0, mesh->triangle_count); copy_psf_splat_stats(psf_stats, stats); return stats.images; @@ -1078,7 +1088,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, + max_magnification, max_cache_psf_flux, psf_relative_tail, psf_min_y, 0, mesh->triangle_count); copy_psf_splat_stats(psf_stats, stats); return stats.images; @@ -1095,13 +1105,14 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, free(private_hdr); const CatalogSplatStats stats = splat_catalog_triangles( mesh, catalog, hdr, width, height, exposure, psf, psf_cache, - max_magnification, max_cache_psf_flux, psf_relative_tail, + max_magnification, max_cache_psf_flux, psf_relative_tail, psf_min_y, 0, mesh->triangle_count); copy_psf_splat_stats(psf_stats, stats); return stats.images; } - size_t images = 0, direct_fallbacks = 0, cached_wing_clipped = 0; + size_t images = 0, direct_fallbacks = 0, cached_wing_clipped = 0, + discarded_below_min_y = 0; #ifdef GR_DEBUG double max_raw_magnification = 0.0; size_t magnification_clamped_triangles = 0; @@ -1109,9 +1120,9 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, /* Keep the ordinary render loop byte-for-byte free of progress checks. */ if (progress != NULL && progress->worker_callback != NULL) { #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) +#pragma omp parallel num_threads((int)worker_count) reduction(+ : images, direct_fallbacks, cached_wing_clipped, discarded_below_min_y, magnification_clamped_triangles) reduction(max : max_raw_magnification) #else -#pragma omp parallel num_threads((int)worker_count) reduction(+ : images, direct_fallbacks, cached_wing_clipped) +#pragma omp parallel num_threads((int)worker_count) reduction(+ : images, direct_fallbacks, cached_wing_clipped, discarded_below_min_y) #endif { const size_t worker = (size_t)omp_get_thread_num(); @@ -1126,10 +1137,12 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, 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); images += stats.images; direct_fallbacks += stats.direct_fallbacks; cached_wing_clipped += stats.cached_wing_clipped; + discarded_below_min_y += stats.discarded_below_min_y; #ifdef GR_DEBUG max_raw_magnification = fmax(max_raw_magnification, stats.max_raw_magnification); magnification_clamped_triangles += stats.magnification_clamped_triangles; @@ -1147,9 +1160,9 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, } } else { #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) +#pragma omp parallel num_threads((int)worker_count) reduction(+ : images, direct_fallbacks, cached_wing_clipped, discarded_below_min_y, magnification_clamped_triangles) reduction(max : max_raw_magnification) #else -#pragma omp parallel num_threads((int)worker_count) reduction(+ : images, direct_fallbacks, cached_wing_clipped) +#pragma omp parallel num_threads((int)worker_count) reduction(+ : images, direct_fallbacks, cached_wing_clipped, discarded_below_min_y) #endif { const size_t worker = (size_t)omp_get_thread_num(); @@ -1158,10 +1171,12 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, 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); images += stats.images; direct_fallbacks += stats.direct_fallbacks; cached_wing_clipped += stats.cached_wing_clipped; + discarded_below_min_y += stats.discarded_below_min_y; #ifdef GR_DEBUG max_raw_magnification = fmax(max_raw_magnification, stats.max_raw_magnification); magnification_clamped_triangles += stats.magnification_clamped_triangles; @@ -1180,6 +1195,7 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, .images = images, .direct_fallbacks = direct_fallbacks, .cached_wing_clipped = cached_wing_clipped, + .discarded_below_min_y = discarded_below_min_y, #ifdef GR_DEBUG .max_raw_magnification = max_raw_magnification, .magnification_clamped_triangles = diff --git a/src/frame.h b/src/frame.h index 5b59cd8..1132bd9 100644 --- a/src/frame.h +++ b/src/frame.h @@ -127,6 +127,7 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, double max_magnification, double max_cache_psf_flux, double psf_relative_tail, + double psf_min_y, int limit_workers_by_memory, int catalog_load_workers, CatalogPrefetchStats *prefetch_stats, diff --git a/src/main.c b/src/main.c index 0fb9bea..19b3775 100644 --- a/src/main.c +++ b/src/main.c @@ -25,6 +25,7 @@ typedef struct { double max_magnification; double max_cache_psf_flux; double psf_relative_tail; + double psf_min_y; double observer_radius; double observer_inward_speed; PointSpreadFunction psf; @@ -188,6 +189,7 @@ static int parse_args(int argc, char **argv, Settings *s, .max_magnification = INFINITY, .max_cache_psf_flux = 1.0, .psf_relative_tail = 1e-8, + .psf_min_y = 0.0, .observer_radius = 30.0, .psf = {2.7, 4.5}, .catalog_path = "assets/sky_grid_5deg.csv", @@ -260,6 +262,8 @@ static int parse_args(int argc, char **argv, Settings *s, } 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-min-y") && i + 1 < argc && + !parse_nonnegative(argv[++i], &s->psf_min_y)) { } else if (!strcmp(argv[i], "--psf-direct")) { s->psf_direct = 1; } else if (!strcmp(argv[i], "--verbose")) { @@ -469,6 +473,7 @@ static int render_observer_frame(const Settings *s, StarCatalog *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, + s->psf_min_y, spacetime_limits_render_workers_by_memory(spacetime), s->catalog_load_workers, &prefetch, &psf_stats, &(FrameSplatProgress){report_splat_progress, @@ -648,6 +653,7 @@ static int render_movie(const Settings *s, StarCatalog *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, + s->psf_min_y, spacetime_limits_render_workers_by_memory(spacetime), s->catalog_load_workers, &prefetch, &psf_stats, &(FrameSplatProgress){report_splat_progress, @@ -693,7 +699,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-relative-tail R] [--psf-min-y Y] " "[--psf-direct] [--verbose] " #ifdef ENABLE_HDR_OUTPUT "[--hdr-output] " diff --git a/src/optics.c b/src/optics.c index 342f8b2..c14ed16 100644 --- a/src/optics.c +++ b/src/optics.c @@ -295,10 +295,14 @@ void psf_kernel_cache_report(const PsfKernelCache *cache, if (stream == NULL) return; (void)cache; - fprintf(stream, "PSF splats: cached %zu, cached wing-clipped %zu, direct fallbacks %zu\n", + fprintf(stream, "PSF splats: cached %zu, cached wing-clipped %zu, direct fallbacks %zu, " + "discarded below min-Y %zu\n", stats == NULL ? 0u : stats->cached_splats, stats == NULL ? 0u : stats->cached_wing_clipped, - stats == NULL ? 0u : stats->direct_fallbacks); + stats == NULL ? 0u : stats->direct_fallbacks, + stats == NULL ? 0u : stats->discarded_below_min_y); + if (stats != NULL && stats->discarded_below_min_y != 0) + fputs("Warning: --psf-min-y discarded one or more PSF events.\n", stream); #ifdef GR_DEBUG if (stats != NULL) fprintf(stream, "Debug: max raw magnification %.6g; magnification-clamped " @@ -325,18 +329,49 @@ static int moffat_row_offset_range(double support_radius, double fx, double fy, return *first_offset_x <= *last_offset_x; } +static double linear_rgb_luminance(LinearRgb color) +{ + return 0.2126729 * color.r + 0.7151522 * color.g + 0.0721750 * color.b; +} + +/* `min_y` is a deliberate rendering cutoff. It compares against the + * continuous Moffat luminance profile rather than adding a second pixel-level + * cache lookup at the boundary. */ +static double moffat_min_y_radius(double alpha, double beta, LinearRgb color, + double flux, double min_y) +{ + if (min_y <= 0.0) + return INFINITY; + const double peak_y = linear_rgb_luminance(color) * flux * + (beta - 1.0) / (pi * alpha * alpha); + if (!(peak_y > min_y) || !isfinite(peak_y)) + return 0.0; + return alpha * sqrt(pow(min_y / peak_y, -1.0 / beta) - 1.0); +} + +static double moffat_effective_support_radius(double alpha, double beta, + LinearRgb color, double flux, + double relative_tail_fraction, + double min_y) +{ + const double strict_radius = moffat_support_radius( + alpha, beta, flux, relative_tail_fraction); + const double min_y_radius = moffat_min_y_radius(alpha, beta, color, flux, min_y); + return fmin(strict_radius, min_y_radius); +} + void splat_moffat_direct(double *hdr, int width, int height, double x, double y, LinearRgb color, double flux, const PointSpreadFunction *psf, - double relative_tail_fraction) + double relative_tail_fraction, double min_y) { 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, relative_tail_fraction); + const double support_radius = moffat_effective_support_radius( + alpha, psf->moffat_beta, color, flux, relative_tail_fraction, min_y); if (!isfinite(support_radius)) return; const int support = (int)ceil(support_radius); @@ -371,20 +406,24 @@ int splat_moffat_cached(double *hdr, int width, int height, double x, double y, const PointSpreadFunction *psf, const PsfKernelCache *cache, double max_cache_psf_flux, - double relative_tail_fraction) + double relative_tail_fraction, double min_y) { 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, relative_tail_fraction); + if (!isfinite(min_y) || min_y < 0.0) + return 1; + const double support_radius = moffat_effective_support_radius( + alpha, psf->moffat_beta, color, flux, relative_tail_fraction, min_y); + if (support_radius == 0.0) + return 3; if (cache == NULL || !cache->ready || cache->fwhm_pixels != psf->fwhm_pixels || 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, - relative_tail_fraction); + relative_tail_fraction, min_y); return 1; } const double base_x = floor(x), base_y = floor(y); @@ -427,10 +466,10 @@ 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, - double relative_tail_fraction) + double relative_tail_fraction, double min_y) { splat_moffat_direct(hdr, width, height, x, y, color, flux, psf, - relative_tail_fraction); + relative_tail_fraction, min_y); } static unsigned char tonemap_channel(double hdr_value) diff --git a/src/optics.h b/src/optics.h index 0c98f58..605a489 100644 --- a/src/optics.h +++ b/src/optics.h @@ -25,6 +25,7 @@ typedef struct { size_t cached_splats; size_t cached_wing_clipped; size_t direct_fallbacks; + size_t discarded_below_min_y; #ifdef GR_DEBUG double max_raw_magnification; size_t magnification_clamped_triangles; @@ -46,11 +47,11 @@ 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, - double relative_tail_fraction); + double relative_tail_fraction, double min_y); void splat_moffat(double *hdr, int width, int height, double x, double y, LinearRgb color, double flux, const PointSpreadFunction *psf, - double relative_tail_fraction); + double relative_tail_fraction, double min_y); /* 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. */ @@ -59,7 +60,7 @@ int splat_moffat_cached(double *hdr, int width, int height, double x, double y, const PointSpreadFunction *psf, const PsfKernelCache *cache, double max_cache_psf_flux, - double relative_tail_fraction); + double relative_tail_fraction, double min_y); /* 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); diff --git a/tests/test_frame.c b/tests/test_frame.c index 2d6f638..6b2f5d3 100644 --- a/tests/test_frame.c +++ b/tests/test_frame.c @@ -79,11 +79,21 @@ int main(void) { const size_t images = frame_splat_catalog(&mesh, &catalog, hdr, width, height, test_exposure, &psf, NULL, INFINITY, 1.0, psf_relative_tail, - 0, 1, NULL, NULL, NULL); + 0.0, 0, 1, NULL, NULL, NULL); if (images != 1 || hdr[3 * (50 * width + 50)] <= 0.0) { fputs("flat-space inverse lens-map regression failed\n", stderr); goto done; } + memset(hdr, 0, (size_t)width * height * 3 * sizeof *hdr); + PsfSplatStats min_y_stats = {0}; + if (frame_splat_catalog(&mesh, &catalog, hdr, width, height, test_exposure, + &psf, NULL, INFINITY, 1.0, psf_relative_tail, + 1e300, 0, 1, NULL, &min_y_stats, NULL) != 1 || + min_y_stats.discarded_below_min_y != 1 || + hdr[3 * (50 * width + 50)] != 0.0) { + fputs("PSF minimum-Y discard regression failed\n", stderr); + goto done; + } /* Private HDR accumulation must preserve the serial splat result. */ double *serial_hdr = calloc((size_t)width * height * 3, sizeof *serial_hdr); double *parallel_hdr = calloc((size_t)width * height * 3, sizeof *parallel_hdr); @@ -97,11 +107,11 @@ int main(void) { omp_set_num_threads(1); const size_t serial_images = frame_splat_catalog( &mesh, &catalog, serial_hdr, width, height, test_exposure, &psf, NULL, - INFINITY, 1.0, psf_relative_tail, 0, 1, NULL, NULL, NULL); + INFINITY, 1.0, psf_relative_tail, 0.0, 0, 1, NULL, NULL, NULL); omp_set_num_threads(4); const size_t parallel_images = frame_splat_catalog( &mesh, &catalog, parallel_hdr, width, height, test_exposure, &psf, NULL, - INFINITY, 1.0, psf_relative_tail, 0, 1, NULL, NULL, NULL); + INFINITY, 1.0, psf_relative_tail, 0.0, 0, 1, NULL, NULL, NULL); omp_set_num_threads(original_threads); for (int value = 0; value < width * height * 3; ++value) if (fabs(serial_hdr[value] - parallel_hdr[value]) > @@ -128,7 +138,7 @@ int main(void) { * beyond the former Gaussian's 3-sigma raster box. */ memset(hdr, 0, (size_t)width * height * 3 * sizeof *hdr); splat_moffat(hdr, width, height, 50.5, 50.5, - (LinearRgb){1.0, 1.0, 1.0}, 1.0, &psf, psf_relative_tail); + (LinearRgb){1.0, 1.0, 1.0}, 1.0, &psf, psf_relative_tail, 0.0); double moffat_flux = 0.0; for (int pixel = 0; pixel < width * height; ++pixel) moffat_flux += hdr[3 * pixel]; @@ -167,14 +177,14 @@ int main(void) { memset(reference_hdr, 0, (size_t)width * height * 3 * sizeof *reference_hdr); if (splat_moffat_cached(cached_hdr, width, height, 50.0 + phases[phase][0], 50.0 + phases[phase][1], (LinearRgb){1.0, 1.0, 1.0}, - 1.0, &psf, &cache, 1.0, psf_relative_tail) != 0) { + 1.0, &psf, &cache, 1.0, psf_relative_tail, 0.0) != 0) { fputs("ordinary PSF cache unexpectedly fell back\n", stderr); free(cached_hdr); free(reference_hdr); psf_kernel_cache_destroy(&cache); goto done; } splat_moffat_direct(reference_hdr, width, height, 50.0 + phases[phase][0], 50.0 + phases[phase][1], - (LinearRgb){1.0, 1.0, 1.0}, 1.0, &psf, psf_relative_tail); + (LinearRgb){1.0, 1.0, 1.0}, 1.0, &psf, psf_relative_tail, 0.0); double peak = 0.0, max_error = 0.0; for (int value = 0; value < width * height * 3; ++value) { peak = fmax(peak, reference_hdr[value]); @@ -186,11 +196,18 @@ int main(void) { goto done; } } + if (splat_moffat_cached(cached_hdr, width, height, 50.5, 50.5, + (LinearRgb){1.0, 1.0, 1.0}, 1.0, &psf, &cache, + 1.0, psf_relative_tail, 1.0) != 3) { + fputs("PSF minimum-Y cached discard regression failed\n", stderr); + free(cached_hdr); free(reference_hdr); psf_kernel_cache_destroy(&cache); + goto done; + } /* A bright event must avoid a cached hard cutoff by selecting the direct * reference path when the requested support exceeds the cache. */ if (splat_moffat_cached(cached_hdr, width, height, 50.5, 50.5, (LinearRgb){1.0, 1.0, 1.0}, 1000.0, &psf, &cache, 1.0, - psf_relative_tail) != 1) { + psf_relative_tail, 0.0) != 1) { fputs("bright PSF direct-fallback regression failed\n", stderr); free(cached_hdr); free(reference_hdr); psf_kernel_cache_destroy(&cache); goto done; @@ -199,13 +216,13 @@ int main(void) { memset(reference_hdr, 0, (size_t)width * height * 3 * sizeof *reference_hdr); if (splat_moffat_cached(cached_hdr, width, height, 50.5, 50.5, (LinearRgb){1.0, 1.0, 1.0}, 1000.0, &psf, &cache, - 1000.0, psf_relative_tail) != 2) { + 1000.0, psf_relative_tail, 0.0) != 2) { fputs("bright PSF cached-wing-clipping regression failed\n", stderr); free(cached_hdr); free(reference_hdr); psf_kernel_cache_destroy(&cache); goto done; } splat_moffat_direct(reference_hdr, width, height, 50.5, 50.5, - (LinearRgb){1.0, 1.0, 1.0}, 1000.0, &psf, psf_relative_tail); + (LinearRgb){1.0, 1.0, 1.0}, 1000.0, &psf, psf_relative_tail, 0.0); const size_t center = 3 * (50 * width + 50); if (cached_hdr[center] <= 0.0 || fabs(cached_hdr[center] - reference_hdr[center]) > @@ -237,7 +254,7 @@ int main(void) { frame_lens_mesh_trace(&fine_mesh, &spacetime, &observer, &trace) || frame_splat_catalog(&fine_mesh, &fine_catalog, hdr, width, height, test_exposure, &psf, NULL, INFINITY, 1.0, - psf_relative_tail, 0, 1, NULL, + psf_relative_tail, 0.0, 0, 1, NULL, NULL, NULL) != 1) { fputs("fine source-triangle containment regression failed\n", stderr); frame_lens_mesh_destroy(&fine_mesh);