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.

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] "

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;

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,