Feat: add initial HIP PSF backend

Doc: update GPU PSF acceleration status
DOC: preserve complete benchmark output
2026-09-05 04:23:57 -04:00 · 2026-09-05 03:47:42 -04:00 · 2026-09-05 03:22:15 -04:00 · 2026-09-05 03:20:40 -04:00 · 2026-09-05 03:06:46 -04:00 · 2026-09-04 20:28:54 -04:00
23 changed files with 1478 additions and 61 deletions
@@ -74,6 +74,8 @@ catalog 内部数据保留 `(direction, temperature, amplitude)`，而非 RGB。

 新增物理、插值或优化时，优先添加能与前一阶段比较的收敛测试或 regression test。尚未由 prototype/convergence test 决定的参数（如时间输出 cadence、时间插值阶数、slab 大小、refinement 阈值、PSF、ODE stepper）不要伪装成既定事实。

+性能 benchmark 记录必须保留完整、可复制的命令及原始终端输出，不能只记录汇总耗时或吞吐量；输出中的 build/cache、输入加载、工作线程、处理数量与 fallback 等统计是后续正确归因性能变化的证据。
+
 ## 修改原则

 - 改动应保持 backend、observer、geodesic、movie mesh 与 optics 的职责分离。
@@ -6,6 +6,10 @@ ENABLE_PNG ?= 1
 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
@@ -36,8 +40,23 @@ 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
@@ -50,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
+HDR_BUILD_TAG := hdr_sink$(PSF_EVENT_SINK)_$(PSF_BACKEND)
 else ifeq ($(ENABLE_HDR),0)
-HDR_BUILD_TAG := standard
+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
@@ -82,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:

@@ -93,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)
@@ -123,3 +166,4 @@ clean:
 	rm -rf build

 -include $(RENDER_DEPS)
+include mk/reference_images.mk
@@ -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.
@@ -54,6 +74,55 @@ synthetic test catalog uses global default exposure `1e-3`; the accelerated
 benchmark explicitly uses `1e-5` because its physical Doppler blue shift
 otherwise clips the later frames.

+## Reuse a completed lens map
+
+Ray tracing and adaptive mesh refinement are independent of catalog lookup,
+PSF evaluation, exposure, and tone mapping.  `--lens-map-output FILE` writes
+the finalized local inverse-lens mesh to one versioned `.grlens` file after all
+ray-trace/refinement generations finish; the same invocation still renders its
+ordinary image.  The file stores every final vertex's image position, camera
+direction, infinity direction, frequency shift, terminal status, and the
+triangle topology.  It also stores the frame dimensions, horizontal FOV, and,
+for a movie, all frame IDs and times.
+
+For example, trace an analytic Schwarzschild frame once and retain the map:
+
+```sh
+make SPACETIME=schwarzschild backend
+mkdir -p output/imgs output/maps
+./build/Release/schwarzschild_sky --catalog assets/sky_grid_5deg.csv \
+  --width 640 --height 360 --coarse-cell-pixels 8 --fov-deg 60 \
+  --refine-max-level 2 --lens-map-output output/maps/schwarzschild.grlens \
+  --output output/imgs/schwarzschild_trace.png
+```
+
+Later, render that saved map with a different catalog, PSF, or exposure without
+constructing an observer, spacetime source, or geodesic rays:
+
+```sh
+./build/Release/schwarzschild_sky --lens-map-input output/maps/schwarzschild.grlens \
+  --catalog assets/sky_grid_5deg.csv --psf-fwhm-pixels 6 --psf-moffat-beta 3 \
+  --exposure 0.002 --output output/imgs/schwarzschild_restyled.png
+```
+
+`--lens-map-input` and `--lens-map-output` are mutually exclusive.  An imported
+map always retains its original pixel width, height, and horizontal FOV; an
+explicit conflicting `--width`, `--height`, or `--fov-deg` is rejected.  This is
+intentional: the stored inverse map and PSF coordinates are in the original
+pixel geometry.  The reader validates the format version, finite values, unit
+directions, triangle indices, and a per-frame CRC before rendering.
+
+A movie export writes every final frame mesh into the same `.grlens` file.
+Export it during the usual movie render, then import it with `--frames-dir` and
+`--frames-prefix`; importing a multi-frame map does not require
+`--observer-track`:
+
+```sh
+./build/Release/minkowski_sky --lens-map-input output/maps/movie.grlens \
+  --catalog assets/sky_grid_5deg.csv --frames-dir output/imgs \
+  --frames-prefix restyled
+```
+
 PNG is the default output and the default build links `libpng`:

 ```sh
@@ -0,0 +1,90 @@
+# 2MASS 全天 1080p CPU PSF event-sink 变更前后记录（2026-09-05）
+
+## 目的
+
+记录引入 CPU `PsfEventSink` 前后的完整渲染耗时。cache-eligible 星像在改动后由每个
+OpenMP render worker 的固定大小块暂存、批量消费；cache/direct/min-Y 判定、HDR 值和
+输出语义不改变。本记录不是 GPU 基准；它以完全没有 event 路径的上一版本为基线，衡量
+event-sink 实现进入生产路径的端到端增量。
+
+## 控制条件
+
+- Host: i7-12700K，16 个 online logical CPUs；未显式设置 `OMP_NUM_THREADS`。
+- Build: Release, C11, `-march=native -O2 -pipe`, OpenMP, PNG；Minkowski。
+- Input: `assets/2mass/processed/all_sky`；9,328 tiles / 18,148,775 stars。
+- Camera: RA 0 deg, Dec -90 deg; 1920 x 1080 px; horizontal FOV 60 deg;
+  coarse cell 16 px; exposure `1e12`。
+- PSF: default FWHM 2.7 px, beta 4.5, relative tail `1e-8`, min-Y disabled,
+  cache/direct threshold 1. 两次均无 direct fallback。
+- 两次均产生 16,554,566 个星像；catalog prefetch 统计一致。
+
+两次运行的完整命令与原始输出如下；保留这些输出以便后续性能变化可结合 cache build、
+catalog prefetch、星像数、fallback 统计与实际 CPU 利用率判断，而非只比较汇总数字。
+
+### stash 后基线（`4953775`）
+
+```console
+time ./build/Release/minkowski_sky \
+  --all-sky-catalog assets/2mass/processed/all_sky \
+  --look-ra-deg 0 --look-dec-deg -90 \
+  --width 1920 --height 1080 \
+  --coarse-cell-pixels 16 --fov-deg 60 \
+  --exposure 1e12 \
+  --output output/imgs/2mass_1080p_60deg_south_pole_1e12_psf_cache_event.png
+PSF cache ready: 64x64 phases, radius 47 px, relative tail 1e-08, tail abs 1e-06, boundary 1e-07, build 0.681 s
+Catalog prefetch: 9328 requested, 9328 newly loaded (18148775 stars), 0 unavailable in 3.942 s; 4 loader workers
+Rendered 16554566 images from 18148775 catalog stars to output/imgs/2mass_1080p_60deg_south_pole_1e12_psf_cache_event.png (ok)
+PSF splats: cached 16554566, cached wing-clipped 0, direct fallbacks 0, discarded below min-Y 0
+./build/Release/minkowski_sky --all-sky-catalog assets/2mass/processed/all_sk  703.03s user 5.78s system 1424% cpu 49.753 total
+```
+
+### stash pop 恢复后的 event-sink 改动（后来提交为 `2a03cbf`）
+
+```console
+time ./build/Release/minkowski_sky \
+  --all-sky-catalog assets/2mass/processed/all_sky \
+  --look-ra-deg 0 --look-dec-deg -90 \
+  --width 1920 --height 1080 \
+  --coarse-cell-pixels 16 --fov-deg 60 \
+  --exposure 1e12 \
+  --output output/imgs/2mass_1080p_60deg_south_pole_1e12_psf_cache_event.png
+PSF cache ready: 64x64 phases, radius 47 px, relative tail 1e-08, tail abs 1e-06, boundary 1e-07, build 0.672 s
+Catalog prefetch: 9328 requested, 9328 newly loaded (18148775 stars), 0 unavailable in 3.926 s; 4 loader workers
+Rendered 16554566 images from 18148775 catalog stars to output/imgs/2mass_1080p_60deg_south_pole_1e12_psf_cache_event.png (ok)
+PSF splats: cached 16554566, cached wing-clipped 0, direct fallbacks 0, discarded below min-Y 0
+./build/Release/minkowski_sky --all-sky-catalog assets/2mass/processed/all_sk  710.02s user 5.90s system 1419% cpu 50.424 total
+```
+
+状态的切换由 Git 完成：先在包含 event-sink 改动的工作树中运行，随后 `git stash`
+暂存该改动、重新编译并运行基线，最后 `git stash pop` 恢复。因而这项测量比较的是
+`4953775` 与后来提交为 `2a03cbf` 的完整源码差异。前者没有 event 路径；后者的差异
+正是 event 准备／消费边界、每 worker 缓冲及其生命周期。因此这不是“开关的微基准”，而是
+event 实现的端到端前后比较。
+
+结果表中的 cache build 与 prefetch 是程序输出；两次均在对应状态下重新编译后运行。
+
+## 结果
+
+| state | revision | cache build | prefetch | user | system | wall | CPU | user/image |
+|---|---|---:|---:|---:|---:|---:|---:|---:|
+| stash 后的基线 | `4953775` | 0.681 s | 3.942 s | 703.03 s | 5.78 s | 49.753 s | 1424% | 42.47 us |
+| stash pop 恢复后的 event-sink 改动 | 后来提交为 `2a03cbf` | 0.672 s | 3.926 s | 710.02 s | 5.90 s | 50.424 s | 1419% | 42.89 us |
+
+恢复改动后的观测差为 +6.99 user s、+0.671 wall s、+0.42 us/star image，即
+**+0.99% user CPU work**。两次的平均并行度近似相同（`user / wall`：14.13 与
+14.08），所以 wall-time 差异并非由明显不同的实际 worker 使用量造成。
+
+因此在这套真实渲染配置上，event-sink 实现的测得开销为约 0.99% user CPU work。该数值与
+固定容量 event 拷贝、额外分支及最终 flush 所预期的小幅 CPU 开销相容。每个状态目前仅有
+一次运行；若要估计波动范围，应在两个 Git 状态下各重复运行多次，而不是把基线改造成一个
+原先不存在的运行时路径。
+
+## Historical comparison boundary
+
+`benchmarks/2mass_all_sky_psf_lookup_2026-08-28.md` recorded 645.24 user s /
+74.48 wall s for lookup, but its analytic render path was restricted to at most
+10 private-HDR workers. The current runs use about 14 active workers on
+average. Therefore the current 49.753--50.424 s wall time is not evidence of
+an event-sink speedup, and the older 645.24 user s cannot be attributed to the
+event sink. 相反，本页的 `4953775..2a03cbf` 前后比较正是 event-sink 实现的性能
+增量记录；它与旧基准的不同仅在于旧基准还混入了已移除的 worker 限制。
@@ -0,0 +1,70 @@
+# PSF event-sink image references
+
+Run date: 2026-09-04
+
+Purpose: fixed pixel-comparison references before replacing the sink consumer
+with HIP.  The renderer used the default CPU `PSF_EVENT_SINK=1` path.  There
+was no adaptive refinement.  The catalog is the repository test catalog, not
+an all-sky tile snapshot.  This is not a performance benchmark.
+
+Environment:
+
+- Git `HEAD`: `ab56e23` (worktree contains the uncommitted event-sink change)
+- Build: `BUILD_TYPE=Release`, `ENABLE_HDR=1`, `ENABLE_PNG=1`, OpenMP
+- Host worker setting: `OMP_NUM_THREADS=16`
+- Images: 640 x 360 px; Moffat FWHM 2.7 px, beta 4.5
+- Observer: radius 30, inward speed 0; ICRS look direction RA 1 deg, Dec 1 deg
+
+The explicit `1e300` magnification cap is the practical command-line spelling
+of the normal unlimited default.  `--draw-mesh` and `--psf-direct` are absent,
+therefore disabled.  All refinement thresholds are included even though level
+zero prevents refinement, so later default changes cannot affect this run.
+
+## Minkowski
+
+Build:
+
+```bash
+make -B SPACETIME=minkowski ENABLE_HDR=1 PSF_EVENT_SINK=1 backend
+```
+
+Render:
+
+```bash
+OMP_NUM_THREADS=16 build/Release/minkowski_sky --catalog assets/sky_grid_5deg.csv --output tests/data/psf_event_sink_reference/minkowski_ra1_dec1_640x360.png --hdr-output --width 640 --height 360 --fov-deg 30 --look-ra-deg 1 --look-dec-deg 1 --exposure 0.1 --observer-radius 30 --observer-inward-speed 0 --psf-fwhm-pixels 2.7 --psf-moffat-beta 4.5 --max-magnification 1e300 --max-cache-psf-flux 1 --psf-relative-tail 1e-8 --psf-min-y 0 --coarse-cell-pixels 16 --refine-max-level 0 --refine-angle-abs-deg 0.001 --refine-angle-rel 0.1 --refine-jacobian-min 1e-3 --refine-min-edge-pixels 0.5 --refine-min-area-pixels2 0.25 --catalog-load-workers 4
+```
+
+Result: 35 images from 5,654 catalog stars; 35 cached, 0 wing-clipped, 0 direct
+fallback, 0 min-Y discarded.
+
+Versioned reference outputs:
+
+- `tests/data/psf_event_sink_reference/minkowski_ra1_dec1_640x360.png`
+- `tests/data/psf_event_sink_reference/minkowski_ra1_dec1_640x360_HDR.fits`
+
+## Analytic Schwarzschild
+
+This reference uses a 60 deg horizontal field of view.  At radius 30, the
+former 30 deg framing mostly contained the shadow and high-order images inside
+the Einstein ring, without enough field to test the first-image star field.
+
+Build:
+
+```bash
+make -B SPACETIME=schwarzschild ENABLE_HDR=1 PSF_EVENT_SINK=1 backend
+```
+
+Render:
+
+```bash
+OMP_NUM_THREADS=16 build/Release/schwarzschild_sky --catalog assets/sky_grid_5deg.csv --output tests/data/psf_event_sink_reference/schwarzschild_ra1_dec1_fov60_640x360.png --hdr-output --width 640 --height 360 --fov-deg 60 --look-ra-deg 1 --look-dec-deg 1 --exposure 0.1 --observer-radius 30 --observer-inward-speed 0 --psf-fwhm-pixels 2.7 --psf-moffat-beta 4.5 --max-magnification 1e300 --max-cache-psf-flux 1 --psf-relative-tail 1e-8 --psf-min-y 0 --coarse-cell-pixels 16 --refine-max-level 0 --refine-angle-abs-deg 0.001 --refine-angle-rel 0.1 --refine-jacobian-min 1e-3 --refine-min-edge-pixels 0.5 --refine-min-area-pixels2 0.25 --catalog-load-workers 4
+```
+
+Result: 5,206 images from 5,654 catalog stars; 5,206 cached, 0 wing-clipped,
+0 direct fallback, 0 min-Y discarded.  This is a visual regression reference,
+not a performance result.
+
+Versioned reference outputs:
+
+- `tests/data/psf_event_sink_reference/schwarzschild_ra1_dec1_fov60_640x360.png`
+- `tests/data/psf_event_sink_reference/schwarzschild_ra1_dec1_fov60_640x360_HDR.fits`
@@ -51,19 +51,26 @@ CPU HDR 与 catalog，仍有可用空间。其可能收益不是减少需要传
 这是一条 ATRI 专用的空间换时间路线，必须以端到端实测证明收益后才保留。它不应成为 ITX、
 Optiplex 或通用 CPU fallback 的前提；在这些机器上仍使用流式块。

-应在 `optics/psf` 层定义一个按块提交的 sink。概念性事件为：
+当前 CPU 已在 `optics` 层建立 cache-eligible 的事件边界：

 ```c
 typedef struct {
-    float x, y;          /* continuous image position */
-    float r, g, b;       /* already premultiplied by flux */
-    float support_radius;
-} PsfEvent;
+    double x, y;
+    LinearRgb color;
+    double flux, support_radius;
+} PsfCachedEvent;
 ```

-CPU worker 产生少量固定上限的事件块，填满后提交；GPU 与 CPU 以双/三缓冲方式重叠执行。
-块的首版大小从 1--16 MiB 试起，由实测决定。GPU 后端仍必须统计：cache splat、wing clip、
-direct fallback 和 min-Y discarded；这些计数在块完成后汇总，而不是在热循环中同步。
+`psf_prepare_cached_event()` 保留既有的 cache/direct/min-Y 判定，状态 0 或 2 产出该
+事件；direct fallback 在 CPU 立即处理，min-Y 在上传前丢弃。`PsfEventSink` 由每个 OpenMP
+worker 初始化一次、贯穿其线程生命周期，固定容量为 16,384 个事件；填满或线程结束时由现有
+CPU cache evaluator 消费。此实现是 GPU producer 的正确性基线，并已通过固定 HDR reference
+和完整 2MASS 前后性能记录验证。
+
+GPU sink 将消费同一 `PsfCachedEvent` 与现有不可变 `PsfKernelCache::weights`，不得重新以
+float 解析采样替代 cache 语义。首版使用有限大小的上传块并复用 device buffer；块大小从
+1--16 MiB 试起，由实测决定。GPU 后端仍必须统计：cache splat、wing clip、direct fallback
+和 min-Y discarded；这些计数在块完成后汇总，而不是在热循环中同步。

 为降低局部热点，CPU 侧应把来自不同 image-plane triangle / worker 的事件块交错提交。
 三角形内部的星像中心必定落在各自不重叠的像平面三角形内；不同三角形的 PSF wings 仍可
@@ -131,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 和参数。

@@ -141,14 +148,36 @@ GPU 后端必须先在小、固定 catalog snapshot 上同 CPU 的 `--psf-direct

 ## 7. 实施顺序

-1. 只做接口：将现有 CPU immediate splat 包装为默认 sink，输出不得变化。
-2. 实现 CPU 事件块 producer，但仍由 CPU sink 消费；验证事件流没有改变 HDR 和统计。
-3. 在工作站实施 HIP direct-atomic 原型，完成正确性和端到端基准。
-4. 在 ATRI 上比较流式块与整帧收集、空间重排/分桶两种输入组织；记录 host RAM 峰值、上传、
-   GPU kernel 和端到端时间。
-5. 根据 profile 决定是否实现 tile-local reduction；不因纸面 TFLOPS 推断需要它。
-6. 实现 Vulkan 设备枚举、feature probe 和 direct-atomic 后端；feature 不足时 CPU fallback。
-7. 在 ITX、Optiplex、工作站分别跑同一小基准和足够长的稳定性试验，记录真实可用矩阵。
-8. 只有当基准显示明显收益时，才为 WX 4100 或 Intel 核显投入 tile-reduction / 专门调优。
+### 当前状态（2026-09-05）
+
+已完成 CPU event 路径：`PsfCachedEvent`、`psf_prepare_cached_event()` 与 per-worker
+`PsfEventSink` 已进入 `frame_splat_catalog()`。固定 HDR reference 与 unit tests 验证该路径
+不改变输出；在 2MASS 全天 1080p 的前后比较中，完整 event 实现相对无 event 路径版本增加
+0.99% user CPU work。该记录含完整命令和原始输出，见
+`benchmarks/2mass_all_sky_psf_event_sink_2026-09-05.md`。
+
+曾有一个 float、解析中心采样的 HIP atomic spike，用于确认 HIP runtime、上传与 float atomic
+可执行；它没有 cache 权重、双精度 HDR 或 production integration。其 `test-hip` target 和测试
+已移除，不能作为 renderer GPU backend 的验证或运行入口。
+
+HIP direct-atomic 的当前实现状态：
+
+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 的
+   整帧收集／重排同样必须以端到端收益证明。
+6. HIP 路径稳定后再评估 Vulkan 的可移植后端及其他机器的实际设备矩阵。

 每一步独立提交，并保持 CPU reference、GPU runtime 层和 shader 资源的提交边界清晰。
@@ -0,0 +1,16 @@
+REFERENCE_DIR := tests/data/psf_event_sink_reference
+REFERENCE_TMP_DIR := /tmp/gr_psf_event_sink_reference
+FLOATDIFF_SCRIPT := scripts/fits_floatdiff.py
+
+.PHONY: test-reference-images
+
+test: test-reference-images
+
+test-reference-images: $(FLOATDIFF_SCRIPT) $(REFERENCE_DIR)/minkowski_ra1_dec1_640x360_HDR.fits $(REFERENCE_DIR)/schwarzschild_ra1_dec1_fov60_640x360_HDR.fits
+	$(MAKE) SPACETIME=minkowski ENABLE_HDR=1 backend
+	mkdir -p $(REFERENCE_TMP_DIR)
+	OMP_NUM_THREADS=16 $(BUILD_DIR)/minkowski_sky --catalog assets/sky_grid_5deg.csv --output $(REFERENCE_TMP_DIR)/minkowski_ra1_dec1_640x360.png --hdr-output --width 640 --height 360 --fov-deg 30 --look-ra-deg 1 --look-dec-deg 1 --exposure 0.1 --observer-radius 30 --observer-inward-speed 0 --psf-fwhm-pixels 2.7 --psf-moffat-beta 4.5 --max-magnification 1e300 --max-cache-psf-flux 1 --psf-relative-tail 1e-8 --psf-min-y 0 --coarse-cell-pixels 16 --refine-max-level 0 --refine-angle-abs-deg 0.001 --refine-angle-rel 0.1 --refine-jacobian-min 1e-3 --refine-min-edge-pixels 0.5 --refine-min-area-pixels2 0.25 --catalog-load-workers 4
+	python3 $(FLOATDIFF_SCRIPT) $(REFERENCE_DIR)/minkowski_ra1_dec1_640x360_HDR.fits $(REFERENCE_TMP_DIR)/minkowski_ra1_dec1_640x360_HDR.fits
+	$(MAKE) SPACETIME=schwarzschild ENABLE_HDR=1 backend
+	OMP_NUM_THREADS=16 $(BUILD_DIR)/schwarzschild_sky --catalog assets/sky_grid_5deg.csv --output $(REFERENCE_TMP_DIR)/schwarzschild_ra1_dec1_fov60_640x360.png --hdr-output --width 640 --height 360 --fov-deg 60 --look-ra-deg 1 --look-dec-deg 1 --exposure 0.1 --observer-radius 30 --observer-inward-speed 0 --psf-fwhm-pixels 2.7 --psf-moffat-beta 4.5 --max-magnification 1e300 --max-cache-psf-flux 1 --psf-relative-tail 1e-8 --psf-min-y 0 --coarse-cell-pixels 16 --refine-max-level 0 --refine-angle-abs-deg 0.001 --refine-angle-rel 0.1 --refine-jacobian-min 1e-3 --refine-min-edge-pixels 0.5 --refine-min-area-pixels2 0.25 --catalog-load-workers 4
+	python3 $(FLOATDIFF_SCRIPT) $(REFERENCE_DIR)/schwarzschild_ra1_dec1_fov60_640x360_HDR.fits $(REFERENCE_TMP_DIR)/schwarzschild_ra1_dec1_fov60_640x360_HDR.fits
@@ -0,0 +1,78 @@
+#!/usr/bin/env python3
+"""Compare primary-image IEEE float values in two simple FITS files.
+
+The renderer writes a single 32-bit RGB primary HDU.  Keeping this dependency-
+free tool in the repository makes HDR reference comparisons available on the
+minimal Gentoo installation as well as in CI.
+"""
+
+import argparse
+import math
+import struct
+import sys
+
+
+def read_fits(path):
+    with open(path, "rb") as stream:
+        cards = []
+        while True:
+            block = stream.read(2880)
+            if len(block) != 2880:
+                raise ValueError(f"{path}: truncated FITS header")
+            cards.extend(block[i:i + 80].decode("ascii") for i in range(0, 2880, 80))
+            if any(card.startswith("END") for card in cards[-36:]):
+                break
+        values = {}
+        for card in cards:
+            if card[8:10] == "= ":
+                values[card[:8].strip()] = card[10:30].strip()
+        if values.get("BITPIX") != "-32" or values.get("NAXIS") != "3":
+            raise ValueError(f"{path}: expected a 32-bit, three-axis primary image")
+        shape = tuple(int(values[f"NAXIS{axis}"]) for axis in (1, 2, 3))
+        count = math.prod(shape)
+        payload = stream.read(count * 4)
+        if len(payload) != count * 4:
+            raise ValueError(f"{path}: truncated image payload")
+        return shape, payload
+
+
+def main():
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("reference")
+    parser.add_argument("candidate")
+    parser.add_argument("--abs-tolerance", type=float, default=0.0)
+    parser.add_argument("--rel-tolerance", type=float, default=0.0)
+    args = parser.parse_args()
+    if args.abs_tolerance < 0 or args.rel_tolerance < 0:
+        parser.error("tolerances must be non-negative")
+    ref_shape, ref_payload = read_fits(args.reference)
+    candidate_shape, candidate_payload = read_fits(args.candidate)
+    if ref_shape != candidate_shape:
+        print(f"shape mismatch: {ref_shape} != {candidate_shape}", file=sys.stderr)
+        return 1
+    maximum_abs = maximum_rel = 0.0
+    mismatch_count = 0
+    worst_index = 0
+    for index, (reference, candidate) in enumerate(
+            zip(struct.iter_unpack(">f", ref_payload),
+                struct.iter_unpack(">f", candidate_payload))):
+        reference, candidate = reference[0], candidate[0]
+        if not math.isfinite(reference) or not math.isfinite(candidate):
+            equal = reference == candidate
+            absolute = relative = math.inf if not equal else 0.0
+        else:
+            absolute = abs(candidate - reference)
+            relative = absolute / max(abs(reference), 1.0e-30)
+            equal = absolute <= args.abs_tolerance or relative <= args.rel_tolerance
+        if absolute > maximum_abs:
+            maximum_abs, worst_index = absolute, index
+        maximum_rel = max(maximum_rel, relative)
+        mismatch_count += not equal
+    print(f"floatdiff: shape={ref_shape} samples={math.prod(ref_shape)} "
+          f"mismatches={mismatch_count} max_abs={maximum_abs:.9g} "
+          f"max_rel={maximum_rel:.9g} worst_sample={worst_index}")
+    return 1 if mismatch_count else 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())
@@ -1,5 +1,9 @@
 #include "frame.h"

+#ifdef PSF_BACKEND_HIP
+#include "hip_psf.h"
+#endif
+
 #include "optics.h"

 #include <limits.h>
@@ -10,6 +14,10 @@
 #include <stdlib.h>
 #include <string.h>

+#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 +123,132 @@ 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;
+#ifdef PSF_BACKEND_HIP
+  HipPsfSink *hip;
+  char hip_message[256];
+#endif
+  int failed;
+} PsfEventSink;
+
+#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) {
+    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 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 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)
+    (void)psf_event_sink_flush(sink);
+  if (sink->failed)
+    return;
+  sink->events[sink->count++] = *event;
+}
+#endif
+
 typedef struct {
  size_t a, b, triangle;
  unsigned int side;
@@ -852,6 +986,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;
@@ -863,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;
@@ -920,11 +1056,29 @@ 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) {
+      /* 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);
+#else
+      splat_prepared_cached_event(context->hdr, context->width, context->height,
+                                  &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;
@@ -951,8 +1105,15 @@ 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) {
+    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) {
    const LensVertex *vertex[3];
    if (!usable_triangle(mesh, &mesh->triangles[t], vertex))
@@ -986,14 +1147,23 @@ 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};
-    if (catalog_visit_source_triangle(catalog, direction, 0, splat_catalog_tile,
-                                      &context) == 0)
+                                    .psf_min_y = psf_min_y,
+                                    .event_sink = event_sink};
+    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))
+    stats.failed = 1;
  return stats;
 }

@@ -1055,13 +1225,30 @@ 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))
    {
      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 +1265,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 +1276,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 +1293,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 +1313,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 +1327,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 +1344,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 +1356,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 +1374,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)
@@ -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
@@ -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;
+}
@@ -0,0 +1,174 @@
+#include "lens_map.h"
+
+#include <errno.h>
+#include <math.h>
+#include <stdint.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+/* All scalar fields are explicitly little-endian; never serialize C structs
+ * because their padding and size_t width are ABI-dependent. */
+static const unsigned char lens_map_magic[8] = {'G', 'R', 'L', 'E', 'N', 'S', 1, 0};
+enum { LENS_MAP_VERSION = 1, LENS_MAP_ENDIAN = 0x01020304u };
+
+static uint32_t crc32_update(uint32_t crc, const void *data, size_t size) {
+  const unsigned char *bytes = data;
+  for (size_t i = 0; i < size; ++i) {
+    crc ^= bytes[i];
+    for (int bit = 0; bit < 8; ++bit)
+      crc = (crc >> 1) ^ (0xedb88320u & (uint32_t)-(int)(crc & 1));
+  }
+  return crc;
+}
+
+static int write_bytes(FILE *file, const void *data, size_t size, uint32_t *crc) {
+  if (fwrite(data, 1, size, file) != size) return -1;
+  if (crc != NULL) *crc = crc32_update(*crc, data, size);
+  return 0;
+}
+static int read_bytes(FILE *file, void *data, size_t size, uint32_t *crc) {
+  if (fread(data, 1, size, file) != size) return -1;
+  if (crc != NULL) *crc = crc32_update(*crc, data, size);
+  return 0;
+}
+static int write_u32(FILE *f, uint32_t v, uint32_t *c) {
+  unsigned char b[4] = {(unsigned char)v, (unsigned char)(v >> 8),
+                        (unsigned char)(v >> 16), (unsigned char)(v >> 24)};
+  return write_bytes(f, b, sizeof b, c);
+}
+static int read_u32(FILE *f, uint32_t *v, uint32_t *c) {
+  unsigned char b[4]; if (read_bytes(f, b, sizeof b, c)) return -1;
+  *v = (uint32_t)b[0] | ((uint32_t)b[1] << 8) | ((uint32_t)b[2] << 16) |
+       ((uint32_t)b[3] << 24); return 0;
+}
+static int write_u64(FILE *f, uint64_t v, uint32_t *c) {
+  unsigned char b[8]; for (int i = 0; i < 8; ++i) b[i] = (unsigned char)(v >> (8*i));
+  return write_bytes(f, b, sizeof b, c);
+}
+static int read_u64(FILE *f, uint64_t *v, uint32_t *c) {
+  unsigned char b[8]; if (read_bytes(f, b, sizeof b, c)) return -1;
+  *v = 0; for (int i = 0; i < 8; ++i) *v |= (uint64_t)b[i] << (8*i); return 0;
+}
+static int write_double(FILE *f, double v, uint32_t *c) {
+  uint64_t bits; memcpy(&bits, &v, sizeof bits); return write_u64(f, bits, c);
+}
+static int read_double(FILE *f, double *v, uint32_t *c) {
+  uint64_t bits; if (read_u64(f, &bits, c)) return -1; memcpy(v, &bits, sizeof bits); return 0;
+}
+static int unit_vector(const double v[3]) {
+  const double n2 = v[0]*v[0] + v[1]*v[1] + v[2]*v[2];
+  return isfinite(n2) && fabs(n2 - 1.0) <= 1e-9;
+}
+static int valid_mesh(const FrameLensMesh *m) {
+  if (m == NULL || m->vertex_count == 0 || m->triangle_count == 0) return 0;
+  for (size_t i = 0; i < m->vertex_count; ++i) {
+    const LensVertex *v = &m->vertices[i];
+    if (!v->traced || v->status < RAY_ENDPOINT_ESCAPED ||
+        v->status > RAY_ENDPOINT_INTEGRATION_FAILURE || !isfinite(v->image_x) ||
+        !isfinite(v->image_y) || !isfinite(v->log_frequency_ratio) ||
+        !unit_vector(v->camera_direction) ||
+        (v->status == RAY_ENDPOINT_ESCAPED && !unit_vector(v->n_infinity))) return 0;
+  }
+  for (size_t i = 0; i < m->triangle_count; ++i)
+    for (int j = 0; j < 3; ++j)
+      if (m->triangles[i].vertex[j] >= m->vertex_count) return 0;
+  return 1;
+}
+
+int lens_map_write(const char *path, int width, int height, double fov,
+                   const LensMapFrame *frames, size_t frame_count) {
+  if (path == NULL || frames == NULL || width <= 0 || height <= 0 ||
+      !isfinite(fov) || fov <= 0.0 || fov >= 179.0 || frame_count == 0 ||
+      frame_count > UINT64_MAX) return -1;
+  for (size_t i = 0; i < frame_count; ++i) if (!valid_mesh(&frames[i].mesh)) return -1;
+  FILE *file = fopen(path, "wb"); if (file == NULL) return -1;
+  int failed = write_bytes(file, lens_map_magic, sizeof lens_map_magic, NULL) ||
+      write_u32(file, LENS_MAP_VERSION, NULL) || write_u32(file, LENS_MAP_ENDIAN, NULL) ||
+      write_u32(file, (uint32_t)width, NULL) || write_u32(file, (uint32_t)height, NULL) ||
+      write_double(file, fov, NULL) || write_u64(file, (uint64_t)frame_count, NULL);
+  for (size_t f = 0; !failed && f < frame_count; ++f) {
+    const FrameLensMesh *m = &frames[f].mesh; uint32_t crc = UINT32_MAX;
+    failed = write_u64(file, frames[f].frame_id, NULL) ||
+        write_double(file, frames[f].coordinate_time, NULL) ||
+        write_double(file, frames[f].proper_time, NULL) ||
+        write_u64(file, (uint64_t)m->vertex_count, NULL) ||
+        write_u64(file, (uint64_t)m->triangle_count, NULL);
+    for (size_t i = 0; !failed && i < m->vertex_count; ++i) {
+      const LensVertex *v = &m->vertices[i];
+      failed = write_double(file, v->image_x, &crc) || write_double(file, v->image_y, &crc);
+      for (int j = 0; !failed && j < 3; ++j) failed = write_double(file, v->camera_direction[j], &crc);
+      for (int j = 0; !failed && j < 3; ++j) failed = write_double(file, v->n_infinity[j], &crc);
+      failed = failed || write_double(file, v->log_frequency_ratio, &crc) ||
+          write_u32(file, (uint32_t)v->status, &crc);
+    }
+    for (size_t i = 0; !failed && i < m->triangle_count; ++i) {
+      for (int j = 0; j < 3; ++j) failed = failed || write_u64(file, m->triangles[i].vertex[j], &crc);
+      failed = failed || write_u32(file, m->triangles[i].level, &crc);
+    }
+    failed = failed || write_u32(file, crc ^ UINT32_MAX, NULL);
+  }
+  if (fclose(file)) failed = 1;
+  return failed ? -1 : 0;
+}
+
+void lens_map_destroy(LensMap *map) {
+  if (map == NULL) return;
+  for (size_t i = 0; i < map->frame_count; ++i) frame_lens_mesh_destroy(&map->frames[i].mesh);
+  free(map->frames); *map = (LensMap){0};
+}
+
+int lens_map_read(const char *path, LensMap *map) {
+  if (path == NULL || map == NULL) return -1;
+  *map = (LensMap){0}; FILE *file = fopen(path, "rb"); if (file == NULL) return -1;
+  unsigned char magic[8]; uint32_t version, endian, width, height; uint64_t count;
+  int failed = read_bytes(file, magic, sizeof magic, NULL) || memcmp(magic, lens_map_magic, sizeof magic) ||
+      read_u32(file, &version, NULL) || read_u32(file, &endian, NULL) ||
+      read_u32(file, &width, NULL) || read_u32(file, &height, NULL) ||
+      read_double(file, &map->horizontal_fov_deg, NULL) || read_u64(file, &count, NULL) ||
+      version != LENS_MAP_VERSION || endian != LENS_MAP_ENDIAN || width == 0 || height == 0 ||
+      width > INT32_MAX || height > INT32_MAX || !isfinite(map->horizontal_fov_deg) ||
+      map->horizontal_fov_deg <= 0.0 || map->horizontal_fov_deg >= 179.0 || count == 0 ||
+      count > SIZE_MAX / sizeof *map->frames;
+  if (failed) goto done;
+  map->width = (int)width; map->height = (int)height; map->frame_count = (size_t)count;
+  map->frames = calloc(map->frame_count, sizeof *map->frames); if (map->frames == NULL) { failed = 1; goto done; }
+  for (size_t f = 0; !failed && f < map->frame_count; ++f) {
+    LensMapFrame *frame = &map->frames[f]; uint64_t vertices, triangles; uint32_t stored_crc, crc = UINT32_MAX;
+    failed = read_u64(file, &frame->frame_id, NULL) || read_double(file, &frame->coordinate_time, NULL) ||
+        read_double(file, &frame->proper_time, NULL) || read_u64(file, &vertices, NULL) || read_u64(file, &triangles, NULL) ||
+        !isfinite(frame->coordinate_time) || !isfinite(frame->proper_time) || vertices == 0 || triangles == 0 ||
+        vertices > SIZE_MAX / sizeof *frame->mesh.vertices || triangles > SIZE_MAX / sizeof *frame->mesh.triangles;
+    if (failed) break;
+    frame->mesh.vertices = calloc((size_t)vertices, sizeof *frame->mesh.vertices);
+    frame->mesh.triangles = calloc((size_t)triangles, sizeof *frame->mesh.triangles);
+    if (frame->mesh.vertices == NULL || frame->mesh.triangles == NULL) { failed = 1; break; }
+    frame->mesh.vertex_count = frame->mesh.vertex_capacity = (size_t)vertices;
+    frame->mesh.triangle_count = frame->mesh.triangle_capacity = (size_t)triangles;
+    for (size_t i = 0; !failed && i < frame->mesh.vertex_count; ++i) {
+      LensVertex *v = &frame->mesh.vertices[i]; uint32_t status;
+      failed = read_double(file, &v->image_x, &crc) || read_double(file, &v->image_y, &crc);
+      for (int j = 0; !failed && j < 3; ++j) failed = read_double(file, &v->camera_direction[j], &crc);
+      for (int j = 0; !failed && j < 3; ++j) failed = read_double(file, &v->n_infinity[j], &crc);
+      failed = failed || read_double(file, &v->log_frequency_ratio, &crc) || read_u32(file, &status, &crc) ||
+          status > RAY_ENDPOINT_INTEGRATION_FAILURE;
+      v->status = (RayEndpointStatus)status; v->traced = 1;
+    }
+    for (size_t i = 0; !failed && i < frame->mesh.triangle_count; ++i) {
+      for (int j = 0; j < 3; ++j) {
+        uint64_t index = 0;
+        if (read_u64(file, &index, &crc) || index > SIZE_MAX) {
+          failed = 1;
+          break;
+        }
+        frame->mesh.triangles[i].vertex[j] = (size_t)index;
+      }
+      failed = failed || read_u32(file, &frame->mesh.triangles[i].level, &crc); frame->mesh.triangles[i].evaluated = 1;
+    }
+    failed = failed || read_u32(file, &stored_crc, NULL) || stored_crc != (crc ^ UINT32_MAX) || !valid_mesh(&frame->mesh);
+  }
+done:
+  if (fclose(file)) failed = 1;
+  if (failed) { lens_map_destroy(map); return -1; }
+  return 0;
+}
@@ -0,0 +1,32 @@
+#ifndef LENS_MAP_H
+#define LENS_MAP_H
+
+#include "frame.h"
+
+#include <stddef.h>
+#include <stdint.h>
+
+/* A finalized, render-only lens map.  It deliberately contains no adaptive
+ * refinement work queues: imported maps may be splatted or drawn, but cannot
+ * be refined without tracing new rays. */
+typedef struct {
+  uint64_t frame_id;
+  double coordinate_time;
+  double proper_time;
+  FrameLensMesh mesh;
+} LensMapFrame;
+
+typedef struct {
+  int width, height;
+  double horizontal_fov_deg;
+  LensMapFrame *frames;
+  size_t frame_count;
+} LensMap;
+
+int lens_map_write(const char *path, int width, int height,
+                   double horizontal_fov_deg, const LensMapFrame *frames,
+                   size_t frame_count);
+int lens_map_read(const char *path, LensMap *map);
+void lens_map_destroy(LensMap *map);
+
+#endif
@@ -1,5 +1,6 @@
 #include "catalog.h"
 #include "frame.h"
+#include "lens_map.h"
 #include "movie.h"
 #include "observer_track.h"
 #include "optics.h"
@@ -33,6 +34,9 @@ typedef struct {
  const char *catalog_path;
  const char *all_sky_catalog_path;
  const char *output_path;
+  const char *lens_map_input_path;
+  const char *lens_map_output_path;
+  int width_specified, height_specified, fov_specified;
 #ifdef ENABLE_HDR_OUTPUT
  int write_hdr_output;
  char hdr_output_path[PATH_MAX];
@@ -216,14 +220,20 @@ static int parse_args(int argc, char **argv, Settings *s,
      s->all_sky_catalog_path = argv[++i];
    else if (!strcmp(argv[i], "--output") && i + 1 < argc)
      s->output_path = argv[++i];
+    else if (!strcmp(argv[i], "--lens-map-input") && i + 1 < argc)
+      s->lens_map_input_path = argv[++i];
+    else if (!strcmp(argv[i], "--lens-map-output") && i + 1 < argc)
+      s->lens_map_output_path = argv[++i];
 #ifdef ENABLE_HDR_OUTPUT
    else if (!strcmp(argv[i], "--hdr-output"))
      s->write_hdr_output = 1;
 #endif
    else if (!strcmp(argv[i], "--width") && i + 1 < argc &&
             !parse_int(argv[++i], &s->width)) {
+      s->width_specified = 1;
    } else if (!strcmp(argv[i], "--height") && i + 1 < argc &&
               !parse_int(argv[++i], &s->height)) {
+      s->height_specified = 1;
    } else if (!strcmp(argv[i], "--coarse-cell-pixels") && i + 1 < argc &&
               !parse_int(argv[++i], &s->coarse_cell_pixels)) {
    } else if (!strcmp(argv[i], "--refine-max-level") && i + 1 < argc &&
@@ -243,6 +253,7 @@ static int parse_args(int argc, char **argv, Settings *s,
      s->draw_mesh = 1;
    } else if (!strcmp(argv[i], "--fov-deg") && i + 1 < argc &&
               !parse_double(argv[++i], &s->horizontal_fov_deg)) {
+      s->fov_specified = 1;
    } else if (!strcmp(argv[i], "--look-ra-deg") && i + 1 < argc &&
               !parse_ra_deg(argv[++i], &s->look_ra_deg)) {
    } else if (!strcmp(argv[i], "--look-dec-deg") && i + 1 < argc &&
@@ -300,6 +311,71 @@ static int parse_args(int argc, char **argv, Settings *s,
  return 0;
 }

+static void print_help(const char *program) {
+#ifdef ENABLE_PNG
+  const char *default_output_path = "output/imgs/minkowski_sky.png";
+#else
+  const char *default_output_path = "output/imgs/minkowski_sky.ppm";
+#endif
+
+  printf("Usage: %s [options]\n\n", program);
+  fputs("Input/output:\n"
+        "  --catalog PATH                 Single CSV catalog (default: assets/sky_grid_5deg.csv)\n"
+        "  --all-sky-catalog DIR          Lazy-loaded all-sky catalog directory (default: disabled)\n"
+        "  --output PATH                  Tonemapped image path (default: ", stdout);
+  printf("%s)\n", default_output_path);
+  fputs("  --lens-map-output FILE         Save finalized ray-traced lens map; keeps rendering normally\n"
+        "  --lens-map-input FILE          Re-render a saved map without observer, spacetime, or ray tracing\n"
+        "                                (map width, height, and FOV are fixed; mutually exclusive with --lens-map-output)\n",
+        stdout);
+  fputs("  --write-catalog PATH           Write the generated test catalog and exit (default: disabled)\n"
+#ifdef ENABLE_HDR_OUTPUT
+        "  --hdr-output                   Also write a single-frame linear HDR FITS file (default: disabled)\n"
+#endif
+        "\nImage and camera:\n"
+        "  --width N                      Image width in pixels (default: 1280)\n"
+        "  --height N                     Image height in pixels (default: 720)\n"
+        "  --fov-deg D                    Horizontal field of view in degrees (default: 30)\n"
+        "  --look-ra-deg D                ICRS look direction right ascension in degrees (default: 90)\n"
+        "  --look-dec-deg D               ICRS look direction declination in degrees (default: -90)\n"
+        "  --exposure E                   Linear exposure multiplier (default: 1e-3)\n"
+        "  --observer-radius R            Observer radius in Schwarzschild units (default: 30)\n"
+        "  --observer-inward-speed V      Inward observer speed as a fraction of c (default: 0)\n"
+        "\nPSF and catalog splatting:\n"
+        "  --psf-fwhm-pixels N            Moffat PSF FWHM in pixels (default: 2.7)\n"
+        "  --psf-moffat-beta N            Moffat PSF beta, greater than 1 (default: 4.5)\n"
+        "  --max-magnification M          Preview magnification cap (default: unlimited)\n"
+        "  --max-cache-psf-flux F         Cache flux limit before direct fallback (default: 1)\n"
+        "  --psf-relative-tail R          Maximum omitted relative PSF tail fraction (default: 1e-8)\n"
+        "  --psf-min-y Y                  Skip events below this linear HDR luminance (default: 0, disabled)\n"
+        "  --psf-direct                   Disable the PSF lookup cache (default: disabled)\n"
+        "  --catalog-load-workers N       All-sky catalog loader workers (default: 4)\n"
+        "\nAdaptive lens mesh:\n"
+        "  --coarse-cell-pixels N         Initial mesh cell size in pixels (default: 16)\n"
+        "  --refine-max-level N           Maximum refinement level (default: 0)\n"
+        "  --refine-angle-abs-deg D       Absolute angular interpolation error limit (default: 0.001)\n"
+        "  --refine-angle-rel R           Relative angular interpolation error limit (default: 0.1)\n"
+        "  --refine-jacobian-min J        Fold-refinement Jacobian threshold (default: 1e-3)\n"
+        "  --refine-min-edge-pixels P     Stop refinement below this edge length (default: 0.5)\n"
+        "  --refine-min-area-pixels2 A    Stop refinement below this triangle area (default: 0.25)\n"
+        "  --draw-mesh                    Draw the final lens mesh overlay (default: disabled)\n"
+        "\nMovie and observer track:\n"
+        "  --observer-track PATH          Observer worldline/tetrad CSV for movie rendering (default: disabled)\n"
+        "  --frames-dir DIR               Write a movie image sequence to this directory (default: disabled)\n"
+        "  --frames-prefix NAME           Movie frame filename prefix (default: frame)\n"
+        "  --start-time T                 Movie start coordinate time (default: 0)\n"
+        "  --duration T                   Movie duration (default: 2)\n"
+        "  --fps N                        Movie frame rate (default: 30)\n"
+        "  --slab-duration T              Metric time-slab duration (default: 64)\n"
+        "  --proper-acceleration A        Minkowski accelerated-track proper acceleration (default: 1.52)\n"
+        "  --write-minkowski-accel-track PATH\n"
+        "                                Write a generated Minkowski acceleration track and exit (default: disabled)\n"
+        "\nOther:\n"
+        "  --verbose                      Print low-frequency rendering progress (default: disabled)\n"
+        "  --help                         Print this help and exit\n",
+        stdout);
+}
+
 typedef struct {
  int verbose;
  size_t frame_id;
@@ -461,6 +537,21 @@ static int render_observer_frame(const Settings *s, StarCatalog *catalog,
            "%zu vertices, %zu triangles.\n",
            omp_get_wtime() - refinement_start, mesh.vertex_count,
            mesh.triangle_count);
+  if (s->lens_map_output_path != NULL) {
+    const LensMapFrame map_frame = {.frame_id = 0,
+                                    .coordinate_time = 0.0,
+                                    .proper_time = 0.0,
+                                    .mesh = mesh};
+    if (lens_map_write(s->lens_map_output_path, s->width, s->height,
+                       s->horizontal_fov_deg, &map_frame, 1)) {
+      fprintf(stderr, "Failed to write lens map: %s\n", s->lens_map_output_path);
+      frame_lens_mesh_destroy(&mesh);
+      free(hdr);
+      return -1;
+    }
+    fprintf(stderr, "Wrote lens map: %s (%zu vertices, %zu triangles)\n",
+            s->lens_map_output_path, mesh.vertex_count, mesh.triangle_count);
+  }
  if (s->verbose)
    fprintf(stderr, "Frame 0: traced %zu lens vertices; starting catalog render.\n",
            mesh.vertex_count);
@@ -481,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
@@ -636,6 +732,25 @@ static int render_movie(const Settings *s, StarCatalog *catalog,
    if (traced == 0)
      break;
  }
+  if (s->lens_map_output_path != NULL) {
+    LensMapFrame *map_frames = calloc(movie.frame_count, sizeof *map_frames);
+    if (map_frames == NULL) goto done;
+    for (size_t i = 0; i < movie.frame_count; ++i)
+      map_frames[i] = (LensMapFrame){.frame_id = movie.frames[i].frame_id,
+                                     .coordinate_time = movie.frames[i].coordinate_time,
+                                     .proper_time = movie.frames[i].proper_time,
+                                     .mesh = movie.frames[i].mesh};
+    const int write_failed = lens_map_write(s->lens_map_output_path, s->width,
+                                            s->height, s->horizontal_fov_deg,
+                                            map_frames, movie.frame_count);
+    free(map_frames);
+    if (write_failed) {
+      fprintf(stderr, "Failed to write lens map: %s\n", s->lens_map_output_path);
+      goto done;
+    }
+    fprintf(stderr, "Wrote lens-map movie: %s (%zu frames)\n",
+            s->lens_map_output_path, movie.frame_count);
+  }
  for (size_t i = 0; i < movie.frame_count; ++i) {
    char output_path[PATH_MAX];
    double *hdr = calloc((size_t)s->width * s->height * 3, sizeof *hdr);
@@ -661,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);
@@ -678,6 +797,79 @@ done:
  return result;
 }

+static int render_lens_map(const Settings *s, StarCatalog *catalog) {
+  LensMap map = {0};
+  if (lens_map_read(s->lens_map_input_path, &map)) {
+    fprintf(stderr, "Failed to read or validate lens map: %s\n", s->lens_map_input_path);
+    return -1;
+  }
+  if ((s->width_specified && s->width != map.width) ||
+      (s->height_specified && s->height != map.height) ||
+      (s->fov_specified && fabs(s->horizontal_fov_deg - map.horizontal_fov_deg) > 1e-12)) {
+    fputs("Imported lens-map geometry conflicts with --width, --height, or --fov-deg.\n", stderr);
+    lens_map_destroy(&map);
+    return -1;
+  }
+  if ((map.frame_count == 1 && s->frames_dir != NULL) ||
+      (map.frame_count > 1 && s->frames_dir == NULL)) {
+    fputs("A single-frame lens map uses --output; a multi-frame lens map requires --frames-dir.\n", stderr);
+    lens_map_destroy(&map);
+    return -1;
+  }
+  int result = 0;
+  for (size_t i = 0; i < map.frame_count; ++i) {
+    const char *output_path = s->output_path;
+    char movie_path[PATH_MAX];
+    if (map.frame_count > 1) {
+      if (frame_output_path(movie_path, s, (size_t)map.frames[i].frame_id)) {
+        fputs("Could not construct imported-map movie output path.\n", stderr);
+        result = -1; break;
+      }
+      output_path = movie_path;
+    }
+    double *hdr = calloc((size_t)map.width * map.height * 3, sizeof *hdr);
+    if (hdr == NULL) { result = -1; break; }
+    CatalogPrefetchStats prefetch = {0};
+    PsfSplatStats psf_stats = {0};
+    RenderProgress progress = {.verbose = s->verbose,
+                               .frame_id = (size_t)map.frames[i].frame_id,
+                               .prefetch = &prefetch,
+                               .catalog_load_workers = s->catalog_load_workers,
+                               .all_sky_catalog = catalog->kind == STAR_CATALOG_ALL_SKY};
+    const size_t images = frame_splat_catalog(
+        &map.frames[i].mesh, catalog, hdr, map.width, map.height, s->exposure,
+        &s->psf, &s->psf_cache, s->max_magnification, s->max_cache_psf_flux,
+        s->psf_relative_tail, s->psf_min_y, 0, s->catalog_load_workers,
+        &prefetch, &psf_stats,
+        &(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
+    if (s->write_hdr_output &&
+        write_hdr_fits(s->hdr_output_path, hdr, map.width, map.height,
+                       map.horizontal_fov_deg)) {
+      fprintf(stderr, "Failed to write HDR FITS image: %s\n", s->hdr_output_path);
+      free(hdr); result = -1; break;
+    }
+#endif
+    const int write_result = write_tonemapped_image(output_path, hdr, map.width, map.height);
+    free(hdr);
+    fprintf(stderr, "Rendered %zu images from %zu catalog stars to %s (%s; imported lens map)\n",
+            images, catalog->count, output_path, write_result == 0 ? "ok" : "write failed");
+    psf_kernel_cache_report(&s->psf_cache, &psf_stats, stderr);
+    if (write_result) { result = -1; break; }
+  }
+  lens_map_destroy(&map);
+  return result;
+}
+
 static int write_minkowski_accel_track(const Settings *s) {
  ObserverTrack track = {0};
  const int result = observer_track_generate_minkowski_acceleration(
@@ -694,10 +886,15 @@ static int write_minkowski_accel_track(const Settings *s) {
 int main(int argc, char **argv) {
  Settings settings;
  const char *write_path;
+  if (argc == 2 && !strcmp(argv[1], "--help")) {
+    print_help(argv[0]);
+    return 0;
+  }
  if (parse_args(argc, argv, &settings, &write_path)) {
    fprintf(stderr,
            "Usage: %s [--catalog PATH | --all-sky-catalog DIR] [--output PATH] [--width N] [--height "
            "N] [--fov-deg D] [--look-ra-deg D] [--look-dec-deg D] "
+            "[--lens-map-input FILE | --lens-map-output FILE] "
            "[--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] "
@@ -729,6 +926,10 @@ int main(int argc, char **argv) {
    return write_minkowski_accel_track(&settings) == 0
               ? 0
               : (perror(settings.write_minkowski_accel_track_path), 1);
+  if (settings.lens_map_input_path != NULL && settings.lens_map_output_path != NULL) {
+    fputs("--lens-map-input and --lens-map-output are mutually exclusive.\n", stderr);
+    return 2;
+  }
 #ifdef ENABLE_HDR_OUTPUT
  if (settings.frames_dir != NULL && settings.write_hdr_output) {
    fputs("--hdr-output is available only for a single-frame render.\n", stderr);
@@ -759,16 +960,23 @@ int main(int argc, char **argv) {
    }
    fprintf(stderr, "Created test catalog: %s\n", settings.catalog_path);
  }
-  SpacetimeSource spacetime = {0};
-  if (spacetime_create_default(&spacetime)) {
-    fputs("Could not create spacetime source\n", stderr);
-    catalog_destroy(&catalog);
-    return 1;
-  }
  if (!settings.psf_direct && psf_kernel_cache_init(&settings.psf_cache, &settings.psf,
                                                     settings.psf_relative_tail))
    fputs("PSF cache construction failed; using direct evaluator.\n", stderr);
  psf_kernel_cache_report_ready(&settings.psf_cache, stderr);
+  if (settings.lens_map_input_path != NULL) {
+    const int result = render_lens_map(&settings, &catalog);
+    catalog_destroy(&catalog);
+    psf_kernel_cache_destroy(&settings.psf_cache);
+    return result == 0 ? 0 : 1;
+  }
+  SpacetimeSource spacetime = {0};
+  if (spacetime_create_default(&spacetime)) {
+    fputs("Could not create spacetime source\n", stderr);
+    catalog_destroy(&catalog);
+    psf_kernel_cache_destroy(&settings.psf_cache);
+    return 1;
+  }
  int result = settings.frames_dir != NULL
                   ? render_movie(&settings, &catalog, &spacetime)
                   : render_frame(&settings, &catalog, &spacetime);
@@ -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,
@@ -32,6 +32,19 @@ 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;
+
+#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);
@@ -61,6 +74,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);
@@ -70,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
@@ -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
@@ -1,4 +1,5 @@
 #include "frame.h"
+#include "lens_map.h"
 #include "optics.h"

 #include <math.h>
@@ -6,6 +7,7 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
+#include <unistd.h>

 static int mesh_has_hanging_vertex(const FrameLensMesh *mesh) {
  for (size_t triangle = 0; triangle < mesh->triangle_count; ++triangle)
@@ -84,6 +86,55 @@ int main(void) {
    fputs("flat-space inverse lens-map regression failed\n", stderr);
    goto done;
  }
+  /* A finalized mesh can be persisted independently of spacetime and then
+   * drive the exact same catalog inverse-map and PSF pass. */
+  const char *lens_map_path = "/tmp/gr_lens_map_test.grlens";
+  const LensMapFrame saved_frame = {.frame_id = 7,
+                                    .coordinate_time = 3.0,
+                                    .proper_time = 2.0,
+                                    .mesh = mesh};
+  LensMap loaded_map = {0};
+  double *roundtrip_hdr = calloc((size_t)width * height * 3, sizeof *roundtrip_hdr);
+  if (roundtrip_hdr == NULL ||
+      lens_map_write(lens_map_path, width, height, 30.0, &saved_frame, 1) ||
+      lens_map_read(lens_map_path, &loaded_map) || loaded_map.frame_count != 1 ||
+      loaded_map.frames[0].frame_id != 7 || loaded_map.width != width ||
+      loaded_map.height != height ||
+      loaded_map.frames[0].mesh.vertex_count != mesh.vertex_count ||
+      frame_splat_catalog(&loaded_map.frames[0].mesh, &catalog, roundtrip_hdr,
+                          width, height, test_exposure, &psf, NULL, INFINITY,
+                          1.0, psf_relative_tail, 0.0, 0, 1, NULL, NULL, NULL) != images) {
+    fputs("lens-map round-trip regression failed\n", stderr);
+    free(roundtrip_hdr); lens_map_destroy(&loaded_map); unlink(lens_map_path);
+    goto done;
+  }
+  for (int value = 0; value < width * height * 3; ++value)
+    if (hdr[value] != roundtrip_hdr[value]) {
+      fputs("lens-map round-trip HDR regression failed\n", stderr);
+      free(roundtrip_hdr); lens_map_destroy(&loaded_map); unlink(lens_map_path);
+      goto done;
+    }
+  free(roundtrip_hdr);
+  lens_map_destroy(&loaded_map);
+  /* A damaged payload must not be mistaken for a reusable physical map. */
+  FILE *damaged = fopen(lens_map_path, "r+b");
+  int damage_failed = damaged == NULL;
+  if (!damage_failed) {
+    if (fseek(damaged, -5L, SEEK_END))
+      damage_failed = 1;
+    const int original = damage_failed ? EOF : fgetc(damaged);
+    if (damage_failed || fseek(damaged, -5L, SEEK_END) || original == EOF ||
+        fputc(original ^ 0xff, damaged) == EOF)
+      damage_failed = 1;
+  }
+  if (damaged != NULL && fclose(damaged))
+    damage_failed = 1;
+  if (damage_failed || !lens_map_read(lens_map_path, &loaded_map)) {
+    fputs("lens-map corruption rejection regression failed\n", stderr);
+    lens_map_destroy(&loaded_map); unlink(lens_map_path);
+    goto done;
+  }
+  unlink(lens_map_path);
  memset(hdr, 0, (size_t)width * height * 3 * sizeof *hdr);
  PsfSplatStats min_y_stats = {0};
  if (frame_splat_catalog(&mesh, &catalog, hdr, width, height, test_exposure,
@@ -171,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);
@@ -0,0 +1,60 @@
+#include "hip_psf.h"
+
+#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+
+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;
+}
Author	SHA1	Message	Date
wyj	7a80086f63	Feat: add initial HIP PSF backend	2026-09-05 04:23:57 -04:00
wyj	e48c490337	Doc: update GPU PSF acceleration status	2026-09-05 03:47:42 -04:00
wyj	e8bc553128	DOC: preserve complete benchmark output	2026-09-05 03:22:15 -04:00
wyj	92d8d294b5	Bench: record PSF event sink overhead	2026-09-05 03:20:40 -04:00
wyj	2a03cbf724	Frame: buffer cached PSF events per worker	2026-09-05 03:06:46 -04:00
wyj	4953775cdd	Test: add fixed PSF image references	2026-09-04 20:28:54 -04:00
wyj	ab56e23fc9	Doc: document reusable lens maps	2026-09-03 20:58:20 -04:00
wyj	d43ad262fc	CLI: document lens map options	2026-09-03 20:55:08 -04:00
wyj	1ee3a8cc26	Feat: add reusable lens map files	2026-09-03 19:46:05 -04:00
wyj	1a54e59eb4	Output: add CLI help text	2026-09-03 19:24:32 -04:00