Compare commits
2
Commits
| Author | SHA1 | Date | |
|---|---|---|---|
 wyj
|
4e34780fa9 | ||
|
wyj
|
9cd933d1f8 |
No files matched your search
@@ -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:
|
||||
|
||||
@@ -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
|
||||
|
||||
+1
-1
@@ -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` 查看完整选项。
|
||||
|
||||
### 示例：紧贴史瓦西视界向外看
|
||||
|
||||
|
||||
@@ -0,0 +1,106 @@
|
||||
# Sensor bloom model benchmark (2026-09-27)
|
||||
|
||||
Isolated benchmark of the optional `src/sensor_bloom.c` limited-response model.
|
||||
It performs no ray tracing, catalog lookup, or PSF splatting: it fills a
|
||||
synthetic HDR framebuffer and times the post-processing model alone. The
|
||||
default frame is `512x288` so a full sweep stays bounded; larger optional sizes
|
||||
are clamped to `1024` and this target is not part of `make test`.
|
||||
|
||||
Environment: Gentoo, `cc` 15.3.0, `16` CPUs.
|
||||
|
||||
## Tested code
|
||||
|
||||
The model under test is committed as
|
||||
`9cd933d1f83cf768249031402cb347a159a70f2d` (the sensor-bloom implementation is
|
||||
the only change it contains). The benchmark was run from that content; the only
|
||||
untracked files at run time were this record and the pre-existing unrelated
|
||||
`nmesh_spacetime_output_format.md`.
|
||||
|
||||
```text
|
||||
$ git rev-parse HEAD
|
||||
9cd933d1f83cf768249031402cb347a159a70f2d
|
||||
|
||||
$ git status --short
|
||||
?? benchmarks/sensor_bloom_2026-09-27.md
|
||||
?? nmesh_spacetime_output_format.md
|
||||
```
|
||||
|
||||
Exact content hashes of the files that define and exercise the model at
|
||||
benchmark time (`nmesh_spacetime_output_format.md` is a pre-existing unrelated
|
||||
untracked file and was not part of this work):
|
||||
|
||||
```text
|
||||
$ sha256sum src/sensor_bloom.c src/sensor_bloom.h tests/test_sensor_bloom.c tests/benchmark_sensor_bloom.c src/main.c Makefile tests/test_camera_cli.py
|
||||
e6f5c8fd3de13756e6804f27b26c27827715667bf0ec46951b4f2bbb8eaabffd src/sensor_bloom.c
|
||||
03e8e58061cfe420630a71b7969350e2e7efecee26d96a57b5a1087a25e4a2ab src/sensor_bloom.h
|
||||
897643aba1631553bdee462a95e6c08c23b16fdf43f0b8d7a4857f369a087f14 tests/test_sensor_bloom.c
|
||||
b1f4a3ed5622ca5c11fc58c87c3aee3bb9149e3035b68ea8c7671f3ce01c1a10 tests/benchmark_sensor_bloom.c
|
||||
fdced6f075f77f17cb795c315fb0ac0590cc11789b556af3792d7ef47fed2629 src/main.c
|
||||
3852dea3c9e9b22a99d03623301b79fc4c0a62cd5808fb97f5f0abb81ee3f379 Makefile
|
||||
5fd9bae9dc1f841105b508ed6b6fce6429357b40e06adf700d5ce333b72ba92a tests/test_camera_cli.py
|
||||
```
|
||||
|
||||
## Commands
|
||||
|
||||
```sh
|
||||
make -B -j4 BUILD_TYPE=Release PSF_BACKEND=cpu sensor-bloom-bench
|
||||
OMP_NUM_THREADS=1 ./build/Release/obj/minkowski/standard_sink1_cpu/benchmark_sensor_bloom
|
||||
OMP_NUM_THREADS=4 ./build/Release/obj/minkowski/standard_sink1_cpu/benchmark_sensor_bloom
|
||||
OMP_NUM_THREADS=16 ./build/Release/obj/minkowski/standard_sink1_cpu/benchmark_sensor_bloom
|
||||
```
|
||||
|
||||
## Raw build output
|
||||
|
||||
```text
|
||||
mkdir -p build/Release/obj/minkowski/standard_sink1_cpu
|
||||
cc -DENABLE_PNG -DFRAME_PSF_EVENT_SINK=1 -DFAST_PSF_FFTW -std=c11 -march=native -pipe -Wall -Wextra -Wpedantic -O2 -DNDEBUG -fopenmp -Isrc tests/benchmark_sensor_bloom.c src/sensor_bloom.c -lm -lpng -lfftw3_omp -lfftw3 -o build/Release/obj/minkowski/standard_sink1_cpu/benchmark_sensor_bloom
|
||||
```
|
||||
|
||||
## Raw benchmark output
|
||||
|
||||
```text
|
||||
### OMP_NUM_THREADS=1
|
||||
sensor-bloom benchmark: size=512x288 threads=1 limit=1.000
|
||||
central e=0.00 size=512x288 saturated=3 predicted=1 actual=1 bbox=0.0001 elapsed=0.001640s checksum=aaec29797d030b5a
|
||||
scattered e=0.00 size=512x288 saturated=192 predicted=1 actual=1 bbox=0.7751 elapsed=0.003504s checksum=0c52973c56f09383
|
||||
random5pct e=0.00 size=512x288 saturated=22223 predicted=1 actual=1 bbox=1.0000 elapsed=0.003353s checksum=d489cdf14f756d3a
|
||||
central e=0.50 size=512x288 saturated=3 predicted=43 actual=15 bbox=0.0065 elapsed=0.001633s checksum=6625d592f76add42
|
||||
scattered e=0.50 size=512x288 saturated=192 predicted=43 actual=15 bbox=0.8433 elapsed=0.020001s checksum=44e39b09a596122f
|
||||
random5pct e=0.50 size=512x288 saturated=22223 predicted=40 actual=17 bbox=1.0000 elapsed=0.027526s checksum=d8c99ec18949278a
|
||||
central e=0.90 size=512x288 saturated=3 predicted=233 actual=43 bbox=0.0513 elapsed=0.002913s checksum=6cf4d8e02d652b62
|
||||
scattered e=0.90 size=512x288 saturated=192 predicted=233 actual=43 bbox=0.9609 elapsed=0.059650s checksum=ef91f6cf7f41249b
|
||||
random5pct e=0.90 size=512x288 saturated=22223 predicted=216 actual=47 bbox=1.0000 elapsed=0.077348s checksum=fd976dae0c08cc92
|
||||
### OMP_NUM_THREADS=4
|
||||
sensor-bloom benchmark: size=512x288 threads=4 limit=1.000
|
||||
central e=0.00 size=512x288 saturated=3 predicted=1 actual=1 bbox=0.0001 elapsed=0.001911s checksum=aaec29797d030b5a
|
||||
scattered e=0.00 size=512x288 saturated=192 predicted=1 actual=1 bbox=0.7751 elapsed=0.002795s checksum=0c52973c56f09383
|
||||
random5pct e=0.00 size=512x288 saturated=22223 predicted=1 actual=1 bbox=1.0000 elapsed=0.001979s checksum=d489cdf14f756d3a
|
||||
central e=0.50 size=512x288 saturated=3 predicted=43 actual=15 bbox=0.0065 elapsed=0.001981s checksum=6625d592f76add42
|
||||
scattered e=0.50 size=512x288 saturated=192 predicted=43 actual=15 bbox=0.8433 elapsed=0.006968s checksum=44e39b09a596122f
|
||||
random5pct e=0.50 size=512x288 saturated=22223 predicted=40 actual=17 bbox=1.0000 elapsed=0.010725s checksum=d8c99ec18949278a
|
||||
central e=0.90 size=512x288 saturated=3 predicted=233 actual=43 bbox=0.0513 elapsed=0.002210s checksum=6cf4d8e02d652b62
|
||||
scattered e=0.90 size=512x288 saturated=192 predicted=233 actual=43 bbox=0.9609 elapsed=0.019092s checksum=ef91f6cf7f41249b
|
||||
random5pct e=0.90 size=512x288 saturated=22223 predicted=216 actual=47 bbox=1.0000 elapsed=0.020793s checksum=fd976dae0c08cc92
|
||||
### OMP_NUM_THREADS=16
|
||||
sensor-bloom benchmark: size=512x288 threads=16 limit=1.000
|
||||
central e=0.00 size=512x288 saturated=3 predicted=1 actual=1 bbox=0.0001 elapsed=0.002165s checksum=aaec29797d030b5a
|
||||
scattered e=0.00 size=512x288 saturated=192 predicted=1 actual=1 bbox=0.7751 elapsed=0.002720s checksum=0c52973c56f09383
|
||||
random5pct e=0.00 size=512x288 saturated=22223 predicted=1 actual=1 bbox=1.0000 elapsed=0.002128s checksum=d489cdf14f756d3a
|
||||
central e=0.50 size=512x288 saturated=3 predicted=43 actual=15 bbox=0.0065 elapsed=0.001764s checksum=6625d592f76add42
|
||||
scattered e=0.50 size=512x288 saturated=192 predicted=43 actual=15 bbox=0.8433 elapsed=0.012356s checksum=44e39b09a596122f
|
||||
random5pct e=0.50 size=512x288 saturated=22223 predicted=40 actual=17 bbox=1.0000 elapsed=0.014079s checksum=d8c99ec18949278a
|
||||
central e=0.90 size=512x288 saturated=3 predicted=233 actual=43 bbox=0.0513 elapsed=0.002201s checksum=6cf4d8e02d652b62
|
||||
scattered e=0.90 size=512x288 saturated=192 predicted=233 actual=43 bbox=0.9609 elapsed=0.030918s checksum=ef91f6cf7f41249b
|
||||
random5pct e=0.90 size=512x288 saturated=22223 predicted=216 actual=47 bbox=1.0000 elapsed=0.020280s checksum=fd976dae0c08cc92
|
||||
```
|
||||
|
||||
- `saturated` is the number of channels initially above `E`; `predicted` is the
|
||||
model's reported `predicted_iterations`: the conservative decay bound plus the
|
||||
fp guard, clamped to the 4096 hard limit; `actual` is the number of synchronous
|
||||
rounds run; `bbox` is the fraction of the frame covered by the final monotonic
|
||||
bounding box. `actual <= predicted` in every row.
|
||||
- The `predicted`/`actual` fields and `checksum` are identical across the three
|
||||
runs above with `OMP_NUM_THREADS=1`, `4`, and `16`; only `elapsed` changes.
|
||||
|
||||
These numbers characterize the model on a bounded synthetic frame only. They do
|
||||
not establish any full-resolution or production-render performance.
|
||||
@@ -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
|
||||
|
||||
@@ -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. 开发顺序
|
||||
|
||||
+67
@@ -5,6 +5,7 @@
|
||||
#include "observer_track.h"
|
||||
#include "optics.h"
|
||||
#include "ray.h"
|
||||
#include "sensor_bloom.h"
|
||||
#include "spacetime.h"
|
||||
|
||||
#include <errno.h>
|
||||
@@ -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] "
|
||||
|
||||
@@ -0,0 +1,320 @@
|
||||
#include "sensor_bloom.h"
|
||||
|
||||
#include <math.h>
|
||||
#include <omp.h>
|
||||
#include <stdint.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
|
||||
/* 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;
|
||||
}
|
||||
@@ -0,0 +1,45 @@
|
||||
#ifndef SENSOR_BLOOM_H
|
||||
#define SENSOR_BLOOM_H
|
||||
|
||||
#include <stddef.h>
|
||||
|
||||
/* 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
|
||||
@@ -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 <math.h>
|
||||
#include <omp.h>
|
||||
#include <stdint.h>
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
|
||||
#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;
|
||||
}
|
||||
@@ -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
|
||||
|
||||
@@ -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 <math.h>
|
||||
#include <omp.h>
|
||||
#include <stdint.h>
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <sys/resource.h>
|
||||
#include <sys/wait.h>
|
||||
#include <unistd.h>
|
||||
|
||||
/* 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;
|
||||
}
|
||||
@@ -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
|
||||
|
||||
Reference in new issue
Block a user