From e47005ea6bfed484105cabf1cab393bcbcfeffdd Mon Sep 17 00:00:00 2001 From: Yingjie Wang Date: Sat, 29 Aug 2026 20:02:45 -0400 Subject: [PATCH] Optics: bound cache-only bright PSF splats --- README.md | 8 ++ src/frame.c | 184 +++++++++++++++++++++++++++++++-------------- src/frame.h | 2 + src/main.c | 28 ++++++- src/optics.c | 18 +++-- src/optics.h | 12 ++- tests/test_frame.c | 36 +++++++-- 7 files changed, 216 insertions(+), 72 deletions(-) diff --git a/README.md b/README.md index c4318f6..b2fc820 100644 --- a/README.md +++ b/README.md @@ -234,6 +234,14 @@ cached/direct-fallback image counts. Pass `--psf-direct` to use the slower whose required HDR-tail support exceeds the cache radius selects that reference path automatically. +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 +true core flux, color, and sub-pixel position, while its Moffat wing is clipped +at the existing cache radius. Images above `F` retain the direct fallback. +`--max-magnification M` (default unlimited) caps the per-triangle rendering +magnification before flux is formed; it is an explicit preview approximation. + The PSF-cache completion line is printed before tracing and catalog splatting begin. For long renders, pass `--verbose` to print catalog-prefetch state, splat-worker local heartbeats (8, 16, 32, ... completed triangles per worker), diff --git a/src/frame.c b/src/frame.c index ce2860b..7c8e368 100644 --- a/src/frame.c +++ b/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); diff --git a/src/frame.h b/src/frame.h index 4c055d6..fd825e2 100644 --- a/src/frame.h +++ b/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, diff --git a/src/main.c b/src/main.c index de9350a..51d86c6 100644 --- a/src/main.c +++ b/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] " diff --git a/src/optics.c b/src/optics.c index 32dccb8..e41517d 100644 --- a/src/optics.c +++ b/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, diff --git a/src/optics.h b/src/optics.h index c8ee4f8..8158db6 100644 --- a/src/optics.h +++ b/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); diff --git a/tests/test_frame.c b/tests/test_frame.c index 537746f..a108914 100644 --- a/tests/test_frame.c +++ b/tests/test_frame.c @@ -27,7 +27,7 @@ int main(void) { goto done; const size_t images = frame_splat_catalog(&mesh, &catalog, hdr, width, height, test_exposure, - &psf, NULL, 0, 1, NULL, NULL, NULL); + &psf, NULL, INFINITY, 1.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; @@ -44,12 +44,12 @@ int main(void) { omp_set_dynamic(0); omp_set_num_threads(1); const size_t serial_images = frame_splat_catalog( - &mesh, &catalog, serial_hdr, width, height, test_exposure, &psf, NULL, 0, - 1, NULL, NULL, NULL); + &mesh, &catalog, serial_hdr, width, height, test_exposure, &psf, NULL, + INFINITY, 1.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, 0, - 1, NULL, NULL, NULL); + &mesh, &catalog, parallel_hdr, width, height, test_exposure, &psf, NULL, + INFINITY, 1.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]) > @@ -105,7 +105,7 @@ 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) != 0) { + 1.0, &psf, &cache, 1.0) != 0) { fputs("ordinary PSF cache unexpectedly fell back\n", stderr); free(cached_hdr); free(reference_hdr); psf_kernel_cache_destroy(&cache); goto done; @@ -127,11 +127,30 @@ int main(void) { /* 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) { + (LinearRgb){1.0, 1.0, 1.0}, 1000.0, &psf, &cache, 1.0) != 1) { fputs("bright PSF direct-fallback regression failed\n", stderr); free(cached_hdr); free(reference_hdr); psf_kernel_cache_destroy(&cache); goto done; } + memset(cached_hdr, 0, (size_t)width * height * 3 * sizeof *cached_hdr); + 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) != 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); + const size_t center = 3 * (50 * width + 50); + if (cached_hdr[center] <= 0.0 || + fabs(cached_hdr[center] - reference_hdr[center]) > + 4e-5 * reference_hdr[center]) { + fputs("cached-wing clipping changed the bright PSF core\n", stderr); + free(cached_hdr); free(reference_hdr); psf_kernel_cache_destroy(&cache); + goto done; + } free(cached_hdr); free(reference_hdr); psf_kernel_cache_destroy(&cache); @@ -154,7 +173,8 @@ int main(void) { if (frame_lens_mesh_build_coarse(&fine_mesh, width, height, 1, 0.1) || frame_lens_mesh_trace(&fine_mesh, &spacetime, &observer, &trace) || frame_splat_catalog(&fine_mesh, &fine_catalog, hdr, width, height, - test_exposure, &psf, NULL, 0, 1, NULL, NULL, NULL) != 1) { + test_exposure, &psf, NULL, INFINITY, 1.0, 0, 1, NULL, + NULL, NULL) != 1) { fputs("fine source-triangle containment regression failed\n", stderr); frame_lens_mesh_destroy(&fine_mesh); goto done;