diff --git a/crates/raw/src/llvc.rs b/crates/raw/src/llvc.rs index 28882c13d..d99acaa64 100644 --- a/crates/raw/src/llvc.rs +++ b/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, +} + +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> { + 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 = 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 { + 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 { + 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 { + 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 { + 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 { + 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 { + 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); + } + } + } } diff --git a/crates/raw/src/vendor/arw6_testenc.rs b/crates/raw/src/vendor/arw6_testenc.rs new file mode 100644 index 000000000..1bb6015ab --- /dev/null +++ b/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, Vec) { + 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 = (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 = (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 } +} diff --git a/crates/raw/src/vendor/mod.rs b/crates/raw/src/vendor/mod.rs index fcab4eabf..25f52cb5e 100644 --- a/crates/raw/src/vendor/mod.rs +++ b/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;