diff --git a/Makefile b/Makefile index 45d7e30..3111d2d 100644 --- a/Makefile +++ b/Makefile @@ -6,6 +6,7 @@ ENABLE_PNG ?= 1 SPACETIME ?= minkowski BUILD_TYPE ?= Release ENABLE_HDR ?= 0 +PSF_EVENT_SINK ?= 1 ifeq ($(BUILD_TYPE),Release) BUILD_CFLAGS := -O2 -DNDEBUG @@ -39,6 +40,11 @@ CATALOG_PREFETCH_TEST_TARGET := $(BUILD_DIR)/test_catalog_prefetch .PHONY: all backend clean run test minkowski schwarzschild FORCE +ifneq ($(filter 0 1,$(PSF_EVENT_SINK)),$(PSF_EVENT_SINK)) +$(error Unknown PSF_EVENT_SINK '$(PSF_EVENT_SINK)'; choose 0 or 1) +endif +BUILD_CPPFLAGS += -DFRAME_PSF_EVENT_SINK=$(PSF_EVENT_SINK) + ifeq ($(SPACETIME),minkowski) BACKEND_CPPFLAGS := -DSPACETIME_MINKOWSKI else ifeq ($(SPACETIME),schwarzschild) @@ -50,9 +56,9 @@ endif ifeq ($(ENABLE_HDR),1) HDR_CPPFLAGS := -DENABLE_HDR_OUTPUT $(shell pkg-config --cflags cfitsio) HDR_LDLIBS := $(shell pkg-config --libs cfitsio) -HDR_BUILD_TAG := hdr +HDR_BUILD_TAG := hdr_sink$(PSF_EVENT_SINK) else ifeq ($(ENABLE_HDR),0) -HDR_BUILD_TAG := standard +HDR_BUILD_TAG := standard_sink$(PSF_EVENT_SINK) else $(error Unknown ENABLE_HDR '$(ENABLE_HDR)'; choose 0 or 1) endif diff --git a/src/frame.c b/src/frame.c index 0f2eaf1..8fdf37c 100644 --- a/src/frame.c +++ b/src/frame.c @@ -10,6 +10,10 @@ #include #include +#ifndef FRAME_PSF_EVENT_SINK +#define FRAME_PSF_EVENT_SINK 1 +#endif + /* Numerical metric backends may reserve substantial memory for slabs and * thread-local evaluators, so they retain this private-HDR allocation budget. * Analytic backends deliberately use all OpenMP render workers instead. */ @@ -115,6 +119,48 @@ int frame_lens_mesh_trace(FrameLensMesh *mesh, const SpacetimeSource *spacetime, return 0; } +enum { PSF_EVENT_SINK_CAPACITY = 16384 }; + +typedef struct { + PsfCachedEvent *events; + size_t count; + double *hdr; + int width, height; + const PsfKernelCache *cache; +} PsfEventSink; + +static void psf_event_sink_init(PsfEventSink *sink, double *hdr, int width, + int height, const PsfKernelCache *cache) { + *sink = (PsfEventSink){.hdr = hdr, .width = width, .height = height, + .cache = cache}; + if (cache != NULL) + sink->events = malloc(PSF_EVENT_SINK_CAPACITY * sizeof *sink->events); +} + +static void psf_event_sink_flush(PsfEventSink *sink) { + for (size_t i = 0; i < sink->count; ++i) + splat_prepared_cached_event(sink->hdr, sink->width, sink->height, + &sink->events[i], sink->cache); + sink->count = 0; +} + +static void psf_event_sink_destroy(PsfEventSink *sink) { + psf_event_sink_flush(sink); + free(sink->events); +} + +#if FRAME_PSF_EVENT_SINK +static void psf_event_sink_emit(PsfEventSink *sink, const PsfCachedEvent *event) { + if (sink->events == NULL) { + splat_prepared_cached_event(sink->hdr, sink->width, sink->height, event, sink->cache); + return; + } + if (sink->count == PSF_EVENT_SINK_CAPACITY) + psf_event_sink_flush(sink); + sink->events[sink->count++] = *event; +} +#endif + typedef struct { size_t a, b, triangle; unsigned int side; @@ -852,6 +898,7 @@ typedef struct { double max_cache_psf_flux; double psf_relative_tail; double psf_min_y; + PsfEventSink *event_sink; size_t images; size_t direct_fallbacks; size_t cached_wing_clipped; @@ -920,11 +967,26 @@ static int splat_catalog_tile(const Star *stars, size_t count, weights[2] * context->vertex[2]->log_frequency_ratio; const LinearRgb color = blackbody_to_linear_rgb(star->temperature_K * exp(log_g)); const double flux = context->exposure * star->amplitude * context->magnification; - const int direct_fallback = splat_moffat_cached( - context->hdr, context->width, context->height, image_x, image_y, color, - flux, - context->psf, context->psf_cache, context->max_cache_psf_flux, - context->psf_relative_tail, context->psf_min_y); + PsfCachedEvent event; + const int direct_fallback = psf_prepare_cached_event( + &event, image_x, image_y, color, flux, context->psf, context->psf_cache, + context->max_cache_psf_flux, context->psf_relative_tail, context->psf_min_y); + if (direct_fallback == 1) { + /* Preserve the old per-event accumulation order exactly while this is + * the CPU reference sink. A future asynchronous GPU sink must instead + * finish the submitted chunk before this CPU reference contribution. */ + psf_event_sink_flush(context->event_sink); + splat_moffat_direct(context->hdr, context->width, context->height, image_x, + image_y, color, flux, context->psf, + context->psf_relative_tail, context->psf_min_y); + } else if (direct_fallback != 3) { +#if FRAME_PSF_EVENT_SINK + psf_event_sink_emit(context->event_sink, &event); +#else + splat_prepared_cached_event(context->hdr, context->width, context->height, + &event, context->psf_cache); +#endif + } context->direct_fallbacks += direct_fallback == 1; context->cached_wing_clipped += direct_fallback == 2; context->discarded_below_min_y += direct_fallback == 3; @@ -951,8 +1013,14 @@ static CatalogSplatStats splat_catalog_triangles( const PsfKernelCache *psf_cache, double max_magnification, double max_cache_psf_flux, double psf_relative_tail, double psf_min_y, - size_t first_triangle, size_t last_triangle) { + size_t first_triangle, size_t last_triangle, PsfEventSink *event_sink) { CatalogSplatStats stats = {0}; + PsfEventSink owned_sink; + const int owns_sink = event_sink == NULL; + if (owns_sink) { + psf_event_sink_init(&owned_sink, hdr, width, height, psf_cache); + event_sink = &owned_sink; + } for (size_t t = first_triangle; t < last_triangle; ++t) { const LensVertex *vertex[3]; if (!usable_triangle(mesh, &mesh->triangles[t], vertex)) @@ -986,7 +1054,8 @@ static CatalogSplatStats splat_catalog_triangles( .psf = psf, .psf_cache = psf_cache, .max_cache_psf_flux = max_cache_psf_flux, .psf_relative_tail = psf_relative_tail, - .psf_min_y = psf_min_y}; + .psf_min_y = psf_min_y, + .event_sink = event_sink}; if (catalog_visit_source_triangle(catalog, direction, 0, splat_catalog_tile, &context) == 0) stats.images += context.images; @@ -994,6 +1063,8 @@ static CatalogSplatStats splat_catalog_triangles( stats.cached_wing_clipped += context.cached_wing_clipped; stats.discarded_below_min_y += context.discarded_below_min_y; } + if (owns_sink) + psf_event_sink_destroy(&owned_sink); return stats; } @@ -1061,7 +1132,7 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, const CatalogSplatStats stats = splat_catalog_triangles( mesh, catalog, hdr, width, height, exposure, psf, psf_cache, max_magnification, max_cache_psf_flux, psf_relative_tail, psf_min_y, - 0, mesh->triangle_count); + 0, mesh->triangle_count, NULL); copy_psf_splat_stats(psf_stats, stats); return stats.images; } @@ -1078,7 +1149,7 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, const CatalogSplatStats stats = splat_catalog_triangles( mesh, catalog, hdr, width, height, exposure, psf, psf_cache, max_magnification, max_cache_psf_flux, psf_relative_tail, psf_min_y, - 0, mesh->triangle_count); + 0, mesh->triangle_count, NULL); copy_psf_splat_stats(psf_stats, stats); return stats.images; } @@ -1089,7 +1160,7 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, const CatalogSplatStats stats = splat_catalog_triangles( mesh, catalog, hdr, width, height, exposure, psf, psf_cache, max_magnification, max_cache_psf_flux, psf_relative_tail, psf_min_y, - 0, mesh->triangle_count); + 0, mesh->triangle_count, NULL); copy_psf_splat_stats(psf_stats, stats); return stats.images; } @@ -1106,7 +1177,7 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, const CatalogSplatStats stats = splat_catalog_triangles( mesh, catalog, hdr, width, height, exposure, psf, psf_cache, max_magnification, max_cache_psf_flux, psf_relative_tail, psf_min_y, - 0, mesh->triangle_count); + 0, mesh->triangle_count, NULL); copy_psf_splat_stats(psf_stats, stats); return stats.images; } @@ -1126,6 +1197,8 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, #endif { const size_t worker = (size_t)omp_get_thread_num(); + PsfEventSink event_sink; + psf_event_sink_init(&event_sink, private_hdr[worker], width, height, psf_cache); size_t local_triangles = 0, next_report = 8; progress->worker_callback(progress->context, worker, worker_count, 0, 0); /* Source density and lens magnification can vary by orders of magnitude @@ -1138,7 +1211,7 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, mesh, catalog, private_hdr[worker], width, height, exposure, psf, psf_cache, max_magnification, max_cache_psf_flux, psf_relative_tail, psf_min_y, - triangle, triangle + 1); + triangle, triangle + 1, &event_sink); images += stats.images; direct_fallbacks += stats.direct_fallbacks; cached_wing_clipped += stats.cached_wing_clipped; @@ -1155,6 +1228,7 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, next_report *= 2; } } + psf_event_sink_destroy(&event_sink); progress->worker_callback(progress->context, worker, worker_count, local_triangles, 1); } @@ -1166,13 +1240,15 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, #endif { const size_t worker = (size_t)omp_get_thread_num(); + PsfEventSink event_sink; + psf_event_sink_init(&event_sink, private_hdr[worker], width, height, psf_cache); #pragma omp for schedule(dynamic, 1) for (size_t triangle = 0; triangle < mesh->triangle_count; ++triangle) { const CatalogSplatStats stats = splat_catalog_triangles( mesh, catalog, private_hdr[worker], width, height, exposure, psf, psf_cache, max_magnification, max_cache_psf_flux, psf_relative_tail, psf_min_y, - triangle, triangle + 1); + triangle, triangle + 1, &event_sink); images += stats.images; direct_fallbacks += stats.direct_fallbacks; cached_wing_clipped += stats.cached_wing_clipped; @@ -1182,6 +1258,7 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, magnification_clamped_triangles += stats.magnification_clamped_triangles; #endif } + psf_event_sink_destroy(&event_sink); } } #pragma omp parallel for schedule(static) diff --git a/src/optics.c b/src/optics.c index c14ed16..a191f42 100644 --- a/src/optics.c +++ b/src/optics.c @@ -401,14 +401,14 @@ void splat_moffat_direct(double *hdr, int width, int height, double x, double y, } } -int splat_moffat_cached(double *hdr, int width, int height, double x, double y, - LinearRgb color, double flux, - const PointSpreadFunction *psf, - const PsfKernelCache *cache, - double max_cache_psf_flux, - double relative_tail_fraction, double min_y) +int psf_prepare_cached_event(PsfCachedEvent *event, double x, double y, + LinearRgb color, double flux, + const PointSpreadFunction *psf, + const PsfKernelCache *cache, + double max_cache_psf_flux, + double relative_tail_fraction, double min_y) { - if (hdr == NULL || width <= 0 || height <= 0 || flux <= 0.0 || !valid_psf(psf)) + if (event == NULL || flux <= 0.0 || !valid_psf(psf)) return 1; const double alpha = moffat_alpha(psf); if (!isfinite(min_y) || min_y < 0.0) @@ -422,24 +422,58 @@ int splat_moffat_cached(double *hdr, int width, int height, double x, double y, cache->relative_tail_fraction != relative_tail_fraction || !isfinite(support_radius) || (support_radius > cache->max_radius_pixels && flux > max_cache_psf_flux)) { - splat_moffat_direct(hdr, width, height, x, y, color, flux, psf, - relative_tail_fraction, min_y); return 1; } - const double base_x = floor(x), base_y = floor(y); - const double fx = x - base_x, fy = y - base_y; + *event = (PsfCachedEvent){.x = x, .y = y, .color = color, .flux = flux, + .support_radius = support_radius}; + return support_radius > cache->max_radius_pixels ? 2 : 0; +} + +int splat_moffat_cached(double *hdr, int width, int height, double x, double y, + LinearRgb color, double flux, + const PointSpreadFunction *psf, + const PsfKernelCache *cache, + double max_cache_psf_flux, + double relative_tail_fraction, double min_y) +{ + if (hdr == NULL || width <= 0 || height <= 0) + return 1; + PsfCachedEvent event; + const int status = psf_prepare_cached_event( + &event, x, y, color, flux, psf, cache, max_cache_psf_flux, + relative_tail_fraction, min_y); + if (status == 1) { + splat_moffat_direct(hdr, width, height, x, y, color, flux, psf, + relative_tail_fraction, min_y); + return status; + } + if (status == 3) + return status; + splat_prepared_cached_event(hdr, width, height, &event, cache); + return status; +} + +void splat_prepared_cached_event(double *hdr, int width, int height, + const PsfCachedEvent *event, + const PsfKernelCache *cache) +{ + if (hdr == NULL || width <= 0 || height <= 0 || event == NULL || cache == NULL || + !cache->ready) + return; + const double base_x = floor(event->x), base_y = floor(event->y); + const double fx = event->x - base_x, fy = event->y - base_y; const double phase_x = fx * cache->phase_resolution; const double phase_y = fy * cache->phase_resolution; const int x0 = (int)floor(phase_x), y0 = (int)floor(phase_y); const int x1 = x0 + 1, y1 = y0 + 1; const double tx = phase_x - x0, ty = phase_y - y0; - const int support = (int)fmin(ceil(support_radius), cache->max_radius_pixels); + const int support = (int)fmin(ceil(event->support_radius), cache->max_radius_pixels); for (int offset_y = -support; offset_y <= support; ++offset_y) { const int py = (int)base_y + offset_y; if (py < 0 || py >= height) continue; int first_offset_x, last_offset_x; - if (!moffat_row_offset_range(support_radius, fx, fy, support, offset_y, + if (!moffat_row_offset_range(event->support_radius, fx, fy, support, offset_y, &first_offset_x, &last_offset_x)) continue; if (first_offset_x < -(int)base_x) @@ -455,12 +489,11 @@ int splat_moffat_cached(double *hdr, int width, int height, double x, double y, const double weight = (1.0 - ty) * ((1.0 - tx) * w00 + tx * w10) + ty * ((1.0 - tx) * w01 + tx * w11); double *pixel = &hdr[3 * (py * width + px)]; - pixel[0] += color.r * flux * weight; - pixel[1] += color.g * flux * weight; - pixel[2] += color.b * flux * weight; + pixel[0] += event->color.r * event->flux * weight; + pixel[1] += event->color.g * event->flux * weight; + pixel[2] += event->color.b * event->flux * weight; } } - return support_radius > cache->max_radius_pixels ? 2 : 0; } void splat_moffat(double *hdr, int width, int height, double x, double y, diff --git a/src/optics.h b/src/optics.h index 605a489..fb62138 100644 --- a/src/optics.h +++ b/src/optics.h @@ -32,6 +32,15 @@ typedef struct { #endif } PsfSplatStats; +/* A cache-eligible star image after all catalog/lens/colour decisions. This + * is the lossless CPU -> PSF-backend boundary; values deliberately remain + * double because the production HDR path is double. */ +typedef struct { + double x, y; + LinearRgb color; + double flux, support_radius; +} PsfCachedEvent; + /* Integrate a Planck spectrum into absolute linear-sRGB spectral radiance * (W m^-2 sr^-1), before catalog amplitude and display exposure. */ LinearRgb blackbody_to_linear_rgb(double temperature_K); @@ -61,6 +70,18 @@ int splat_moffat_cached(double *hdr, int width, int height, double x, double y, const PsfKernelCache *cache, double max_cache_psf_flux, double relative_tail_fraction, double min_y); +/* Applies precisely the eligibility rules used by splat_moffat_cached(). + * Returns its status code; status 0 and 2 populate event for cache backends. */ +int psf_prepare_cached_event(PsfCachedEvent *event, double x, double y, + LinearRgb color, double flux, + const PointSpreadFunction *psf, + const PsfKernelCache *cache, + double max_cache_psf_flux, + double relative_tail_fraction, double min_y); +/* Consumes an event already accepted by psf_prepare_cached_event(). */ +void splat_prepared_cached_event(double *hdr, int width, int height, + const PsfCachedEvent *event, + const PsfKernelCache *cache); /* Writes PNG when built with libpng; non-libpng builds use PPM fallback. */ int write_tonemapped_image(const char *path, const double *hdr, int width, int height); diff --git a/tests/test_frame.c b/tests/test_frame.c index b4b052c..e472bd0 100644 --- a/tests/test_frame.c +++ b/tests/test_frame.c @@ -222,6 +222,24 @@ int main(void) { goto done; } psf_kernel_cache_destroy(&loose_tail_cache); + /* This is the exact eligibility split that a future event sink exposes to + * HIP: cache event, CPU direct fallback, or min-Y discard. */ + PsfCachedEvent prepared = {0}; + if (psf_prepare_cached_event(&prepared, 12.25, 14.75, + (LinearRgb){1.0, 1.0, 1.0}, 1.0, &psf, &cache, + 1.0, psf_relative_tail, 0.0) != 0 || + prepared.x != 12.25 || prepared.y != 14.75 || + !(prepared.support_radius > 0.0) || + psf_prepare_cached_event(&prepared, 12.25, 14.75, + (LinearRgb){1.0, 1.0, 1.0}, 1000.0, &psf, &cache, + 1.0, psf_relative_tail, 0.0) != 1 || + psf_prepare_cached_event(&prepared, 12.25, 14.75, + (LinearRgb){1.0, 1.0, 1.0}, 1.0, &psf, &cache, + 1.0, psf_relative_tail, 1.0) != 3) { + fputs("PSF event eligibility regression failed\n", stderr); + free(cached_hdr); free(reference_hdr); psf_kernel_cache_destroy(&cache); + goto done; + } const double phases[][2] = {{0.01, 0.99}, {0.499, 0.501}, {0.999, 0.001}}; for (size_t phase = 0; phase < sizeof phases / sizeof *phases; ++phase) { memset(cached_hdr, 0, (size_t)width * height * 3 * sizeof *cached_hdr);