Feat: Add optional three-channel sensor bloom model

Add an opt-in, post-processing limited-response model applied to the
finished linear HDR before tone mapping. Each RGB channel is processed
independently and isotropically: overflow above E spreads to the eight
neighbours with a fixed 9-point stencil, while the rest is absorbed or
lost at the image boundary. The synchronous ping-pong update uses a
monotonic bounding box and a row-parallel, deterministic reduction; the
conservative round bound reserves fp guard rounds inside a 4096 hard
limit and fails before touching HDR when exceeded.

Expose --sensor-bloom-limit E and --sensor-bloom-transfer e (both
required together, default disabled), validate them before expensive
initialization, and route every output path through the same hook in
write_frame_outputs: raw FITS first, bloom, tone-mapped PNG/PPM, then the
mesh overlay. The raw --hdr-output FITS therefore stays pre-bloom.

Add a standalone unit test (stencil, boundary loss, cascade reference,
symmetry, thread determinism, convergence limits, validation, allocation
failure), CLI integration and regression coverage, an isolated
sensor-bloom-bench target, and document the model in the design, usage,
README, and build docs.

src/main.c

@@ -5,6 +5,7 @@
 #include "observer_track.h"
 #include "optics.h"
 #include "ray.h"
+#include "sensor_bloom.h"
 #include "spacetime.h"
 
 #include <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] "

src/sensor_bloom.c

@@ -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;
+}

src/sensor_bloom.h

@@ -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