LR-IMP-FORMATS: a float raw with NaN or infinite samples no longer panics the demosaic

Cause: a floating-point CFA DNG is read as stored (no check of its float
samples). AHD and PPG clamp each green estimate between its two neighbouring
greens with `est.clamp(a.min(b), a.max(b))`; when both greens are NaN the
bounds are NaN and `f32::clamp` panics (ppg.rs:45, ahd.rs green pass). Any
float DNG with two adjacent NaN greens crashed decode + develop. Infinities
alone did not panic but turned into NaN in the demosaic arithmetic.

Fix, at two layers:
- demosaic: a `clamp_between` helper replaces the three clamps. With finite
  or infinite bounds it is the same `f32::clamp` call (bit-identical
  results); with NaN bounds it returns the estimate unclamped, so the panic
  is gone whatever puts NaN into the mosaic.
- DNG decode: non-finite float samples are replaced where the white level is
  known. NaN and -inf read as black (0), +inf as the white level (a clipped
  highlight). No sensor records NaN or infinity, so a valid file is never
  touched; finite samples, negative or above white included, stay as stored.

Tests: `demosaic::tests::non_finite_samples_do_not_panic` (NaN block, +/-inf
and an all-NaN mosaic through Bilinear, PPG and AHD; plus `clamp_between`
equal to the old clamp bit for bit over finite, signed-zero and infinite
values) and `non_finite_float_samples_decode_and_develop` (a float CFA DNG,
uncompressed and Deflate, with NaN and +/-inf samples: decoded values
checked, develop with all three methods gives a finite image). Both panicked
before the fix at ppg.rs:45 / ahd.rs.

Co-Authored-By: Claude <noreply@anthropic.com>

crates/raw/src/demosaic/ahd.rs

@@ -77,10 +77,10 @@ fn tile(m: &Mosaic, x0: usize, y0: usize, tw: usize, th: usize) -> Vec<[f32; 3]>
             }
             let (l, r) = (mos[cl(x - 1, y)], mos[cl(x + 1, y)]);
             let est = (l + r) * 0.5 + (2.0 * c - mos[cl(x - 2, y)] - mos[cl(x + 2, y)]) * 0.25;
-            gh[i] = est.clamp(l.min(r), l.max(r));
+            gh[i] = super::clamp_between(est, l, r);
             let (u, d) = (mos[cl(x, y - 1)], mos[cl(x, y + 1)]);
             let est = (u + d) * 0.5 + (2.0 * c - mos[cl(x, y - 2)] - mos[cl(x, y + 2)]) * 0.25;
-            gv[i] = est.clamp(u.min(d), u.max(d));
+            gv[i] = super::clamp_between(est, u, d);
         }
     }
     // 2. full candidates

crates/raw/src/demosaic/mod.rs

@@ -52,6 +52,14 @@ pub(crate) fn reflect(i: isize, n: usize) -> usize {
     i as usize
 }
 
+/// `v` clamped between `a` and `b`: `v.clamp(a.min(b), a.max(b))`, except that it returns `v` as is when both
+/// bounds are NaN (NaN samples in a float raw), where `f32::clamp` panics.
+#[inline]
+pub(crate) fn clamp_between(v: f32, a: f32, b: f32) -> f32 {
+    let (lo, hi) = (a.min(b), a.max(b));
+    if lo <= hi { v.clamp(lo, hi) } else { v }
+}
+
 impl Mosaic<'_> {
     #[inline]
     pub fn get(&self, x: isize, y: isize) -> f32 {
@@ -222,6 +230,37 @@ mod tests {
         assert_eq!(demosaic(&n, Method::Ahd).data, vec![[0.1, 0.2, 0.3]]);
     }
 
+    #[test]
+    fn non_finite_samples_do_not_panic() {
+        // regression: NaN samples (a float DNG) gave AHD and PPG NaN bounds for `f32::clamp`, which panics
+        let cfa = Cfa::bayer("RGGB").unwrap();
+        let mut patchy = mosaic_from_rgb(&edgy_scene(24, 20), &cfa);
+        for (i, v) in patchy.data.iter_mut().enumerate() {
+            let (x, y) = (i % 24, i / 24);
+            if (8..12).contains(&x) && (6..10).contains(&y) {
+                *v = f32::NAN;
+            } else if i % 7 == 0 {
+                *v = if i % 2 == 0 { f32::INFINITY } else { f32::NEG_INFINITY };
+            }
+        }
+        let all_nan = Normalized { width: 9, height: 7, cpp: 1, data: vec![f32::NAN; 63], cfa: Some(cfa) };
+        for n in [&patchy, &all_nan] {
+            for m in [Method::Bilinear, Method::Ppg, Method::Ahd] {
+                let out = demosaic(n, m);
+                assert_eq!((out.width, out.height), (n.width, n.height), "{m:?}");
+            }
+        }
+        // without NaN bounds it is the `f32::clamp` it replaced, bit for bit (finite results stay identical)
+        let vals = [f32::NEG_INFINITY, -3e38, -2.5, -0.0, 0.0, 1e-40, 0.5, 1.0, 3e38, f32::INFINITY];
+        for v in vals {
+            for a in vals {
+                for b in vals {
+                    assert_eq!(clamp_between(v, a, b).to_bits(), v.clamp(a.min(b), a.max(b)).to_bits(), "{v} {a} {b}");
+                }
+            }
+        }
+    }
+
     #[test]
     fn reflect_indices() {
         assert_eq!(reflect(-1, 5), 1);

crates/raw/src/demosaic/ppg.rs

@@ -42,7 +42,7 @@ pub(crate) fn green_plane(m: &Mosaic) -> Vec<f32> {
                 _ => (ge, gw, ce),
             };
             let est = (near * 3.0 + far + c - cfar) * 0.25;
-            *o = est.clamp(near.min(far), near.max(far));
+            *o = super::clamp_between(est, near, far);
         }
     });
     g

crates/raw/src/dng.rs

@@ -234,6 +234,16 @@ pub(crate) fn decode_as(bytes: &[u8], mode: Mode, format: RawFormat) -> Result<R
     if white.is_empty() {
         white = vec![default_white];
     }
+    // No sensor records NaN or infinity, so such a float sample is damage. Left in, it turns into NaN in demosaic
+    // arithmetic and spreads; NaN and −∞ read as black (0), +∞ as the white level (a clipped highlight). Finite
+    // samples, negative or above white included, are kept as stored.
+    if let RawData::F32(v) = &mut data {
+        for (i, s) in v.iter_mut().enumerate() {
+            if !s.is_finite() {
+                *s = if *s > 0.0 { white.get(i % cpp).or(white.first()).copied().unwrap_or(default_white) } else { 0.0 };
+            }
+        }
+    }
 
     let dim = raw.u64s(t::BLACK_LEVEL_REPEAT_DIM).unwrap_or_default();
     let (br, bc) = match dim.as_slice() {

crates/raw/tests/dng.rs

@@ -378,6 +378,43 @@ fn floating_point_dngs() {
     assert_eq!(decode(&dng(raw, ByteOrder::Big)).unwrap().data, RawData::F32(vals));
 }
 
+#[test]
+fn non_finite_float_samples_decode_and_develop() {
+    // regression: NaN samples in a float CFA DNG panicked AHD and PPG (`f32::clamp` with NaN bounds)
+    let (w, h) = (24, 20);
+    let mut raw = synthetic(w, h, Some(Cfa::bayer("RGGB").unwrap()), 1);
+    let mut vals: Vec<f32> = (0..w * h).map(|i| ((i * 37) % 200) as f32 / 256.0).collect();
+    for (i, v) in vals.iter_mut().enumerate() {
+        let (x, y) = (i % w, i / w);
+        if (8..14).contains(&x) && (6..12).contains(&y) {
+            *v = f32::NAN;
+        } else if i % 11 == 0 {
+            *v = if i % 2 == 0 { f32::INFINITY } else { f32::NEG_INFINITY };
+        }
+    }
+    raw.data = RawData::F32(vals.clone());
+    raw.black = BlackLevel::uniform(0.0);
+    raw.white = vec![0.75];
+    for compression in [DngCompression::Uncompressed, DngCompression::Deflate { tile: 16, half: false }] {
+        let bytes = write_dng(&raw, &DngWriteOptions { compression, ..Default::default() }).unwrap();
+        let img = decode(&bytes).unwrap();
+        // NaN and −∞ read as black, +∞ as the white level; finite samples are kept bit-exactly
+        let RawData::F32(got) = &img.data else { panic!("float data expected") };
+        for (a, b) in vals.iter().zip(got) {
+            let want = match *a {
+                f32::INFINITY => 0.75,
+                v if v.is_finite() => v,
+                _ => 0.0,
+            };
+            assert_eq!(want.to_bits(), b.to_bits(), "{a} read as {b}");
+        }
+        for m in [Method::Bilinear, Method::Ppg, Method::Ahd] {
+            let rgb = img.develop(m).unwrap();
+            assert!(rgb.data.iter().all(|p| p.iter().all(|v| v.is_finite())), "{m:?}");
+        }
+    }
+}
+
 #[test]
 fn lj92_strips_and_levels_and_linearization() {
     let (w, h) = (20, 12);