From 7a80086f63e8b1e5f27d9524ba0455f2d6e49be5 Mon Sep 17 00:00:00 2001 From: Yingjie Wang Date: Sat, 5 Sep 2026 04:23:57 -0400 Subject: [PATCH] Feat: add initial HIP PSF backend --- Makefile | 45 ++++++- README.md | 20 ++++ gpu_psf_acceleration_plan.md | 26 ++-- src/frame.c | 148 ++++++++++++++++++++--- src/hip_psf.h | 44 +++++++ src/hip_psf.hip | 173 +++++++++++++++++++++++++++ src/main.c | 14 +++ src/optics.h | 8 ++ tests/data/hip_psf_mixed_catalog.csv | 4 + tests/test_hip_psf.hip | 60 ++++++++++ 10 files changed, 512 insertions(+), 30 deletions(-) create mode 100644 src/hip_psf.h create mode 100644 src/hip_psf.hip create mode 100644 tests/data/hip_psf_mixed_catalog.csv create mode 100644 tests/test_hip_psf.hip diff --git a/Makefile b/Makefile index 3111d2d..853e7fe 100644 --- a/Makefile +++ b/Makefile @@ -7,6 +7,9 @@ SPACETIME ?= minkowski BUILD_TYPE ?= Release ENABLE_HDR ?= 0 PSF_EVENT_SINK ?= 1 +PSF_BACKEND ?= cpu +HIPCC ?= hipcc +HIP_CXXFLAGS ?= -std=c++17 -O2 -Wall -Wextra -Wpedantic ifeq ($(BUILD_TYPE),Release) BUILD_CFLAGS := -O2 -DNDEBUG @@ -37,14 +40,24 @@ FRAME_TEST_TARGET := $(BUILD_DIR)/test_frame SCHWARZSCHILD_TEST_TARGET := $(BUILD_DIR)/test_schwarzschild OBSERVER_TRACK_TEST_TARGET := $(BUILD_DIR)/test_observer_track CATALOG_PREFETCH_TEST_TARGET := $(BUILD_DIR)/test_catalog_prefetch +HIP_PSF_TEST_TARGET := $(BUILD_DIR)/test_hip_psf -.PHONY: all backend clean run test minkowski schwarzschild FORCE +.PHONY: all backend clean run test hip-psf-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 ($(PSF_BACKEND),cpu) +RENDER_LINKER := $(CC) +else ifeq ($(PSF_BACKEND),hip) +RENDER_LINKER := $(HIPCC) +BUILD_CPPFLAGS += -DPSF_BACKEND_HIP +else +$(error Unknown PSF_BACKEND '$(PSF_BACKEND)'; choose cpu or hip) +endif + ifeq ($(SPACETIME),minkowski) BACKEND_CPPFLAGS := -DSPACETIME_MINKOWSKI else ifeq ($(SPACETIME),schwarzschild) @@ -56,16 +69,24 @@ 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_sink$(PSF_EVENT_SINK) +HDR_BUILD_TAG := hdr_sink$(PSF_EVENT_SINK)_$(PSF_BACKEND) else ifeq ($(ENABLE_HDR),0) -HDR_BUILD_TAG := standard_sink$(PSF_EVENT_SINK) +HDR_BUILD_TAG := standard_sink$(PSF_EVENT_SINK)_$(PSF_BACKEND) else $(error Unknown ENABLE_HDR '$(ENABLE_HDR)'; choose 0 or 1) endif +ifeq ($(PSF_BACKEND),hip) +TARGET := $(BUILD_DIR)/$(TARGET_BASENAME)_hip +else TARGET := $(BUILD_DIR)/$(TARGET_BASENAME) +endif RENDER_SOURCES := $(COMMON_SOURCES) $(PROVIDER_SOURCE) src/main.c RENDER_OBJECTS := $(patsubst %.c,$(OBJECT_DIR)/$(HDR_BUILD_TAG)/%.o,$(RENDER_SOURCES)) +ifeq ($(PSF_BACKEND),hip) +HIP_PSF_OBJECT := $(OBJECT_DIR)/$(HDR_BUILD_TAG)/src/hip_psf.o +RENDER_OBJECTS += $(HIP_PSF_OBJECT) +endif RENDER_DEPS := $(RENDER_OBJECTS:.o=.d) # With no explicit backend choice, all means every currently supported @@ -88,7 +109,7 @@ schwarzschild: backend: $(TARGET) $(TARGET): $(RENDER_OBJECTS) FORCE | $(BUILD_DIR) - $(CC) $(BUILD_CFLAGS) $(OPENMP_FLAGS) $(RENDER_OBJECTS) $(LDLIBS) $(HDR_LDLIBS) -o $@ + $(RENDER_LINKER) $(BUILD_CFLAGS) $(OPENMP_FLAGS) $(RENDER_OBJECTS) $(LDLIBS) $(HDR_LDLIBS) -o $@ FORCE: @@ -99,6 +120,22 @@ $(OBJECT_DIR)/$(HDR_BUILD_TAG)/%.o: %.c $(BUILD_DIR): mkdir -p $@ +$(OBJECT_DIR)/$(HDR_BUILD_TAG)/src/hip_psf.o: src/hip_psf.hip src/hip_psf.h src/optics.h + @mkdir -p $(dir $@) + $(HIPCC) $(HIP_CXXFLAGS) $(BUILD_CPPFLAGS) -Isrc -c $< -o $@ + +# This is a production-sink regression, not the removed float analytic spike. +# It requires PSF_BACKEND=hip and an actual HIP GPU agent at execution time. +ifeq ($(PSF_BACKEND),hip) +hip-psf-test: $(HIP_PSF_TEST_TARGET) + +$(HIP_PSF_TEST_TARGET): tests/test_hip_psf.hip $(HIP_PSF_OBJECT) $(OBJECT_DIR)/$(HDR_BUILD_TAG)/src/optics.o | $(BUILD_DIR) + $(HIPCC) $(HIP_CXXFLAGS) $(OPENMP_FLAGS) -Isrc -x hip $< -x none $(filter-out $<,$^) $(LDLIBS) -o $@ +else +hip-psf-test: + $(error hip-psf-test requires PSF_BACKEND=hip) +endif + run: $(TARGET) mkdir -p output/imgs ./$(TARGET) --catalog assets/sky_grid_5deg.csv --output output/imgs/$(SPACETIME)_sky.$(IMAGE_EXT) diff --git a/README.md b/README.md index 9a71843..1d17583 100644 --- a/README.md +++ b/README.md @@ -23,6 +23,26 @@ equator to the northern hemisphere, so boundary stars have a deterministic color. The default output is PNG at `output/imgs/minkowski_sky.png`. +### Optional HIP PSF backend + +The default `PSF_BACKEND=cpu` uses the established OpenMP/private-HDR path. +Build the direct-atomic HIP PSF backend explicitly with +`PSF_BACKEND=hip`; it requires HIP/ROCm and a usable GPU agent, and writes a +separate `_hip` binary so it cannot overwrite the CPU renderer: + +```sh +make PSF_BACKEND=hip SPACETIME=minkowski backend +./build/Release/minkowski_sky_hip --catalog assets/sky_grid_5deg.csv \ + --output output/imgs/minkowski_sky_hip.png +``` + +The HIP backend accelerates only cache-eligible PSF events. Existing direct +fallbacks remain CPU reference evaluations, with an ordered HDR transfer before +and after each fallback. A HIP initialization, upload, kernel, or download +error terminates the render; it never switches to the CPU backend silently. +Use `make PSF_BACKEND=hip SPACETIME=minkowski hip-psf-test` for the small +GPU-vs-CPU cache-HDR regression. + ## Movie PNG sequence (Phase A) Movie mode consumes a canonical observer-track CSV rather than a fixed camera. diff --git a/gpu_psf_acceleration_plan.md b/gpu_psf_acceleration_plan.md index d34263a..df40426 100644 --- a/gpu_psf_acceleration_plan.md +++ b/gpu_psf_acceleration_plan.md @@ -138,8 +138,8 @@ GPU 后端必须先在小、固定 catalog snapshot 上同 CPU 的 `--psf-direct 2. 记录每通道绝对/相对误差、最大像素误差、总 RGB/Y 通量差和事件统计。 3. 单星（不同 phase、亮/暗、边缘和大 support）、稀疏场、密集重叠场各有测试。 4. 验证 `--psf-relative-tail`、`--psf-min-y`、cache/direct fallback 的统计及图像语义一致。 -5. GPU 不可用、feature 不足、初始化失败或运行错误时，明确报告并回退 CPU；不得默默输出 - 部分帧。 +5. 当用户显式选择 GPU 后端时，GPU 不可用、feature 不足、初始化失败或运行错误必须明确报错 + 并终止该次渲染；不得回退 CPU 或默默输出部分帧。未选择 GPU 时，CPU 路径仍是独立基线。 6. 以端到端 wall time 为准，同时记录 CPU 线程数、GPU、driver/runtime、块大小、分辨率、 catalog snapshot、Git hash 和参数。 @@ -160,15 +160,21 @@ GPU 后端必须先在小、固定 catalog snapshot 上同 CPU 的 `--psf-direct 可执行；它没有 cache 权重、双精度 HDR 或 production integration。其 `test-hip` target 和测试 已移除，不能作为 renderer GPU backend 的验证或运行入口。 -下一步是仅将现有 `PsfCachedEvent` 块与 `PsfKernelCache::weights` 接入 HIP： +HIP direct-atomic 的当前实现状态： -1. 定义 production HIP sink 的 C ABI，上传 double event 与 immutable cache 权重，并复用 - device event/HDR buffer。 -2. 保持 CPU 对 direct fallback 的处理和提交顺序；每个 cache 事件仅在 HIP 端按现有 bilinear - phase cache、圆形支持和 wing clip 规则累积到 double HDR。 -3. 用固定 catalog reference 比较 HDR RGB/Y、最大像素误差、总通量与统计；GPU atomic 的 - 非确定加法顺序使用事先记录的容差，而非 bytewise 比较。 -4. 记录端到端基准的完整命令和原始输出：Git hash、CPU threads、GPU/runtime、块大小、上传、 +1. 已定义 production HIP sink 的 C ABI，上传 double event 与 immutable cache 权重，并复用 + device/HDR buffer。正式构建选项为 `PSF_BACKEND=cpu|hip`；HIP 选择会把该层链接到 renderer， + 而 CPU 选择不要求 HIP 工具链或 runtime。 + 在 RX 9070（gfx1201）上，`make -B PSF_BACKEND=hip SPACETIME=minkowski hip-psf-test` + 的三事件重叠 cache 对照给出 `max_abs=3.4694469519536142e-18`、 + `max_rel=2.7555746842446215e-16`。 +2. 已接入 `frame` 的流式 sink。`PSF_BACKEND=hip` 使用单一 GPU double HDR framebuffer； + cache 事件按 16,384 项块提交，仍使用既有 bilinear phase cache、圆形支持和 wing clip 规则。 + direct fallback 会完成前序 GPU 工作、在 CPU 计算、再上传 HDR 后继续，因而保持提交顺序。 +3. RX 9070 上的 640×360 固定 Minkowski catalog HDR 对照逐样本一致；含 2 个 cache 事件与 + 1 个 direct fallback 的 `tests/data/hip_psf_mixed_catalog.csv` 场同样逐样本一致。GPU failure + 会打印阶段与 HIP 错误并令该帧返回失败，不会切换到 CPU backend。 +4. 待记录端到端基准的完整命令和原始输出：Git hash、CPU threads、GPU/runtime、块大小、上传、 kernel、download、分辨率、catalog、PSF 及所有 fallback 统计。 5. 只在 direct atomic profile 显示热点竞争是主要瓶颈时，比较 tile-local reduction；ATRI 的 整帧收集／重排同样必须以端到端收益证明。 diff --git a/src/frame.c b/src/frame.c index 8fdf37c..201a793 100644 --- a/src/frame.c +++ b/src/frame.c @@ -1,5 +1,9 @@ #include "frame.h" +#ifdef PSF_BACKEND_HIP +#include "hip_psf.h" +#endif + #include "optics.h" #include @@ -127,36 +131,120 @@ typedef struct { double *hdr; int width, height; const PsfKernelCache *cache; +#ifdef PSF_BACKEND_HIP + HipPsfSink *hip; + char hip_message[256]; +#endif + int failed; } PsfEventSink; -static void psf_event_sink_init(PsfEventSink *sink, double *hdr, int width, - int height, const PsfKernelCache *cache) { +#ifdef PSF_BACKEND_HIP +static void psf_event_sink_mark_failed(PsfEventSink *sink, const char *stage) { + sink->failed = 1; + fprintf(stderr, "HIP PSF backend failed during %s: %s\n", stage, + sink->hip_message[0] == '\0' ? "unknown HIP error" : sink->hip_message); +} +#endif + +static int 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) + if (cache != NULL) { sink->events = malloc(PSF_EVENT_SINK_CAPACITY * sizeof *sink->events); + if (sink->events == NULL) { +#ifdef PSF_BACKEND_HIP + snprintf(sink->hip_message, sizeof sink->hip_message, "host event allocation failed"); + psf_event_sink_mark_failed(sink, "initialization"); + return -1; +#else + return 0; +#endif + } +#ifdef PSF_BACKEND_HIP + if (hip_psf_sink_create(&sink->hip, width, height, cache, + PSF_EVENT_SINK_CAPACITY, sink->hip_message, + sizeof sink->hip_message)) { + psf_event_sink_mark_failed(sink, "initialization"); + free(sink->events); + sink->events = NULL; + return -1; + } +#endif + } + return 0; } -static void psf_event_sink_flush(PsfEventSink *sink) { +static int psf_event_sink_flush(PsfEventSink *sink) { + if (sink->failed) + return -1; +#ifdef PSF_BACKEND_HIP + if (sink->hip != NULL) { + if (hip_psf_sink_submit(sink->hip, sink->events, sink->count, + sink->hip_message, sizeof sink->hip_message)) { + psf_event_sink_mark_failed(sink, "event submission"); + return -1; + } + sink->count = 0; + return 0; + } +#endif 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; + return 0; } -static void psf_event_sink_destroy(PsfEventSink *sink) { - psf_event_sink_flush(sink); +static int psf_event_sink_finish_for_cpu(PsfEventSink *sink) { + if (psf_event_sink_flush(sink)) + return -1; +#ifdef PSF_BACKEND_HIP + if (sink->hip != NULL && hip_psf_sink_finish(sink->hip, sink->hdr, + sink->hip_message, + sizeof sink->hip_message)) { + psf_event_sink_mark_failed(sink, "HDR download"); + return -1; + } +#endif + return 0; +} + +static int psf_event_sink_resume_gpu(PsfEventSink *sink) { +#ifdef PSF_BACKEND_HIP + if (sink->hip != NULL && hip_psf_sink_load_hdr(sink->hip, sink->hdr, + sink->hip_message, + sizeof sink->hip_message)) { + psf_event_sink_mark_failed(sink, "HDR upload"); + return -1; + } +#else + (void)sink; +#endif + return 0; +} + +static int psf_event_sink_destroy(PsfEventSink *sink) { + int result = psf_event_sink_finish_for_cpu(sink); +#ifdef PSF_BACKEND_HIP + hip_psf_sink_destroy(sink->hip); +#endif free(sink->events); + return result; } #if FRAME_PSF_EVENT_SINK static void psf_event_sink_emit(PsfEventSink *sink, const PsfCachedEvent *event) { + if (sink->failed) + return; 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); + (void)psf_event_sink_flush(sink); + if (sink->failed) + return; sink->events[sink->count++] = *event; } #endif @@ -910,6 +998,7 @@ typedef struct { size_t direct_fallbacks; size_t cached_wing_clipped; size_t discarded_below_min_y; + int failed; #ifdef GR_DEBUG double max_raw_magnification; size_t magnification_clamped_triangles; @@ -972,13 +1061,14 @@ static int splat_catalog_tile(const Star *stars, size_t count, &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); + /* A direct fallback forms an explicit ordered CPU/GPU boundary. */ + if (psf_event_sink_finish_for_cpu(context->event_sink)) + return -1; 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); + if (psf_event_sink_resume_gpu(context->event_sink)) + return -1; } else if (direct_fallback != 3) { #if FRAME_PSF_EVENT_SINK psf_event_sink_emit(context->event_sink, &event); @@ -987,6 +1077,8 @@ static int splat_catalog_tile(const Star *stars, size_t count, &event, context->psf_cache); #endif } + if (context->event_sink->failed) + return -1; context->direct_fallbacks += direct_fallback == 1; context->cached_wing_clipped += direct_fallback == 2; context->discarded_below_min_y += direct_fallback == 3; @@ -1018,7 +1110,8 @@ static CatalogSplatStats splat_catalog_triangles( PsfEventSink owned_sink; const int owns_sink = event_sink == NULL; if (owns_sink) { - psf_event_sink_init(&owned_sink, hdr, width, height, psf_cache); + if (psf_event_sink_init(&owned_sink, hdr, width, height, psf_cache)) + return (CatalogSplatStats){.failed = 1}; event_sink = &owned_sink; } for (size_t t = first_triangle; t < last_triangle; ++t) { @@ -1056,15 +1149,21 @@ static CatalogSplatStats splat_catalog_triangles( .psf_relative_tail = psf_relative_tail, .psf_min_y = psf_min_y, .event_sink = event_sink}; - if (catalog_visit_source_triangle(catalog, direction, 0, splat_catalog_tile, - &context) == 0) + const int visit_result = catalog_visit_source_triangle( + catalog, direction, 0, splat_catalog_tile, &context); + if (visit_result == 0) stats.images += context.images; + else { + stats.failed = event_sink->failed; + if (stats.failed) + break; + } stats.direct_fallbacks += context.direct_fallbacks; 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); + if (owns_sink && psf_event_sink_destroy(&owned_sink)) + stats.failed = 1; return stats; } @@ -1126,6 +1225,23 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, progress->callback(progress->context, FRAME_SPLAT_PROGRESS_BEGIN, 0, mesh->triangle_count); +#ifdef PSF_BACKEND_HIP + /* A single GPU HDR framebuffer owns all cached events. Keep catalog/lens + * work serial for this first direct-atomic integration; the CPU parallel + * private-HDR path remains the PSF_BACKEND=cpu implementation. */ + const CatalogSplatStats hip_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, NULL); + copy_psf_splat_stats(psf_stats, hip_stats); + if (hip_stats.failed) + return SIZE_MAX; + if (progress != NULL && progress->callback != NULL) + progress->callback(progress->context, FRAME_SPLAT_PROGRESS_END, + mesh->triangle_count, mesh->triangle_count); + return hip_stats.images; +#endif + const size_t pixel_count = (size_t)width * height * 3; if (pixel_count > SIZE_MAX / sizeof(double)) { diff --git a/src/hip_psf.h b/src/hip_psf.h new file mode 100644 index 0000000..59b2da5 --- /dev/null +++ b/src/hip_psf.h @@ -0,0 +1,44 @@ +#ifndef HIP_PSF_H +#define HIP_PSF_H + +#include "optics.h" + +#include + +#ifdef __cplusplus +extern "C" { +#endif + +/* Opaque HIP state for one HDR framebuffer. It owns reusable device buffers + * for PsfCachedEvent chunks, immutable cache weights, and double RGB HDR. */ +typedef struct HipPsfSink HipPsfSink; + +int hip_psf_available(char *message, size_t message_size); + +/* Creates a zeroed GPU HDR framebuffer. event_capacity is the reusable upload + * chunk capacity, not a full-frame event limit. cache remains caller-owned. */ +int hip_psf_sink_create(HipPsfSink **sink, int width, int height, + const PsfKernelCache *cache, size_t event_capacity, + char *message, size_t message_size); + +/* Completes preceding GPU work, then replaces device HDR with hdr. This is + * the ordered boundary used before resuming GPU cache splats after a CPU + * direct fallback. */ +int hip_psf_sink_load_hdr(HipPsfSink *sink, const double *hdr, + char *message, size_t message_size); + +/* Queues a cached-event chunk. Direct fallbacks and min-Y discards stay with + * the caller; event_count must not exceed the creation capacity. */ +int hip_psf_sink_submit(HipPsfSink *sink, const PsfCachedEvent *events, + size_t event_count, char *message, size_t message_size); + +/* Completes all work and overwrites hdr with double linear RGB device HDR. */ +int hip_psf_sink_finish(HipPsfSink *sink, double *hdr, char *message, + size_t message_size); +void hip_psf_sink_destroy(HipPsfSink *sink); + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/src/hip_psf.hip b/src/hip_psf.hip new file mode 100644 index 0000000..c0bf5e2 --- /dev/null +++ b/src/hip_psf.hip @@ -0,0 +1,173 @@ +#include "hip_psf.h" + +#include + +#include +#include +#include + +struct HipPsfSink { + PsfCachedEvent *events = nullptr; + float *weights = nullptr; + double *hdr = nullptr; + hipStream_t stream = nullptr; + size_t event_capacity = 0, hdr_values = 0; + int width = 0, height = 0, phase_resolution = 0, radius_pixels = 0; + double max_radius_pixels = 0.0; +}; + +static int report(hipError_t status, char *message, size_t message_size) { + if (status == hipSuccess) return 0; + if (message && message_size) std::snprintf(message, message_size, "%s", hipGetErrorString(status)); + return -1; +} +static void ok(char *message, size_t message_size) { + if (message && message_size) std::snprintf(message, message_size, "ok"); +} + +__device__ static size_t weight_index(int phase_resolution, int radius_pixels, + int phase_x, int phase_y, int offset_x, int offset_y) { + const size_t nodes = (size_t)phase_resolution + 1; + const size_t side = (size_t)radius_pixels * 2 + 1; + return (((size_t)phase_y * nodes + phase_x) * side + (size_t)(offset_y + radius_pixels)) * side + + (size_t)(offset_x + radius_pixels); +} + +__device__ static int row_range(double radius, double fx, double fy, int support, + int offset_y, int *first, int *last) { + const double dy = offset_y + 0.5 - fy; + const double remaining = radius * radius - dy * dy; + if (remaining < 0.0) return 0; + const double half_span = sqrt(remaining); + const int low = (int)ceil(fx - 0.5 - half_span); + const int high = (int)floor(fx - 0.5 + half_span); + *first = low > -support ? low : -support; + *last = high < support ? high : support; + return *first <= *last; +} + +__global__ static void splat_kernel(const PsfCachedEvent *events, size_t count, + int width, int height, const float *weights, + int phase_resolution, int radius_pixels, + double max_radius_pixels, double *hdr) { + const size_t index = (size_t)blockIdx.x * blockDim.x + threadIdx.x; + if (index >= count) return; + const PsfCachedEvent event = events[index]; + const int base_x = (int)floor(event.x), base_y = (int)floor(event.y); + const double fx = event.x - base_x, fy = event.y - base_y; + const double phase_x = fx * phase_resolution, phase_y = fy * 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(event.support_radius), max_radius_pixels); + for (int offset_y = -support; offset_y <= support; ++offset_y) { + const int py = base_y + offset_y; + if (py < 0 || py >= height) continue; + int first, last; + if (!row_range(event.support_radius, fx, fy, support, offset_y, &first, &last)) continue; + if (first < -base_x) first = -base_x; + if (last >= width - base_x) last = width - base_x - 1; + for (int offset_x = first; offset_x <= last; ++offset_x) { + const int px = base_x + offset_x; + const double w00 = weights[weight_index(phase_resolution, radius_pixels, x0, y0, offset_x, offset_y)]; + const double w10 = weights[weight_index(phase_resolution, radius_pixels, x1, y0, offset_x, offset_y)]; + const double w01 = weights[weight_index(phase_resolution, radius_pixels, x0, y1, offset_x, offset_y)]; + const double w11 = weights[weight_index(phase_resolution, radius_pixels, x1, y1, offset_x, offset_y)]; + const double weight = (1.0 - ty) * ((1.0 - tx) * w00 + tx * w10) + + ty * ((1.0 - tx) * w01 + tx * w11); + double *pixel = &hdr[3 * ((size_t)py * width + px)]; + atomicAdd(&pixel[0], event.color.r * event.flux * weight); + atomicAdd(&pixel[1], event.color.g * event.flux * weight); + atomicAdd(&pixel[2], event.color.b * event.flux * weight); + } + } +} + +extern "C" int hip_psf_available(char *message, size_t message_size) { + int count = 0; + if (report(hipGetDeviceCount(&count), message, message_size)) return -1; + if (count < 1) { + if (message && message_size) std::snprintf(message, message_size, "no HIP GPU agent found"); + return -1; + } + if (message && message_size) std::snprintf(message, message_size, "HIP device 0 of %d available", count); + return 0; +} + +extern "C" int hip_psf_sink_create(HipPsfSink **out, int width, int height, + const PsfKernelCache *cache, size_t event_capacity, + char *message, size_t message_size) { + if (!out || width <= 0 || height <= 0 || !cache || !cache->ready || !cache->weights || + cache->phase_resolution <= 0 || cache->radius_pixels < 0 || event_capacity == 0) { + if (message && message_size) std::snprintf(message, message_size, "invalid HIP PSF sink arguments"); + return -1; + } + const size_t nodes = (size_t)cache->phase_resolution + 1; + const size_t side = (size_t)cache->radius_pixels * 2 + 1; + if (nodes > std::numeric_limits::max() / nodes || nodes * nodes > std::numeric_limits::max() / side || + nodes * nodes * side > std::numeric_limits::max() / side || (size_t)width > std::numeric_limits::max() / (size_t)height || + (size_t)width * (size_t)height > std::numeric_limits::max() / 3) { + if (message && message_size) std::snprintf(message, message_size, "HIP PSF sink size overflow"); + return -1; + } + *out = nullptr; + HipPsfSink *sink = new HipPsfSink; + sink->event_capacity = event_capacity; + sink->hdr_values = (size_t)width * (size_t)height * 3; + sink->width = width; sink->height = height; sink->phase_resolution = cache->phase_resolution; + sink->radius_pixels = cache->radius_pixels; sink->max_radius_pixels = cache->max_radius_pixels; + const size_t weight_count = nodes * nodes * side * side; + if (report(hipStreamCreate(&sink->stream), message, message_size) || + report(hipMalloc(&sink->events, event_capacity * sizeof *sink->events), message, message_size) || + report(hipMalloc(&sink->weights, weight_count * sizeof *sink->weights), message, message_size) || + report(hipMalloc(&sink->hdr, sink->hdr_values * sizeof *sink->hdr), message, message_size) || + report(hipMemcpyAsync(sink->weights, cache->weights, weight_count * sizeof *sink->weights, hipMemcpyHostToDevice, sink->stream), message, message_size) || + report(hipMemsetAsync(sink->hdr, 0, sink->hdr_values * sizeof *sink->hdr, sink->stream), message, message_size)) { + hip_psf_sink_destroy(sink); return -1; + } + *out = sink; ok(message, message_size); return 0; +} + +extern "C" int hip_psf_sink_submit(HipPsfSink *sink, const PsfCachedEvent *events, + size_t event_count, char *message, size_t message_size) { + if (!sink || (event_count && !events) || event_count > sink->event_capacity) { + if (message && message_size) std::snprintf(message, message_size, "invalid HIP PSF event chunk"); + return -1; + } + if (!event_count) { ok(message, message_size); return 0; } + const unsigned int threads = 128; + const size_t blocks = (event_count + threads - 1) / threads; + if (blocks > std::numeric_limits::max() || + report(hipMemcpyAsync(sink->events, events, event_count * sizeof *events, hipMemcpyHostToDevice, sink->stream), message, message_size)) + return -1; + hipLaunchKernelGGL(splat_kernel, dim3((unsigned int)blocks), dim3(threads), 0, sink->stream, + sink->events, event_count, sink->width, sink->height, sink->weights, + sink->phase_resolution, sink->radius_pixels, sink->max_radius_pixels, sink->hdr); + if (report(hipGetLastError(), message, message_size)) return -1; + ok(message, message_size); return 0; +} + +extern "C" int hip_psf_sink_load_hdr(HipPsfSink *sink, const double *hdr, + char *message, size_t message_size) { + if (!sink || !hdr) { + if (message && message_size) std::snprintf(message, message_size, "invalid HIP PSF HDR upload"); + return -1; + } + if (report(hipMemcpyAsync(sink->hdr, hdr, sink->hdr_values * sizeof *hdr, + hipMemcpyHostToDevice, sink->stream), message, message_size)) + return -1; + ok(message, message_size); return 0; +} + +extern "C" int hip_psf_sink_finish(HipPsfSink *sink, double *hdr, char *message, size_t message_size) { + if (!sink || !hdr) { if (message && message_size) std::snprintf(message, message_size, "invalid HIP PSF HDR download"); return -1; } + if (report(hipMemcpyAsync(hdr, sink->hdr, sink->hdr_values * sizeof *hdr, hipMemcpyDeviceToHost, sink->stream), message, message_size) || + report(hipStreamSynchronize(sink->stream), message, message_size)) return -1; + ok(message, message_size); return 0; +} + +extern "C" void hip_psf_sink_destroy(HipPsfSink *sink) { + if (!sink) return; + (void)hipFree(sink->events); (void)hipFree(sink->weights); (void)hipFree(sink->hdr); + (void)hipStreamDestroy(sink->stream); delete sink; +} diff --git a/src/main.c b/src/main.c index 3f84332..aa35b57 100644 --- a/src/main.c +++ b/src/main.c @@ -572,6 +572,11 @@ static int render_observer_frame(const Settings *s, StarCatalog *catalog, &(FrameSplatProgress){report_splat_progress, s->verbose ? report_splat_worker_progress : NULL, &progress}); + if (images == SIZE_MAX) { + frame_lens_mesh_destroy(&mesh); + free(hdr); + return -1; + } if (s->draw_mesh) frame_draw_mesh(&mesh, hdr, s->width, s->height, 0.5, 0.5); #ifdef ENABLE_HDR_OUTPUT @@ -771,6 +776,10 @@ static int render_movie(const Settings *s, StarCatalog *catalog, &(FrameSplatProgress){report_splat_progress, s->verbose ? report_splat_worker_progress : NULL, &progress}); + if (images == SIZE_MAX) { + free(hdr); + goto done; + } if (s->draw_mesh) frame_draw_mesh(&movie.frames[i].mesh, hdr, s->width, s->height, 0.5, 0.5); const int write_result = write_tonemapped_image(output_path, hdr, s->width, s->height); @@ -835,6 +844,11 @@ static int render_lens_map(const Settings *s, StarCatalog *catalog) { &(FrameSplatProgress){report_splat_progress, s->verbose ? report_splat_worker_progress : NULL, &progress}); + if (images == SIZE_MAX) { + free(hdr); + result = -1; + break; + } if (s->draw_mesh) frame_draw_mesh(&map.frames[i].mesh, hdr, map.width, map.height, 0.5, 0.5); #ifdef ENABLE_HDR_OUTPUT diff --git a/src/optics.h b/src/optics.h index fb62138..2a6b2f7 100644 --- a/src/optics.h +++ b/src/optics.h @@ -41,6 +41,10 @@ typedef struct { double flux, support_radius; } PsfCachedEvent; +#ifdef __cplusplus +extern "C" { +#endif + /* 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); @@ -91,4 +95,8 @@ int write_hdr_fits(const char *path, const double *hdr, int width, int height, double horizontal_fov_deg); #endif +#ifdef __cplusplus +} +#endif + #endif diff --git a/tests/data/hip_psf_mixed_catalog.csv b/tests/data/hip_psf_mixed_catalog.csv new file mode 100644 index 0000000..0eec936 --- /dev/null +++ b/tests/data/hip_psf_mixed_catalog.csv @@ -0,0 +1,4 @@ +longitude_deg,latitude_deg,temperature_K,amplitude +0,0,3000,1 +2,0,12000,0.00141095580387 +4,0,12000,0.00141095580387 diff --git a/tests/test_hip_psf.hip b/tests/test_hip_psf.hip new file mode 100644 index 0000000..e6a8880 --- /dev/null +++ b/tests/test_hip_psf.hip @@ -0,0 +1,60 @@ +#include "hip_psf.h" + +#include +#include +#include + +int main() { + constexpr int width = 64, height = 48; + constexpr size_t values = (size_t)width * height * 3; + const PointSpreadFunction psf = {2.7, 4.5}; + PsfKernelCache cache = {}; + if (psf_kernel_cache_init(&cache, &psf, 1e-8)) { + std::fputs("could not build PSF cache\n", stderr); + return 1; + } + PsfCachedEvent events[3] = {}; + const LinearRgb colors[3] = {{1.0, 0.2, 0.6}, {0.1, 0.8, 0.3}, {0.7, 0.4, 1.0}}; + const double positions[3][2] = {{23.25, 19.75}, {24.50, 20.125}, {23.875, 21.25}}; + for (size_t i = 0; i < 3; ++i) { + if (psf_prepare_cached_event(&events[i], positions[i][0], positions[i][1], + colors[i], 0.1 + 0.03 * i, &psf, &cache, + 1.0, 1e-8, 0.0) != 0) { + std::fputs("test event unexpectedly missed the cache\n", stderr); + psf_kernel_cache_destroy(&cache); + return 1; + } + } + double *cpu_hdr = (double *)std::calloc(values, sizeof *cpu_hdr); + double *gpu_hdr = (double *)std::calloc(values, sizeof *gpu_hdr); + if (!cpu_hdr || !gpu_hdr) { + std::fputs("HDR allocation failed\n", stderr); + std::free(cpu_hdr); std::free(gpu_hdr); psf_kernel_cache_destroy(&cache); + return 1; + } + for (const PsfCachedEvent &event : events) + splat_prepared_cached_event(cpu_hdr, width, height, &event, &cache); + + char message[256] = {}; + HipPsfSink *sink = nullptr; + if (hip_psf_available(message, sizeof message) || + hip_psf_sink_create(&sink, width, height, &cache, 3, message, sizeof message) || + hip_psf_sink_submit(sink, events, 3, message, sizeof message) || + hip_psf_sink_finish(sink, gpu_hdr, message, sizeof message)) { + std::fprintf(stderr, "HIP PSF test failed: %s\n", message); + hip_psf_sink_destroy(sink); + std::free(cpu_hdr); std::free(gpu_hdr); psf_kernel_cache_destroy(&cache); + return 1; + } + double max_abs = 0.0, max_rel = 0.0; + for (size_t i = 0; i < values; ++i) { + const double absolute = std::fabs(cpu_hdr[i] - gpu_hdr[i]); + max_abs = std::fmax(max_abs, absolute); + if (std::fabs(cpu_hdr[i]) > 1e-30) + max_rel = std::fmax(max_rel, absolute / std::fabs(cpu_hdr[i])); + } + std::printf("HIP PSF cache comparison: max_abs=%.17g max_rel=%.17g\n", max_abs, max_rel); + hip_psf_sink_destroy(sink); + std::free(cpu_hdr); std::free(gpu_hdr); psf_kernel_cache_destroy(&cache); + return max_abs <= 1e-12 && max_rel <= 1e-12 ? 0 : 1; +}