Misc: add a test catalog for testing PSF

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.
2026-09-27 02:48:52 -04:00 · 2026-09-27 00:44:19 -04:00
11 changed files with 441 additions and 30 deletions
@@ -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:
@@ -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

@@ -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 自由落体相机轨迹

@@ -0,0 +1,82 @@
+ra_deg,dec_deg,temperature_K,amplitude_sr
+-0.500000,0.500000,2500,1.288768695e-13
+-0.375000,0.500000,2500,5.130680588e-14
+-0.250000,0.500000,2500,2.042560732e-14
+-0.125000,0.500000,2500,8.131580736e-15
+0.000000,0.500000,2500,3.237240599e-15
+0.125000,0.500000,2500,1.288768695e-15
+0.250000,0.500000,2500,5.130680588e-16
+0.375000,0.500000,2500,2.042560732e-16
+0.500000,0.500000,2500,8.131580736e-17
+-0.500000,0.375000,3044.7082,5.529309491e-14
+-0.375000,0.375000,3044.7082,2.201257756e-14
+-0.250000,0.375000,3044.7082,8.763364971e-15
+-0.125000,0.375000,3044.7082,3.488758433e-15
+0.000000,0.375000,3044.7082,1.388899748e-15
+0.125000,0.375000,3044.7082,5.529309491e-16
+0.250000,0.375000,3044.7082,2.201257756e-16
+0.375000,0.375000,3044.7082,8.763364971e-17
+0.500000,0.375000,3044.7082,3.488758433e-17
+-0.500000,0.250000,3708.0992,2.727501648e-14
+-0.375000,0.250000,3708.0992,1.085837964e-14
+-0.250000,0.250000,3708.0992,4.322798794e-15
+-0.125000,0.250000,3708.0992,1.720937197e-15
+0.000000,0.250000,3708.0992,6.851174381e-16
+0.125000,0.250000,3708.0992,2.727501648e-16
+0.250000,0.250000,3708.0992,1.085837964e-16
+0.375000,0.250000,3708.0992,4.322798794e-17
+0.500000,0.250000,3708.0992,1.720937197e-17
+-0.500000,0.125000,4516.0321,1.501787774e-14
+-0.375000,0.125000,4516.0321,5.978724815e-15
+-0.250000,0.125000,4516.0321,2.380173219e-15
+-0.125000,0.125000,4516.0321,9.475640258e-16
+0.000000,0.125000,4516.0321,3.772320332e-16
+0.125000,0.125000,4516.0321,1.501787774e-16
+0.250000,0.125000,4516.0321,5.978724815e-17
+0.375000,0.125000,4516.0321,2.380173219e-17
+0.500000,0.125000,4516.0321,9.475640258e-18
+-0.500000,0.000000,5500,9.010857559e-15
+-0.375000,0.000000,5500,3.587287007e-15
+-0.250000,0.000000,5500,1.428124680e-15
+-0.125000,0.000000,5500,5.685466757e-16
+0.000000,0.000000,5500,2.263425084e-16
+0.125000,0.000000,5500,9.010857559e-17
+0.250000,0.000000,5500,3.587287007e-17
+0.375000,0.000000,5500,1.428124680e-17
+0.500000,0.000000,5500,5.685466757e-18
+-0.500000,-0.125000,7067.9722,5.170563867e-15
+-0.375000,-0.125000,7067.9722,2.058438551e-15
+-0.250000,-0.125000,7067.9722,8.194791473e-16
+-0.125000,-0.125000,7067.9722,3.262405247e-16
+0.000000,-0.125000,7067.9722,1.298786922e-16
+0.125000,-0.125000,7067.9722,5.170563867e-17
+0.250000,-0.125000,7067.9722,2.058438551e-17
+0.375000,-0.125000,7067.9722,8.194791473e-18
+0.500000,-0.125000,7067.9722,3.262405247e-18
+-0.500000,-0.250000,9082.9511,3.209818520e-15
+-0.375000,-0.250000,9082.9511,1.277851769e-15
+-0.250000,-0.250000,9082.9511,5.087219521e-16
+-0.125000,-0.250000,9082.9511,2.025258569e-16
+0.000000,-0.250000,9082.9511,8.062699587e-17
+0.125000,-0.250000,9082.9511,3.209818520e-17
+0.250000,-0.250000,9082.9511,1.277851769e-17
+0.375000,-0.250000,9082.9511,5.087219521e-18
+0.500000,-0.250000,9082.9511,2.025258569e-18
+-0.500000,-0.375000,11672.372,2.109932559e-15
+-0.375000,-0.375000,11672.372,8.399792812e-16
+-0.250000,-0.375000,11672.372,3.344017750e-16
+-0.125000,-0.375000,11672.372,1.331277445e-16
+0.000000,-0.375000,11672.372,5.299910967e-17
+0.125000,-0.375000,11672.372,2.109932559e-17
+0.250000,-0.375000,11672.372,8.399792812e-18
+0.375000,-0.375000,11672.372,3.344017750e-18
+0.500000,-0.375000,11672.372,1.331277445e-18
+-0.500000,-0.500000,15000,1.446233663e-15
+-0.375000,-0.500000,15000,5.757559917e-16
+-0.250000,-0.500000,15000,2.292125888e-16
+-0.125000,-0.500000,15000,9.125117517e-17
+0.000000,-0.500000,15000,3.632774716e-17
+0.125000,-0.500000,15000,1.446233663e-17
+0.250000,-0.500000,15000,5.757559917e-18
+0.375000,-0.500000,15000,2.292125888e-18
+0.500000,-0.500000,15000,9.125117517e-19
@@ -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)` 或其他更方便的量。

@@ -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] "
@@ -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;
@@ -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,
@@ -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'),
@@ -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 <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+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;
+}
@@ -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