LR-EDIT-OPTICS-PROFILE: share the radial warp fit between RW2 and Sony

Move the forward-map inversion, least-squares fit and solve4 out of the RW2
decoder into vendor::warp_from_forward_map; RW2 calls it in DNG-normalised units.

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

crates/raw/src/vendor/mod.rs

@@ -111,6 +111,99 @@ pub(crate) fn white_from_data(data: &[u16], bits: u32) -> f32 {
     }
 }
 
+/// Output radii sampled to fit the inverse mapping.
+const DISTORTION_SAMPLES: usize = 64;
+
+/// Least-squares radial warp from a camera's forward distortion map.
+///
+/// `corrected` maps a source radius to the output radius it lands at; it must rise monotonically up to `t_max`, the
+/// output radius of the image corners. Radii use the caller's unit; pass `WarpRectilinear`'s normalisation (distance
+/// from the centre to the **farthest** corner = 1) so the result is directly a `WarpRectilinear` plane. The map is
+/// scanned at 4096 points up to a source radius of `3 * t_max`: it must reach `t_max` within that range and rise at
+/// every scan point (monotonicity is checked only there), otherwise `None`; a NaN at a scan point is rejected the
+/// same way, other non-finite values are caught by the fit-error gate. The map is inverted by bisection and
+/// `source = u·(k0 + k1u² + k2u⁴ + k3u⁶)` is fitted over output radii `u·t_max`, `u` in [0, 1]: the returned
+/// `[k0, k1, k2, k3]` give the source radius in the caller's unit as a function of `u`. `None` too when the worst
+/// fit error exceeds `max_fit_error`.
+pub(crate) fn warp_from_forward_map(corrected: &dyn Fn(f64) -> f64, t_max: f64, max_fit_error: f64) -> Option<[f64; 4]> {
+    if !(t_max.is_finite() && t_max > 0.0) {
+        return None;
+    }
+    let mut top = None;
+    let mut prev = 0.0;
+    for i in 1..=4096 {
+        let r = 3.0 * t_max * i as f64 / 4096.0;
+        let v = corrected(r);
+        if v.is_nan() || v <= prev {
+            return None;
+        }
+        prev = v;
+        if v >= t_max {
+            top = Some(r);
+            break;
+        }
+    }
+    let top = top?;
+    let source = |t: f64| {
+        let (mut lo, mut hi) = (0.0, top);
+        for _ in 0..60 {
+            let m = 0.5 * (lo + hi);
+            if corrected(m) < t { lo = m } else { hi = m }
+        }
+        0.5 * (lo + hi)
+    };
+    let samples: Vec<(f64, f64)> = (0..=DISTORTION_SAMPLES)
+        .map(|i| {
+            let u = i as f64 / DISTORTION_SAMPLES as f64;
+            (u, source(u * t_max))
+        })
+        .collect();
+    let basis = |u: f64| {
+        let u2 = u * u;
+        [u, u * u2, u * u2 * u2, u * u2 * u2 * u2]
+    };
+    let mut a = [[0.0; 4]; 4];
+    let mut rhs = [0.0; 4];
+    for &(u, v) in &samples {
+        let f = basis(u);
+        for (i, fi) in f.iter().enumerate() {
+            for (j, fj) in f.iter().enumerate() {
+                a[i][j] += fi * fj;
+            }
+            rhs[i] += fi * v;
+        }
+    }
+    let k = solve4(a, rhs)?;
+    let fit = |u: f64| basis(u).iter().zip(&k).map(|(f, c)| f * c).sum::<f64>();
+    let worst = samples.iter().map(|&(u, v)| (fit(u) - v).abs()).fold(0.0, f64::max);
+    (worst <= max_fit_error).then_some(k)
+}
+
+/// Solve the 4 × 4 system `a·x = b` (Gaussian elimination with partial pivoting); `None` when singular or not finite.
+fn solve4(mut a: [[f64; 4]; 4], mut b: [f64; 4]) -> Option<[f64; 4]> {
+    for c in 0..4 {
+        let p = (c..4).max_by(|&i, &j| a[i][c].abs().total_cmp(&a[j][c].abs()))?;
+        if a[p][c].is_nan() || a[p][c].abs() <= 1e-300 {
+            return None;
+        }
+        a.swap(c, p);
+        b.swap(c, p);
+        for r in c + 1..4 {
+            let f = a[r][c] / a[c][c];
+            for k in c..4 {
+                a[r][k] -= f * a[c][k];
+            }
+            b[r] -= f * b[c];
+        }
+    }
+    let mut x = [0.0; 4];
+    for c in (0..4).rev() {
+        let s: f64 = (c + 1..4).map(|k| a[c][k] * x[k]).sum();
+        x[c] = (b[c] - s) / a[c][c];
+    }
+    x.iter().all(|v| v.is_finite()).then_some(x)
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
@@ -160,4 +253,32 @@ mod tests {
         let w = white_from_data(&d, 14);
         assert!((13500.0..=13583.0).contains(&w), "{w}");
     }
+
+    #[test]
+    fn warp_from_forward_map_recovers_identity_and_rejects_folds() {
+        let k = warp_from_forward_map(&|r| r, 1.0, 1e-9).unwrap();
+        assert!((k[0] - 1.0).abs() < 1e-9 && k[1].abs() < 1e-9 && k[2].abs() < 1e-9 && k[3].abs() < 1e-9);
+        let barrel = |r: f64| r * (1.0 + 0.05 * r * r); // mild barrel: source < output radius
+        let t_max = 0.9;
+        let k = warp_from_forward_map(&barrel, t_max, 1e-3).unwrap();
+        assert!(k[0] < t_max && k[0] > 0.8 && k[1] < 0.0); // source radius at u = 1 is below t_max
+        for i in 0..=20 {
+            let u = i as f64 / 20.0;
+            let s = u * (k[0] + u * u * (k[1] + u * u * (k[2] + u * u * k[3])));
+            assert!((barrel(s) - u * t_max).abs() < 1e-5, "u {u}");
+        }
+        assert!(warp_from_forward_map(&|r| if r < 0.5 { r } else { 0.5 - (r - 0.5) }, 1.0, 1e-3).is_none()); // fold, never reaches t_max
+        // a tiny dip (two scan points wide, 1e-4 deep) that still reaches t_max and fits well: only the monotonicity check rejects it
+        let dip = |r: f64| {
+            if r < 0.5 {
+                r
+            } else if r < 0.5015 {
+                0.5 - 0.05 * (r - 0.5)
+            } else {
+                r - 0.001575
+            }
+        };
+        assert!(warp_from_forward_map(&dip, 1.0, 1e-2).is_none());
+        assert!(warp_from_forward_map(&|_| f64::NAN, 1.0, 1e-3).is_none());
+    }
 }

crates/raw/src/vendor/rw2.rs

@@ -145,8 +145,6 @@ const FORMAT8_FIXED: [(u16, &[u16]); 9] = [
 const DISTORTION_INFO: u16 = 0x0119;
 /// Fixed-point one of the distortion coefficients and zoom.
 const DISTORTION_ONE: f64 = 32768.0;
-/// Output radii sampled to fit the inverse mapping.
-const DISTORTION_SAMPLES: usize = 64;
 /// Largest error of the fitted inverse (relative to the half diagonal) accepted: real lenses fit within 1.2e-4.
 const DISTORTION_MAX_FIT_ERROR: f64 = 1e-3;
 /// Most strips seen in a file is 3 (the 12008 × 8008 high-resolution mode).
@@ -650,87 +648,12 @@ fn distortion(ifd: &Ifd, active: Rect) -> Option<Opcode> {
         let r2 = r * r;
         r * (1.0 + r2 * (k2 + r2 * (k4 + r2 * k6)))
     };
-    // the output corners show corrected radius `t_max`; the map must rise monotonically to it (a fold is no lens)
-    let t_max = half_diag * zoom / radius;
-    let mut top = None;
-    let mut prev = 0.0;
-    for i in 1..=4096 {
-        let r = 3.0 * t_max * i as f64 / 4096.0;
-        let v = corrected(r);
-        if v.is_nan() || v <= prev {
-            return None;
-        }
-        prev = v;
-        if v >= t_max {
-            top = Some(r);
-            break;
-        }
-    }
-    let top = top?;
-    let source = |t: f64| {
-        let (mut lo, mut hi) = (0.0, top);
-        for _ in 0..60 {
-            let m = 0.5 * (lo + hi);
-            if corrected(m) < t { lo = m } else { hi = m }
-        }
-        0.5 * (lo + hi)
-    };
-    // least-squares fit of source radius = u·(kr0 + kr1·u² + kr2·u⁴ + kr3·u⁶) over the output radii u ∈ [0, 1]
-    let samples: Vec<(f64, f64)> = (0..=DISTORTION_SAMPLES)
-        .map(|i| {
-            let u = i as f64 / DISTORTION_SAMPLES as f64;
-            (u, source(u * t_max) * radius / half_diag)
-        })
-        .collect();
-    let basis = |u: f64| {
-        let u2 = u * u;
-        [u, u * u2, u * u2 * u2, u * u2 * u2 * u2]
-    };
-    let mut a = [[0.0; 4]; 4];
-    let mut rhs = [0.0; 4];
-    for &(u, v) in &samples {
-        let f = basis(u);
-        for (i, fi) in f.iter().enumerate() {
-            for (j, fj) in f.iter().enumerate() {
-                a[i][j] += fi * fj;
-            }
-            rhs[i] += fi * v;
-        }
-    }
-    let k = solve4(a, rhs)?;
-    let fit = |u: f64| basis(u).iter().zip(&k).map(|(f, c)| f * c).sum::<f64>();
-    let worst = samples.iter().map(|&(u, v)| (fit(u) - v).abs()).fold(0.0, f64::max);
-    if worst.is_nan() || worst > DISTORTION_MAX_FIT_ERROR {
-        return None;
-    }
+    // work in DNG-normalised units (corner distance = 1), where the output corners sit at radius `zoom`
+    let scale = radius / half_diag;
+    let k = super::warp_from_forward_map(&|s| corrected(s / scale) * scale, zoom, DISTORTION_MAX_FIT_ERROR)?;
     Some(Opcode::WarpRectilinear { planes: vec![[k[0], k[1], k[2], k[3], 0.0, 0.0]], center: [0.5, 0.5] })
 }
 
-/// Solve the 4 × 4 system `a·x = b` (Gaussian elimination with partial pivoting); `None` when singular or not finite.
-fn solve4(mut a: [[f64; 4]; 4], mut b: [f64; 4]) -> Option<[f64; 4]> {
-    for c in 0..4 {
-        let p = (c..4).max_by(|&i, &j| a[i][c].abs().total_cmp(&a[j][c].abs()))?;
-        if a[p][c].is_nan() || a[p][c].abs() <= 1e-300 {
-            return None;
-        }
-        a.swap(c, p);
-        b.swap(c, p);
-        for r in c + 1..4 {
-            let f = a[r][c] / a[c][c];
-            for k in c..4 {
-                a[r][k] -= f * a[c][k];
-            }
-            b[r] -= f * b[c];
-        }
-    }
-    let mut x = [0.0; 4];
-    for c in (0..4).rev() {
-        let s: f64 = (c + 1..4).map(|k| a[c][k] * x[k]).sum();
-        x[c] = (b[c] - s) / a[c][c];
-    }
-    x.iter().all(|v| v.is_finite()).then_some(x)
-}
-
 #[cfg(test)]
 mod tests {
     use super::*;