From 9cd933d1f83cf768249031402cb347a159a70f2d Mon Sep 17 00:00:00 2001 From: Yingjie Wang Date: Sun, 27 Sep 2026 04:39:38 -0400 Subject: [PATCH] Feat: Add optional three-channel sensor bloom model Add an opt-in, post-processing limited-response model applied to the finished linear HDR before tone mapping. Each RGB channel is processed independently and isotropically: overflow above E spreads to the eight neighbours with a fixed 9-point stencil, while the rest is absorbed or lost at the image boundary. The synchronous ping-pong update uses a monotonic bounding box and a row-parallel, deterministic reduction; the conservative round bound reserves fp guard rounds inside a 4096 hard limit and fails before touching HDR when exceeded. Expose --sensor-bloom-limit E and --sensor-bloom-transfer e (both required together, default disabled), validate them before expensive initialization, and route every output path through the same hook in write_frame_outputs: raw FITS first, bloom, tone-mapped PNG/PPM, then the mesh overlay. The raw --hdr-output FITS therefore stays pre-bloom. Add a standalone unit test (stencil, boundary loss, cascade reference, symmetry, thread determinism, convergence limits, validation, allocation failure), CLI integration and regression coverage, an isolated sensor-bloom-bench target, and document the model in the design, usage, README, and build docs. --- Makefile | 21 +- README.md | 6 +- README.zh-CN.md | 2 +- build.md | 6 +- nr_spacetime_movie_renderer_design.md | 12 + src/main.c | 67 +++ src/sensor_bloom.c | 320 +++++++++++++++ src/sensor_bloom.h | 45 ++ tests/benchmark_sensor_bloom.c | 121 ++++++ tests/test_camera_cli.py | 89 ++++ tests/test_sensor_bloom.c | 568 ++++++++++++++++++++++++++ usage.md | 30 +- 12 files changed, 1280 insertions(+), 7 deletions(-) create mode 100644 src/sensor_bloom.c create mode 100644 src/sensor_bloom.h create mode 100644 tests/benchmark_sensor_bloom.c create mode 100644 tests/test_sensor_bloom.c diff --git a/Makefile b/Makefile index 066710d..cd00a00 100644 --- a/Makefile +++ b/Makefile @@ -36,7 +36,7 @@ TARGET_BASENAME := $(SPACETIME)_sky OBJECT_DIR := $(BUILD_DIR)/obj/$(SPACETIME) CORE_MINKOWSKI_SOURCES := $(COMMON_SOURCES) src/spacetime_minkowski.c -.PHONY: all backend clean run test tone-map-test hip-psf-test hip-psf-bench fast-psf-fftw-bench minkowski schwarzschild FORCE +.PHONY: all backend clean run test tone-map-test sensor-bloom-test sensor-bloom-bench hip-psf-test hip-psf-bench fast-psf-fftw-bench 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) @@ -114,6 +114,8 @@ CAMERA_TEST_TARGETS := $(TEST_OUT_DIR)/test_observer_minkowski $(TEST_OUT_DIR)/t FAST_PSF_FFTW_TEST_TARGET := $(TEST_OUT_DIR)/test_fast_psf_fftw FAST_PSF_FFTW_BENCH_TARGET := $(TEST_OUT_DIR)/benchmark_fast_psf_fftw TONE_MAP_TEST_TARGET := $(TEST_OUT_DIR)/test_tone_map +SENSOR_BLOOM_TEST_TARGET := $(TEST_OUT_DIR)/test_sensor_bloom +SENSOR_BLOOM_BENCH_TARGET := $(TEST_OUT_DIR)/benchmark_sensor_bloom ifeq ($(PSF_BACKEND),hip) TARGET := $(BUILD_DIR)/$(TARGET_BASENAME)_hip @@ -224,6 +226,15 @@ $(TONE_MAP_TEST_TARGET): tests/test_tone_map.c src/optics.c src/optics.h $(CPU_F $(FAST_PSF_FFTW_BENCH_TARGET): tests/benchmark_fast_psf_fftw.c $(CORE_MINKOWSKI_SOURCES) $(CPU_FFTW_SOURCES) | $(TEST_OUT_DIR) $(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc $^ $(LDLIBS) -o $@ +# The sensor-bloom checks link only the standalone model, so they need neither +# a catalog nor ray tracing, and build in every ENABLE_HDR/PSF_BACKEND +# configuration. +$(SENSOR_BLOOM_TEST_TARGET): tests/test_sensor_bloom.c src/sensor_bloom.c src/sensor_bloom.h | $(TEST_OUT_DIR) + $(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc tests/test_sensor_bloom.c src/sensor_bloom.c $(LDLIBS) -o $@ + +$(SENSOR_BLOOM_BENCH_TARGET): tests/benchmark_sensor_bloom.c src/sensor_bloom.c src/sensor_bloom.h | $(TEST_OUT_DIR) + $(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc tests/benchmark_sensor_bloom.c src/sensor_bloom.c $(LDLIBS) -o $@ + # The FFTW-vs-spatial test is meaningful only in the CPU PSF build. ifneq ($(CPU_FFTW_SOURCES),) FAST_PSF_FFTW_TEST_DEP := $(FAST_PSF_FFTW_TEST_TARGET) @@ -233,7 +244,7 @@ FAST_PSF_FFTW_TEST_DEP := FAST_PSF_FFTW_TEST_RUN := endif -test: $(CAMERA_TEST_TARGETS) $(TEST_TARGET) $(FRAME_TEST_TARGET) $(SCHWARZSCHILD_TEST_TARGET) $(OBSERVER_TRACK_TEST_TARGET) $(CATALOG_PREFETCH_TEST_TARGET) $(FAST_PSF_FFTW_TEST_DEP) $(TONE_MAP_TEST_TARGET) +test: $(CAMERA_TEST_TARGETS) $(TEST_TARGET) $(FRAME_TEST_TARGET) $(SCHWARZSCHILD_TEST_TARGET) $(OBSERVER_TRACK_TEST_TARGET) $(CATALOG_PREFETCH_TEST_TARGET) $(FAST_PSF_FFTW_TEST_DEP) $(TONE_MAP_TEST_TARGET) $(SENSOR_BLOOM_TEST_TARGET) $(TEST_OUT_DIR)/test_observer_minkowski $(TEST_OUT_DIR)/test_observer_schwarzschild $(TEST_TARGET) @@ -243,11 +254,17 @@ test: $(CAMERA_TEST_TARGETS) $(TEST_TARGET) $(FRAME_TEST_TARGET) $(SCHWARZSCHILD $(CATALOG_PREFETCH_TEST_TARGET) $(FAST_PSF_FFTW_TEST_RUN) $(TONE_MAP_TEST_TARGET) + $(SENSOR_BLOOM_TEST_TARGET) python3 tests/test_camera_cli.py $(BUILD_DIR) $(TEST_OUT_DIR) tone-map-test: $(TONE_MAP_TEST_TARGET) $(TONE_MAP_TEST_TARGET) +sensor-bloom-test: $(SENSOR_BLOOM_TEST_TARGET) + $(SENSOR_BLOOM_TEST_TARGET) + +sensor-bloom-bench: $(SENSOR_BLOOM_BENCH_TARGET) + fast-psf-fftw-bench: $(FAST_PSF_FFTW_BENCH_TARGET) clean: diff --git a/README.md b/README.md index 11d2fed..e976ded 100644 --- a/README.md +++ b/README.md @@ -138,8 +138,10 @@ be configured for the desired image accuracy; it is disabled by default. Tone-mapped PNG/PPM output uses a per-channel soft-clip display transform by default (`--tone-map softclip --tone-map-p 2`); `--tone-map reinhard` restores the historical curve used by the reference images above. The linear HDR FITS -output never applies it. See [tone mapping and -display](usage.md#tone-mapping-and-display). +output never applies it. An optional limited-response sensor bloom +(`--sensor-bloom-limit E --sensor-bloom-transfer e`, default disabled) can be +applied to the linear HDR before display; it never changes the FITS output. See +[tone mapping and display](usage.md#tone-mapping-and-display). Camera controls, movie sequences, lens-map reuse, PSF settings, and HDR output are described in [usage.md](usage.md). For fast single-frame previews, `--fast-mode` replaces per-event PSF splats with a supersampled delta deposit diff --git a/README.zh-CN.md b/README.zh-CN.md index ebbf81b..caa1250 100644 --- a/README.zh-CN.md +++ b/README.zh-CN.md @@ -86,7 +86,7 @@ mkdir -p output/imgs 这里使用更新后的参考图像的渲染设置。`--max-cache-psf-flux 1e8` 允许将明亮 PSF 的翼部截断在缓存半径内，但本次运行没有发生翼部裁剪或 direct fallback。[CPU/HIP benchmark 记录](benchmarks/2mass_galactic_center_cpu_hip_2026-09-07.md)保留了命令、终端输出，以及 CPU/HIP PNG 比特级一致的结果；上述命令省略了可选的 HDR 和 lens-map 导出。 -上述 Schwarzschild 示例省略位置，因此按指向与半径推导出朝向黑洞的相机。自适应细分默认关闭，应根据所需图像精度配置。tone-mapped PNG/PPM 默认使用逐通道 soft-clip 显示变换（`--tone-map softclip --tone-map-p 2`）；`--tone-map reinhard` 恢复上方参考图像使用的历史曲线，线性 HDR FITS 输出不受其影响。相机控制、图像序列、透镜映射复用、PSF 设置及 HDR 输出参见 [usage.md](usage.md)（英文）。两个可执行文件都可通过 `--help` 查看完整选项。 +上述 Schwarzschild 示例省略位置，因此按指向与半径推导出朝向黑洞的相机。自适应细分默认关闭，应根据所需图像精度配置。tone-mapped PNG/PPM 默认使用逐通道 soft-clip 显示变换（`--tone-map softclip --tone-map-p 2`）；`--tone-map reinhard` 恢复上方参考图像使用的历史曲线，线性 HDR FITS 输出不受其影响。可选的有限响应 sensor bloom（`--sensor-bloom-limit E --sensor-bloom-transfer e`，默认关闭）在显示前作用于线性 HDR，且不会改变 FITS 输出。相机控制、图像序列、透镜映射复用、PSF 设置及 HDR 输出参见 [usage.md](usage.md)（英文）。两个可执行文件都可通过 `--help` 查看完整选项。 ### 示例：紧贴史瓦西视界向外看 diff --git a/build.md b/build.md index 1ffc96e..1b42fcb 100644 --- a/build.md +++ b/build.md @@ -173,7 +173,11 @@ make test ``` The CPU suite includes both coordinate-camera builders, CLI validation, -small PNG renders, and single-frame/movie lens-map agreement. The reference +small PNG renders, single-frame/movie lens-map agreement, the tone-map +regression, and the standalone sensor-bloom model regression +(`make sensor-bloom-test`). An isolated `make sensor-bloom-bench` target times +only the bloom model on synthetic frames; it is not part of `make test`. The +reference image checks require Python 3 and CFITSIO. The original Schwarzschild HDR fixture is retained with a float32 relative tolerance of `2^-23` (zero absolute tolerance): replacing the specialized static tetrad with metric-based diff --git a/nr_spacetime_movie_renderer_design.md b/nr_spacetime_movie_renderer_design.md index c260314..71e27e5 100644 --- a/nr_spacetime_movie_renderer_design.md +++ b/nr_spacetime_movie_renderer_design.md @@ -1467,6 +1467,18 @@ void rays_trace_generation( pre-tone-map 线性 RGB。这里的 operator 只是当前简化显示管线，不是最终相机/ 传感器模型。 +可选的后处理式 sensor bloom（`--sensor-bloom-limit E +--sensor-bloom-transfer e`，默认关闭）在 PSF 累积完成的线性 HDR 上、display +operator 之前运行，CPU 与 HIP、标准与 fast-mode、直接 trace 与 imported lens +map 共用同一入口。每个 RGB 通道独立、各向同性地把超过有限响应上限 \(E\) 的 +溢出按固定比例 \(e\) 传播到 8 邻域，其余 \((1-e)D\) 被 drain 吸收，指向图像外的 +部分在边界丢失。它是唯象的有限响应近似，不模拟特定 CCD/CMOS/Bayer 结构； +\(E\) 使用 exposure 之后的 renderer-scale 线性 HDR；\(e\) 同时决定每轮保留传播的 +比例和有效传播距离；模型允许信号损失，不守恒。原始 `--hdr-output` FITS 在模型 +运行前写出，因此始终是 bloom 前的 PSF HDR；tone-mapped PNG/PPM 与视频帧在模型 +之后写出。视觉式多尺度 bloom 不在当前范围内。mesh overlay 在模型之后绘制， +不参与溢出传播。 + --- # 28. 开发顺序 diff --git a/src/main.c b/src/main.c index a152178..253f95a 100644 --- a/src/main.c +++ b/src/main.c @@ -5,6 +5,7 @@ #include "observer_track.h" #include "optics.h" #include "ray.h" +#include "sensor_bloom.h" #include "spacetime.h" #include @@ -60,6 +61,11 @@ typedef struct { const char *blackbody_table_path; RefinementConfig refinement; ToneMapSettings tone_map; + int sensor_bloom_enabled; + int sensor_bloom_limit_specified; + int sensor_bloom_transfer_specified; + double sensor_bloom_limit; + double sensor_bloom_transfer; } Settings; static int parse_int(const char *text, int *value) { @@ -145,6 +151,17 @@ static int parse_finite_positive(const char *text, double *value) { return errno || end == text || *end || !isfinite(*value) || *value <= 0.0 ? -1 : 0; } +/* Sensor-bloom transfer coefficient: finite and in [0, 1). */ +static int parse_sensor_bloom_transfer(const char *text, double *value) { + char *end; + errno = 0; + *value = strtod(text, &end); + return errno || end == text || *end || !isfinite(*value) || *value < 0.0 || + *value >= 1.0 + ? -1 + : 0; +} + static int validate_tonemapped_output_path(const char *path) { const size_t path_length = strlen(path); #ifdef ENABLE_PNG @@ -241,6 +258,30 @@ static int write_frame_outputs(const Settings *s, const FrameLensMesh *mesh, (void)s; (void)fov_deg; #endif + if (s->sensor_bloom_enabled) { + const SensorBloomSettings bloom = { + .response_limit = s->sensor_bloom_limit, + .transfer = s->sensor_bloom_transfer}; + SensorBloomStats bloom_stats; + if (sensor_bloom_apply(hdr, width, height, &bloom, &bloom_stats)) { + fputs("Sensor bloom failed (invalid settings, non-finite HDR, allocation " + "failure, or a round bound above the 4096 hard limit; lower " + "--sensor-bloom-transfer, raise --sensor-bloom-limit, or reduce " + "--exposure); aborting tone-mapped output.\n", + stderr); + return -1; + } + fprintf(stderr, + "Sensor bloom: saturated=%zu clamped=%zu iterations=%zu/%zu " + "peak=%.6g max_overflow=%.6g absorbed=%.6g boundary=%.6g " + "residual=%.6g elapsed=%.6fs\n", + bloom_stats.initially_saturated_channels, + bloom_stats.final_clamped_channels, bloom_stats.iterations, + bloom_stats.predicted_iterations, bloom_stats.peak_input, + bloom_stats.initial_max_overflow, bloom_stats.absorbed_signal, + bloom_stats.boundary_loss, bloom_stats.residual_clamp_loss, + bloom_stats.elapsed_seconds); + } const int write_result = write_tonemapped_image(paths->output_path, hdr, width, height, &s->tone_map); @@ -369,6 +410,13 @@ static int parse_args(int argc, char **argv, Settings *s, } else if (!strcmp(argv[i], "--tone-map-p") && i + 1 < argc && !parse_at_least_one(argv[++i], &s->tone_map.p)) { tone_map_p_specified = 1; + } else if (!strcmp(argv[i], "--sensor-bloom-limit") && i + 1 < argc && + !parse_finite_positive(argv[++i], &s->sensor_bloom_limit)) { + s->sensor_bloom_limit_specified = 1; + } else if (!strcmp(argv[i], "--sensor-bloom-transfer") && i + 1 < argc && + !parse_sensor_bloom_transfer(argv[++i], + &s->sensor_bloom_transfer)) { + s->sensor_bloom_transfer_specified = 1; } else if ((!strcmp(argv[i], "--observer-position") || !strcmp(argv[i], "--observer-velocity")) && i + 3 < argc) { const int position = !strcmp(argv[i], "--observer-position"); @@ -447,6 +495,14 @@ static int parse_args(int argc, char **argv, Settings *s, fputs("--tone-map-p applies only to --tone-map softclip.\n", stderr); return -1; } + if (s->sensor_bloom_limit_specified != + s->sensor_bloom_transfer_specified) { + fputs("--sensor-bloom-limit and --sensor-bloom-transfer must be " + "specified together.\n", + stderr); + return -1; + } + s->sensor_bloom_enabled = s->sensor_bloom_limit_specified; return 0; } @@ -481,6 +537,16 @@ static void print_help(const char *program) { " --tone-map MODE Display transform: softclip or reinhard\n" " (default: softclip)\n" " --tone-map-p P Softclip hardness P >= 1 (default: 2)\n" + " --sensor-bloom-limit E Enable the optional three-channel sensor overflow\n" + " model; E is the finite response limit on the\n" + " post-exposure linear HDR renderer scale\n" + " (default: disabled)\n" + " --sensor-bloom-transfer e Limited-response overflow transfer ratio in [0,1)\n" + " (e=0 clamps only). Both bloom options must appear\n" + " together. The isotropic, non-conservative model is a\n" + " phenomenological approximation, not a specific\n" + " CCD/CMOS/Bayer structure; the raw --hdr-output FITS\n" + " never includes it.\n" " --observer-position X Y Z Coordinate position; alone implies looking at the origin\n" " --observer-radius R Infer position = -R * look direction (default R: 30); conflicts with position\n" " --observer-velocity VX VY VZ Coordinate dx/dt, dy/dt, dz/dt (default: 0 0 0); must be timelike\n" @@ -1180,6 +1246,7 @@ int main(int argc, char **argv) { "N] [--fov-deg D] [--look-ra-deg D] [--look-dec-deg D] " "[--lens-map-input FILE | --lens-map-output FILE] " "[--exposure E] [--tone-map softclip|reinhard] [--tone-map-p P] " + "[--sensor-bloom-limit E --sensor-bloom-transfer e] " "[--observer-radius R | --observer-position X Y Z] " "[--observer-velocity VX VY VZ] [--camera-roll-deg ANGLE] " "[--psf-fwhm-pixels N] [--psf-moffat-beta N] " diff --git a/src/sensor_bloom.c b/src/sensor_bloom.c new file mode 100644 index 0000000..b06f901 --- /dev/null +++ b/src/sensor_bloom.c @@ -0,0 +1,320 @@ +#include "sensor_bloom.h" + +#include +#include +#include +#include +#include + +/* Fixed 9-point isotropic stencil: the four axial neighbours each carry 4/20 + * and the four diagonal neighbours each carry 1/20, so an interior pixel + * spreads exactly its full overflow. Edge and corner pixels are not + * renormalized: the missing weight leaves the image and is counted as loss. */ +#define SB_AXIAL (4.0 / 20.0) +#define SB_DIAGONAL (1.0 / 20.0) +#define SB_INTERNAL_TOLERANCE (1e-9) +#define SB_MAX_ROUNDS ((size_t)4096) +#define SB_GUARD_ROUNDS ((size_t)8) + +/* Sequential sum in the same order the gather loop accumulates it. The + * real-valued stencil sums to 1, but the double additions round up by one ULP + * to 1.0000000000000002; using that upper bound lets the round prediction stay + * conservative instead of assuming a decay of exactly `transfer`. */ +static double stencil_weight_total(void) { + double total = 0.0; + total += SB_AXIAL; + total += SB_AXIAL; + total += SB_AXIAL; + total += SB_AXIAL; + total += SB_DIAGONAL; + total += SB_DIAGONAL; + total += SB_DIAGONAL; + total += SB_DIAGONAL; + return total; +} + +static double overflow_at(const double *src, size_t index, double limit) { + const double difference = src[index] - limit; + return difference > 0.0 ? difference : 0.0; +} + +/* Total stencil weight leaving the image at (x, y), in [0, 1]. A location + * fully inside the image has zero boundary weight. */ +static double boundary_weight(int x, int y, int width, int height) { + const int has_left = x > 0; + const int has_right = x + 1 < width; + const int has_up = y > 0; + const int has_down = y + 1 < height; + double weight = 0.0; + + if (!has_left) weight += SB_AXIAL; + if (!has_right) weight += SB_AXIAL; + if (!has_up) weight += SB_AXIAL; + if (!has_down) weight += SB_AXIAL; + if (!has_left || !has_up) weight += SB_DIAGONAL; + if (!has_right || !has_up) weight += SB_DIAGONAL; + if (!has_left || !has_down) weight += SB_DIAGONAL; + if (!has_right || !has_down) weight += SB_DIAGONAL; + return weight; +} + +/* Conservative round bound: an interior pixel transfers at most + * `transfer * stencil_weight_total()` of its overflow, so + * D_max^(n+1) <= decay * D_max^(n) with decay = transfer * stencil total. + * Returns -1 if the bound is not finite or exceeds the hard limit, so the + * caller can fail before touching the framebuffer. */ +static int round_bound(double limit, double transfer, double max_overflow, + size_t *rounds) { + const double tolerance = SB_INTERNAL_TOLERANCE * limit; + if (max_overflow <= tolerance) { + *rounds = 0; + return 0; + } + const double decay = transfer * stencil_weight_total(); + if (decay == 0.0) { + *rounds = 1; + return 0; + } + const double estimate = log(tolerance / max_overflow) / log(decay); + if (!isfinite(estimate) || estimate < 0.0 || + estimate > (double)SB_MAX_ROUNDS) + return -1; + const size_t count = (size_t)ceil(estimate); + if (count > SB_MAX_ROUNDS) + return -1; + *rounds = count; + return 0; +} + +int sensor_bloom_apply(double *hdr, int width, int height, + const SensorBloomSettings *settings, + SensorBloomStats *stats) { + SensorBloomStats local = {0}; + if (stats == NULL) + stats = &local; + *stats = (SensorBloomStats){0}; + + if (hdr == NULL || settings == NULL || width <= 0 || height <= 0) + return -1; + const double limit = settings->response_limit; + const double transfer = settings->transfer; + if (!isfinite(limit) || limit <= 0.0 || !isfinite(transfer) || + transfer < 0.0 || transfer >= 1.0) + return -1; + + const size_t size_width = (size_t)width; + const size_t size_height = (size_t)height; + if (size_width > SIZE_MAX / size_height) + return -1; + const size_t pixels = size_width * size_height; + if (pixels > SIZE_MAX / 3) + return -1; + const size_t channels = pixels * 3; + if (channels > SIZE_MAX / sizeof(double)) + return -1; + + const double start = omp_get_wtime(); + + /* Initial scan: finite-input validation, peak, initial saturation count and + * overflow, and the bounding box of every saturated channel. */ + size_t saturated = 0; + double peak_input = -INFINITY; + double initial_max_overflow = 0.0; + int bbox_left = width, bbox_right = -1, bbox_top = height, bbox_bottom = -1; + for (size_t i = 0; i < channels; ++i) { + const double value = hdr[i]; + if (!isfinite(value)) + return -1; + if (value > peak_input) + peak_input = value; + const double difference = value - limit; + if (difference > 0.0) { + ++saturated; + if (difference > initial_max_overflow) + initial_max_overflow = difference; + const size_t pixel = i / 3; + const int x = (int)(pixel % size_width); + const int y = (int)(pixel / size_width); + if (x < bbox_left) bbox_left = x; + if (x > bbox_right) bbox_right = x; + if (y < bbox_top) bbox_top = y; + if (y > bbox_bottom) bbox_bottom = y; + } + } + stats->peak_input = peak_input; + stats->initial_max_overflow = initial_max_overflow; + stats->initially_saturated_channels = saturated; + + if (saturated == 0) { + /* Nothing exceeds E: leave the framebuffer untouched and allocate nothing. */ + stats->elapsed_seconds = omp_get_wtime() - start; + return 0; + } + + size_t predicted = 0; + if (round_bound(limit, transfer, initial_max_overflow, &predicted)) { + stats->elapsed_seconds = omp_get_wtime() - start; + return -1; + } + /* Reserve a few fp guard rounds inside the hard limit, and report the cap + * actually used so `iterations <= predicted_iterations` always holds. At the + * boundary the guard is clamped, never added on top of SB_MAX_ROUNDS. */ + size_t round_cap = predicted; + if (round_cap > 0 && transfer > 0.0) { + round_cap += SB_GUARD_ROUNDS; + if (round_cap > SB_MAX_ROUNDS) + round_cap = SB_MAX_ROUNDS; + } + stats->predicted_iterations = round_cap; + + const double tolerance = SB_INTERNAL_TOLERANCE * limit; + + double *scratch = malloc(channels * sizeof *scratch); + double *row_accumulator = calloc((size_t)height * 3, sizeof *row_accumulator); + if (scratch == NULL || row_accumulator == NULL) { + free(scratch); + free(row_accumulator); + stats->elapsed_seconds = omp_get_wtime() - start; + return -1; + } + double *row_absorbed = row_accumulator; + double *row_boundary = row_accumulator + (size_t)height; + double *row_max = row_accumulator + (size_t)height * 2; + memcpy(scratch, hdr, channels * sizeof *scratch); + + double *source = hdr; + double *destination = scratch; + size_t iterations = 0; + int final_left = bbox_left, final_right = bbox_right; + int final_top = bbox_top, final_bottom = bbox_bottom; + + while (iterations < round_cap) { + /* The bounding box only ever grows, and every round widens it by one pixel + * first so it keeps covering every pixel any earlier round could have + * modified and every pixel this round can reach. It must not shrink with + * the current overflow, or the ping-pong buffers would restore stale + * values. */ + if (bbox_left > 0) --bbox_left; + if (bbox_top > 0) --bbox_top; + if (bbox_right + 1 < width) ++bbox_right; + if (bbox_bottom + 1 < height) ++bbox_bottom; + + const int x0 = bbox_left, x1 = bbox_right; + const int y0 = bbox_top, y1 = bbox_bottom; + final_left = x0; + final_right = x1; + final_top = y0; + final_bottom = y1; + + memset(row_accumulator, 0, (size_t)height * 3 * sizeof *row_accumulator); + +#pragma omp parallel for schedule(static) + for (int y = y0; y <= y1; ++y) { + const size_t row = (size_t)y * size_width; + double absorbed = 0.0, boundary = 0.0, maximum = 0.0; + for (int x = x0; x <= x1; ++x) { + const size_t base = (row + (size_t)x) * 3; + const double edge_weight = boundary_weight(x, y, width, height); + for (int c = 0; c < 3; ++c) { + const size_t index = base + (size_t)c; + const double value = source[index]; + const double overflow = value > limit ? value - limit : 0.0; + absorbed += (1.0 - transfer) * overflow; + boundary += transfer * overflow * edge_weight; + + double incoming = 0.0; + if (x > 0) + incoming += SB_AXIAL * overflow_at(source, (row + x - 1) * 3 + c, limit); + if (x + 1 < width) + incoming += SB_AXIAL * overflow_at(source, (row + x + 1) * 3 + c, limit); + if (y > 0) + incoming += SB_AXIAL * overflow_at(source, (row - size_width + x) * 3 + c, limit); + if (y + 1 < height) + incoming += SB_AXIAL * overflow_at(source, (row + size_width + x) * 3 + c, limit); + if (x > 0 && y > 0) + incoming += SB_DIAGONAL * + overflow_at(source, (row - size_width + x - 1) * 3 + c, limit); + if (x + 1 < width && y > 0) + incoming += SB_DIAGONAL * + overflow_at(source, (row - size_width + x + 1) * 3 + c, limit); + if (x > 0 && y + 1 < height) + incoming += SB_DIAGONAL * + overflow_at(source, (row + size_width + x - 1) * 3 + c, limit); + if (x + 1 < width && y + 1 < height) + incoming += SB_DIAGONAL * + overflow_at(source, (row + size_width + x + 1) * 3 + c, limit); + + const double next = (value < limit ? value : limit) + transfer * incoming; + destination[index] = next; + const double next_overflow = next - limit; + if (next_overflow > maximum) + maximum = next_overflow; + } + } + row_absorbed[y] = absorbed; + row_boundary[y] = boundary; + row_max[y] = maximum; + } + + /* Row-ordered reduction keeps the reported totals independent of the + * worker count. Pixel values are already deterministic because each + * target reads its eight source neighbours in a fixed order. */ + double round_absorbed = 0.0, round_boundary = 0.0, next_max = 0.0; + for (int y = y0; y <= y1; ++y) { + round_absorbed += row_absorbed[y]; + round_boundary += row_boundary[y]; + if (row_max[y] > next_max) + next_max = row_max[y]; + } + stats->absorbed_signal += round_absorbed; + stats->boundary_loss += round_boundary; + ++iterations; + + double *swap_source = destination; + destination = source; + source = swap_source; + + if (next_max <= tolerance) + break; + } + stats->iterations = iterations; + + /* Residual values still inside (E, E + tolerance] are snapped to E so every + * successfully returned channel is at most the response limit. */ + for (int y = final_top; y <= final_bottom; ++y) { + const size_t row = (size_t)y * size_width; + for (int x = final_left; x <= final_right; ++x) { + const size_t base = (row + (size_t)x) * 3; + for (int c = 0; c < 3; ++c) { + const size_t index = base + (size_t)c; + double *value = &source[index]; + if (*value > limit) { + stats->residual_clamp_loss += *value - limit; + *value = limit; + ++stats->final_clamped_channels; + } + } + } + } + + /* `source` now holds the final state; the other framebuffer may still hold an + * older state inside the final box. Copy the box back into the caller's + * buffer, whose pixels outside the box were never modified. */ + if (source != hdr) { + for (int y = final_top; y <= final_bottom; ++y) { + const size_t row = (size_t)y * size_width; + memcpy(hdr + (row + (size_t)final_left) * 3, + source + (row + (size_t)final_left) * 3, + (size_t)(final_right - final_left + 1) * 3 * sizeof *hdr); + } + } + + const double box_width = (double)(final_right - final_left + 1); + const double box_height = (double)(final_bottom - final_top + 1); + stats->bbox_coverage = (box_width * box_height) / (double)pixels; + + free(scratch); + free(row_accumulator); + stats->elapsed_seconds = omp_get_wtime() - start; + return 0; +} diff --git a/src/sensor_bloom.h b/src/sensor_bloom.h new file mode 100644 index 0000000..7be4252 --- /dev/null +++ b/src/sensor_bloom.h @@ -0,0 +1,45 @@ +#ifndef SENSOR_BLOOM_H +#define SENSOR_BLOOM_H + +#include + +/* Optional post-PSF three-channel sensor saturation/overflow model. + * + * Each RGB channel is processed independently. A channel value above the + * finite response limit E contributes an overflow D = max(H - E, 0); a fixed + * fraction `transfer` of that overflow is spread to the eight neighbours with + * a 9-point stencil (4 axial weights 4/20, 4 diagonal weights 1/20), while the + * remainder (1 - transfer) * D is absorbed. Overflow directed off the image + * is lost at the boundary. The response is isotropic, non-conservative, and a + * phenomenological approximation, not a specific CCD/CMOS/Bayer structure. */ +typedef struct { + double response_limit; /* E, in post-exposure linear HDR renderer scale */ + double transfer; /* e, in [0, 1): per-round retained transfer ratio */ +} SensorBloomSettings; + +typedef struct { + size_t initially_saturated_channels; + size_t final_clamped_channels; + size_t iterations; + size_t predicted_iterations; + double bbox_coverage; + double peak_input; + double initial_max_overflow; + double absorbed_signal; + double boundary_loss; + double residual_clamp_loss; + double elapsed_seconds; +} SensorBloomStats; + +/* Applies the model in place to an interleaved RGB HDR framebuffer of + * width * height * 3 doubles. Returns 0 on success (including the trivial + * case of no channel above E, where the buffer is left byte-for-byte + * unchanged and no scratch is allocated), and -1 for invalid settings, + * dimensions, allocation failure, a non-finite input sample, or a converged + * round bound above the internal hard limit. On failure the buffer is not + * modified. `stats` may be NULL. */ +int sensor_bloom_apply(double *hdr, int width, int height, + const SensorBloomSettings *settings, + SensorBloomStats *stats); + +#endif diff --git a/tests/benchmark_sensor_bloom.c b/tests/benchmark_sensor_bloom.c new file mode 100644 index 0000000..1157d0e --- /dev/null +++ b/tests/benchmark_sensor_bloom.c @@ -0,0 +1,121 @@ +/* Isolated benchmark for the sensor-bloom model in src/sensor_bloom.c. + * + * This performs no ray tracing, catalog lookup, or PSF splatting; it fills a + * synthetic HDR framebuffer and times the post-processing model alone. The + * default 512x288 size keeps a full sweep bounded, and the optional size + * arguments are clamped so the benchmark cannot be turned into a production + * render. It is a non-default target and is not part of `make test`. */ +#include "sensor_bloom.h" + +#include +#include +#include +#include +#include +#include + +#define DEFAULT_WIDTH 512 +#define DEFAULT_HEIGHT 288 +#define MAX_DIMENSION 1024 +#define RESPONSE_LIMIT 1.0 + +static uint32_t next_random(uint32_t *state) { + *state = *state * 1664525u + 1013904223u; + return *state; +} + +static uint64_t checksum(const double *hdr, size_t count) { + const unsigned char *bytes = (const unsigned char *)hdr; + uint64_t hash = 1469598103934665603ull; + for (size_t i = 0; i < count * sizeof(double); ++i) { + hash ^= bytes[i]; + hash *= 1099511628211ull; + } + return hash; +} + +static void fill_central(double *hdr, int width, int height) { + const size_t index = ((size_t)(height / 2) * width + width / 2) * 3; + hdr[index + 0] = 20.0; + hdr[index + 1] = 20.0; + hdr[index + 2] = 20.0; +} + +static void fill_scattered(double *hdr, int width, int height) { + int placed = 0; + for (int y = 0; y < height && placed < 64; y += height > 8 ? height / 8 : 1) { + for (int x = 0; x < width && placed < 64; x += width > 8 ? width / 8 : 1) { + const size_t index = ((size_t)y * width + x) * 3; + hdr[index + 0] = 20.0; + hdr[index + 1] = 16.0; + hdr[index + 2] = 12.0; + ++placed; + } + } +} + +static void fill_random_fraction(double *hdr, int width, int height) { + const size_t count = (size_t)width * height * 3; + uint32_t state = 20260927u; + for (size_t i = 0; i < count; ++i) + hdr[i] = 1.5 + (double)(next_random(&state) % 2500u) / 1000.0; + for (size_t i = 0; i < count; ++i) + if (next_random(&state) % 100u >= 5u) + hdr[i] = 0.0; +} + +typedef void (*FillFunction)(double *, int, int); + +static void run_case(const char *name, FillFunction fill, int width, int height, + double transfer) { + const size_t count = (size_t)width * height * 3; + double *hdr = calloc(count, sizeof *hdr); + if (hdr == NULL) { + fprintf(stderr, "allocation failed for %s\n", name); + exit(EXIT_FAILURE); + } + fill(hdr, width, height); + + const SensorBloomSettings settings = {RESPONSE_LIMIT, transfer}; + SensorBloomStats stats; + if (sensor_bloom_apply(hdr, width, height, &settings, &stats)) { + fprintf(stderr, "%s e=%.2f failed\n", name, transfer); + free(hdr); + exit(EXIT_FAILURE); + } + /* Report the model's own guard-inclusive conservative cap, not a duplicate + * of the bound formula. */ + printf("%-10s e=%.2f size=%dx%d saturated=%zu predicted=%zu actual=%zu " + "bbox=%.4f elapsed=%.6fs checksum=%016llx\n", + name, transfer, width, height, stats.initially_saturated_channels, + stats.predicted_iterations, stats.iterations, stats.bbox_coverage, + stats.elapsed_seconds, (unsigned long long)checksum(hdr, count)); + free(hdr); +} + +int main(int argc, char **argv) { + int width = DEFAULT_WIDTH, height = DEFAULT_HEIGHT; + if (argc >= 2) + width = atoi(argv[1]); + if (argc >= 3) + height = atoi(argv[2]); + if (width < 1) + width = 1; + if (height < 1) + height = 1; + if (width > MAX_DIMENSION) + width = MAX_DIMENSION; + if (height > MAX_DIMENSION) + height = MAX_DIMENSION; + + const int workers = omp_get_max_threads(); + printf("sensor-bloom benchmark: size=%dx%d threads=%d limit=%.3f\n", + width, height, workers, RESPONSE_LIMIT); + const double transfers[] = {0.0, 0.5, 0.9}; + for (size_t e = 0; e < sizeof transfers / sizeof transfers[0]; ++e) { + run_case("central", fill_central, width, height, transfers[e]); + run_case("scattered", fill_scattered, width, height, transfers[e]); + run_case("random5pct", fill_random_fraction, width, height, transfers[e]); + } + return EXIT_SUCCESS; +} diff --git a/tests/test_camera_cli.py b/tests/test_camera_cli.py index 503024c..616a5a3 100644 --- a/tests/test_camera_cli.py +++ b/tests/test_camera_cli.py @@ -8,6 +8,28 @@ import sys import tempfile import zlib + +def fits_max(path): + """Largest positive sample in the renderer's three-plane float FITS.""" + data = path.read_bytes() + cards, offset = [], 0 + while True: + block = data[offset:offset + 2880] + assert len(block) == 2880, f'truncated FITS header: {path}' + offset += 2880 + cards.extend(block[i:i + 80] for i in range(0, 2880, 80)) + if any(card.startswith(b'END') for card in cards[-36:]): + break + values = {} + for card in cards: + if card[8:10] == b'= ': + values[card[:8].decode().strip()] = card[10:30].decode().strip() + shape = tuple(int(values[f'NAXIS{axis}']) for axis in (1, 2, 3)) + count = shape[0] * shape[1] * shape[2] + payload = data[offset:offset + count * 4] + assert len(payload) == count * 4, f'truncated FITS payload: {path}' + return max(struct.unpack(f'>{count}f', payload)) + BUILD = Path(sys.argv[1] if len(sys.argv) > 1 else 'build/Release').resolve() TESTDIR = Path(sys.argv[2]).resolve() if len(sys.argv) > 2 else BUILD ENV = dict(os.environ, OMP_NUM_THREADS='4') @@ -89,6 +111,8 @@ with tempfile.TemporaryDirectory(prefix='gr-camera-cli-') as directory: for option in ('--observer-position', '--observer-velocity', '--camera-roll-deg'): assert option in help_text assert '--tone-map' in help_text and '--tone-map-p' in help_text + assert '--sensor-bloom-limit' in help_text + assert '--sensor-bloom-transfer' in help_text assert '--observer-inward-speed' not in help_text assert '_mesh.' in help_text, help_text # The synthetic grid is calibrated for the renderer's default exposure. @@ -205,6 +229,17 @@ with tempfile.TemporaryDirectory(prefix='gr-camera-cli-') as directory: (['--tone-map-p', 'nan'], None), (['--tone-map-p', 'inf'], None), (['--tone-map', 'reinhard', '--tone-map-p', 2], 'applies only'), + (['--sensor-bloom-limit', 1], 'specified together'), + (['--sensor-bloom-transfer', 0.5], 'specified together'), + (['--sensor-bloom-limit', 0], None), + (['--sensor-bloom-limit', -1], None), + (['--sensor-bloom-limit', 'nan'], None), + (['--sensor-bloom-limit', 'inf'], None), + (['--sensor-bloom-limit', 1, '--sensor-bloom-transfer', -0.1], None), + (['--sensor-bloom-limit', 1, '--sensor-bloom-transfer', 1], None), + (['--sensor-bloom-limit', 1, '--sensor-bloom-transfer', 1.5], None), + (['--sensor-bloom-limit', 1, '--sensor-bloom-transfer', 'nan'], None), + (['--sensor-bloom-limit', 1, '--sensor-bloom-transfer', 'inf'], None), ] if backend == 'schwarzschild': errors += [(['--observer-position', 1.5, 0, 0, '--observer-velocity', -0.5, 0, 0], 'capture cutoff'), @@ -279,6 +314,60 @@ with tempfile.TemporaryDirectory(prefix='gr-camera-cli-') as directory: assert image_payload(hdr_mesh_output) == baseline assert image_payload(tmp / f'{backend}_hdr_mesh_mesh.{ext}') != baseline + # Sensor-bloom integration runs in every build. With linear HDR the + # response limit is derived from the rendered peak; without it the + # fixture's calibrated exposure puts the display shoulder near 1.0, + # which is guaranteed to saturate this scene, so the default ENABLE_HDR=0 + # suite still exercises the output-pipeline hook. The limit is placed + # well below the peak so the clamped region falls in the tone map's + # sensitive range and the display comparison stays discriminating. + if hdr_available: + base_peak = fits_max(base_fits) + assert base_peak > 0.0 + bloom_limit = base_peak / 100.0 + else: + bloom_limit = 1.0 + hdr_args = ['--hdr-output'] if hdr_available else [] + bloom_output = tmp / f'{backend}_bloom.{ext}' + bloom_run = run(binary, *common, *hdr_args, + '--sensor-bloom-limit', bloom_limit, + '--sensor-bloom-transfer', 0.5, + '--output', bloom_output) + assert image_payload(bloom_output) != baseline + assert 'Sensor bloom:' in bloom_run.stderr, bloom_run.stderr + report = bloom_run.stderr.split('Sensor bloom:', 1)[1].splitlines()[0] + fields = dict(token.split('=', 1) for token in report.split() if '=' in token) + assert int(fields['saturated']) > 0, report + assert int(fields['iterations'].split('/')[0]) >= 1, report + assert float(fields['peak']) >= bloom_limit, report + if hdr_available: + # The raw FITS is written before the model runs, so it stays + # byte-identical to the baseline even though the display changes. + bloom_fits = tmp / f'{backend}_bloom_HDR.fits' + assert bloom_fits.exists(), bloom_run.stderr + diff = subprocess.run([sys.executable, str(FITSDIFF), str(base_fits), + str(bloom_fits)], capture_output=True, text=True) + assert diff.returncode == 0, diff.stdout + diff.stderr + assert 'mismatches=0 max_abs=0 max_rel=0' in diff.stdout, diff.stdout + + # The mesh overlay is drawn on the already-bloomed frame, so the primary + # image is unchanged by --draw-mesh and the diagnostic lines do not feed + # back into the overflow model. + bloom_mesh_output = tmp / f'{backend}_bloom_mesh.{ext}' + run(binary, *common, *hdr_args, '--sensor-bloom-limit', bloom_limit, + '--sensor-bloom-transfer', 0.5, '--draw-mesh', + '--output', bloom_mesh_output) + assert image_payload(bloom_mesh_output) == image_payload(bloom_output) + bloom_mesh_sibling = tmp / f'{backend}_bloom_mesh_mesh.{ext}' + assert bloom_mesh_sibling.exists() + assert image_payload(bloom_mesh_sibling) != image_payload(bloom_output) + if hdr_available: + bloom_mesh_fits = tmp / f'{backend}_bloom_mesh_HDR.fits' + diff = subprocess.run([sys.executable, str(FITSDIFF), str(base_fits), + str(bloom_mesh_fits)], capture_output=True, text=True) + assert diff.returncode == 0, diff.stdout + diff.stderr + assert 'mismatches=0 max_abs=0 max_rel=0' in diff.stdout, diff.stdout + if backend == 'schwarzschild': # Two inward-looking free-fall samples at r=6.2696 and r=3.1593. # At 16:9 the latter frame finishes in generation 0, while the diff --git a/tests/test_sensor_bloom.c b/tests/test_sensor_bloom.c new file mode 100644 index 0000000..101e96e --- /dev/null +++ b/tests/test_sensor_bloom.c @@ -0,0 +1,568 @@ +/* Regression tests for the standalone sensor-bloom model in + * src/sensor_bloom.c. + * + * These link only the model: no catalog, ray tracing, image writer, FFTW, or + * GPU is involved, so they build in every ENABLE_HDR/PSF_BACKEND + * configuration. They call the production entry point and compare it against + * an independently written dense synchronous reference where a closed form is + * not simpler. */ +#include "sensor_bloom.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +/* The scratch malloc-failure path is covered by the ordinary non-ASan run of + * `allocation_failure_is_reported()` below. AddressSanitizer intercepts + * allocation and its runtime cannot mmap under a forced RLIMIT_AS, so that + * probe is skipped in ASan builds rather than being made to pass; the + * dimension-overflow validation is a separate check and does not substitute + * for it. */ +#if defined(__SANITIZE_ADDRESS__) +#define SB_HAVE_ASAN 1 +#elif defined(__has_feature) +#if __has_feature(address_sanitizer) +#define SB_HAVE_ASAN 1 +#endif +#endif + +static int failures = 0; + +static void check(int condition, const char *message) { + if (!condition) { + fprintf(stderr, "FAIL: %s\n", message); + ++failures; + } +} + +static int close_to(double value, double expected, double tolerance) { + return fabs(value - expected) <= tolerance; +} + +static double *alloc_rgb(int width, int height) { + return calloc((size_t)width * height * 3, sizeof(double)); +} + +static size_t at(int width, int x, int y, int c) { + return ((size_t)y * width + x) * 3 + (size_t)c; +} + +static double sample(const double *rgb, int width, int x, int y, int c) { + return rgb[at(width, x, y, c)]; +} + +static void set(double *rgb, int width, int x, int y, int c, double value) { + rgb[at(width, x, y, c)] = value; +} + +static int all_finite_and_bounded(const double *rgb, size_t count, + double limit) { + for (size_t i = 0; i < count; ++i) + if (!isfinite(rgb[i]) || rgb[i] > limit) + return 0; + return 1; +} + +/* Independently written dense synchronous reference: every round reads the + * whole previous state and writes the whole next state, with the same 9-point + * weights. Only used to cross-check cascades. */ +static double reference_overflow(const double *src, size_t index, double limit) { + const double difference = src[index] - limit; + return difference > 0.0 ? difference : 0.0; +} + +static void reference_bloom(double *hdr, int width, int height, double limit, + double transfer) { + const double axial = 4.0 / 20.0, diagonal = 1.0 / 20.0; + const double tolerance = 1e-9 * limit; + const size_t count = (size_t)width * height * 3; + double *first = malloc(count * sizeof *first); + double *second = malloc(count * sizeof *second); + memcpy(first, hdr, count * sizeof *first); + memcpy(second, hdr, count * sizeof *second); + double *src = first, *dst = second; + for (int round = 0; round < 20000; ++round) { + double next_max = 0.0; + for (int y = 0; y < height; ++y) { + for (int x = 0; x < width; ++x) { + for (int c = 0; c < 3; ++c) { + const size_t index = at(width, x, y, c); + const double value = src[index]; + double incoming = 0.0; + if (x > 0) + incoming += axial * reference_overflow(src, at(width, x - 1, y, c), limit); + if (x + 1 < width) + incoming += axial * reference_overflow(src, at(width, x + 1, y, c), limit); + if (y > 0) + incoming += axial * reference_overflow(src, at(width, x, y - 1, c), limit); + if (y + 1 < height) + incoming += axial * reference_overflow(src, at(width, x, y + 1, c), limit); + if (x > 0 && y > 0) + incoming += diagonal * reference_overflow(src, at(width, x - 1, y - 1, c), limit); + if (x + 1 < width && y > 0) + incoming += diagonal * reference_overflow(src, at(width, x + 1, y - 1, c), limit); + if (x > 0 && y + 1 < height) + incoming += diagonal * reference_overflow(src, at(width, x - 1, y + 1, c), limit); + if (x + 1 < width && y + 1 < height) + incoming += diagonal * reference_overflow(src, at(width, x + 1, y + 1, c), limit); + const double next = (value < limit ? value : limit) + transfer * incoming; + dst[index] = next; + const double overflow = next - limit; + if (overflow > next_max) + next_max = overflow; + } + } + } + double *swap = src; + src = dst; + dst = swap; + if (next_max <= tolerance) + break; + } + for (size_t i = 0; i < count; ++i) + if (src[i] > limit) + src[i] = limit; + memcpy(hdr, src, count * sizeof *hdr); + free(first); + free(second); +} + +static void test_no_saturation(void) { + const int width = 8, height = 6; + const size_t count = (size_t)width * height * 3; + double *hdr = alloc_rgb(width, height); + for (size_t i = 0; i < count; ++i) + hdr[i] = 0.25 * (double)(i % 5); + double *original = malloc(count * sizeof *original); + memcpy(original, hdr, count * sizeof *original); + + const SensorBloomSettings settings = {1.0, 0.5}; + SensorBloomStats stats; + check(sensor_bloom_apply(hdr, width, height, &settings, &stats) == 0, + "no-saturation apply succeeds"); + check(memcmp(hdr, original, count * sizeof *hdr) == 0, + "no-saturation buffer is byte-for-byte unchanged"); + check(stats.iterations == 0, "no-saturation uses zero iterations"); + check(stats.predicted_iterations == 0, + "no-saturation predicts zero iterations"); + check(stats.initially_saturated_channels == 0, + "no-saturation reports no saturated channel"); + check(stats.final_clamped_channels == 0, + "no-saturation clamps nothing"); + free(original); + free(hdr); +} + +static void test_zero_transfer(void) { + const int width = 5, height = 5; + double *hdr = alloc_rgb(width, height); + set(hdr, width, 2, 2, 0, 1.7); + const SensorBloomSettings settings = {1.0, 0.0}; + SensorBloomStats stats; + check(sensor_bloom_apply(hdr, width, height, &settings, &stats) == 0, + "zero-transfer apply succeeds"); + check(close_to(sample(hdr, width, 2, 2, 0), 1.0, 1e-15), + "zero-transfer clamps the saturated channel to E"); + check(sample(hdr, width, 1, 2, 0) == 0.0 && sample(hdr, width, 3, 2, 0) == 0.0, + "zero-transfer leaves axial neighbours unchanged"); + check(sample(hdr, width, 1, 1, 0) == 0.0 && sample(hdr, width, 3, 3, 0) == 0.0, + "zero-transfer leaves diagonal neighbours unchanged"); + check(close_to(stats.absorbed_signal, 0.7, 1e-15), + "zero-transfer absorbs the full overflow"); + check(stats.boundary_loss == 0.0, "zero-transfer loses nothing at the boundary"); + check(stats.iterations >= 1, "zero-transfer iterates at least once"); + free(hdr); +} + +static void test_single_impulse_stencil(void) { + const int width = 5, height = 5; + const double limit = 1.0, transfer = 0.5; + double *hdr = alloc_rgb(width, height); + set(hdr, width, 2, 2, 0, 1.8); + set(hdr, width, 2, 2, 1, 1.5); + set(hdr, width, 2, 2, 2, 1.2); + const SensorBloomSettings settings = {limit, transfer}; + SensorBloomStats stats; + check(sensor_bloom_apply(hdr, width, height, &settings, &stats) == 0, + "single-impulse apply succeeds"); + check(sample(hdr, width, 2, 2, 0) == limit && + sample(hdr, width, 2, 2, 1) == limit && + sample(hdr, width, 2, 2, 2) == limit, + "single-impulse centre sits exactly at E"); + const double red_axis = transfer * 0.8 * (4.0 / 20.0); + const double red_diagonal = transfer * 0.8 * (1.0 / 20.0); + const double green_axis = transfer * 0.5 * (4.0 / 20.0); + const double blue_axis = transfer * 0.2 * (4.0 / 20.0); + check(close_to(sample(hdr, width, 1, 2, 0), red_axis, 1e-15) && + close_to(sample(hdr, width, 3, 2, 0), red_axis, 1e-15) && + close_to(sample(hdr, width, 2, 1, 0), red_axis, 1e-15) && + close_to(sample(hdr, width, 2, 3, 0), red_axis, 1e-15), + "axial neighbours receive e * D * 4/20"); + check(close_to(sample(hdr, width, 1, 1, 0), red_diagonal, 1e-15) && + close_to(sample(hdr, width, 3, 3, 0), red_diagonal, 1e-15), + "diagonal neighbours receive e * D * 1/20"); + check(close_to(sample(hdr, width, 1, 2, 1), green_axis, 1e-15) && + close_to(sample(hdr, width, 1, 2, 2), blue_axis, 1e-15), + "RGB channels transfer independently"); + check(sample(hdr, width, 0, 0, 0) == 0.0 && sample(hdr, width, 4, 4, 2) == 0.0, + "far pixels stay untouched"); + free(hdr); +} + +static void test_boundary_loss(void) { + const int width = 5, height = 5; + const double limit = 1.0, transfer = 0.5; + double *hdr = alloc_rgb(width, height); + set(hdr, width, 0, 0, 0, 1.8); + const SensorBloomSettings settings = {limit, transfer}; + SensorBloomStats stats; + check(sensor_bloom_apply(hdr, width, height, &settings, &stats) == 0, + "corner apply succeeds"); + const double overflow = 0.8; + check(close_to(sample(hdr, width, 1, 0, 0), transfer * overflow * 0.2, 1e-15) && + close_to(sample(hdr, width, 0, 1, 0), transfer * overflow * 0.2, 1e-15), + "in-image axial propagation at a corner uses the fixed weight"); + check(close_to(sample(hdr, width, 1, 1, 0), transfer * overflow * 0.05, 1e-15), + "in-image diagonal propagation at a corner uses the fixed weight"); + check(close_to(stats.boundary_loss, + transfer * overflow * (2.0 * 0.2 + 3.0 * 0.05), 1e-15), + "off-image weight is counted as boundary loss without renormalization"); + check(close_to(stats.absorbed_signal, (1.0 - transfer) * overflow, 1e-15), + "absorbed signal is (1 - e) * D"); + free(hdr); +} + +static void test_cascade_against_reference(void) { + const int width = 9, height = 9; + const double limit = 0.5, transfer = 0.6; + const size_t count = (size_t)width * height * 3; + double *hdr = alloc_rgb(width, height); + set(hdr, width, 4, 4, 0, 10.0); + set(hdr, width, 4, 4, 1, 4.0); + set(hdr, width, 2, 6, 2, 3.0); + set(hdr, width, 7, 1, 0, 1.5); + double *expected = malloc(count * sizeof *expected); + memcpy(expected, hdr, count * sizeof *expected); + reference_bloom(expected, width, height, limit, transfer); + + const SensorBloomSettings settings = {limit, transfer}; + SensorBloomStats stats; + check(sensor_bloom_apply(hdr, width, height, &settings, &stats) == 0, + "cascade apply succeeds"); + check(stats.iterations >= 2, "cascade needs at least two rounds"); + double max_difference = 0.0; + for (size_t i = 0; i < count; ++i) + max_difference = fmax(max_difference, fabs(hdr[i] - expected[i])); + check(max_difference < 1e-12, "cascade matches the dense synchronous reference"); + check(all_finite_and_bounded(hdr, count, limit), + "cascade output is finite and at most E"); + free(expected); + free(hdr); +} + +static void rotate_90(const double *src, double *dst, int size) { + for (int y = 0; y < size; ++y) + for (int x = 0; x < size; ++x) + for (int c = 0; c < 3; ++c) + set(dst, size, size - 1 - y, x, c, sample(src, size, x, y, c)); +} + +static void test_symmetry(void) { + const int size = 9; + const double limit = 1.0, transfer = 0.5; + double *hdr = alloc_rgb(size, size); + set(hdr, size, 4, 4, 0, 10.0); + const SensorBloomSettings settings = {limit, transfer}; + SensorBloomStats stats; + check(sensor_bloom_apply(hdr, size, size, &settings, &stats) == 0, + "symmetry apply succeeds"); + check(sample(hdr, size, 3, 4, 0) == sample(hdr, size, 5, 4, 0) && + sample(hdr, size, 3, 4, 0) == sample(hdr, size, 4, 3, 0) && + sample(hdr, size, 3, 4, 0) == sample(hdr, size, 4, 5, 0), + "four axial neighbours of a centred point are equal"); + check(sample(hdr, size, 3, 3, 0) == sample(hdr, size, 5, 3, 0) && + sample(hdr, size, 3, 3, 0) == sample(hdr, size, 3, 5, 0) && + sample(hdr, size, 3, 3, 0) == sample(hdr, size, 5, 5, 0), + "four diagonal neighbours of a centred point are equal"); + + /* 90-degree rotational covariance: rotating the input then blooming must + * equal blooming then rotating the output. This starts from a fresh, + * asymmetric and still-saturated input; the already-bloomed centre-impulse + * buffer above has no overflow and would only exercise the no-op path. */ + const size_t count = (size_t)size * size * 3; + double *input = alloc_rgb(size, size); + set(input, size, 2, 5, 0, 30.0); + set(input, size, 6, 3, 1, 8.0); + set(input, size, 4, 1, 2, 5.0); + set(input, size, 7, 7, 0, 4.0); + double *bloomed_input = malloc(count * sizeof *bloomed_input); + double *rotated_input = alloc_rgb(size, size); + rotate_90(input, rotated_input, size); + double *bloomed_rotated = malloc(count * sizeof *bloomed_rotated); + memcpy(bloomed_input, input, count * sizeof *bloomed_input); + memcpy(bloomed_rotated, rotated_input, count * sizeof *bloomed_rotated); + check(sensor_bloom_apply(bloomed_input, size, size, &settings, &stats) == 0 && + stats.initially_saturated_channels > 0, + "covariance input is saturated and blooms"); + check(sensor_bloom_apply(bloomed_rotated, size, size, &settings, &stats) == 0, + "rotated covariance input blooms"); + double *expected = alloc_rgb(size, size); + rotate_90(bloomed_input, expected, size); + double max_difference = 0.0; + for (size_t i = 0; i < count; ++i) + max_difference = fmax(max_difference, fabs(bloomed_rotated[i] - expected[i])); + check(max_difference < 1e-9, "model is 90-degree rotation covariant"); + free(expected); + free(bloomed_rotated); + free(rotated_input); + free(bloomed_input); + free(input); + free(hdr); +} + +static uint32_t next_random(uint32_t *state) { + *state = *state * 1664525u + 1013904223u; + return *state; +} + +static void test_thread_determinism(void) { + const int width = 40, height = 30; + const size_t count = (size_t)width * height * 3; + double *hdr = alloc_rgb(width, height); + uint32_t state = 12345u; + for (size_t i = 0; i < count; ++i) + hdr[i] = (double)(next_random(&state) % 3000u) / 1000.0; + double *first = malloc(count * sizeof *first); + double *second = malloc(count * sizeof *second); + memcpy(first, hdr, count * sizeof *first); + memcpy(second, hdr, count * sizeof *second); + + const SensorBloomSettings settings = {1.5, 0.5}; + SensorBloomStats stats_one, stats_four; + omp_set_num_threads(1); + check(sensor_bloom_apply(first, width, height, &settings, &stats_one) == 0, + "single-thread apply succeeds"); + omp_set_num_threads(4); + check(sensor_bloom_apply(second, width, height, &settings, &stats_four) == 0, + "four-thread apply succeeds"); + + check(memcmp(first, second, count * sizeof *first) == 0, + "final HDR is byte-identical for 1 and 4 threads"); + check(stats_one.absorbed_signal == stats_four.absorbed_signal && + stats_one.boundary_loss == stats_four.boundary_loss && + stats_one.iterations == stats_four.iterations && + stats_one.residual_clamp_loss == stats_four.residual_clamp_loss, + "reported statistics are identical for 1 and 4 threads"); + free(first); + free(second); + free(hdr); +} + +static void test_convergence_and_hard_limit(void) { + const int width = 9, height = 9; + const size_t count = (size_t)width * height * 3; + double *hdr = alloc_rgb(width, height); + set(hdr, width, 4, 4, 0, 100.0); + const SensorBloomSettings settings = {1.0, 0.5}; + SensorBloomStats stats; + check(sensor_bloom_apply(hdr, width, height, &settings, &stats) == 0, + "converging apply succeeds"); + check(all_finite_and_bounded(hdr, count, 1.0), + "converged output is finite and at most E"); + check(stats.iterations <= stats.predicted_iterations, + "actual iterations do not exceed the conservative bound"); + + /* A round bound just inside the hard limit must still converge. The reported + * cap includes the fp guard but must never exceed 4096, and it must bound the + * actual number of rounds. */ + const int near_width = 5, near_height = 5; + const size_t near_count = (size_t)near_width * near_height * 3; + double *near = alloc_rgb(near_width, near_height); + set(near, near_width, 2, 2, 0, 1.0 + 7.4e-9); + const SensorBloomSettings near_settings = {1.0, 0.9995}; + SensorBloomStats near_stats; + check(sensor_bloom_apply(near, near_width, near_height, &near_settings, + &near_stats) == 0, + "near-limit round bound still converges"); + check(near_stats.predicted_iterations <= (size_t)4096, + "reported round cap never exceeds the 4096 hard limit"); + check(near_stats.iterations <= near_stats.predicted_iterations, + "actual rounds never exceed the reported guard-inclusive cap"); + check(all_finite_and_bounded(near, near_count, 1.0), + "near-limit output is finite and at most E"); + free(near); + + /* A bound one round above the hard limit (4097 with these parameters) must + * fail in the pre-check, before any HDR sample is touched. */ + double *near_fail = alloc_rgb(near_width, near_height); + set(near_fail, near_width, 2, 2, 0, 1.0 + 7.758e-9); + double *near_original = malloc(near_count * sizeof *near_original); + memcpy(near_original, near_fail, near_count * sizeof *near_original); + const SensorBloomSettings near_hard = {1.0, 0.9995}; + SensorBloomStats near_hard_stats; + check(sensor_bloom_apply(near_fail, near_width, near_height, &near_hard, + &near_hard_stats) == -1, + "round bound just above the hard limit fails"); + check(memcmp(near_fail, near_original, near_count * sizeof *near_fail) == 0, + "near-limit failure leaves the buffer byte-for-byte unchanged"); + free(near_original); + free(near_fail); + + /* An unbounded-overflow parameter set (e = 1 - 1e-6) needs far more than the + * 4096-round hard limit and must fail before modifying the input. */ + double *failing = alloc_rgb(width, height); + set(failing, width, 4, 4, 0, 2.0); + double *original = malloc(count * sizeof *original); + memcpy(original, failing, count * sizeof *original); + const SensorBloomSettings hard = {1.0, 0.999999}; + SensorBloomStats hard_stats; + check(sensor_bloom_apply(failing, width, height, &hard, &hard_stats) == -1, + "round bound above the hard limit fails"); + check(memcmp(failing, original, count * sizeof *failing) == 0, + "hard-limit failure leaves the buffer byte-for-byte unchanged"); + check(hard_stats.iterations == 0, + "hard-limit failure reports no iterations"); + check(hard_stats.predicted_iterations == 0 && + hard_stats.final_clamped_channels == 0, + "hard-limit failure reports no predicted rounds or clamps"); + free(original); + free(failing); + free(hdr); +} + +static void test_input_validation(void) { + const int width = 4, height = 3; + const size_t count = (size_t)width * height * 3; + double *hdr = alloc_rgb(width, height); + double *original = malloc(count * sizeof *original); + const SensorBloomSettings valid = {1.0, 0.5}; + SensorBloomStats stats; + + check(sensor_bloom_apply(NULL, width, height, &valid, &stats) == -1, + "NULL framebuffer is rejected"); + check(sensor_bloom_apply(hdr, width, height, NULL, &stats) == -1, + "NULL settings are rejected"); + check(sensor_bloom_apply(hdr, 0, height, &valid, &stats) == -1, + "zero width is rejected"); + check(sensor_bloom_apply(hdr, width, -1, &valid, &stats) == -1, + "negative height is rejected"); + check(sensor_bloom_apply(hdr, width, height, &valid, NULL) == 0, + "NULL stats are permitted for an unsaturated buffer"); + + double dummy = 0.0; + check(sensor_bloom_apply(&dummy, INT32_MAX, INT32_MAX, &valid, &stats) == -1, + "multiplication-overflowing dimensions are rejected before scanning"); + + const double bad_limits[] = {0.0, -1.0, NAN, INFINITY, -INFINITY}; + for (size_t i = 0; i < sizeof bad_limits / sizeof bad_limits[0]; ++i) { + const SensorBloomSettings bad = {bad_limits[i], 0.5}; + check(sensor_bloom_apply(hdr, width, height, &bad, &stats) == -1, + "invalid response limit is rejected"); + } + const double bad_transfers[] = {-0.1, 1.0, 1.5, NAN, INFINITY, -INFINITY}; + for (size_t i = 0; i < sizeof bad_transfers / sizeof bad_transfers[0]; ++i) { + const SensorBloomSettings bad = {1.0, bad_transfers[i]}; + check(sensor_bloom_apply(hdr, width, height, &bad, &stats) == -1, + "invalid transfer is rejected"); + } + + const double non_finite[] = {NAN, INFINITY, -INFINITY}; + for (size_t i = 0; i < sizeof non_finite / sizeof non_finite[0]; ++i) { + for (size_t j = 0; j < count; ++j) + hdr[j] = 0.5; + hdr[count / 2] = non_finite[i]; + memcpy(original, hdr, count * sizeof *original); + check(sensor_bloom_apply(hdr, width, height, &valid, &stats) == -1, + "non-finite HDR sample is rejected"); + check(memcmp(hdr, original, count * sizeof *hdr) == 0, + "non-finite rejection leaves the buffer unchanged"); + } + free(original); + free(hdr); +} + +#ifndef SB_HAVE_ASAN +/* Forces the scratch allocation to fail by lowering RLIMIT_AS in a child after + * the input buffer is already mapped, then checks that the model reports the + * allocation failure instead of touching the framebuffer. */ +static int allocation_failure_is_reported(void) { + const pid_t pid = fork(); + if (pid < 0) + return 0; + if (pid == 0) { + const int width = 64, height = 64; + const size_t count = (size_t)width * height * 3; + double *hdr = calloc(count, sizeof *hdr); + if (hdr == NULL) + _exit(2); + hdr[0] = 2.0; + struct rlimit existing; + if (getrlimit(RLIMIT_AS, &existing)) + _exit(2); + const struct rlimit tiny = {.rlim_cur = 1, .rlim_max = existing.rlim_max}; + if (setrlimit(RLIMIT_AS, &tiny)) + _exit(2); + const SensorBloomSettings settings = {1.0, 0.5}; + SensorBloomStats stats; + _exit(sensor_bloom_apply(hdr, width, height, &settings, &stats) == -1 ? 0 + : 1); + } + int status = 0; + if (waitpid(pid, &status, 0) < 0) + return 0; + return WIFEXITED(status) && WEXITSTATUS(status) == 0; +} +#endif + +static void test_nonnegative_output(void) { + const int width = 12, height = 7; + const size_t count = (size_t)width * height * 3; + double *hdr = alloc_rgb(width, height); + uint32_t state = 99u; + for (size_t i = 0; i < count; ++i) + hdr[i] = (double)(next_random(&state) % 4000u) / 1000.0; + const SensorBloomSettings settings = {1.2, 0.5}; + SensorBloomStats stats; + check(sensor_bloom_apply(hdr, width, height, &settings, &stats) == 0, + "nonnegative apply succeeds"); + int nonnegative = 1; + for (size_t i = 0; i < count; ++i) + if (!(hdr[i] >= 0.0) || !isfinite(hdr[i])) + nonnegative = 0; + check(nonnegative, "nonnegative input produces nonnegative finite output"); + check(all_finite_and_bounded(hdr, count, 1.2), + "nonnegative output respects the response limit"); + free(hdr); +} + +int main(void) { +#ifndef SB_HAVE_ASAN + /* Run the fork-based allocation-failure check before any OpenMP region has + * been created in the parent. */ + check(allocation_failure_is_reported(), + "scratch allocation failure is reported without touching the input"); +#endif + test_no_saturation(); + test_zero_transfer(); + test_single_impulse_stencil(); + test_boundary_loss(); + test_cascade_against_reference(); + test_symmetry(); + test_thread_determinism(); + test_convergence_and_hard_limit(); + test_input_validation(); + test_nonnegative_output(); + + if (failures != 0) { + fprintf(stderr, "%d sensor-bloom assertion(s) failed\n", failures); + return EXIT_FAILURE; + } + puts("sensor-bloom tests passed"); + return EXIT_SUCCESS; +} diff --git a/usage.md b/usage.md index 6b9b3a2..e5d9086 100644 --- a/usage.md +++ b/usage.md @@ -307,6 +307,33 @@ Commands written before this change that omit `--tone-map` actually used Reinhard; add `--tone-map reinhard` explicitly to reproduce their PNG/PPM output. Historical FITS/HDR benchmarks are unaffected. +### Sensor bloom (optional) + +`--sensor-bloom-limit E --sensor-bloom-transfer e` (both required together, +default disabled) enable an optional post-PSF sensor saturation model applied to +the finished linear HDR framebuffer before the display operator. Each RGB +channel is processed independently and isotropically. A channel value above the +finite response limit \(E\) contributes overflow \(D=\max(H-E,0)\). A fraction +\(e\in[0,1)\) of that overflow is spread to the eight neighbours with a fixed +9-point stencil (four axial weights \(4/20\), four diagonal weights \(1/20\)), +while the remaining \((1-e)D\) is absorbed. Overflow directed outside the image +is lost at the boundary and the stencil is not renormalized there. `e=0` clamps +every over-limit channel to \(E\) without spreading to neighbours. + +\(E\) is expressed in the post-exposure linear HDR renderer scale, so it scales +with `--exposure`. The model is a phenomenological limited-response +approximation, not a specific CCD/CMOS/Bayer structure. It is non-conservative: +signal is lost both to the drain \((1-e)D\) and to the image boundary, and \(e\) +controls both the per-round retained fraction and the effective propagation +distance. The mesh overlay is drawn after the model, so the diagnostic lines +neither bloom nor feed back into overflow propagation. Visual multi-scale bloom +is not implemented. + +When enabled, each frame prints one `Sensor bloom:` line with the initial +saturated and final clamped channel counts, the actual and conservative round +counts, the peak and maximum initial overflow, absorbed and boundary loss, +residual clamp loss, and elapsed time. + ### Fast preview mode `--fast-mode` replaces the per-event PSF splat with a two-stage approximation: @@ -413,7 +440,8 @@ ordinary output, replacing its extension with `_HDR.fits`; for example, framebuffer as a three-plane, 32-bit float FITS image. Values remain linear HDR at the renderer's arbitrary scale; no tone mapping or per-frame normalization is applied, and the `--tone-map` operator and `--tone-map-p` value never affect -this file. The ordinary +this file. The `--sensor-bloom-*` model runs only after this file is written, so +enabling bloom never changes the stored FITS. The ordinary binaries do not contain this option or writer. The FITS header describes a synthetic 8640-by-5760, 36-by-24 mm full-frame