From 5f2f013c8eb32e5a02a75d3071b82a52027efb24 Mon Sep 17 00:00:00 2001 From: dev-null-undefined Date: Mon, 12 Oct 2026 18:34:52 +0200 Subject: [PATCH 1/5] Perf: a full-size raw decode no longer copies its 64 MP image to "fit" it MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `load_bytes_now` ran every development through `resample::fit`, which returns a clone when the image already fits. At full size (1:1 zoom, export) that cloned the whole picture: 797 MB and ~0.3 s for a 9984 × 6656 Sony ILCE-7RM6 raw. It now fits only when the image is larger, as the other preview paths in files.rs already do. Output unchanged (same pixels, by construction). Co-Authored-By: Claude Opus 5.5 --- crates/engine/src/files.rs | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/crates/engine/src/files.rs b/crates/engine/src/files.rs index dc79b8b31..feefb21a0 100644 --- a/crates/engine/src/files.rs +++ b/crates/engine/src/files.rs @@ -450,7 +450,8 @@ fn load_bytes_now( stages.push(("highlights", t0.elapsed())); colour.apply(&mut img); stages.push(("colour", t0.elapsed())); - let img = fit(&img, max_edge, max_edge, Filter::Box); + // (`fit` copies an image that already fits: 797 MB and 0.3 s for a full-size 64 MP decode) + let img = if img.width.max(img.height) > max_edge { fit(&img, max_edge, max_edge, Filter::Box) } else { img }; stages.push(("fit", t0.elapsed())); let img = img.into_oriented(raw.orientation); stages.push(("orient", t0.elapsed())); From f4078135690352e8591e08a6a7fd2f0b1fa896f4 Mon Sep 17 00:00:00 2001 From: dev-null-undefined Date: Mon, 12 Oct 2026 18:35:00 +0200 Subject: [PATCH 2/5] =?UTF-8?q?Perf:=20AHD=20demosaic=201.7=C3=97=20faster?= =?UTF-8?q?=20and=200.8=20GB=20lighter=20on=20a=2064=20MP=20raw,=20same=20?= =?UTF-8?q?pixels?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Three changes, each leaving the output bit for bit as it was: - Each step runs only where the tile's output depends on it (choice ±0, homogeneity ±1, Lab ±2, greens ±3, mosaic ±5 inside the 6 px margin), with direct neighbour indices instead of an index clamped on every read; the clamped values at the buffer edge never reached the output. - A candidate pixel is built in registers before its Lab conversion: writing its channels through an index and reading the pixel back stalled on store forwarding (the hottest instructions of the profile). - Tiles are written into the output one band at a time, instead of all 64 MP of tiles being collected and then copied (797 MB held twice). Tiles inside the image copy their mosaic rows instead of reflecting and taking a modulo per sample. Sony ILCE-7RM6 9984 × 6656 Compressed HQ raw, demosaic only, best of 5 on 4 cores: 1.85 s → 1.05 s wall, 6.3 s → 4.0 s CPU. Verified identical (f32 bits) on 7 real raws (ILCE-7RM6 FF/APS-C, ARW6 HQ/C and lossless, ILCE-7RM4A) and, against the previous code, on random mosaics of all four Bayer layouts at 1 × 1 … 300 × 7, 257 × 130 (tile edges, reflection). Co-Authored-By: Claude Opus 5.5 --- crates/raw/src/demosaic/ahd.rs | 193 ++++++++++++++++----------------- 1 file changed, 94 insertions(+), 99 deletions(-) diff --git a/crates/raw/src/demosaic/ahd.rs b/crates/raw/src/demosaic/ahd.rs index 3f17413ff..a75c395aa 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; @@ -19,20 +22,22 @@ const MARGIN: usize = 6; pub(crate) fn ahd(m: &Mosaic) -> Rgb32f { let (w, h) = (m.w, m.h); - let tiles: Vec<(usize, usize)> = (0..h.div_ceil(TILE)).flat_map(|ty| (0..w.div_ceil(TILE)).map(move |tx| (tx * TILE, ty * TILE))).collect(); - let results: Vec<(usize, usize, usize, usize, Vec<[f32; 3]>)> = tiles - .par_iter() - .map(|&(x0, y0)| { - let (tw, th) = (TILE.min(w - x0), TILE.min(h - y0)); - (x0, y0, tw, th, tile(m, x0, y0, tw, th)) - }) - .collect(); let mut out = Rgb32f::new(w, h); - for (x0, y0, tw, th, px) in results { - for y in 0..th { - out.data[(y0 + y) * w + x0..(y0 + y) * w + x0 + tw].copy_from_slice(&px[y * tw..(y + 1) * tw]); + 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); + } } - } + }); out } @@ -56,110 +61,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 } From 79ab54e167fe4cdd12425cfb0f1ea46faeb2c009 Mon Sep 17 00:00:00 2001 From: dev-null-undefined Date: Mon, 12 Oct 2026 18:59:47 +0200 Subject: [PATCH 3/5] Perf: ARW6 decode 20-33 % faster, same samples MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - The vertical inverse 5/3 lifts whole rows at once (`inverse_lift` with samples `w` values wide) instead of gathering each column through a strided iterator and scattering it back: the same integer arithmetic, sequential memory. The horizontal pass is the same function with w = 1. - `0^n 1` codes are counted with one `leading_zeros` over the next MAX_DPT + 1 bits instead of bit by bit (same bits consumed, same errors); each line reserves its values up front. - Colour reconstruction reads plane rows as slices and runs over rows in parallel (it was one thread per tile: 4 threads on a full-frame file whatever the machine). Best of 5 on 4 CPUs, previous and new decoder alternating: ILCE-7RM6 FF HQ 0.56 → 0.46 s, FF Compressed 0.54 → 0.38 s, APS-C HQ 0.34 → 0.23 s, APS-C Compressed 0.35 → 0.24 s (docs/arw6-compression.md → Timings). The decoded mosaic is identical (hash of every sample) on the 8 ARW6 corpus files and 2 private full-frame photos (ISO 100, ISO 12800); the pinned corpus sums still match. Co-Authored-By: Claude Opus 5.5 --- crates/raw/src/llvc.rs | 105 +++++++++++++++++++--------------- crates/raw/src/vendor/arw6.rs | 46 +++++++-------- docs/arw6-compression.md | 11 +++- 3 files changed, 93 insertions(+), 69 deletions(-) diff --git a/crates/raw/src/llvc.rs b/crates/raw/src/llvc.rs index 0e2cc39cc..ad386157c 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,56 @@ 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: `w = 1` for one line, a band's width for all its columns at once (row +/// by row, as each column alone would give). +fn inverse_lift(l: &[i32], h: &[i32], w: usize, phase: u8) -> Result> { + 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]; + for (j, (e, lv)) in even.chunks_exact_mut(w).zip(l.chunks_exact(w)).enumerate() { + 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, hi) = (usize::from(phase), 1 - usize::from(phase)); + for (o, e) in out.chunks_exact_mut(w).skip(lo).step_by(2).zip(even.chunks_exact(w)) { + o.copy_from_slice(e); + } + for (k, (o, hv)) in out.chunks_exact_mut(w).skip(hi).step_by(2).zip(h.chunks_exact(w)).enumerate() { + 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 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,17 +244,16 @@ 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) }; + 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)? }); + } 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) { + let x = inverse_lift(low.row(i).unwrap_or_default(), high.row(i).unwrap_or_default(), 1, phase)?; + if let Some(r) = out.data.get_mut(i * w..(i + 1) * w) { r.copy_from_slice(&x); } } diff --git a/crates/raw/src/vendor/arw6.rs b/crates/raw/src/vendor/arw6.rs index f6f535bd6..cb903a6f3 100644 --- a/crates/raw/src/vendor/arw6.rs +++ b/crates/raw/src/vendor/arw6.rs @@ -468,36 +468,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..a45f8365f 100644 --- a/docs/arw6-compression.md +++ b/docs/arw6-compression.md @@ -171,4 +171,13 @@ 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 | From bc1f138ee3e8beade33bd60e039a1f8f418d453b Mon Sep 17 00:00:00 2001 From: dev-null-undefined Date: Tue, 13 Oct 2026 01:22:34 +0200 Subject: [PATCH 4/5] Perf: ARW6 decode another 20-34 % faster, same samples MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - The two halves of each stream are independent bit streams: they now decode in parallel (`decode_stream`) instead of one after the other inside each of a tile's 13 streams. - The vertical inverse 5/3 runs its rows in parallel. The horizontal pass lifts each row in place (`lift_row`, scratch space reused per thread), rows in parallel, instead of allocating two vectors per row and copying the result back; a test checks it against the vertical pass, row by row. Best of 3 on 4 CPUs, previous and new decoder alternating: ILCE-7RM6 FF HQ 0.60 → 0.48 s (process CPU 1.14 → 0.87 s), FF Compressed 0.62 → 0.48 s, APS-C HQ 0.34 → 0.23 s, APS-C Compressed 0.31 → 0.24 s (docs/arw6-compression.md → Timings). The decoded mosaic is identical (hash of every sample) on the 8 ARW6 files of the RawDB ILCE-7RM6 set and 3 private A7R VI photos. Co-Authored-By: Claude --- crates/raw/src/llvc.rs | 82 ++++++++++++++++++++++++++++------- crates/raw/src/vendor/arw6.rs | 5 ++- docs/arw6-compression.md | 9 ++++ 3 files changed, 78 insertions(+), 18 deletions(-) diff --git a/crates/raw/src/llvc.rs b/crates/raw/src/llvc.rs index ad386157c..17408c760 100644 --- a/crates/raw/src/llvc.rs +++ b/crates/raw/src/llvc.rs @@ -200,9 +200,10 @@ fn band_sample<'a>(v: &'a [i32], w: usize, n: usize, i: isize, zeros: &'a [i32]) } /// 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: `w = 1` for one line, a band's width for all its columns at once (row -/// by row, as each column alone would give). +/// 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)); @@ -211,26 +212,49 @@ fn inverse_lift(l: &[i32], h: &[i32], w: usize, phase: u8) -> Result> { } let zeros = vec![0; w]; let mut even = vec![0; nl * w]; - for (j, (e, lv)) in even.chunks_exact_mut(w).zip(l.chunks_exact(w)).enumerate() { + 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, hi) = (usize::from(phase), 1 - usize::from(phase)); - for (o, e) in out.chunks_exact_mut(w).skip(lo).step_by(2).zip(even.chunks_exact(w)) { - o.copy_from_slice(e); - } - for (k, (o, hv)) in out.chunks_exact_mut(w).skip(hi).step_by(2).zip(h.chunks_exact(w)).enumerate() { + 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 { @@ -250,13 +274,16 @@ fn lines(low: &Plane, high: &Plane, phase: u8, vertical: bool) -> Result if vertical { return Ok(Plane { width: w, height: h, data: inverse_lift(&low.data, &high.data, w, phase)? }); } - let mut out = Plane::zeros(w, h); - for i in 0..across_l { - let x = inverse_lift(low.row(i).unwrap_or_default(), high.row(i).unwrap_or_default(), 1, phase)?; - if let Some(r) = out.data.get_mut(i * w..(i + 1) * w) { - r.copy_from_slice(&x); - } + 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) } @@ -463,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/vendor/arw6.rs b/crates/raw/src/vendor/arw6.rs index cb903a6f3..70f438190 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"))?; diff --git a/docs/arw6-compression.md b/docs/arw6-compression.md index a45f8365f..334b651f1 100644 --- a/docs/arw6-compression.md +++ b/docs/arw6-compression.md @@ -181,3 +181,12 @@ The 1.0 s budget for the full-frame file is met. Parallelised: rayon over the ti | 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 | From 6467d5432552ce86106773d5438e9adc72f0467f Mon Sep 17 00:00:00 2001 From: dev-null-undefined Date: Mon, 12 Oct 2026 20:43:12 +0200 Subject: [PATCH 5/5] Perf: camera-look colour step 1.5-1.9x faster on full-size raws, same pixels The hue/saturation table of a fitted camera look (ILCE-7RM4A, APS-C ILCE-7RM6, ...) ran HsvTable::apply pixel by pixel: branches on the data (largest channel, hue sector), saturating float-to-usize casts and bounds-checked table reads, ~550 instructions and 1.2 branch misses per pixel. CameraColourModel::apply now runs it 16 pixels at a time (HueSat::apply_lanes -> HsvTable::apply_lanes, structure of arrays): per pixel the same IEEE operations in the same order, with the branches turned into selects and the table indices computed as exact f32 integers, so the compiler vectorises everything but the table reads, powf and the rare rem_euclid, which run lane by lane. DNG profile tables and the other colour models keep the per-pixel path. Measured on a Core Ultra 9 185H P-core (busy machine, min of N): - HsvTable alone, 72x5x5 fitted-look-shaped table, random colours, one thread: 243 -> 106 cycles (57 -> 26 ns) per pixel. - Colour stage of a full-size load, one thread: 60 MP ILCE-7RM4A 3.80 -> 1.99 s, 26 MP APS-C ILCE-7RM6 HQ 1.88 -> 0.98 s; 4 CPUs: 1.02 -> 0.67 s, 0.56 -> 0.30 s. Exactness: new tests compare the lanes with the one-by-one path bit for bit (to_bits) on tens of thousands of colours (every combination of zero, -0, tiny, negative, around 1, huge, NaN and +-inf components; grey; hues at the 0/360 degree wrap; near-black noise in the restore's blend) for fitted, value-scaled, 2D, many-turn and degenerate tables, 16 and 5 lanes, and on an image whose rows aren't a whole number of lanes. Stage fingerprints (decode, AHD, 2560 px and full-size loads) of 7 Sony raws are identical to origin/main. Sources: own design. Co-Authored-By: Claude Opus 5.5 --- crates/engine/src/camera_preview.rs | 81 ++++++++ crates/engine/src/files.rs | 75 ++++++-- crates/raw/src/profile.rs | 277 ++++++++++++++++++++++++++++ 3 files changed, 420 insertions(+), 13 deletions(-) diff --git a/crates/engine/src/camera_preview.rs b/crates/engine/src/camera_preview.rs index 741fc391d..a3d13bfdc 100644 --- a/crates/engine/src/camera_preview.rs +++ b/crates/engine/src/camera_preview.rs @@ -1180,6 +1180,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. @@ -1951,6 +1987,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 feefb21a0..0940f3125 100644 --- a/crates/engine/src/files.rs +++ b/crates/engine/src/files.rs @@ -266,19 +266,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)) } } @@ -1647,6 +1671,31 @@ 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/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) }