Feat: add initial HIP PSF backend

src/frame.c

@@ -1,5 +1,9 @@
 #include "frame.h"
 
+#ifdef PSF_BACKEND_HIP
+#include "hip_psf.h"
+#endif
+
 #include "optics.h"
 
 #include <limits.h>
@@ -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))
     {

src/hip_psf.h

@@ -0,0 +1,44 @@
+#ifndef HIP_PSF_H
+#define HIP_PSF_H
+
+#include "optics.h"
+
+#include <stddef.h>
+
+#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

src/hip_psf.hip

@@ -0,0 +1,173 @@
+#include "hip_psf.h"
+
+#include <hip/hip_runtime.h>
+
+#include <cmath>
+#include <cstdio>
+#include <limits>
+
+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<size_t>::max() / nodes || nodes * nodes > std::numeric_limits<size_t>::max() / side ||
+      nodes * nodes * side > std::numeric_limits<size_t>::max() / side || (size_t)width > std::numeric_limits<size_t>::max() / (size_t)height ||
+      (size_t)width * (size_t)height > std::numeric_limits<size_t>::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<unsigned int>::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;
+}

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

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