Merge PR #858

crates/engine/src/camera_preview.rs

@@ -833,6 +833,22 @@ fn to_prophoto() -> Mat3 {
     PROPHOTO.from_xyz().mul(&bradford(D65, D50)).mul(&REC2020.to_xyz())
 }
 
+/// [`HueSat::apply`] restores the luminance in full only if the table kept at least this share of it. A colour
+/// the table darkens to [`NO_LUMINANCE_KEPT`] or less (only near-black noise, whose channels go negative) keeps
+/// its colour; in between the two are blended, so noisy shadows show no step.
+const MIN_LUMINANCE_KEPT: f32 = 0.25;
+const NO_LUMINANCE_KEPT: f32 = 0.15;
+
+/// How much of the table's colour a pixel gets when the table kept `share` of its luminance: 0 up to
+/// [`NO_LUMINANCE_KEPT`], 1 from [`MIN_LUMINANCE_KEPT`], a smoothstep between.
+fn restore_weight(share: f32) -> f32 {
+    if !share.is_finite() {
+        return 0.0;
+    }
+    let t = ((share - NO_LUMINANCE_KEPT) / (MIN_LUMINANCE_KEPT - NO_LUMINANCE_KEPT)).clamp(0.0, 1.0);
+    t * t * (3.0 - 2.0 * t)
+}
+
 /// A fitted [`CameraLook::hue_sat`] table, ready to apply to linear Rec.2020 pixels. It changes
 /// hue and saturation only: luminance is restored, the camera tone curve owns it.
 pub(crate) struct HueSat<'a> {
@@ -852,7 +868,19 @@ impl<'a> HueSat<'a> {
         let mul = |m: &[[f32; 3]; 3], v: [f32; 3]| -> [f32; 3] { std::array::from_fn(|i| m[i][0] * v[0] + m[i][1] * v[1] + m[i][2] * v[2]) };
         let out = mul(&self.from, self.table.apply(mul(&self.to, rgb)));
         let (before, after) = (luminance_2020(rgb), luminance_2020(out));
-        if before > 0.0 && after > 0.0 && out.iter().all(|v| v.is_finite()) { out.map(|v| v * before / after) } else { rgb }
+        // restoring the luminance divides by the table's: near black, where noise makes some channels
+        // negative, the table can bring it close to 0, and the gain would turn the pixel white (#820)
+        if before <= 0.0 || after <= 0.0 || !out.iter().all(|v| v.is_finite()) {
+            return rgb;
+        }
+        let w = restore_weight(after / before);
+        if w >= 1.0 {
+            return out.map(|v| v * before / after);
+        }
+        if w <= 0.0 {
+            return rgb;
+        }
+        std::array::from_fn(|i| rgb[i] + (out[i] * before / after - rgb[i]) * w)
     }
 }
 
@@ -1439,6 +1467,61 @@ mod tests {
         assert!(fit_hue_sat(&[], &Mat3::IDENTITY).is_none(), "too few samples");
     }
 
+    #[test]
+    fn hue_sat_never_turns_near_black_noise_white() {
+        // a saturating table, as fitted to high-ISO Sony files; near-black noise with mixed-sign channels
+        let table = HsvTable { hue_divisions: 4, sat_divisions: 2, val_divisions: 1, data: vec![[10.0, 1.5, 1.0]; 8], srgb_value: false };
+        let hue_sat = HueSat::new(&table).unwrap();
+        let mut seed = 0x9e37_79b9_7f4a_7c15u64;
+        let mut next = || {
+            seed ^= seed << 13;
+            seed ^= seed >> 7;
+            seed ^= seed << 17;
+            (seed >> 11) as f32 / (1u64 << 53) as f32 * 2.0 - 1.0
+        };
+        let mut worst = 0.0f32;
+        for _ in 0..200_000 {
+            let rgb = [0.002 + 0.004 * next(), 0.002 + 0.004 * next(), 0.002 + 0.004 * next()];
+            if luminance_2020(rgb) <= 0.0 {
+                continue;
+            }
+            let out = hue_sat.apply(rgb);
+            let scale = rgb.iter().fold(0.0f32, |m, v| m.max(v.abs()));
+            worst = worst.max(out.iter().fold(0.0f32, |m, v| m.max(v.abs())) / scale);
+        }
+        assert!(worst < 20.0, "a near-black pixel came out {worst}× its largest channel");
+    }
+
+    #[test]
+    fn hue_sat_restore_weight_blends_smoothly_between_the_thresholds() {
+        assert_eq!(restore_weight(0.0), 0.0);
+        assert_eq!(restore_weight(NO_LUMINANCE_KEPT), 0.0);
+        assert_eq!(restore_weight(MIN_LUMINANCE_KEPT), 1.0);
+        assert_eq!(restore_weight(1.0), 1.0);
+        assert_eq!(restore_weight(f32::NAN), 0.0);
+        // monotonic and without a step: neighbouring shares 0.001 apart differ by little
+        let mut last = 0.0f32;
+        for i in 0..=400 {
+            let w = restore_weight(i as f32 * 0.001);
+            assert!(w >= last && w - last < 0.02, "step at {}: {last} → {w}", i as f32 * 0.001);
+            last = w;
+        }
+    }
+
+    #[test]
+    fn hue_sat_keeps_colours_the_table_leaves_bright_exactly_as_before() {
+        // the blend only touches colours whose luminance the table cut below MIN_LUMINANCE_KEPT
+        let table = HsvTable { hue_divisions: 4, sat_divisions: 2, val_divisions: 1, data: vec![[10.0, 1.5, 1.0]; 8], srgb_value: false };
+        let hue_sat = HueSat::new(&table).unwrap();
+        let mul = |m: &[[f32; 3]; 3], v: [f32; 3]| -> [f32; 3] { std::array::from_fn(|i| m[i][0] * v[0] + m[i][1] * v[1] + m[i][2] * v[2]) };
+        for rgb in [[0.4, 0.2, 0.1], [0.05, 0.3, 0.2], [0.18, 0.18, 0.18], [0.6, 0.1, 0.5]] {
+            let out = mul(&hue_sat.from, table.apply(mul(&hue_sat.to, rgb)));
+            let (before, after) = (luminance_2020(rgb), luminance_2020(out));
+            assert!(after >= MIN_LUMINANCE_KEPT * before);
+            assert_eq!(hue_sat.apply(rgb), out.map(|v| v * before / after));
+        }
+    }
+
     fn raw_of(format: RawFormat, make: &str, model: &str, color: lightcraft_raw::ColorData) -> RawImage {
         let area = lightcraft_raw::Rect::new(0, 0, 2, 2);
         RawImage {