LR-IMP-FORMATS: ARW6 parity dequantiser and 5/3 inverse lifting with phases

crates/raw/src/llvc.rs

@@ -145,6 +145,138 @@ pub(crate) fn vld_decode_line(bits: &mut Bits<'_>, source_bits: usize, width: us
     }
 }
 
+/// Parity dequantiser, Phase 0 (oracle-verified): 0 -> 0, qi = 0 -> q, else
+/// `sign(q) * (((2|q| + 1) << (qi - 1)) - (|q| & 1))`. Saturates instead of overflowing on hostile input.
+pub(crate) fn dequant(q: i32, qi: u32) -> i32 {
+    if q == 0 || qi == 0 {
+        return q;
+    }
+    let a = i64::from(q).abs();
+    let m = (((2 * a + 1) << (qi - 1).min(30)) - (a & 1)).min(i64::from(i32::MAX));
+    (if q < 0 { -m } else { m }) as i32
+}
+
+/// A row-major plane of integer samples.
+#[derive(Clone, Debug, PartialEq, Eq)]
+pub(crate) struct Plane {
+    pub width: usize,
+    pub height: usize,
+    pub data: Vec<i32>,
+}
+
+impl Plane {
+    pub fn zeros(width: usize, height: usize) -> Plane {
+        Plane { width, height, data: vec![0; width.saturating_mul(height)] }
+    }
+    pub fn row(&self, y: usize) -> Option<&[i32]> {
+        let start = y.checked_mul(self.width)?;
+        self.data.get(start..start.checked_add(self.width)?)
+    }
+    fn ok(&self) -> bool {
+        self.width.checked_mul(self.height) == Some(self.data.len())
+    }
+}
+
+/// `(low, high)` sample counts of an `n`-sample signal whose low-pass samples sit at indices = `phase` (mod 2).
+pub(crate) fn band_rows(n: usize, phase: u8) -> (usize, usize) {
+    let low = (n + 1).saturating_sub(usize::from(phase)) / 2;
+    (low, n - low)
+}
+
+/// `v[i]` with `i` clamped into range (the whole-sample symmetric extension of the band); 0 for an empty band.
+fn at(v: &[i32], i: isize) -> i32 {
+    let last = v.len() as isize - 1;
+    if last < 0 {
+        return 0;
+    }
+    v.get(i.clamp(0, last) as usize).copied().unwrap_or(0)
+}
+
+/// 1-D inverse 5/3 lifting, RDD 34 §6.5, with X[-1] := X[1], X[n] := X[n-2] (see `band_rows` for the sizes).
+fn inverse_1d(l: &[i32], h: &[i32], phase: u8) -> Result<Vec<i32>> {
+    let (n, p) = (l.len() + h.len(), isize::from(phase));
+    if band_rows(n, phase) != (l.len(), h.len()) || phase > 1 {
+        return Err(RawError::Corrupt("ARW6 wavelet: band sizes do not match".into()));
+    }
+    let even: Vec<i32> = l.iter().enumerate().map(|(j, &lv)| lv - ((at(h, j as isize - 1 + p) + at(h, j as isize + p) + 2) >> 2)).collect();
+    let mut x = vec![0; n];
+    let (lo, hi) = (usize::from(phase), 1 - usize::from(phase));
+    for (xv, &e) in x.iter_mut().skip(lo).step_by(2).zip(&even) {
+        *xv = e;
+    }
+    for (k, (xv, &hv)) in x.iter_mut().skip(hi).step_by(2).zip(h).enumerate() {
+        *xv = hv + ((at(&even, k as isize - p) + at(&even, k as isize + 1 - p)) >> 1);
+    }
+    Ok(x)
+}
+
+/// Gathers the lines of `a` along one axis (columns when `vertical`), applies `f`, and scatters back.
+fn lines(low: &Plane, high: &Plane, phase: u8, vertical: bool) -> Result<Plane> {
+    if !low.ok() || !high.ok() {
+        return Err(RawError::Corrupt("ARW6 wavelet: plane size mismatch".into()));
+    }
+    let (across_l, across_h) = if vertical { (low.width, high.width) } else { (low.height, high.height) };
+    if across_l != across_h {
+        return Err(RawError::Corrupt("ARW6 wavelet: band sizes do not match".into()));
+    }
+    let (nl, nh) = if vertical { (low.height, high.height) } else { (low.width, high.width) };
+    let n = nl + nh;
+    let (w, h) = if vertical { (across_l, n) } else { (n, across_l) };
+    let mut out = Plane::zeros(w, h);
+    let get = |p: &Plane, i: usize| -> Vec<i32> {
+        if vertical { p.data.iter().skip(i).step_by(p.width.max(1)).copied().collect() } else { p.row(i).map(<[i32]>::to_vec).unwrap_or_default() }
+    };
+    for i in 0..across_l {
+        let x = inverse_1d(&get(low, i), &get(high, i), phase)?;
+        if vertical {
+            for (o, v) in out.data.iter_mut().skip(i).step_by(w.max(1)).zip(x) {
+                *o = v;
+            }
+        } else if let Some(r) = out.data.get_mut(i * w..(i + 1) * w) {
+            r.copy_from_slice(&x);
+        }
+    }
+    Ok(out)
+}
+
+/// Vertical inverse 5/3 (along rows of samples): `low`/`high` are stacked bands of equal width.
+pub(crate) fn inverse_53_1d_rows(low: &Plane, high: &Plane, phase: u8) -> Result<Plane> {
+    lines(low, high, phase, true)
+}
+
+/// Horizontal inverse 5/3: `low`/`high` are side-by-side bands of equal height.
+pub(crate) fn inverse_53_1d_cols(low: &Plane, high: &Plane, phase: u8) -> Result<Plane> {
+    lines(low, high, phase, false)
+}
+
+/// 2-D inverse, RDD 34 §6.5 order: vertical IWT (phase `vphase`) on (ll, lh) and (hl, hh), then horizontal (phase 0).
+pub(crate) fn inverse_53_2d(ll: &Plane, lh: &Plane, hl: &Plane, hh: &Plane, vphase: u8) -> Result<Plane> {
+    let lo = inverse_53_1d_rows(ll, lh, vphase)?;
+    let hi = inverse_53_1d_rows(hl, hh, vphase)?;
+    inverse_53_1d_cols(&lo, &hi, 0)
+}
+
+/// The ten bands of a 3-level decomposition (LL only at the deepest level).
+pub(crate) struct Bands3 {
+    pub ll3: Plane,
+    pub hl3: Plane,
+    pub lh3: Plane,
+    pub hh3: Plane,
+    pub hl2: Plane,
+    pub lh2: Plane,
+    pub hh2: Plane,
+    pub hl1: Plane,
+    pub lh1: Plane,
+    pub hh1: Plane,
+}
+
+/// Full inverse: level 3 uses `phases[2]`, level 2 `phases[1]`, level 1 `phases[0]`.
+pub(crate) fn reconstruct3(b: &Bands3, phases: [u8; 3]) -> Result<Plane> {
+    let ll2 = inverse_53_2d(&b.ll3, &b.lh3, &b.hl3, &b.hh3, phases[2])?;
+    let ll1 = inverse_53_2d(&ll2, &b.lh2, &b.hl2, &b.hh2, phases[1])?;
+    inverse_53_2d(&ll1, &b.lh1, &b.hl1, &b.hh1, phases[0])
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
@@ -229,4 +361,68 @@ mod tests {
             let _ = vld_decode_line(&mut b, src.len() * 8, width);
         }
     }
+
+    fn col(v: &[i32]) -> Plane {
+        Plane { width: 1, height: v.len(), data: v.to_vec() }
+    }
+    #[test]
+    fn dequant_parity_rule() {
+        for (qi, table) in
+            [(0u32, [-3, -2, -1, 0, 1, 2, 3]), (1, [-6, -5, -2, 0, 2, 5, 6]), (2, [-13, -10, -5, 0, 5, 10, 13]), (3, [-27, -20, -11, 0, 11, 20, 27])]
+        {
+            for (q, v) in (-3..=3).zip(table) {
+                assert_eq!(dequant(q, qi), v, "q={q} qi={qi}");
+            }
+        }
+    }
+    #[test]
+    fn quantize_inverts_dequant() {
+        for qi in 0..6 {
+            for q in -50..=50 {
+                assert_eq!(crate::vendor::arw6_testenc::quantize(dequant(q, qi), qi), q, "q={q} qi={qi}");
+            }
+        }
+    }
+    #[test]
+    fn dequant_does_not_overflow() {
+        assert_eq!(dequant(i32::MAX, 31), i32::MAX);
+        assert_eq!(dequant(i32::MIN, 40), -i32::MAX);
+    }
+    #[test]
+    fn band_rows_counts() {
+        assert_eq!(band_rows(1093, 0), (547, 546));
+        assert_eq!(band_rows(1093, 1), (546, 547));
+        assert_eq!(band_rows(1668, 0), (834, 834));
+        assert_eq!(band_rows(5, 1), (2, 3));
+        assert_eq!(band_rows(1, 1), (0, 1));
+    }
+    #[test]
+    fn inverse_53_hand_vectors() {
+        // computed by hand from RDD 34 §6.5 with X[-1] := X[1], X[n] := X[n-2]; checked against the Phase 0 model
+        assert_eq!(inverse_53_1d_rows(&col(&[10, 20]), &col(&[3, -4]), 0).unwrap().data, vec![8, 17, 20, 16]);
+        assert_eq!(inverse_53_1d_rows(&col(&[10, 20]), &col(&[3, -4]), 1).unwrap().data, vec![13, 10, 12, 22]);
+        assert_eq!(inverse_53_1d_rows(&col(&[10, 20, 30]), &col(&[3, -4]), 0).unwrap().data, vec![8, 17, 20, 22, 32]);
+        assert_eq!(inverse_53_1d_rows(&col(&[3, -4]), &col(&[10, 20, 30]), 1).unwrap().data, vec![5, -5, 9, -17, 13]);
+    }
+    #[test]
+    fn inverse_53_rejects_mismatched_bands() {
+        assert!(inverse_53_1d_rows(&col(&[1, 2, 3]), &col(&[1]), 0).is_err());
+        assert!(inverse_53_1d_rows(&col(&[1, 2]), &Plane::zeros(2, 2), 0).is_err());
+        assert!(inverse_53_1d_cols(&col(&[1, 2]), &col(&[1, 2, 3]), 0).is_err());
+    }
+    proptest::proptest! {
+        #[test]
+        fn forward_then_inverse_is_identity(w in 1usize..24, h in 1usize..40, p1 in 0u8..2, p2 in 0u8..2, p3 in 0u8..2, seed in proptest::num::u64::ANY) {
+            use crate::vendor::arw6_testenc::{forward3, forward_53_2d};
+            use proptest::prop_assert_eq;
+            let mut s = seed | 1;
+            let mut next = || { s ^= s << 13; s ^= s >> 7; s ^= s << 17; (s % 8191) as i32 - 4095 };
+            let x = Plane { width: w, height: h, data: (0..w * h).map(|_| next()).collect() };
+            let (ll, lh, hl, hh) = forward_53_2d(&x, p1);
+            prop_assert_eq!(inverse_53_2d(&ll, &lh, &hl, &hh, p1).unwrap(), x.clone());
+            if w >= 8 && h >= 8 {
+                prop_assert_eq!(reconstruct3(&forward3(&x, [p1, p2, p3]), [p1, p2, p3]).unwrap(), x);
+            }
+        }
+    }
 }

crates/raw/src/vendor/arw6_testenc.rs

@@ -0,0 +1,79 @@
+//! Test-only ARW6 synthetic encoder pieces: forward 5/3 transforms (exact inverse of `llvc::inverse_*`) and `quantize`.
+
+use crate::llvc::{Bands3, Plane, band_rows};
+
+/// Exact inverse of `dequant` on representable values: `sign(v) * (|v| >> qi)`.
+pub(crate) fn quantize(v: i32, qi: u32) -> i32 {
+    v.signum() * (v.abs() >> qi)
+}
+
+fn at(v: &[i32], i: isize) -> i32 {
+    if v.is_empty() { 0 } else { v[i.clamp(0, v.len() as isize - 1) as usize] }
+}
+
+/// Forward 1-D 5/3 with whole-sample symmetric extension (`plan/arw6/scratch/fwd.py::fwd1d`).
+fn forward_1d(x: &[i32], phase: u8) -> (Vec<i32>, Vec<i32>) {
+    let n = x.len() as isize;
+    let p = isize::from(phase);
+    let xs = |i: isize| -> i32 {
+        let r = if i < 0 {
+            -i
+        } else if i >= n {
+            2 * (n - 1) - i
+        } else {
+            i
+        };
+        if (0..n).contains(&r) { x[r as usize] } else { 0 }
+    };
+    let (nl, nh) = band_rows(x.len(), phase);
+    let h: Vec<i32> = (0..nh as isize)
+        .map(|k| {
+            let pos = 2 * k + 1 - p;
+            xs(pos) - ((xs(pos - 1) + xs(pos + 1)) >> 1)
+        })
+        .collect();
+    let l = (0..nl as isize).map(|j| xs(2 * j + p) + ((at(&h, j - 1 + p) + at(&h, j + p) + 2) >> 2)).collect();
+    (l, h)
+}
+
+fn split(x: &Plane, phase: u8, vertical: bool) -> (Plane, Plane) {
+    let (across, n) = if vertical { (x.width, x.height) } else { (x.height, x.width) };
+    let (nl, nh) = band_rows(n, phase);
+    let dims = |m| if vertical { (across, m) } else { (m, across) };
+    let (mut lo, mut hi) = (Plane::zeros(dims(nl).0, dims(nl).1), Plane::zeros(dims(nh).0, dims(nh).1));
+    for i in 0..across {
+        let line: Vec<i32> = (0..n).map(|k| if vertical { x.data[k * x.width + i] } else { x.data[i * x.width + k] }).collect();
+        let (l, h) = forward_1d(&line, phase);
+        for (p, v) in [(&mut lo, l), (&mut hi, h)] {
+            for (k, v) in v.into_iter().enumerate() {
+                let idx = if vertical { k * p.width + i } else { i * p.width + k };
+                p.data[idx] = v;
+            }
+        }
+    }
+    (lo, hi)
+}
+
+pub(crate) fn forward_53_1d_rows(x: &Plane, phase: u8) -> (Plane, Plane) {
+    split(x, phase, true)
+}
+
+pub(crate) fn forward_53_1d_cols(x: &Plane, phase: u8) -> (Plane, Plane) {
+    split(x, phase, false)
+}
+
+/// Horizontal first (phase 0), then vertical; returns `(ll, lh, hl, hh)`, the exact inverse of `inverse_53_2d`.
+pub(crate) fn forward_53_2d(x: &Plane, vphase: u8) -> (Plane, Plane, Plane, Plane) {
+    let (lo, hi) = forward_53_1d_cols(x, 0);
+    let (ll, lh) = forward_53_1d_rows(&lo, vphase);
+    let (hl, hh) = forward_53_1d_rows(&hi, vphase);
+    (ll, lh, hl, hh)
+}
+
+/// 3-level forward, the exact inverse of `reconstruct3` with the same `phases`.
+pub(crate) fn forward3(x: &Plane, phases: [u8; 3]) -> Bands3 {
+    let (ll1, lh1, hl1, hh1) = forward_53_2d(x, phases[0]);
+    let (ll2, lh2, hl2, hh2) = forward_53_2d(&ll1, phases[1]);
+    let (ll3, lh3, hl3, hh3) = forward_53_2d(&ll2, phases[2]);
+    Bands3 { ll3, hl3, lh3, hh3, hl2, lh2, hh2, hl1, lh1, hh1 }
+}

crates/raw/src/vendor/mod.rs

@@ -1,6 +1,8 @@
 //! Vendor raw formats (TIFF-based) and helpers shared by them.
 
 pub mod arw;
+#[cfg(test)]
+pub(crate) mod arw6_testenc;
 mod arw_lens;
 pub mod cr2;
 pub mod cr3;