From e31e1ed27ef6cc5505028d49ce637f067181894e Mon Sep 17 00:00:00 2001 From: Yingjie Wang Date: Sun, 27 Sep 2026 00:44:19 -0400 Subject: [PATCH] Feat: Add soft-clip tone mapping with legacy Reinhard Replace the per-channel Reinhard display transform with a parameterized soft clip T_p(x) = tanh(x^p)^(1/p), default softclip p=2, exposed through --tone-map and --tone-map-p. Keep --tone-map reinhard bit-compatible with the previous x/(1+x) curve for existing images and reject combining it with an explicit --tone-map-p. Route the primary image and the mesh overlay through the same ToneMapSettings; the linear HDR FITS writer stays pre-tone-map. Add a focused optics-linked tone-map test target, CLI success/error coverage, and document the display operator versus the Moffat effective PSF. --- Makefile | 14 ++- README.md | 23 +++-- README.zh-CN.md | 14 +-- nr_spacetime_movie_renderer_design.md | 13 ++- src/main.c | 32 ++++++- src/optics.c | 73 ++++++++++++--- src/optics.h | 27 +++++- tests/test_camera_cli.py | 25 +++++- tests/test_tone_map.c | 125 ++++++++++++++++++++++++++ usage.md | 43 ++++++++- 10 files changed, 359 insertions(+), 30 deletions(-) create mode 100644 tests/test_tone_map.c diff --git a/Makefile b/Makefile index e83769b..066710d 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 hip-psf-test hip-psf-bench fast-psf-fftw-bench minkowski schwarzschild FORCE +.PHONY: all backend clean run test tone-map-test 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) @@ -113,6 +113,7 @@ HIP_PSF_BENCH_TARGET := $(TEST_OUT_DIR)/benchmark_hip_psf CAMERA_TEST_TARGETS := $(TEST_OUT_DIR)/test_observer_minkowski $(TEST_OUT_DIR)/test_observer_schwarzschild 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 ifeq ($(PSF_BACKEND),hip) TARGET := $(BUILD_DIR)/$(TARGET_BASENAME)_hip @@ -215,6 +216,11 @@ $(TEST_OUT_DIR)/test_observer_schwarzschild: tests/test_observer.c $(COMMON_SOUR $(FAST_PSF_FFTW_TEST_TARGET): tests/test_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 tone-map regression links production optics.c only, so it neither needs a +# catalog nor touches ray tracing, the mesh, or any image writer at runtime. +$(TONE_MAP_TEST_TARGET): tests/test_tone_map.c src/optics.c src/optics.h $(CPU_FFTW_SOURCES) | $(TEST_OUT_DIR) + $(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc tests/test_tone_map.c src/optics.c $(CPU_FFTW_SOURCES) $(LDLIBS) -o $@ + $(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 $@ @@ -227,7 +233,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) +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_OUT_DIR)/test_observer_minkowski $(TEST_OUT_DIR)/test_observer_schwarzschild $(TEST_TARGET) @@ -236,8 +242,12 @@ test: $(CAMERA_TEST_TARGETS) $(TEST_TARGET) $(FRAME_TEST_TARGET) $(SCHWARZSCHILD $(OBSERVER_TRACK_TEST_TARGET) $(CATALOG_PREFETCH_TEST_TARGET) $(FAST_PSF_FFTW_TEST_RUN) + $(TONE_MAP_TEST_TARGET) python3 tests/test_camera_cli.py $(BUILD_DIR) $(TEST_OUT_DIR) +tone-map-test: $(TONE_MAP_TEST_TARGET) + $(TONE_MAP_TEST_TARGET) + fast-psf-fftw-bench: $(FAST_PSF_FFTW_BENCH_TARGET) clean: diff --git a/README.md b/README.md index d4fd931..11d2fed 100644 --- a/README.md +++ b/README.md @@ -10,7 +10,8 @@ four-dimensional spacetime. [![Schwarzschild black-hole lensing of the 2MASS Galactic-center star field](assets/images/schwarzschild_galactic_center.png)](assets/images/schwarzschild_galactic_center.png) *The 2MASS Galactic-center star field seen through Schwarzschild spacetime, -from a camera at radius 100 M. Click the image for the full 4K render; +from a camera at radius 100 M. This reference image uses the legacy Reinhard +tone map (`--tone-map reinhard`). Click the image for the full 4K render; [rendering command below](#example-galactic-center-field-through-schwarzschild-spacetime).* Stars are individual catalog point sources with direction, temperature, and @@ -92,12 +93,14 @@ mkdir -p output/imgs --all-sky-catalog assets/2mass/processed/all_sky \ --look-ra-deg 296 --look-dec-deg 27 --fov-deg 72 \ --width 3840 --height 2160 --exposure 1e12 \ + --tone-map reinhard \ --output output/imgs/summer_triangle.png ``` [![Summer Triangle in flat spacetime, exposure 1e12](assets/images/summer_triangle.png)](assets/images/summer_triangle.png) -*4K reference image with exposure `1e12`. Click to view at full resolution.* +*4K reference image with exposure `1e12`; it uses the legacy Reinhard tone +map (`--tone-map reinhard`). Click to view at full resolution.* Angles are in degrees; `--fov-deg` is the horizontal field of view. Exposure is an adjustable display multiplier. Use `--verbose` for progress during long @@ -119,6 +122,7 @@ mkdir -p output/imgs --psf-fwhm-pixels 2.7 --psf-moffat-beta 4.5 \ --psf-relative-tail 1e-8 --psf-min-y 0 --max-cache-psf-flux 1e8 \ --catalog-load-workers 4 \ + --tone-map reinhard \ --output output/imgs/schwarzschild_galactic_center.png ``` @@ -131,6 +135,11 @@ PNGs. The command above omits the optional HDR and lens-map exports. This example infers a position facing the hole from its look direction and radius. Adaptive refinement should 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). 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 @@ -158,13 +167,15 @@ mkdir -p output/imgs --width 3840 --height 2160 --fov-deg 90 \ --coarse-cell-pixels 16 --refine-max-level 3 \ --exposure 0.01 \ + --tone-map reinhard \ --output output/imgs/schwarzschild_near_horizon_outward_R2.1_0.01_refine3.png ``` [![Synthetic stellar sky seen outward by a static observer at r=2.1M](assets/images/schwarzschild_near_horizon_outward_R2.1_0.01_refine3.png)](assets/images/schwarzschild_near_horizon_outward_R2.1_0.01_refine3.png) *4K render with a 90° horizontal field of view, exposure `0.01`, and maximum -refinement level 3. Click to view at full resolution.* +refinement level 3; it uses the legacy Reinhard tone map (`--tone-map +reinhard`). Click to view at full resolution.* The distant sky occupies a bounded angular region around the outward direction, with repeated images crowded near its edge. Strong gravitational @@ -190,13 +201,15 @@ mkdir -p output/imgs --coarse-cell-pixels 32 \ --observer-radius 100 --exposure 0.2 \ --psf-fwhm-pixels 2.7 --psf-moffat-beta 4.5 --draw-mesh \ + --tone-map reinhard \ --output output/imgs/schwarzschild_test_grid.png ``` [![Synthetic stellar grid lensed by a Schwarzschild black hole, with adaptive mesh overlay](assets/images/schwarzschild_test_grid_mesh.png)](assets/images/schwarzschild_test_grid_mesh.png) -*4K test-grid reference image showing the `_mesh.png` overlay. Click to view at -full resolution.* +*4K test-grid reference image showing the `_mesh.png` overlay; it uses the +legacy Reinhard tone map (`--tone-map reinhard`). Click to view at full +resolution.* ### Freely falling Schwarzschild movie camera diff --git a/README.zh-CN.md b/README.zh-CN.md index 74b0154..ebbf81b 100644 --- a/README.zh-CN.md +++ b/README.zh-CN.md @@ -6,7 +6,7 @@ [![Schwarzschild 黑洞对 2MASS 银心方向星场的引力透镜效果](assets/images/schwarzschild_galactic_center.png)](assets/images/schwarzschild_galactic_center.png) -*从半径 100 M 处的相机观察 Schwarzschild 时空中的 2MASS 银心方向星场。点击图片查看完整 4K 图像；[渲染命令见下文](#示例schwarzschild-时空中的银心方向星场)。* +*从半径 100 M 处的相机观察 Schwarzschild 时空中的 2MASS 银心方向星场。该参考图像使用 legacy Reinhard tone map（`--tone-map reinhard`）。点击图片查看完整 4K 图像；[渲染命令见下文](#示例schwarzschild-时空中的银心方向星场)。* 恒星以独立的星表点源表示，保留方向、温度和振幅。渲染器将它们映射到相机图像中，处理多像、引力透镜放大和频移，再将保留亚像素位置的点扩散函数（PSF）累积到 HDR 图像。电影中的运动相机由世界线与四标架（tetrad）轨迹描述。 @@ -55,12 +55,13 @@ mkdir -p output/imgs --all-sky-catalog assets/2mass/processed/all_sky \ --look-ra-deg 296 --look-dec-deg 27 --fov-deg 72 \ --width 3840 --height 2160 --exposure 1e12 \ + --tone-map reinhard \ --output output/imgs/summer_triangle.png ``` [![平直时空中的夏季大三角，曝光为 1e12](assets/images/summer_triangle.png)](assets/images/summer_triangle.png) -*曝光为 `1e12` 的 4K 参考图像。点击查看完整分辨率。* +*曝光为 `1e12` 的 4K 参考图像；使用 legacy Reinhard tone map（`--tone-map reinhard`）。点击查看完整分辨率。* 角度单位为度，`--fov-deg` 表示水平视场角。曝光是可调节的显示倍率。长时间渲染时可使用 `--verbose` 查看进度。 @@ -79,12 +80,13 @@ mkdir -p output/imgs --psf-fwhm-pixels 2.7 --psf-moffat-beta 4.5 \ --psf-relative-tail 1e-8 --psf-min-y 0 --max-cache-psf-flux 1e8 \ --catalog-load-workers 4 \ + --tone-map reinhard \ --output output/imgs/schwarzschild_galactic_center.png ``` 这里使用更新后的参考图像的渲染设置。`--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 示例省略位置，因此按指向与半径推导出朝向黑洞的相机。自适应细分默认关闭，应根据所需图像精度配置。相机控制、图像序列、透镜映射复用、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 输出不受其影响。相机控制、图像序列、透镜映射复用、PSF 设置及 HDR 输出参见 [usage.md](usage.md)（英文）。两个可执行文件都可通过 `--help` 查看完整选项。 ### 示例：紧贴史瓦西视界向外看 @@ -103,12 +105,13 @@ mkdir -p output/imgs --width 3840 --height 2160 --fov-deg 90 \ --coarse-cell-pixels 16 --refine-max-level 3 \ --exposure 0.01 \ + --tone-map reinhard \ --output output/imgs/schwarzschild_near_horizon_outward_R2.1_0.01_refine3.png ``` [![r=2.1M 处的静态观者向外观察合成恒星天空](assets/images/schwarzschild_near_horizon_outward_R2.1_0.01_refine3.png)](assets/images/schwarzschild_near_horizon_outward_R2.1_0.01_refine3.png) -*水平视场角 90°、曝光 `0.01`、最大细分级别 3 的 4K 图像。点击查看完整分辨率。* +*水平视场角 90°、曝光 `0.01`、最大细分级别 3 的 4K 图像；使用 legacy Reinhard tone map（`--tone-map reinhard`）。点击查看完整分辨率。* 远方天空集中在朝外方向的有限角域中，多级成像在边缘附近密集堆叠。 强烈的引力蓝移使 3000 K 和 12000 K 的测试恒星都被推向蓝白色。 @@ -128,12 +131,13 @@ mkdir -p output/imgs --coarse-cell-pixels 32 \ --observer-radius 100 --exposure 0.2 \ --psf-fwhm-pixels 2.7 --psf-moffat-beta 4.5 --draw-mesh \ + --tone-map reinhard \ --output output/imgs/schwarzschild_test_grid.png ``` [![Schwarzschild 黑洞对合成恒星网格的透镜效果，叠加自适应网格](assets/images/schwarzschild_test_grid_mesh.png)](assets/images/schwarzschild_test_grid_mesh.png) -*4K 测试网格参考图像，展示 `_mesh.png` 网格叠加结果。点击查看完整分辨率。* +*4K 测试网格参考图像，展示 `_mesh.png` 网格叠加结果；使用 legacy Reinhard tone map（`--tone-map reinhard`）。点击查看完整分辨率。* ### Schwarzschild 自由落体相机轨迹 diff --git a/nr_spacetime_movie_renderer_design.md b/nr_spacetime_movie_renderer_design.md index 9aea84d..c260314 100644 --- a/nr_spacetime_movie_renderer_design.md +++ b/nr_spacetime_movie_renderer_design.md @@ -1454,10 +1454,19 @@ void rays_trace_generation( - HDR frame； - exposure； +- 线性 HDR → 可选 display operator → sRGB encoding； - tone mapping； - PNG/EXR； - 最终视频编码接口。 +当前实现决策：预览用的 tone-mapped PNG/PPM 路径为 +线性 RGB → display operator → sRGB 传输曲线 → 8-bit。默认 display operator +是逐通道 soft clip \(T_2(x)=\sqrt{\tanh(x^2)}\)（`--tone-map softclip +--tone-map-p 2`）；历史 Reinhard \(x/(1+x)\) 以 `--tone-map reinhard` +保留，用于复现旧输出。FITS `--hdr-output` 绕过 tone mapping，始终写 +pre-tone-map 线性 RGB。这里的 operator 只是当前简化显示管线，不是最终相机/ +传感器模型。 + --- # 28. 开发顺序 @@ -1553,7 +1562,9 @@ renderer 顶层架构原则上不应为 BBH 重新设计。 - PSF 模型； - ODE integrator； - observer trajectory 标准； -- tone mapping / exposure 规则； +- 当前预览/PNG 默认采用逐通道 \(T_2(x)=\sqrt{\tanh(x^2)}\) soft clip， + Reinhard 作为兼容模式保留；最终 production color management、传感器模型、 + 曝光标定和 HDR 视频编码规则仍未决定； - 是否需要 diffuse Milky Way background； - 是否将 ray redshift 变量定义为 `log(alpha p^0)` 或其他更方便的量。 diff --git a/src/main.c b/src/main.c index 411876b..a152178 100644 --- a/src/main.c +++ b/src/main.c @@ -59,6 +59,7 @@ typedef struct { int catalog_load_workers; const char *blackbody_table_path; RefinementConfig refinement; + ToneMapSettings tone_map; } Settings; static int parse_int(const char *text, int *value) { @@ -241,7 +242,8 @@ static int write_frame_outputs(const Settings *s, const FrameLensMesh *mesh, (void)fov_deg; #endif const int write_result = - write_tonemapped_image(paths->output_path, hdr, width, height); + write_tonemapped_image(paths->output_path, hdr, width, height, + &s->tone_map); fprintf(stderr, "Rendered %zu images from %zu catalog stars to %s (%s%s)\n", images, stars, paths->output_path, write_result == 0 ? "ok" : "write failed", note); @@ -249,7 +251,8 @@ static int write_frame_outputs(const Settings *s, const FrameLensMesh *mesh, return -1; if (paths->draw_mesh) { frame_draw_mesh(mesh, hdr, width, height, 0.5, 0.5); - if (write_tonemapped_image(paths->mesh_path, hdr, width, height)) { + if (write_tonemapped_image(paths->mesh_path, hdr, width, height, + &s->tone_map)) { fprintf(stderr, "Failed to write mesh overlay image: %s\n", paths->mesh_path); return -1; @@ -302,8 +305,10 @@ static int parse_args(int argc, char **argv, Settings *s, .angle_relative = 0.1, .jacobian_minimum = 1e-3, .min_edge_pixels = 0.5, - .min_area_pixels2 = 0.25}}; + .min_area_pixels2 = 0.25}, + .tone_map = {.op = TONE_MAP_SOFTCLIP, .p = 2.0}}; *write_path = NULL; + int tone_map_p_specified = 0; for (int i = 1; i < argc; ++i) { if (!strcmp(argv[i], "--catalog") && i + 1 < argc) s->catalog_path = argv[++i]; @@ -353,6 +358,17 @@ static int parse_args(int argc, char **argv, Settings *s, s->look_specified = 1; } else if (!strcmp(argv[i], "--exposure") && i + 1 < argc && !parse_positive(argv[++i], &s->exposure)) { + } else if (!strcmp(argv[i], "--tone-map") && i + 1 < argc) { + const char *mode = argv[++i]; + if (!strcmp(mode, "softclip")) + s->tone_map.op = TONE_MAP_SOFTCLIP; + else if (!strcmp(mode, "reinhard")) + s->tone_map.op = TONE_MAP_REINHARD; + else + return -1; + } 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], "--observer-position") || !strcmp(argv[i], "--observer-velocity")) && i + 3 < argc) { const int position = !strcmp(argv[i], "--observer-position"); @@ -427,6 +443,10 @@ static int parse_args(int argc, char **argv, Settings *s, else return -1; } + if (tone_map_p_specified && s->tone_map.op == TONE_MAP_REINHARD) { + fputs("--tone-map-p applies only to --tone-map softclip.\n", stderr); + return -1; + } return 0; } @@ -458,6 +478,9 @@ static void print_help(const char *program) { " --look-ra-deg D ICRS look direction right ascension in degrees (default: 90)\n" " --look-dec-deg D ICRS look direction declination in degrees (default: -90)\n" " --exposure E Linear exposure multiplier (default: 1e-3)\n" + " --tone-map MODE Display transform: softclip or reinhard\n" + " (default: softclip)\n" + " --tone-map-p P Softclip hardness P >= 1 (default: 2)\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" @@ -1156,7 +1179,8 @@ int main(int argc, char **argv) { "Usage: %s [--catalog PATH | --all-sky-catalog DIR] [--output PATH] [--width N] [--height " "N] [--fov-deg D] [--look-ra-deg D] [--look-dec-deg D] " "[--lens-map-input FILE | --lens-map-output FILE] " - "[--exposure E] [--observer-radius R | --observer-position X Y Z] " + "[--exposure E] [--tone-map softclip|reinhard] [--tone-map-p P] " + "[--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] " "[--max-magnification M] [--max-cache-psf-flux F] " diff --git a/src/optics.c b/src/optics.c index a84608f..487783a 100644 --- a/src/optics.c +++ b/src/optics.c @@ -757,10 +757,62 @@ void fast_psf_accumulator_report(const FastPsfAccumulator *accumulator, } -static unsigned char tonemap_channel(double hdr_value) +/* ToneMapSettings is always owned by the caller; invalid or missing settings + * fall back to the production default so the pure functions remain total. */ +static ToneMapSettings resolved_tone_map(const ToneMapSettings *settings) { - /* Reinhard tone mapping followed by the sRGB display transfer curve. */ - const double linear = hdr_value / (1.0 + hdr_value); + ToneMapSettings resolved = {.op = TONE_MAP_SOFTCLIP, .p = 2.0}; + if (settings == NULL) + return resolved; + if (settings->op == TONE_MAP_SOFTCLIP || + settings->op == TONE_MAP_REINHARD) + resolved.op = settings->op; + if (isfinite(settings->p) && settings->p >= 1.0) + resolved.p = settings->p; + return resolved; +} + +/* For z = x^p with small z, T_p(x) = x - x^(2p+1)/(3p) + ..., so the relative + * deviation from x is z^2/(3p). Requiring that to stay below the double + * rounding floor gives z < sqrt(3 p eps), i.e. about 2.6e-8 for the worst case + * p = 1. Below 1e-8 the exact identity branch is therefore accurate to double + * precision and avoids a subnormal x^p that would otherwise collapse to x = 0. + * This threshold follows from double precision, not image appearance. */ +static const double tone_map_small_z = 1e-8; + +/* tanh(z) rounds to 1.0 in double once 1 - tanh(z) ~ 2 exp(-2z) < eps/2, i.e. + * z > ln(4/eps)/2 = 18.7. At or above 20 the result is exactly 1, so the + * tanh/pow pair only adds rounding noise and overflowing x^p must short-circuit. */ +static const double tone_map_saturated_z = 20.0; + +double tone_map_linear_channel(double hdr_value, const ToneMapSettings *settings) +{ + const ToneMapSettings resolved = resolved_tone_map(settings); + /* Zero, negatives and NaN all map to black; only a positive finite value can + * acquire brightness, so squaring a negative softclip input cannot create + * light and a negative Reinhard input cannot become white. +Infinity is the + * only saturating input, and handling it here keeps lround() away from NaN. */ + if (!(hdr_value > 0.0)) + return 0.0; + if (!isfinite(hdr_value)) + return 1.0; + if (resolved.op == TONE_MAP_REINHARD) + return clamp(hdr_value / (1.0 + hdr_value), 0.0, 1.0); + const double p = resolved.p; + const double z = pow(hdr_value, p); + if (z == 0.0 || z < tone_map_small_z) + return clamp(hdr_value, 0.0, 1.0); + if (!isfinite(z) || z >= tone_map_saturated_z) + return 1.0; + return clamp(pow(tanh(z), 1.0 / p), 0.0, 1.0); +} + +unsigned char tone_map_srgb8_channel(double hdr_value, + const ToneMapSettings *settings) +{ + /* The sRGB display transfer curve and 8-bit quantization are unchanged from + * the original Reinhard-only writer. */ + const double linear = tone_map_linear_channel(hdr_value, settings); const double display = linear <= 0.0031308 ? 12.92 * linear : 1.055 * pow(linear, 1.0 / 2.4) - 0.055; return (unsigned char)lround(255.0 * clamp(display, 0.0, 1.0)); @@ -768,13 +820,13 @@ static unsigned char tonemap_channel(double hdr_value) #ifndef ENABLE_PNG static int write_tonemapped_ppm(const char *path, const double *hdr, int width, - int height) + int height, const ToneMapSettings *settings) { FILE *file = fopen(path, "wb"); if (file == NULL) return -1; fprintf(file, "P6\n%d %d\n255\n", width, height); for (int i = 0; i < width * height * 3; ++i) { - const unsigned char value = tonemap_channel(hdr[i]); + const unsigned char value = tone_map_srgb8_channel(hdr[i], settings); if (fwrite(&value, 1, 1, file) != 1) { fclose(file); return -1; } } return fclose(file) == 0 ? 0 : -1; @@ -783,7 +835,7 @@ static int write_tonemapped_ppm(const char *path, const double *hdr, int width, #ifdef ENABLE_PNG static int write_tonemapped_png(const char *path, const double *hdr, int width, - int height) + int height, const ToneMapSettings *settings) { FILE *file = fopen(path, "wb"); png_structp png = NULL; @@ -798,7 +850,7 @@ static int write_tonemapped_png(const char *path, const double *hdr, int width, pixels = malloc((size_t)width * height * 3); if (pixels == NULL) goto done; for (int i = 0; i < width * height * 3; ++i) - pixels[i] = tonemap_channel(hdr[i]); + pixels[i] = tone_map_srgb8_channel(hdr[i], settings); png_init_io(png, file); png_set_IHDR(png, info, (png_uint_32)width, (png_uint_32)height, 8, PNG_COLOR_TYPE_RGB, PNG_INTERLACE_NONE, @@ -816,17 +868,18 @@ done: } #endif -int write_tonemapped_image(const char *path, const double *hdr, int width, int height) +int write_tonemapped_image(const char *path, const double *hdr, int width, + int height, const ToneMapSettings *settings) { const size_t path_length = strlen(path); #ifdef ENABLE_PNG if (path_length >= 4 && strcmp(path + path_length - 4, ".png") == 0) - return write_tonemapped_png(path, hdr, width, height); + return write_tonemapped_png(path, hdr, width, height, settings); fputs("PNG output is enabled; use a .png output path.\n", stderr); return -1; #else if (path_length >= 4 && strcmp(path + path_length - 4, ".ppm") == 0) - return write_tonemapped_ppm(path, hdr, width, height); + return write_tonemapped_ppm(path, hdr, width, height, settings); fputs("PNG output is unavailable; use a .ppm output path or rebuild with libpng.\n", stderr); return -1; diff --git a/src/optics.h b/src/optics.h index ffdb25e..d3fc57e 100644 --- a/src/optics.h +++ b/src/optics.h @@ -10,6 +10,21 @@ typedef struct { double moffat_beta; } PointSpreadFunction; +/* Per-channel display transform applied to the linear HDR framebuffer before + * sRGB encoding. TONE_MAP_SOFTCLIP is the production preview default + * T_p(x) = tanh(x^p)^(1/p) with p >= 1; TONE_MAP_REINHARD keeps the historical + * x / (1 + x) curve so pre-soft-clip images remain reproducible. The FITS HDR + * writer never sees these settings. */ +typedef enum { + TONE_MAP_SOFTCLIP, + TONE_MAP_REINHARD +} ToneMapOperator; + +typedef struct { + ToneMapOperator op; + double p; +} ToneMapSettings; + /* One immutable process-wide kernel for the one PSF parameter/tail-policy pair * accepted by the current renderer. Its storage remains private to optics.c. */ typedef struct { @@ -87,6 +102,16 @@ typedef struct { extern "C" { #endif +/* Pure, testable display-transform pieces. `settings` may be NULL, in which + * case the production default (softclip, p = 2) is used. tone_map_linear_channel + * maps a linear HDR channel value to tone-mapped linear RGB in [0, 1]; + * tone_map_srgb8_channel applies the unchanged sRGB transfer curve and 8-bit + * quantization used by the PNG and PPM writers. */ +double tone_map_linear_channel(double hdr_value, + const ToneMapSettings *settings); +unsigned char tone_map_srgb8_channel(double hdr_value, + const ToneMapSettings *settings); + /* Integrate a Planck spectrum into absolute linear-sRGB spectral radiance * (W m^-2 sr^-1), before catalog amplitude and display exposure. */ LinearRgb blackbody_to_linear_rgb(double temperature_K); @@ -167,7 +192,7 @@ void splat_prepared_cached_event(double *hdr, int width, int height, const PsfKernelCache *cache); /* Writes PNG when built with libpng; non-libpng builds use PPM fallback. */ int write_tonemapped_image(const char *path, const double *hdr, int width, - int height); + int height, const ToneMapSettings *settings); #ifdef ENABLE_HDR_OUTPUT /* Writes the pre-tone-mapping framebuffer as a 32-bit RGB FITS image. */ int write_hdr_fits(const char *path, const double *hdr, int width, int height, diff --git a/tests/test_camera_cli.py b/tests/test_camera_cli.py index 2ae5062..503024c 100644 --- a/tests/test_camera_cli.py +++ b/tests/test_camera_cli.py @@ -88,10 +88,15 @@ with tempfile.TemporaryDirectory(prefix='gr-camera-cli-') as directory: hdr_available = '--hdr-output' in help_text 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 '--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. + # Using it (rather than the former 0.1) keeps the soft-clip default from + # saturating the whole frame, so the image comparisons below stay + # discriminating and the mesh overlay remains visible. common = ['--catalog', 'assets/sky_grid_5deg.csv', '--width', 64, - '--height', 48, '--fov-deg', 80, '--exposure', 0.1, + '--height', 48, '--fov-deg', 80, '--exposure', 1e-3, '--coarse-cell-pixels', 8, '--refine-max-level', 0, '--psf-relative-tail', 1e-4] def render(name, *options): @@ -129,6 +134,16 @@ with tempfile.TemporaryDirectory(prefix='gr-camera-cli-') as directory: assert render('partial', partial, value) == render( 'complete', '--look-ra-deg', ra, '--look-dec-deg', dec) + # Tone-map selection: the default must equal explicit softclip p=2, the + # hardness must be adjustable, and legacy Reinhard must remain available + # while producing a visibly different image. + softclip2 = render('tonemap_softclip2', '--tone-map', 'softclip', + '--tone-map-p', 2) + assert softclip2 == default + assert softclip2 == render('tonemap_softclip', '--tone-map', 'softclip') + assert render('tonemap_p1', '--tone-map', 'softclip', '--tone-map-p', 1) + assert render('tonemap_reinhard', '--tone-map', 'reinhard') != softclip2 + # --draw-mesh must leave the primary image untouched and only add a # mesh-overlay sibling. baseline = render('mesh_base') @@ -182,6 +197,14 @@ with tempfile.TemporaryDirectory(prefix='gr-camera-cli-') as directory: (['--observer-track', 'missing.csv', '--observer-velocity', 0, 0, 0], 'cannot be combined'), (['--frames-dir', tmp, '--look-ra-deg', 0], 'cannot be combined'), (['--lens-map-input', 'missing.grlens', '--camera-roll-deg', 0], 'cannot be combined'), + (['--tone-map', 'unknown'], None), + (['--tone-map'], None), + (['--tone-map-p'], None), + (['--tone-map-p', 0], None), + (['--tone-map-p', 0.5], None), + (['--tone-map-p', 'nan'], None), + (['--tone-map-p', 'inf'], None), + (['--tone-map', 'reinhard', '--tone-map-p', 2], 'applies only'), ] if backend == 'schwarzschild': errors += [(['--observer-position', 1.5, 0, 0, '--observer-velocity', -0.5, 0, 0], 'capture cutoff'), diff --git a/tests/test_tone_map.c b/tests/test_tone_map.c new file mode 100644 index 0000000..394def8 --- /dev/null +++ b/tests/test_tone_map.c @@ -0,0 +1,125 @@ +/* Regression tests for the display tone mapping in src/optics.c. + * + * These call the production public functions so the test cannot drift from the + * implementation; no formula is duplicated here. The test needs neither a + * catalog, ray tracing, a GPU, nor image files, and it builds in both the + * ENABLE_PNG=1 and ENABLE_PNG=0 configurations. */ +#include "optics.h" + +#include +#include +#include + +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; +} + +int main(void) +{ + const double grid[] = {0.0, 1e-12, 1e-9, 1e-6, 1e-3, 0.01, 0.1, 0.25, + 0.5, 0.75, 1.0, 1.5, 2.0, 4.0, 10.0, 1e6}; + const size_t grid_count = sizeof grid / sizeof grid[0]; + const double hardness[] = {1.0, 2.0, 4.0}; + const double small[] = {1e-12, 1e-9, 1e-6}; + const double negative[] = {-1e-12, -0.5, -1.0, -1e6}; + const ToneMapSettings softclip2 = {TONE_MAP_SOFTCLIP, 2.0}; + const ToneMapSettings reinhard = {TONE_MAP_REINHARD, 2.0}; + + /* Basic properties for every supported softclip hardness. */ + for (size_t h = 0; h < 3; ++h) { + const ToneMapSettings settings = {TONE_MAP_SOFTCLIP, hardness[h]}; + check(tone_map_linear_channel(0.0, &settings) == 0.0, + "softclip T(0) == 0"); + check(tone_map_linear_channel(INFINITY, &settings) == 1.0, + "softclip T(+infinity) == 1"); + check(tone_map_linear_channel(1e300, &settings) == 1.0, + "softclip overflow saturates at 1"); + for (size_t n = 0; n < sizeof negative / sizeof negative[0]; ++n) + check(tone_map_linear_channel(negative[n], &settings) == 0.0, + "softclip negative input maps to 0"); + double previous = tone_map_linear_channel(grid[0], &settings); + for (size_t i = 1; i < grid_count; ++i) { + const double current = tone_map_linear_channel(grid[i], &settings); + check(isfinite(current) && current >= 0.0 && current <= 1.0, + "softclip output is finite and bounded"); + check(current >= previous, "softclip is monotonically non-decreasing"); + previous = current; + } + for (size_t i = 0; i < sizeof small / sizeof small[0]; ++i) { + const double ratio = + tone_map_linear_channel(small[i], &settings) / small[i]; + check(close_to(ratio, 1.0, 1e-9), + "softclip T(x)/x tends to 1 for small x"); + } + } + + /* Known double-precision values of T_2(x) = sqrt(tanh(x^2)). */ + check(close_to(tone_map_linear_channel(0.5, &softclip2), + 0.49489257663023107, 1e-12), "T2(0.5) value"); + check(close_to(tone_map_linear_channel(1.0, &softclip2), + 0.8726936208978296, 1e-12), "T2(1.0) value"); + check(close_to(tone_map_linear_channel(2.0, &softclip2), + 0.9996645936208139, 1e-12), "T2(2.0) value"); + + /* NULL settings select the production default softclip, p = 2. */ + check(tone_map_linear_channel(0.5, NULL) == + tone_map_linear_channel(0.5, &softclip2), + "NULL settings default to softclip p=2"); + + /* Both operators must sanitize non-finite and nonpositive inputs: the + * public contract is a finite result in [0, 1] and lround() must never see + * NaN. This includes Reinhard, whose raw formula would map x < -1 to white + * and produce NaN for +Infinity or NaN. */ + const ToneMapSettings operators[] = {softclip2, reinhard}; + for (size_t i = 0; i < sizeof operators / sizeof operators[0]; ++i) { + check(tone_map_linear_channel(NAN, &operators[i]) == 0.0, + "NaN maps to 0 for every operator"); + check(tone_map_linear_channel(INFINITY, &operators[i]) == 1.0, + "+infinity maps to 1 for every operator"); + check(tone_map_linear_channel(-2.0, &operators[i]) == 0.0, + "negative input maps to 0 for every operator"); + check(tone_map_linear_channel(-0.5, &operators[i]) == 0.0, + "small negative input maps to 0 for every operator"); + check(tone_map_srgb8_channel(NAN, &operators[i]) == 0, + "NaN 8-bit output is black"); + check(tone_map_srgb8_channel(INFINITY, &operators[i]) == 255, + "+infinity 8-bit output is white"); + } + + /* Reinhard must reproduce the pre-soft-clip formula bit for bit. */ + check(tone_map_linear_channel(0.5, &reinhard) == 0.5 / (1.0 + 0.5), + "reinhard x=0.5 -> 1/3"); + check(tone_map_linear_channel(1.0, &reinhard) == 1.0 / (1.0 + 1.0), + "reinhard x=1.0 -> 1/2"); + check(tone_map_linear_channel(2.0, &reinhard) == 2.0 / (1.0 + 2.0), + "reinhard x=2.0 -> 2/3"); + + /* The 8-bit path is what the PNG and PPM writers actually emit. */ + check(tone_map_srgb8_channel(0.0, &softclip2) == 0, + "softclip p=2 x=0 is black"); + check(tone_map_srgb8_channel(2.0, &softclip2) == 255, + "softclip p=2 x=2 saturates the white channel"); + check(tone_map_srgb8_channel(2.0, &reinhard) != 255, + "reinhard x=2 is not fully saturated"); + check(tone_map_srgb8_channel(0.5, &softclip2) > + tone_map_srgb8_channel(0.5, &reinhard), + "softclip is brighter than reinhard at x=0.5"); + + if (failures != 0) { + fprintf(stderr, "%d tone-map assertion(s) failed\n", failures); + return EXIT_FAILURE; + } + puts("tone-map tests passed"); + return EXIT_SUCCESS; +} diff --git a/usage.md b/usage.md index 2eaded1..6b9b3a2 100644 --- a/usage.md +++ b/usage.md @@ -267,6 +267,46 @@ at the existing cache radius. Images above `F` retain the direct fallback. `--max-magnification M` (default unlimited) caps the per-triangle rendering magnification before flux is formed; it is an explicit preview approximation. +### Tone mapping and display + +`--exposure` is a linear multiplier applied while the HDR framebuffer and its +PSF splats are accumulated. The Moffat profile is the empirical effective +optical PSF, absorbing seeing, jitter, and residual aberrations; it is not +changed by any display choice. + +The tone map is a later camera/display response applied to the finished linear +HDR image, so it does not alter the linear HDR PSF. The default operator is a +parameterized, per-channel soft clip + +\[ +T_p(x) = \left[\tanh(x^p)\right]^{1/p}, \qquad p \ge 1, +\] + +selected with `--tone-map softclip --tone-map-p 2`, i.e. + +\[ +T_2(x) = \sqrt{\tanh(x^2)}. +\] + +It is close to the identity in the dark, then rolls smoothly into a saturated +shoulder above \(x \sim 1\). Because it is applied independently to each linear +RGB channel, a bright star transitions from its colored Moffat wings into a +white saturated core. Larger `p` approaches \(\min(x,1)\). This per-channel +response is what gives the current simplified preview pipeline the crisp white +cores seen in real astrophotography, rather than the gray, soft cores a +Reinhard curve leaves on bright stars. + +`--tone-map reinhard` selects the historical \(x/(1+x)\) display transform and +is retained only to reproduce output rendered before the soft-clip operator was +introduced. FITS output always remains the pre-tone-map linear RGB framebuffer, +independent of operator and `p`; use it to measure physical flux, peak values, +and FWHM. Tone-mapped PNG/PPM data are display-referred and must not be used for +those measurements. + +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. + ### Fast preview mode `--fast-mode` replaces the per-event PSF splat with a two-stage approximation: @@ -372,7 +412,8 @@ ordinary output, replacing its extension with `_HDR.fits`; for example, `output/imgs/ring_HDR.fits`. It writes the pre-tone-mapping RGB 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. The ordinary +is applied, and the `--tone-map` operator and `--tone-map-p` value never affect +this file. 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