diff --git a/crates/engine/src/camera_preview.rs b/crates/engine/src/camera_preview.rs index c1efa7506..790ae4aab 100644 --- a/crates/engine/src/camera_preview.rs +++ b/crates/engine/src/camera_preview.rs @@ -1425,6 +1425,42 @@ impl<'a> HueSat<'a> { } std::array::from_fn(|i| rgb[i] + (out[i] * before / after - rgb[i]) * w) } + + /// [`Self::apply`] to `N` pixels at once (`rgb[channel][lane]`), bit for bit, through + /// [`HsvTable::apply_lanes`], with the guards as selects. + #[inline] + pub fn apply_lanes(&self, rgb: [[f32; N]; 3]) -> [[f32; N]; 3] { + let mul = |m: &[[f32; 3]; 3], v: &[[f32; N]; 3]| { + let mut out = [[0.0f32; N]; 3]; + for (row, out) in m.iter().zip(&mut out) { + for l in 0..N { + out[l] = row[0] * v[0][l] + row[1] * v[1][l] + row[2] * v[2][l]; + } + } + out + }; + let out = mul(&self.from, &self.table.apply_lanes(mul(&self.to, &rgb))); + let mut restored = rgb; + for l in 0..N { + let (input, out) = ([rgb[0][l], rgb[1][l], rgb[2][l]], [out[0][l], out[1][l], out[2][l]]); + let (before, after) = (luminance_2020(input), luminance_2020(out)); + let untouched = (before <= 0.0) | (after <= 0.0) | !(out[0].is_finite() & out[1].is_finite() & out[2].is_finite()); + let w = restore_weight(after / before); + for (c, restored) in restored.iter_mut().enumerate() { + let full = out[c] * before / after; + let blend = input[c] + (full - input[c]) * w; + // (`w >= 1` and `w <= 0` exclude each other) + restored[l] = if untouched | (w <= 0.0) { + input[c] + } else if w >= 1.0 { + full + } else { + blend + }; + } + } + restored + } } /// HSV hue (degrees) and saturation; `None` for black or non-finite colours. @@ -2197,6 +2233,51 @@ mod tests { assert!(worst < 20.0, "a near-black pixel came out {worst}× its largest channel"); } + /// `apply_lanes` gives `apply`'s bits for any colour: random, near-black noise (the restore's blend), grey, + /// negative, tiny, above 1, huge, NaN and infinite components; on a fitted look's shape and a saturating table. + #[test] + fn hue_sat_lanes_match_one_by_one_bit_for_bit() { + let mut seed = 0x2545_f491_4f6c_dd1du64; + let mut next = move || { + seed ^= seed << 13; + seed ^= seed >> 7; + seed ^= seed << 17; + seed + }; + let unit = |r: u64, i: usize| ((r >> (16 * i)) & 0xffff) as f32 / 65535.0; + let fitted = (0..72 * 5 * 5) + .map(|_| { + let r = next(); + [unit(r, 0) * 60.0 - 30.0, 0.5 + unit(r, 1), 1.0] + }) + .collect(); + let fitted = HsvTable { hue_divisions: 72, sat_divisions: 5, val_divisions: 5, data: fitted, srgb_value: true }; + let saturating = HsvTable { hue_divisions: 4, sat_divisions: 2, val_divisions: 1, data: vec![[10.0, 1.5, 1.0]; 8], srgb_value: false }; + let special = [0.0, -0.0, 1e-30, -1e-6, 0.002, 0.5, 1.0, 3.0, 1e30, f32::INFINITY, f32::NEG_INFINITY, f32::NAN]; + let mut colours = Vec::new(); + for a in special { + for b in special { + colours.extend(special.map(|c| [a, b, c])); + } + } + for _ in 0..60_000 { + let r = next(); + let (offset, scale) = [(0.15, 1.3), (0.15e-3, 1.3e-3), (0.6, 5.0), (-0.002, 0.008)][(r >> 62) as usize]; + colours.push(std::array::from_fn(|i| unit(r, i) * scale - offset)); + } + for table in [fitted, saturating] { + let hue_sat = HueSat::new(&table).unwrap(); + for chunk in colours.chunks(16) { + let rgb: [[f32; 16]; 3] = std::array::from_fn(|c| std::array::from_fn(|l| chunk.get(l).map_or(0.0, |p| p[c]))); + let out = hue_sat.apply_lanes(rgb); + for (l, p) in chunk.iter().enumerate() { + let (got, want) = ([out[0][l], out[1][l], out[2][l]], hue_sat.apply(*p)); + assert_eq!(got.map(f32::to_bits), want.map(f32::to_bits), "{p:?}: {got:?} vs {want:?}"); + } + } + } + } + #[test] fn hue_sat_restore_weight_blends_smoothly_between_the_thresholds() { assert_eq!(restore_weight(0.0), 0.0); diff --git a/crates/engine/src/files.rs b/crates/engine/src/files.rs index 92bd0aa77..350f0e12d 100644 --- a/crates/engine/src/files.rs +++ b/crates/engine/src/files.rs @@ -267,19 +267,43 @@ impl CameraColourModel<'_> { /// White balance, matrix, tables and baseline exposure, in place. pub fn apply(&self, img: &mut Rgb32f) { let (wb, m, gain) = (self.wb, self.matrix.to_f32(), self.gain); - img.map_in_place(|p| { - let c = [p[0] * wb[0], p[1] * wb[1], p[2] * wb[2]]; - let rgb = [ - m[0][0] * c[0] + m[0][1] * c[1] + m[0][2] * c[2], - m[1][0] * c[0] + m[1][1] * c[1] + m[1][2] * c[2], - m[2][0] * c[0] + m[2][1] * c[1] + m[2][2] * c[2], - ]; - let rgb = match self.tables { - Some(tables) => tables.apply(rgb, gain), - None => rgb.map(|v| v * gain), - }; - self.hue_sat.map_or(rgb, |h| h.apply(rgb)).map(|v| v.max(0.0)) - }); + if let (None, Some(hue_sat)) = (self.tables, self.hue_sat) { + // a fitted camera look: `pixel`, `LANES` pixels at a time so the table step vectorises (same bits) + const LANES: usize = 16; + lightcraft_raster::par_rows(&mut img.data, img.width, |_, row| { + for chunk in row.chunks_mut(LANES) { + let mut rgb = [[0.0f32; LANES]; 3]; + for (l, p) in chunk.iter().enumerate() { + let c = [p[0] * wb[0], p[1] * wb[1], p[2] * wb[2]]; + for (i, rgb) in rgb.iter_mut().enumerate() { + rgb[l] = (m[i][0] * c[0] + m[i][1] * c[1] + m[i][2] * c[2]) * gain; + } + } + let out = hue_sat.apply_lanes(rgb); + for (l, p) in chunk.iter_mut().enumerate() { + *p = [out[0][l], out[1][l], out[2][l]].map(|v| v.max(0.0)); + } + } + }); + return; + } + img.map_in_place(|p| self.pixel(&m, p)); + } + + /// [`Self::apply`] to one pixel; `m` is the matrix as `f32`. + fn pixel(&self, m: &[[f32; 3]; 3], p: [f32; 3]) -> [f32; 3] { + let (wb, gain) = (self.wb, self.gain); + let c = [p[0] * wb[0], p[1] * wb[1], p[2] * wb[2]]; + let rgb = [ + m[0][0] * c[0] + m[0][1] * c[1] + m[0][2] * c[2], + m[1][0] * c[0] + m[1][1] * c[1] + m[1][2] * c[2], + m[2][0] * c[0] + m[2][1] * c[1] + m[2][2] * c[2], + ]; + let rgb = match self.tables { + Some(tables) => tables.apply(rgb, gain), + None => rgb.map(|v| v * gain), + }; + self.hue_sat.map_or(rgb, |h| h.apply(rgb)).map(|v| v.max(0.0)) } } @@ -1809,6 +1833,31 @@ pub(crate) mod tests { lightcraft_color::Mat3([[0.42, 0.548, 0.04], [0.246, 0.728, 0.033], [-0.247, 0.679, 0.447]]) } + /// A fitted camera look applied in lanes, on rows that aren't a whole number of them, gives the bits of the + /// pixel-by-pixel path, awkward pixels (black, negative, huge, NaN, infinite) included. + #[test] + fn camera_look_in_lanes_matches_pixel_by_pixel() { + let data = (0..72 * 5 * 5).map(|i| [((i * 37) % 61) as f32 - 30.0, 0.6 + ((i * 11) % 13) as f32 / 13.0, 1.0]).collect(); + let table = lightcraft_raw::profile::HsvTable { hue_divisions: 72, sat_divisions: 5, val_divisions: 5, data, srgb_value: true }; + let hue_sat = crate::camera_preview::HueSat::new(&table).unwrap(); + let colour = CameraColourModel { wb: [2.1, 1.0, 1.5], matrix: yellowing_look(), gain: 1.3, tables: None, hue_sat: Some(&hue_sat) }; + let special = [0.0, -0.0, -0.01, 1e-30, 0.3, 1.0, 5.0, 1e30, f32::INFINITY, f32::NAN]; + let img = Rgb32f::from_fn(37, 5, |x, y| { + let k = y * 37 + x; + if k % 5 == 0 { + [special[k % 10], special[k / 10 % 10], special[k / 100 % 10]] + } else { + [(x as f32 * 0.37).sin().abs(), ((x + y) as f32 * 0.11).cos().abs(), y as f32 / 4.0 - 0.05] + } + }); + let mut lanes = img.clone(); + colour.apply(&mut lanes); + let m = colour.matrix.to_f32(); + for (p, got) in img.data.iter().zip(&lanes.data) { + assert_eq!(got.map(f32::to_bits), colour.pixel(&m, *p).map(f32::to_bits), "{p:?}"); + } + } + #[test] fn clipped_white_renders_neutral_through_the_colour_model() { use lightcraft_color::{LUMA_2020, Mat3}; diff --git a/crates/raw/src/demosaic/ahd.rs b/crates/raw/src/demosaic/ahd.rs index d97b56215..d9fc27d76 100644 --- a/crates/raw/src/demosaic/ahd.rs +++ b/crates/raw/src/demosaic/ahd.rs @@ -8,7 +8,10 @@ //! adaptive luminance / chrominance tolerances (ε_L, ε_C from the paper) — the homogeneity. //! 4. Sum homogeneity over a 3×3 window and pick the more homogeneous candidate (average on ties). //! -//! Work is done in overlapping tiles (bounded memory, rayon-parallel). +//! Work is done in overlapping tiles (bounded memory, rayon-parallel), one band of tiles at a time written +//! straight into the output. Each step runs only where the tile's output depends on it: the choice reads the +//! homogeneity 1 pixel around the tile, which reads Lab 2 pixels around, the candidates read the greens 3 +//! pixels around, and those read the mosaic 5 pixels around, inside the `MARGIN`. use super::Mosaic; use crate::Rgb32f; @@ -23,14 +26,19 @@ pub(crate) fn ahd(m: &Mosaic) -> Rgb32f { if w == 0 || h == 0 { return out; } - // one row of tiles per task, written straight into the picture (collecting the tiles first held a second copy) - out.data.par_chunks_mut(w * TILE).enumerate().for_each(|(ty, band)| { - let (y0, th) = (ty * TILE, band.len() / w); - for x0 in (0..w).step_by(TILE) { - let tw = TILE.min(w - x0); - let px = tile(m, x0, y0, tw, th); - for y in 0..th { - band[y * w + x0..y * w + x0 + tw].copy_from_slice(&px[y * tw..(y + 1) * tw]); + // one row of tiles per task, its tiles in parallel too, written straight into the picture + out.data.par_chunks_mut((w * TILE).max(1)).enumerate().for_each(|(ty, band)| { + let (y0, th) = (ty * TILE, band.len() / w.max(1)); + let tiles: Vec<(usize, usize, Vec<[f32; 3]>)> = (0..w.div_ceil(TILE)) + .into_par_iter() + .map(|tx| { + let (x0, tw) = (tx * TILE, TILE.min(w - tx * TILE)); + (x0, tw, tile(m, x0, y0, tw, th)) + }) + .collect(); + for (x0, tw, px) in tiles { + for (row, src) in band.chunks_exact_mut(w).zip(px.chunks_exact(tw)) { + row[x0..x0 + tw].copy_from_slice(src); } } }); @@ -57,110 +65,100 @@ fn tile(m: &Mosaic, x0: usize, y0: usize, tw: usize, th: usize) -> Vec<[f32; 3]> let (ox, oy) = (x0 as isize - MARGIN as isize, y0 as isize - MARGIN as isize); let mut mos = vec![0f32; n]; let mut col = vec![0u8; n]; - for y in 0..bh { - for x in 0..bw { - mos[y * bw + x] = m.get(ox + x as isize, oy + y as isize); - col[y * bw + x] = m.color(ox + x as isize, oy + y as isize); + for (y, (mrow, crow)) in mos.chunks_exact_mut(bw).zip(col.chunks_exact_mut(bw)).enumerate() { + let ry = super::reflect(oy + y as isize, m.h); + if ox >= 0 && ox as usize + bw <= m.w { + // inside the image: the row as it is, and its colours repeating with the pattern + let x0 = ry * m.w + ox as usize; + mrow.copy_from_slice(&m.data[x0..x0 + bw]); + let period: Vec = (0..m.cfa.width).map(|k| m.cfa.color_at(ox as usize + k, ry)).collect(); + for (c, &p) in crow.iter_mut().zip(period.iter().cycle()) { + *c = p; + } + } else { + for x in 0..bw { + mrow[x] = m.get(ox + x as isize, oy + y as isize); + crow[x] = m.color(ox + x as isize, oy + y as isize); + } } } - let cl = |x: isize, y: isize| -> usize { y.clamp(0, bh as isize - 1) as usize * bw + x.clamp(0, bw as isize - 1) as usize }; + // flat indices of the pixels `d` or more from the buffer's edge, row by row + let inner = |d: usize| (d..bh - d).flat_map(move |y| y * bw + d..y * bw + bw - d); // 1. directional greens let mut gh = vec![0f32; n]; let mut gv = vec![0f32; n]; - for y in 0..bh as isize { - for x in 0..bw as isize { - let i = cl(x, y); - let c = mos[i]; - if col[i] == 1 { - gh[i] = c; - gv[i] = c; - continue; - } - 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] = 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] = super::clamp_between(est, u, d); + for i in inner(MARGIN - 3) { + let c = mos[i]; + if col[i] == 1 { + gh[i] = c; + gv[i] = c; + continue; } + let (l, r) = (mos[i - 1], mos[i + 1]); + let est = (l + r) * 0.5 + (2.0 * c - mos[i - 2] - mos[i + 2]) * 0.25; + gh[i] = super::clamp_between(est, l, r); + let (u, d) = (mos[i - bw], mos[i + bw]); + let est = (u + d) * 0.5 + (2.0 * c - mos[i - 2 * bw] - mos[i + 2 * bw]) * 0.25; + gv[i] = super::clamp_between(est, u, d); } - // 2. full candidates - let complete = |g: &[f32]| -> Vec<[f32; 3]> { + // 2. full candidates, and 3. their Lab + // (`p` with channel `c` set to `v`, in registers: writing channels through an index and reading the pixel + // back for its Lab stalled on store forwarding) + let set = |p: [f32; 3], c: u8, v: f32| [if c == 0 { v } else { p[0] }, if c == 1 { v } else { p[1] }, if c >= 2 { v } else { p[2] }]; + let complete = |g: &[f32]| -> (Vec<[f32; 3]>, Vec<[f32; 3]>) { let mut rgb = vec![[0f32; 3]; n]; - for y in 0..bh as isize { - for x in 0..bw as isize { - let i = cl(x, y); - let own = col[i] as usize; - let gvv = g[i]; - let d = |dx: isize, dy: isize| { - let j = cl(x + dx, y + dy); - mos[j] - g[j] - }; - let p = &mut rgb[i]; - p[1] = gvv; - if own == 1 { - let ch = col[cl(x + 1, y)] as usize; - let cv = col[cl(x, y + 1)] as usize; - p[ch] = gvv + (d(-1, 0) + d(1, 0)) * 0.5; - p[cv] = gvv + (d(0, -1) + d(0, 1)) * 0.5; - } else { - p[own] = mos[i]; - p[2 - own] = gvv + (d(-1, -1) + d(1, -1) + d(-1, 1) + d(1, 1)) * 0.25; - } - } + let mut lab_of = vec![[0f32; 3]; n]; + for i in inner(MARGIN - 2) { + let own = col[i]; + let gvv = g[i]; + let d = |j: usize| mos[j] - g[j]; + let p = if own == 1 { + let p = set([0.0, gvv, 0.0], col[i + 1], gvv + (d(i - 1) + d(i + 1)) * 0.5); + set(p, col[i + bw], gvv + (d(i - bw) + d(i + bw)) * 0.5) + } else { + let p = set([0.0, gvv, 0.0], own, mos[i]); + set(p, 2 - own, gvv + (d(i - bw - 1) + d(i - bw + 1) + d(i + bw - 1) + d(i + bw + 1)) * 0.25) + }; + rgb[i] = p; + lab_of[i] = lab(p); } - rgb + (rgb, lab_of) }; - let rh = complete(&gh); - let rv = complete(&gv); - // 3. Lab + homogeneity - let lh: Vec<[f32; 3]> = rh.iter().map(|&p| lab(p)).collect(); - let lv: Vec<[f32; 3]> = rv.iter().map(|&p| lab(p)).collect(); + let (rh, lh) = complete(&gh); + let (rv, lv) = complete(&gv); + // homogeneity let dl = |a: &[f32; 3], b: &[f32; 3]| (a[0] - b[0]).abs(); let dc = |a: &[f32; 3], b: &[f32; 3]| (a[1] - b[1]).powi(2) + (a[2] - b[2]).powi(2); let mut hh = vec![0u8; n]; let mut hv = vec![0u8; n]; - for y in 1..bh as isize - 1 { - for x in 1..bw as isize - 1 { - let i = cl(x, y); - let (l, r, u, d) = (cl(x - 1, y), cl(x + 1, y), cl(x, y - 1), cl(x, y + 1)); - let eps_l = dl(&lh[i], &lh[l]).max(dl(&lh[i], &lh[r])).min(dl(&lv[i], &lv[u]).max(dl(&lv[i], &lv[d]))); - let eps_c = dc(&lh[i], &lh[l]).max(dc(&lh[i], &lh[r])).min(dc(&lv[i], &lv[u]).max(dc(&lv[i], &lv[d]))); - let mut ch = 0; - let mut cv = 0; - for j in [l, r, u, d] { - if dl(&lh[i], &lh[j]) <= eps_l && dc(&lh[i], &lh[j]) <= eps_c { - ch += 1; - } - if dl(&lv[i], &lv[j]) <= eps_l && dc(&lv[i], &lv[j]) <= eps_c { - cv += 1; - } - } - hh[i] = ch; - hv[i] = cv; - } + for i in inner(MARGIN - 1) { + // left, right, up, down + let (l, r, u, d) = (i - 1, i + 1, i - bw, i + bw); + let (a, b) = (lh[i], lv[i]); + let dlh = [dl(&a, &lh[l]), dl(&a, &lh[r]), dl(&a, &lh[u]), dl(&a, &lh[d])]; + let dch = [dc(&a, &lh[l]), dc(&a, &lh[r]), dc(&a, &lh[u]), dc(&a, &lh[d])]; + let dlv = [dl(&b, &lv[l]), dl(&b, &lv[r]), dl(&b, &lv[u]), dl(&b, &lv[d])]; + let dcv = [dc(&b, &lv[l]), dc(&b, &lv[r]), dc(&b, &lv[u]), dc(&b, &lv[d])]; + let eps_l = dlh[0].max(dlh[1]).min(dlv[2].max(dlv[3])); + let eps_c = dch[0].max(dch[1]).min(dcv[2].max(dcv[3])); + hh[i] = (0..4).map(|k| u8::from(dlh[k] <= eps_l && dch[k] <= eps_c)).sum(); + hv[i] = (0..4).map(|k| u8::from(dlv[k] <= eps_l && dcv[k] <= eps_c)).sum(); } // 4. choose let mut out = Vec::with_capacity(tw * th); - for y in MARGIN..MARGIN + th { - for x in MARGIN..MARGIN + tw { - let (mut sh, mut sv) = (0u32, 0u32); - for dy in -1isize..=1 { - for dx in -1isize..=1 { - let j = cl(x as isize + dx, y as isize + dy); - sh += hh[j] as u32; - sv += hv[j] as u32; - } - } - let i = y * bw + x; - out.push(if sh > sv { - rh[i] - } else if sv > sh { - rv[i] - } else { - [(rh[i][0] + rv[i][0]) * 0.5, (rh[i][1] + rv[i][1]) * 0.5, (rh[i][2] + rv[i][2]) * 0.5] - }); + for i in inner(MARGIN) { + let (mut sh, mut sv) = (0u32, 0u32); + for j in [i - bw - 1, i - bw, i - bw + 1, i - 1, i, i + 1, i + bw - 1, i + bw, i + bw + 1] { + sh += hh[j] as u32; + sv += hv[j] as u32; } + out.push(if sh > sv { + rh[i] + } else if sv > sh { + rv[i] + } else { + [(rh[i][0] + rv[i][0]) * 0.5, (rh[i][1] + rv[i][1]) * 0.5, (rh[i][2] + rv[i][2]) * 0.5] + }); } out } diff --git a/crates/raw/src/llvc.rs b/crates/raw/src/llvc.rs index 0e2cc39cc..17408c760 100644 --- a/crates/raw/src/llvc.rs +++ b/crates/raw/src/llvc.rs @@ -14,6 +14,7 @@ fn corrupt(why: &str) -> RawError { } /// One bit, failing when the reader has run past the `source_bits` of its half. +#[inline(always)] fn bit(bits: &mut Bits<'_>, source_bits: usize) -> Result { let b = bits.get(1); if bits.consumed_bits() > source_bits { @@ -24,12 +25,14 @@ fn bit(bits: &mut Bits<'_>, source_bits: usize) -> Result { /// Count zero bits up to and including the terminating 1 (the `0^n 1` codes), bounded by `MAX_DPT`. fn zeros_then_one(bits: &mut Bits<'_>, source_bits: usize) -> Result { - let mut n = 0; - while bit(bits, source_bits)? == 0 { - n += 1; - if n > MAX_DPT { - return Err(corrupt("DPT code exceeds MAX_DPT")); - } + // the next MAX_DPT + 1 bits at once (more zeros than that is corrupt either way) + let n = bits.peek(MAX_DPT + 1).leading_zeros() - (32 - (MAX_DPT + 1)); + bits.skip((n + 1).min(MAX_DPT + 1)); + if bits.consumed_bits() > source_bits { + return Err(corrupt("past the end of its half")); + } + if n > MAX_DPT { + return Err(corrupt("DPT code exceeds MAX_DPT")); } Ok(n) } @@ -42,7 +45,7 @@ pub(crate) fn vld_decode_line(bits: &mut Bits<'_>, source_bits: usize, width: us let sets = width.div_ceil(4); let (mut abs, mut dpt, mut cnt, mut dpts) = ([0u32; 4], 0u32, sets, 0u8); let (mut abs_v, mut abs_e, mut dpt_v, mut dpt_e, mut zr) = (false, false, true, false, 0u64); - let mut out: Vec = Vec::new(); + let mut out: Vec = Vec::with_capacity(sets * 4); loop { // DPT value decoding (Figure 6.6). if cnt == 0 { @@ -180,43 +183,80 @@ pub(crate) fn band_rows(n: usize, phase: u8) -> (usize, usize) { (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) -} - /// Saturating narrowing: hostile coefficients clamp instead of overflowing. fn sat(v: i64) -> i32 { v.clamp(i64::from(i32::MIN), i64::from(i32::MAX)) as i32 } -/// 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())); +/// Sample `i` (`w` values wide) of a band of `n` samples, clamped into range (the whole-sample symmetric extension); +/// `zeros` for an empty band. +#[inline] +fn band_sample<'a>(v: &'a [i32], w: usize, n: usize, i: isize, zeros: &'a [i32]) -> &'a [i32] { + if n == 0 { + return zeros; } - let even: Vec = l - .iter() - .enumerate() - .map(|(j, &lv)| sat(i64::from(lv) - ((i64::from(at(h, j as isize - 1 + p)) + i64::from(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 = sat(i64::from(hv) + ((i64::from(at(&even, k as isize - p)) + i64::from(at(&even, k as isize + 1 - p))) >> 1)); - } - Ok(x) + let i = i.clamp(0, n as isize - 1) as usize; + v.get(i * w..(i + 1) * w).unwrap_or(zeros) } -/// Gathers the lines of the `low`/`high` pair along one axis (columns when `vertical`), transforms them, and scatters back. +/// Inverse 5/3 lifting, RDD 34 §6.5, with X[-1] := X[1], X[n] := X[n-2] (see `band_rows` for the sizes), of a +/// signal whose samples are `w` values wide: a band's width for all its columns at once (row by row, as each column +/// alone would give; rows in parallel). One line (`w = 1`) goes through [`lift_row`]. +fn inverse_lift(l: &[i32], h: &[i32], w: usize, phase: u8) -> Result> { + use rayon::prelude::*; + let w = w.max(1); + let (nl, nh) = (l.len() / w, h.len() / w); + let (n, p) = (nl + nh, isize::from(phase)); + if band_rows(n, phase) != (nl, nh) || phase > 1 { + return Err(RawError::Corrupt("ARW6 wavelet: band sizes do not match".into())); + } + let zeros = vec![0; w]; + let mut even = vec![0; nl * w]; + even.par_chunks_exact_mut(w).zip(l.par_chunks_exact(w)).enumerate().for_each(|(j, (e, lv))| { + let (a, b) = (band_sample(h, w, nh, j as isize - 1 + p, &zeros), band_sample(h, w, nh, j as isize + p, &zeros)); + for (((e, &lv), &a), &b) in e.iter_mut().zip(lv).zip(a).zip(b) { + *e = sat(i64::from(lv) - ((i64::from(a) + i64::from(b) + 2) >> 2)); + } + }); + let mut out = vec![0; n * w]; + let lo = usize::from(phase); + out.par_chunks_exact_mut(w).enumerate().for_each(|(t, o)| { + if t % 2 == lo { + if let Some(e) = even.get((t - lo) / 2 * w..((t - lo) / 2 + 1) * w) { + o.copy_from_slice(e); + } + return; + } + let k = (t - (1 - lo)) / 2; + let Some(hv) = h.get(k * w..(k + 1) * w) else { return }; + let (a, b) = (band_sample(&even, w, nl, k as isize - p, &zeros), band_sample(&even, w, nl, k as isize + 1 - p, &zeros)); + for (((o, &hv), &a), &b) in o.iter_mut().zip(hv).zip(a).zip(b) { + *o = sat(i64::from(hv) + ((i64::from(a) + i64::from(b)) >> 1)); + } + }); + Ok(out) +} + +/// [`inverse_lift`] of one line (`w = 1`) into `out` (`l.len() + h.len()` values), without allocating: `even` is +/// scratch space. The caller has checked the band sizes. +fn lift_row(l: &[i32], h: &[i32], phase: u8, even: &mut Vec, out: &mut [i32]) { + let (nl, nh, p) = (l.len(), h.len(), isize::from(phase)); + let hs = |i: isize| if nh == 0 { 0 } else { h.get(i.clamp(0, nh as isize - 1) as usize).map_or(0, |&v| i64::from(v)) }; + even.clear(); + even.extend(l.iter().enumerate().map(|(j, &lv)| sat(i64::from(lv) - ((hs(j as isize - 1 + p) + hs(j as isize + p) + 2) >> 2)))); + let even = &*even; + let es = |i: isize| if nl == 0 { 0 } else { even.get(i.clamp(0, nl as isize - 1) as usize).map_or(0, |&v| i64::from(v)) }; + let (lo, hi) = (usize::from(phase), 1 - usize::from(phase)); + for (o, &e) in out.iter_mut().skip(lo).step_by(2).zip(even.iter()) { + *o = e; + } + for (k, (o, &hv)) in out.iter_mut().skip(hi).step_by(2).zip(h).enumerate() { + *o = sat(i64::from(hv) + ((es(k as isize - p) + es(k as isize + 1 - p)) >> 1)); + } +} + +/// Inverse transform of the `low`/`high` pair along one axis: down the columns when `vertical` (all at once, row by +/// row), else along each row. 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())); @@ -228,20 +268,22 @@ fn lines(low: &Plane, high: &Plane, phase: u8, vertical: bool) -> Result 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); - } + if across_l == 0 { + return Ok(Plane::zeros(w, h)); } + if vertical { + return Ok(Plane { width: w, height: h, data: inverse_lift(&low.data, &high.data, w, phase)? }); + } + if band_rows(n, phase) != (nl, nh) || phase > 1 { + return Err(RawError::Corrupt("ARW6 wavelet: band sizes do not match".into())); + } + use rayon::prelude::*; + let mut out = Plane::zeros(w, h); + out.data.par_chunks_exact_mut(w.max(1)).enumerate().for_each_init(Vec::new, |even, (i, r)| { + if let (Some(l), Some(hi)) = (low.row(i), high.row(i)) { + lift_row(l, hi, phase, even, r); + } + }); Ok(out) } @@ -448,6 +490,29 @@ mod tests { assert_eq!(inverse_53_1d_rows(&col(&[3, -4]), &col(&[10, 20, 30]), 1).unwrap().data, vec![5, -5, 9, -17, 13]); } #[test] + fn horizontal_pass_gives_the_vertical_passs_samples() { + // two implementations of one lifting: `lift_row` along rows, `inverse_lift` down columns + let mut s = 0x9e37_79b9_7f4a_7c15u64; + let mut next = || { + s ^= s << 13; + s ^= s >> 7; + s ^= s << 17; + (s % 8191) as i32 - 4095 + }; + for n in 1..12usize { + for phase in 0..2u8 { + let ((nl, nh), rows) = (band_rows(n, phase), 3); + let low = Plane { width: nl, height: rows, data: (0..nl * rows).map(|_| next()).collect() }; + let high = Plane { width: nh, height: rows, data: (0..nh * rows).map(|_| next()).collect() }; + let across = inverse_53_1d_cols(&low, &high, phase).unwrap(); + for r in 0..rows { + let down = inverse_53_1d_rows(&col(low.row(r).unwrap()), &col(high.row(r).unwrap()), phase).unwrap(); + assert_eq!(across.row(r).unwrap(), &down.data[..], "n {n}, phase {phase}, row {r}"); + } + } + } + } + #[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()); diff --git a/crates/raw/src/profile.rs b/crates/raw/src/profile.rs index 566c2b1c3..7ca6b0e38 100644 --- a/crates/raw/src/profile.rs +++ b/crates/raw/src/profile.rs @@ -123,6 +123,177 @@ impl HsvTable { }; hsv_to_rgb(hue + shift, s2, v2) } + + /// [`Self::apply`] to `N` colours at once (`rgb[channel][lane]`), bit for bit: the same operations, with the + /// branches on the data (largest channel, hue sector, wrap-around) turned into selects so the compiler + /// vectorises them. The table reads, `powf` and the rare `rem_euclid` run lane by lane. + pub fn apply_lanes(&self, rgb: [[f32; N]; 3]) -> [[f32; N]; 3] { + let [r, g, b] = rgb; + if !self.fits_lanes() { + let mut out = rgb; + for l in 0..N { + [out[0][l], out[1][l], out[2][l]] = self.apply([r[l], g[l], b[l]]); + } + return out; + } + // `active`: the colours the table applies to (`apply` returns black, NaN and infinite ones as they are) + let (mut max, mut active, mut first) = ([0.0f32; N], [true; N], [false; N]); + let (mut d, mut ratio, mut second, mut third) = ([0.0f32; N], [0.0f32; N], [0.0f32; N], [0.0f32; N]); + for l in 0..N { + let m = r[l].max(g[l]).max(b[l]); + active[l] = !(m.is_nan() | (m <= 0.0) | !m.is_finite()); + d[l] = m - r[l].min(g[l]).min(b[l]); + // hue_of, for each largest channel + ratio[l] = (g[l] - b[l]) / d[l]; + second[l] = (b[l] - r[l]) / d[l] + 2.0; + third[l] = (r[l] - g[l]) / d[l] + 4.0; + first[l] = active[l] & (d[l] > 0.0) & (m == r[l]); + max[l] = m; + } + let ratio = wrap_lanes(ratio, 6.0, first); + let (mut hue, mut sat) = ([0.0f32; N], [0.0f32; N]); + for l in 0..N { + let h = if max[l] == r[l] { + ratio[l] + } else if max[l] == g[l] { + second[l] + } else { + third[l] + }; + let h = h * 60.0; + hue[l] = if d[l] <= 0.0 { 0.0 } else { h }; + sat[l] = d[l] / max[l]; + } + let mut v_index = max; + if self.srgb_value { + for l in 0..N { + if active[l] { + v_index[l] = linear_to_srgb(max[l].min(1.0)); + } + } + } + let [shift, ss, vs] = self.lookup_lanes(&hue, &sat, &v_index, active); + let (mut s2, mut v2, mut curve) = ([0.0f32; N], [0.0f32; N], [false; N]); + for l in 0..N { + s2[l] = limited_select(sat[l], ss[l]); + v2[l] = if self.srgb_value { max[l] } else { limited_select(max[l], vs[l]) }; + curve[l] = self.srgb_value & active[l] & !((max[l] <= 1.0) & (vs[l] == 1.0)); + } + for l in 0..N { + if curve[l] { + v2[l] = if max[l] <= 1.0 { srgb_to_linear((v_index[l] * vs[l]).min(1.0)) } else { max[l] * srgb_to_linear(vs[l].min(1.0)).min(1.0) }; + } + } + // hsv_to_rgb + let hue = wrap_lanes(std::array::from_fn(|l| hue[l] + shift[l]), 360.0, active); + let mut out = rgb; + for l in 0..N { + let hp = hue[l] / 60.0; + // `(hp as u32).min(5) as f32` for every `hp`, NaN included + let step = |k: f32| if hp >= k { 1.0 } else { 0.0 }; + let sector = step(1.0) + step(2.0) + step(3.0) + step(4.0) + step(5.0); + let f = hp - sector; + let (s, v) = (s2[l], v2[l]); + let (p, q, t) = (v * (1.0 - s), v * (1.0 - s * f), v * (1.0 - s * (1.0 - f))); + let red = if (sector == 0.0) | (sector == 5.0) { + v + } else if sector == 1.0 { + q + } else if sector == 4.0 { + t + } else { + p + }; + let green = if sector == 0.0 { + t + } else if sector <= 2.0 { + v + } else if sector == 3.0 { + q + } else { + p + }; + let blue = if sector <= 1.0 { + p + } else if sector == 2.0 { + t + } else if sector == 5.0 { + q + } else { + v + }; + out[0][l] = if active[l] { red } else { r[l] }; + out[1][l] = if active[l] { green } else { g[l] }; + out[2][l] = if active[l] { blue } else { b[l] }; + } + out + } + + /// Whether [`Self::lookup_lanes`] can do this table's index arithmetic in `f32`: every index below 2²³. + fn fits_lanes(&self) -> bool { + let entries = self.hue_divisions.checked_mul(self.sat_divisions).and_then(|n| n.checked_mul(self.val_divisions.max(1))); + self.hue_divisions >= 1 && self.sat_divisions >= 1 && entries.is_some_and(|n| n <= 1 << 23) + } + + /// [`Self::lookup`] at `N` points (`(shift, saturation scale, value scale)[lane]`), bit for bit for the lanes + /// `wanted`, on a table that [`Self::fits_lanes`]. The indices are computed as `f32` (exact integers below 2²³: + /// `x as usize` is a saturating conversion that doesn't vectorise). + fn lookup_lanes(&self, h: &[f32; N], s: &[f32; N], v: &[f32; N], wanted: [bool; N]) -> [[f32; N]; 3] { + const EXACT: f32 = 8_388_608.0; // 2²³: ulp 1 from here + // `x as usize` for `0 ≤ x < 2²³`: rounded to an integer, minus 1 if it rounded up; NaN is 0 + let floor = |x: f32| { + let nearest = (x + EXACT) - EXACT; + let floor = if nearest > x { nearest - 1.0 } else { nearest }; + if x.is_nan() { 0.0 } else { floor } + }; + let (hn, sn, vn) = (self.hue_divisions as f32, self.sat_divisions as f32, self.val_divisions as f32); + let axis = |x: f32, n: f32| { + let f = x.clamp(0.0, 1.0) * (n - 1.0); + let i0 = floor(f).min((n - 2.0).max(0.0)); + let i1 = (i0 + 1.0).min(n - 1.0); + (i0, i1, (f - i0).clamp(0.0, 1.0)) + }; + let h = wrap_lanes(*h, 360.0, wanted); + // the corners' indices (v, h, s), s fastest, and the weights between them + let (mut at, mut t) = ([[0u32; N]; 8], [[0.0f32; N]; 3]); + for l in 0..N { + let hf = (h[l] / 360.0 * hn).clamp(0.0, hn); + let h0 = floor(hf).min(hn - 1.0); + // (h0 + 1) % hn + let h1 = if h0 + 1.0 == hn { 0.0 } else { h0 + 1.0 }; + let (s0, s1, ts) = axis(s[l], sn); + let (v0, v1, tv) = axis(v[l], vn); + let (v0, v1, tv) = if self.val_divisions > 1 { (v0, v1, tv) } else { (0.0, 0.0, 0.0) }; + for (k, (vi, hi)) in [(v0, h0), (v0, h1), (v1, h0), (v1, h1)].into_iter().enumerate() { + let row = (vi * hn + hi) * sn; + // an exact integer below 2²³ plus 2²³ has it as its mantissa + at[2 * k][l] = (row + s0 + EXACT).to_bits().wrapping_sub(EXACT.to_bits()); + at[2 * k + 1][l] = (row + s1 + EXACT).to_bits().wrapping_sub(EXACT.to_bits()); + } + [t[0][l], t[1][l], t[2][l]] = [(hf - h0).clamp(0.0, 1.0), ts, tv]; + } + let mut corner = [[[0.0f32; N]; 3]; 8]; + for l in 0..N { + for (corner, at) in corner.iter_mut().zip(&at) { + [corner[0][l], corner[1][l], corner[2][l]] = self.data.get(at[l] as usize).copied().unwrap_or([0.0, 1.0, 1.0]); + } + } + let lerp = |a: &[f32; N], b: &[f32; N], t: &[f32; N]| { + let mut out = [0.0f32; N]; + for l in 0..N { + out[l] = a[l] + (b[l] - a[l]) * t[l]; + } + out + }; + let [th, ts, tv] = &t; + let mut out = [[0.0f32; N]; 3]; + for (c, out) in out.iter_mut().enumerate() { + let plane = |k: usize| lerp(&lerp(&corner[k][c], &corner[k + 1][c], ts), &lerp(&corner[k + 2][c], &corner[k + 3][c], ts), th); + // `lookup`'s `v0 == v1` for every lane: `axis` gives `v1 = v0 + 1` whenever there are two values or more + *out = if self.val_divisions > 1 { lerp(&plane(0), &plane(4), tv) } else { plane(0) }; + } + out + } } /// `x · k` clipped to 1, for `x ≤ 1`; above 1, only scale factors below 1 apply (continuous at 1). @@ -145,6 +316,37 @@ fn wrap(x: f32, p: f32) -> f32 { } } +/// [`wrap`] of `N` values, bit for bit for the lanes `wanted`: branch-free but for the `rem_euclid` fallback. +#[inline] +fn wrap_lanes(x: [f32; N], p: f32, wanted: [bool; N]) -> [f32; N] { + let (mut out, mut slow) = ([0.0f32; N], [false; N]); + for l in 0..N { + let (first, below, second) = ((x[l] >= 0.0) & (x[l] < p), (x[l] > -p) & (x[l] < 0.0), (x[l] >= p) & (x[l] < 2.0 * p)); + let (up, down) = (x[l] + p, x[l] - p); + out[l] = if first { + x[l] + } else if below { + up + } else { + down + }; + slow[l] = wanted[l] & !(first | below | second); + } + for l in 0..N { + if slow[l] { + out[l] = x[l].rem_euclid(p); + } + } + out +} + +/// [`limited`] without a branch. +#[inline] +fn limited_select(x: f32, k: f32) -> f32 { + let (below, above) = ((x * k).min(1.0), x * k.min(1.0)); + if x <= 1.0 { below } else { above } +} + fn hue_of(c: [f32; 3], max: f32, d: f32) -> f32 { if d <= 0.0 { return 0.0; @@ -330,6 +532,81 @@ mod tests { } } + /// Every combination of awkward components (zero, negative, tiny, around 1, huge, NaN, ±infinity), and random + /// colours from black to above 1: grey, each channel largest, hues either side of the wrap-around at 0°. + fn awkward_colours() -> Vec<[f32; 3]> { + let special = [0.0, -0.0, 1e-30, -1e-30, 1e-45, 0.003_130_8, 0.5, -0.5, 1.0, 1.0f32.next_up(), 2.0, 1e30, -1e30]; + let special = [&special[..], &[f32::MAX, f32::INFINITY, f32::NEG_INFINITY, f32::NAN]].concat(); + let mut out = Vec::new(); + for a in &special { + for b in &special { + out.extend(special.iter().map(|c| [*a, *b, *c])); + } + } + let mut seed = 0x2545_f491_4f6c_dd1du64; + for _ in 0..40_000 { + seed ^= seed << 13; + seed ^= seed >> 7; + seed ^= seed << 17; + let scale = [1.0, 1e-4, 8.0, 1e-38][(seed >> 56) as usize & 3]; + let c: [f32; 3] = std::array::from_fn(|i| (((seed >> (16 * i)) & 0xffff) as f32 / 65535.0 * 1.6 - 0.2) * scale); + out.push(match seed >> 60 { + 0 => [c[0]; 3], + 1 => [c[0].abs() + 0.1, c[1], c[1].next_up()], + 2 => [c[0].abs() + 0.1, c[1], c[1].next_down()], + _ => c, + }); + } + out + } + + /// `apply_lanes` against `apply`, bit for bit, `N` colours at a time (the last batch padded). + fn assert_lanes_match(t: &HsvTable, colours: &[[f32; 3]]) { + for chunk in colours.chunks(N) { + let rgb: [[f32; N]; 3] = std::array::from_fn(|c| std::array::from_fn(|l| chunk.get(l).map_or(0.25, |p| p[c]))); + let out = t.apply_lanes(rgb); + for (l, p) in chunk.iter().enumerate() { + let (got, want) = ([out[0][l], out[1][l], out[2][l]], t.apply(*p)); + assert_eq!(got.map(f32::to_bits), want.map(f32::to_bits), "{p:?}: {got:?} vs {want:?}"); + } + } + } + + #[test] + fn lanes_match_one_by_one_bit_for_bit() { + let mut seed = 0x9e37_79b9_7f4a_7c15u64; + let mut table = |h: usize, s: usize, v: usize, srgb_value: bool, shift: f32, value_scale: bool| { + let data = (0..h * s * v) + .map(|_| { + seed ^= seed << 13; + seed ^= seed >> 7; + seed ^= seed << 17; + let x = |k: u64| ((seed >> k) & 0xffff) as f32 / 65535.0; + [(x(0) * 2.0 - 1.0) * shift, 0.5 + x(16), if value_scale { 0.5 + x(32) } else { 1.0 }] + }) + .collect(); + HsvTable { hue_divisions: h, sat_divisions: s, val_divisions: v, data, srgb_value } + }; + let colours = awkward_colours(); + // a fitted camera look's shape; with value scales; a 2D DNG map turning hues by up to 5 turns; one hue and + // one saturation + for t in [ + table(72, 5, 5, true, 30.0, false), + table(72, 5, 5, true, 30.0, true), + table(90, 30, 1, false, 1800.0, true), + table(1, 1, 3, true, 30.0, true), + ] { + assert_lanes_match::<16>(&t, &colours); + assert_lanes_match::<5>(&t, &colours); + } + // entries missing (the identity fallback); too big for the lanes' index arithmetic (one by one) + let short = HsvTable { data: vec![[10.0, 1.5, 0.8]; 7], ..identity(4, 3, 2) }; + let huge = HsvTable { data: Vec::new(), ..identity(4096, 4096, 1) }; + for t in [short, huge] { + assert_lanes_match::<16>(&t, &colours); + } + } + fn close(a: [f32; 3], b: [f32; 3], tol: f32) -> bool { a.iter().zip(&b).all(|(x, y)| (x - y).abs() <= tol) } diff --git a/crates/raw/src/vendor/arw6.rs b/crates/raw/src/vendor/arw6.rs index 46af909dd..9837f3e7b 100644 --- a/crates/raw/src/vendor/arw6.rs +++ b/crates/raw/src/vendor/arw6.rs @@ -329,9 +329,10 @@ type HalfLines = (Half, Vec>); fn decode_stream(tile: &[u8], st: &StreamRef, ntu: usize, w2: usize, layouts: &[HalfRows]) -> Result> { let halves = stream_halves(tile, st, ntu)?; let width = line_width(st.group, w2); + // halves are independent bit streams: decode them in parallel halves - .into_iter() - .zip(layouts) + .into_par_iter() + .zip(layouts.par_iter()) .map(|(half, rows)| { let end = half.start.checked_add(half.len).ok_or_else(|| corrupt("half position"))?; let bytes = tile.get(half.start..end).ok_or_else(|| corrupt("half outside the tile"))?; @@ -468,36 +469,36 @@ pub(crate) fn colour(p: &TilePlanes) -> Result> { } let mw = w2.checked_mul(2).ok_or_else(|| corrupt("tile size"))?; let total = len.checked_mul(4).ok_or_else(|| corrupt("tile size"))?; - let at = |q: &Plane, j: usize, x: usize| -> i64 { q.data.get(j * w2 + x).copied().unwrap_or(0) as i64 }; + // (rows as slices: `at` is 0 past a row's end, as never happens once the sizes above match) + let at = |r: &[i32], x: usize| i64::from(r.get(x).copied().unwrap_or(0)); let mut g1 = vec![0i64; len]; - for j in 0..h2 { - let jp = j.saturating_sub(1); - for x in 0..w2 { + g1.par_chunks_mut(w2.max(1)).enumerate().for_each(|(j, row)| { + let (Some(m), Some(rp), Some(rj)) = (p.m.row(j), p.res.row(j.saturating_sub(1)), p.res.row(j)) else { return }; + for (x, g) in row.iter_mut().enumerate() { let xn = (x + 1).min(w2 - 1); - let r = at(&p.res, jp, x) + at(&p.res, jp, xn) + at(&p.res, j, x) + at(&p.res, j, xn); - if let Some(g) = g1.get_mut(j * w2 + x) { - *g = at(&p.m, j, x) - ((r + 4) >> 3); - } + let r = at(rp, x) + at(rp, xn) + at(rj, x) + at(rj, xn); + *g = at(m, x) - ((r + 4) >> 3); } - } - let g = |j: usize, x: usize| g1.get(j * w2 + x).copied().unwrap_or(0); + }); + let g1_row = |j: usize| g1.get(j * w2..(j + 1) * w2).unwrap_or(&[]); + let g = |r: &[i64], x: usize| r.get(x).copied().unwrap_or(0); let mut out = vec![0u16; total]; - for j in 0..h2 { - let jn = (j + 1).min(h2 - 1); - for x in 0..w2 { + // each plane row gives two mosaic rows + out.par_chunks_mut((2 * mw).max(1)).enumerate().for_each(|(j, rows)| { + let (gj, gn) = (g1_row(j), g1_row((j + 1).min(h2 - 1))); + let (Some(res), Some(c1), Some(c2), Some((top, bottom))) = (p.res.row(j), p.c1.row(j), p.c2.row(j), rows.split_at_mut_checked(mw)) else { + return; + }; + for (x, (t, b)) in top.as_chunks_mut::<2>().0.iter_mut().zip(bottom.as_chunks_mut::<2>().0.iter_mut()).enumerate() { let xp = x.saturating_sub(1); - let g2 = ((g(j, xp) + g(j, x) + g(jn, xp) + g(jn, x)) >> 2) + at(&p.res, j, x); - let mean = (g(j, x).clamp(0, 4095) + g2.clamp(0, 4095)) >> 1; - let r = 2 * at(&p.c1, j, x) + mean; - let b = 2 * at(&p.c2, j, x) + mean; - let cells = [(0, r), (1, g(j, x)), (mw, g2), (mw + 1, b)]; - for (off, v) in cells { - if let Some(o) = out.get_mut(2 * j * mw + 2 * x + off) { - *o = curve(v.clamp(0, 4095) as i32); - } + let g2 = ((g(gj, xp) + g(gj, x) + g(gn, xp) + g(gn, x)) >> 2) + at(res, x); + let mean = (g(gj, x).clamp(0, 4095) + g2.clamp(0, 4095)) >> 1; + let cells = [2 * at(c1, x) + mean, g(gj, x), g2, 2 * at(c2, x) + mean]; + for (o, v) in t.iter_mut().chain(b.iter_mut()).zip(cells) { + *o = curve(v.clamp(0, 4095) as i32); } } - } + }); Ok(out) } diff --git a/docs/arw6-compression.md b/docs/arw6-compression.md index 3b65c1d37..334b651f1 100644 --- a/docs/arw6-compression.md +++ b/docs/arw6-compression.md @@ -171,4 +171,22 @@ Raw decode only, 4 CPUs, release build, best of 3 wall-clock (2026-10-10). Metho | ILCE-7RM6 APS-C, Compressed (HQ) | 29.7 MB | 6592x4372 | 0.35 s | | ILCE-7RM6 APS-C, Compressed | 30.9 MB | 6592x4372 | 0.34 s | -The 1.0 s budget for the full-frame file is met. Parallelised: rayon over the tiles, over the 13 streams of each tile and over its 3 components. +The 1.0 s budget for the full-frame file is met. Parallelised: rayon over the tiles, over the 13 streams of each tile and over its 3 components, and the colour reconstruction over rows. + +2026-10-12 performance pass (output unchanged, bit for bit): the vertical wavelet pass lifts whole rows instead of gathering each column, `0^n 1` codes are counted with one `leading_zeros`, and the colour reconstruction reads rows as slices, in parallel. Same method, best of 5, the previous and the new decoder run alternately on the same machine (4 CPUs, which was busier than on 2026-10-10), best of two rounds: + +| File | Before | After | +|---|---|---| +| ILCE-7RM6 FF, Compressed (HQ) | 0.56 s | 0.46 s | +| ILCE-7RM6 FF, Compressed | 0.54 s | 0.38 s | +| ILCE-7RM6 APS-C, Compressed (HQ) | 0.34 s | 0.23 s | +| ILCE-7RM6 APS-C, Compressed | 0.35 s | 0.24 s | + +2026-10-13 (output unchanged, bit for bit): the two halves of each stream, which are independent bit streams, decode in parallel; the vertical wavelet pass runs its rows in parallel; the horizontal pass lifts each row in place (scratch space reused per thread), rows in parallel, instead of allocating two vectors per row. Same method, best of 3, the previous and the new decoder run alternately, on another machine (Ryzen 7 9800X3D, a 4-CPU cgroup quota shared with other builds); CPU is the process CPU time of the decode (all threads): + +| File | Before | After | CPU before | CPU after | +|---|---|---|---|---| +| ILCE-7RM6 FF, Compressed (HQ) | 0.60 s | 0.48 s | 1.14 s | 0.87 s | +| ILCE-7RM6 FF, Compressed | 0.62 s | 0.48 s | 1.14 s | 0.89 s | +| ILCE-7RM6 APS-C, Compressed (HQ) | 0.34 s | 0.23 s | 0.50 s | 0.40 s | +| ILCE-7RM6 APS-C, Compressed | 0.31 s | 0.24 s | 0.51 s | 0.42 s |