From ec6637f238e3f738dd79b05423821dde291b9cbc Mon Sep 17 00:00:00 2001 From: Brandon Thomas Date: Thu, 1 Oct 2026 09:01:53 -0400 Subject: [PATCH] =?UTF-8?q?M5.4b:=20faster=20raw=20path=20=E2=80=94=20bin?= =?UTF-8?q?=20previews/thumbnails=20straight=20from=20the=20mosaic?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - lightcraft-raw: RawImage::develop_binned(k) — k×k superpixel binning of the CFA (Bayer k even, X-Trans k multiple of 3), black/white normalised row by row, default crop only, clipping preserved for highlight reconstruction. 2560 px previews of a 24 MP raw bin 2×, thumbnails 4–8×; exports / 1:1 still demosaic (AHD) at full size. - highlight::reconstruct and its pyramid fill are rayon-parallel (bit-exact). - X-Trans directional demosaic: precomputed per-phase neighbour offsets and no reflection in the interior, no per-pixel allocation (bit-exact, tested against the reference formulation): 2.9 s → 0.75 s CPU at 24 MP. - files::load_bytes runs on a rayon worker so its short parallel loops don't each wait for the pool to wake (dominant cost on a loaded machine); LIGHTCRAFT_PROFILE prints per-stage source timings; files::bin_factor picks the block size. - render_bench ONLY=source: per-file decode / normalise / demosaic / binning / load. - Fix: u64 overflow on hostile IFD sizes in DNG/ARW/NEF raw-IFD selection (proptest). load_bytes(2560) CPU time, ARW a7 III 24 MP: 2347 → 240 ms; X-T20: 3120 → 1220 ms. Co-Authored-By: Claude Opus 5.5 --- Cargo.lock | 1 + crates/engine/Cargo.toml | 1 + crates/engine/examples/render_bench.rs | 37 +++ crates/engine/src/files.rs | 60 ++++- crates/pipeline/src/lib.rs | 3 +- crates/raw/src/binned.rs | 194 +++++++++++++++ crates/raw/src/demosaic/xtrans.rs | 236 +++++++++++++++++-- crates/raw/src/dng.rs | 2 +- crates/raw/src/highlight.rs | 87 +++---- crates/raw/src/lib.rs | 4 +- crates/raw/src/vendor/arw.rs | 2 +- crates/raw/src/vendor/nef.rs | 2 +- crates/raw/tests/robust.proptest-regressions | 7 + 13 files changed, 561 insertions(+), 75 deletions(-) create mode 100644 crates/raw/src/binned.rs create mode 100644 crates/raw/tests/robust.proptest-regressions diff --git a/Cargo.lock b/Cargo.lock index 5bcb8ca4..111bb6e1 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -2093,6 +2093,7 @@ dependencies = [ "lightcraft-raw", "lightcraft-scenes", "log", + "rayon", "serde", "serde_json", "thiserror 2.0.21", diff --git a/crates/engine/Cargo.toml b/crates/engine/Cargo.toml index 59922a08..dcaab969 100644 --- a/crates/engine/Cargo.toml +++ b/crates/engine/Cargo.toml @@ -25,6 +25,7 @@ serde = { workspace = true } serde_json = { workspace = true } thiserror = { workspace = true } log = { workspace = true } +rayon = { workspace = true } ab_glyph = { workspace = true } [dev-dependencies] diff --git a/crates/engine/examples/render_bench.rs b/crates/engine/examples/render_bench.rs index 8e0027f8..c1cb9040 100644 --- a/crates/engine/examples/render_bench.rs +++ b/crates/engine/examples/render_bench.rs @@ -88,6 +88,43 @@ fn main() { let n: usize = std::env::var("N").ok().and_then(|v| v.parse().ok()).unwrap_or(5); let only = std::env::var("ONLY").unwrap_or_default(); let run = |name: &str| only.is_empty() || only.split(',').any(|o| name.contains(o)); + if only == "source" { + // opening a photo: where the time goes before the first sharp loupe (per file) + for path in &args { + let bytes = std::fs::read(path).expect("read"); + println!("{path}"); + let raw = lightcraft_raw::decode(&bytes).expect("decode"); + println!(" {}×{} cpp {} cfa {:?}", raw.width, raw.height, raw.cpp, raw.cfa.as_ref().map(|c| c.name())); + println!(" {:<36} {}", "raw decode:", best(n, || drop(lightcraft_raw::decode(&bytes).expect("decode")))); + println!(" {:<36} {}", "normalize:", best(n, || drop(raw.normalized().expect("norm")))); + let norm = raw.normalized().expect("norm"); + println!(" {:<36} {}", "demosaic AHD:", best(n, || drop(lightcraft_raw::demosaic(&norm, lightcraft_raw::Method::Ahd)))); + if let Some(k) = lightcraft_engine::files::bin_factor(&raw, 2560) { + println!(" {:<36} {}", format!("binned ×{k} (2560 preview):"), best(n, || drop(raw.develop_binned(k, 0.99)))); + let img = raw.develop_binned(k, 0.99).expect("bin").expect("binnable"); + let t = lightcraft_raw::color::camera_transform(&raw, lightcraft_raw::color::as_shot_white_xy(&raw)); + println!( + " {:<36} {}", + "highlight reconstruct (binned):", + best(n, || { + let mut i = img.clone(); + lightcraft_raw::highlight::reconstruct(&mut i, t.wb, 0.99); + }) + ); + println!(" {:<36} {}", "clone (binned):", best(n, || drop(img.clone()))); + println!(" {:<36} {}", "fit 2560 (binned):", best(n, || drop(fit(&img, 2560, 2560, Filter::Box)))); + } + println!(" {:<36} {}", "embedded preview (2560):", best(n, || drop(lightcraft_engine::files::load_embedded_preview(&bytes, 2560)))); + for edge in [512usize, 2560, usize::MAX] { + println!( + " {:<36} {}", + format!("load_bytes({edge}):"), + best(n, || drop(lightcraft_engine::files::load_bytes(&bytes, edge).expect("load"))) + ); + } + } + return; + } let (full, info) = match args.first() { Some(path) => { let bytes = std::fs::read(path).expect("read"); diff --git a/crates/engine/src/files.rs b/crates/engine/src/files.rs index f9d08522..cd1cc4b5 100644 --- a/crates/engine/src/files.rs +++ b/crates/engine/src/files.rs @@ -167,8 +167,29 @@ pub fn probe_bytes(name: &str, bytes: &[u8]) -> Result { }) } +/// Sensor clip level (normalised) for highlight reconstruction. +const HIGHLIGHT_CLIP: f32 = 0.99; + +/// The largest block size to bin a raw's mosaic by for a source of at most `max_edge` pixels: +/// the binned image must keep at least 90 % of `max_edge` (a 16 MP sensor still bins 2× for the +/// 2560 px preview). `None` = demosaic at full size. +pub fn bin_factor(raw: &lightcraft_raw::RawImage, max_edge: usize) -> Option { + let c = raw.crop.clipped(raw.active_area.width, raw.active_area.height); + let long = if c.width > 1 && c.height > 1 { c.width.max(c.height) } else { raw.active_area.width.max(raw.active_area.height) }; + let need = max_edge.saturating_mul(9) / 10; + [8usize, 6, 4, 3, 2].into_iter().find(|&k| long / k >= need.max(1) && raw.can_bin(k)) +} + /// Decode a file into a linear Rec.2020 image no larger than `max_edge`, oriented. +/// +/// Runs on a rayon worker: its many short parallel loops then start on the worker's own queue +/// instead of each one waking the pool from outside and waiting for it (which costs more than the +/// loops themselves when the machine is busy). pub fn load_bytes(bytes: &[u8], max_edge: usize) -> Result<(Rgb32f, SourceInfo), String> { + rayon::scope(|_| load_bytes_now(bytes, max_edge)) +} + +fn load_bytes_now(bytes: &[u8], max_edge: usize) -> Result<(Rgb32f, SourceInfo), String> { if lightcraft_raw::probe(bytes).is_some() { let mut raw = match lightcraft_raw::decode(bytes) { Ok(r) => r, @@ -181,11 +202,26 @@ pub fn load_bytes(bytes: &[u8], max_edge: usize) -> Result<(Rgb32f, SourceInfo), // Embedded lens corrections are applied by the pipeline ("Enable Profile Corrections"), not baked in. let lens = embedded_lens(&raw); raw.opcodes.list3.retain(|op| !is_lens_opcode(op)); - let method = if max_edge <= 600 { lightcraft_raw::Method::Bilinear } else { lightcraft_raw::Method::Ahd }; - let mut img = raw.develop(method).map_err(|e| e.to_string())?; + // Previews and thumbnails bin the mosaic straight to (about) the size they need; only + // larger levels (exports, 1:1) demosaic the whole sensor. + let t0 = std::time::Instant::now(); + let binned = match bin_factor(&raw, max_edge) { + Some(k) => raw.develop_binned(k, HIGHLIGHT_CLIP).map_err(|e| e.to_string())?, + None => None, + }; + let mut img = match binned { + Some(img) => img, + None => { + let method = if max_edge <= 600 { lightcraft_raw::Method::Bilinear } else { lightcraft_raw::Method::Ahd }; + raw.develop(method).map_err(|e| e.to_string())? + } + }; + let mut stages = vec![("develop", t0.elapsed())]; let xy = lightcraft_raw::color::as_shot_white_xy(&raw); let t = lightcraft_raw::color::camera_transform(&raw, xy); - lightcraft_raw::highlight::reconstruct(&mut img, t.wb, 0.99); + stages.push(("transform", t0.elapsed())); + lightcraft_raw::highlight::reconstruct(&mut img, t.wb, HIGHLIGHT_CLIP); + stages.push(("highlights", t0.elapsed())); let m = t.matrix.to_f32(); let gain = 2f32.powf(t.baseline_exposure as f32); let wb = t.wb; @@ -197,7 +233,23 @@ pub fn load_bytes(bytes: &[u8], max_edge: usize) -> Result<(Rgb32f, SourceInfo), ((m[2][0] * c[0] + m[2][1] * c[1] + m[2][2] * c[2]) * gain).max(0.0), ] }); - let img = fit(&img, max_edge, max_edge, Filter::Box).oriented(raw.orientation); + stages.push(("colour", t0.elapsed())); + let img = fit(&img, max_edge, max_edge, Filter::Box); + stages.push(("fit", t0.elapsed())); + let img = img.oriented(raw.orientation); + stages.push(("orient", t0.elapsed())); + if lightcraft_pipeline::profiling() { + let mut prev = std::time::Duration::ZERO; + let parts: Vec = stages + .iter() + .map(|(n, t)| { + let d = *t - prev; + prev = *t; + format!("{n} {:.1}", d.as_secs_f64() * 1e3) + }) + .collect(); + eprintln!("[profile] raw source {}×{} (max {max_edge}, ms after decode): {}", img.width, img.height, parts.join(", ")); + } let (temp, tint) = xy_to_temp_tint(xy); return Ok((img, SourceInfo { raw: true, as_shot_temp: temp.round(), as_shot_tint: tint.round(), lens })); } diff --git a/crates/pipeline/src/lib.rs b/crates/pipeline/src/lib.rs index 3ccee51b..4ba0d886 100644 --- a/crates/pipeline/src/lib.rs +++ b/crates/pipeline/src/lib.rs @@ -351,9 +351,8 @@ pub(crate) fn is_bw(s: &DevelopSettings) -> bool { s.treatment == Treatment::Bw || s.profile.id == "lc.mono" } -/// Parallel map over output rows with index. /// `LIGHTCRAFT_PROFILE` is set: print per-stage timings to stderr. -pub(crate) fn profiling() -> bool { +pub fn profiling() -> bool { static ON: std::sync::OnceLock = std::sync::OnceLock::new(); *ON.get_or_init(|| std::env::var_os("LIGHTCRAFT_PROFILE").is_some()) } diff --git a/crates/raw/src/binned.rs b/crates/raw/src/binned.rs new file mode 100644 index 00000000..fde6c2c0 --- /dev/null +++ b/crates/raw/src/binned.rs @@ -0,0 +1,194 @@ +//! Reduced-size development straight from the mosaic ("superpixel" binning). +//! +//! When the caller needs an image at most half (or a third, a quarter…) of the sensor size — a +//! 2560 px loupe preview or a grid thumbnail of a 24 MP raw — demosaicing the full mosaic and then +//! downscaling throws most of the work away. Instead every `k × k` block of CFA samples becomes one +//! output pixel whose channels are the means of the block's samples of that colour. Every block must +//! contain all three colours (Bayer: `k` even; X-Trans: `k` a multiple of 3). There is no +//! interpolation, so there are no demosaicing artefacts; the colour planes are offset by less than +//! one output pixel, which is invisible at these scales. +//! +//! Clipping is preserved for [`crate::highlight::reconstruct`]: a channel whose block contains a +//! sample at or above `clip` takes that block's maximum sample of the colour (instead of a mean that +//! would fall just below the clip level and escape reconstruction). + +use crate::{Normalized, RawData, RawImage, Result, Rgb32f}; +use rayon::prelude::*; + +impl RawImage { + /// Can `k × k` blocks be binned (single-plane CFA data, every `k × k` window containing R, G + /// and B; Bayer: `k` even, so that every block has the same layout)? + pub fn can_bin(&self, k: usize) -> bool { + let Some(cfa) = &self.cfa else { return false }; + if k < 2 || self.cpp != 1 || !cfa.valid() || (cfa.is_bayer() && !k.is_multiple_of(2)) { + return false; + } + // every window position relative to the pattern + (0..cfa.height).all(|j| (0..cfa.width).all(|i| block_layout(cfa, i, j, k).1.iter().all(|&n| n > 0))) + } + + /// Like [`RawImage::develop`] at `1/k` of the size (camera RGB, white = 1.0, default crop applied, + /// not oriented). `None` when the data cannot be binned ([`RawImage::can_bin`]) or carries an + /// `OpcodeList3` (whose operations are defined at full resolution); callers then demosaic. + pub fn develop_binned(&self, k: usize, clip: f32) -> Result> { + self.validate()?; + if !self.can_bin(k) || !self.opcodes.list3.is_empty() { + return Ok(None); + } + let Some(cfa) = &self.cfa else { return Ok(None) }; + let a = self.active_area; + // bin the default crop (relative to the active area) only + let c = self.crop.clipped(a.width, a.height); + let c = if c.width < k || c.height < k { crate::Rect::new(0, 0, a.width, a.height) } else { c }; + let (bw, bh) = (c.width / k, c.height / k); + if bw == 0 || bh == 0 { + return Ok(None); + } + // Opcode lists 1/2 need the normalised copy; plain files are normalised row by row. + let norm: Option = if self.opcodes.list1.is_empty() && self.opcodes.list2.is_empty() { None } else { Some(self.normalized()?) }; + let cfa = cfa.shifted(a.x, a.y); + let range = self.white_at(0) - self.black.mean(); + let scale = if range > 0.0 { 1.0 / range } else { 1.0 }; + let uniform_black = + self.black.delta_h.is_empty() && self.black.delta_v.is_empty() && self.black.values.iter().all(|&v| v == self.black.values[0]); + let black0 = self.black.at(0, 0, 0, 1); + // colour of every sample of a block, per pattern phase of the block's corner + let (pw, ph) = (cfa.width, cfa.height); + let layouts: Vec<(Vec, [f32; 3])> = (0..ph) + .flat_map(|j| (0..pw).map(move |i| (i, j))) + .map(|(i, j)| { + let (l, n) = block_layout(&cfa, i, j, k); + (l, n.map(|n| 1.0 / n.max(1) as f32)) + }) + .collect(); + let wlen = bw * k; + let mut out = Rgb32f::new(bw, bh); + out.data.par_chunks_mut(bw).enumerate().for_each(|(by, row)| { + // the block row's k sample rows, normalised + let mut rows = vec![0f32; k * wlen]; + for dy in 0..k { + let y = c.y + by * k + dy; + let dst = &mut rows[dy * wlen..(dy + 1) * wlen]; + match (&norm, &self.data) { + (Some(nm), _) => dst.copy_from_slice(&nm.data[y * nm.width + c.x..][..wlen]), + (None, RawData::U16(d)) if uniform_black => { + let src = &d[(a.y + y) * self.width + a.x + c.x..][..wlen]; + for (o, &v) in dst.iter_mut().zip(src) { + *o = (v as f32 - black0) * scale; + } + } + (None, data) => { + let base = (a.y + y) * self.width + a.x + c.x; + for (x, o) in dst.iter_mut().enumerate() { + *o = (data.get(base + x) - self.black.at(c.x + x, y, 0, 1)) * scale; + } + } + } + } + let py = (c.y + by * k) % ph; + for (bx, px) in row.iter_mut().enumerate() { + let (layout, inv) = &layouts[py * pw + (c.x + bx * k) % pw]; + let (mut sum, mut max) = ([0f32; 3], [f32::MIN; 3]); + for dy in 0..k { + let src = &rows[dy * wlen + bx * k..][..k]; + for (&v, &col) in src.iter().zip(&layout[dy * k..(dy + 1) * k]) { + let col = col as usize; + sum[col] += v; + max[col] = max[col].max(v); + } + } + for ch in 0..3 { + px[ch] = if max[ch] >= clip { max[ch] } else { sum[ch] * inv[ch] }; + } + } + }); + Ok(Some(out)) + } +} + +/// Colours of the `k × k` samples of a block whose corner sits at pattern position `(i, j)` +/// (row-major) and the count of each colour. +fn block_layout(cfa: &crate::Cfa, i: usize, j: usize, k: usize) -> (Vec, [u32; 3]) { + let mut n = [0u32; 3]; + let l: Vec = (0..k * k) + .map(|t| { + let c = cfa.color_at(i + t % k, j + t / k).min(2); + n[c as usize] += 1; + c + }) + .collect(); + (l, n) +} + +#[cfg(test)] +mod tests { + use crate::demosaic::mosaic_from_rgb; + use crate::*; + + fn raw_from(n: &Normalized, black: f32, white: f32) -> RawImage { + let data = n.data.iter().map(|&v| (v * (white - black) + black).round() as u16).collect(); + RawImage { + format: RawFormat::Dng, + width: n.width, + height: n.height, + cpp: 1, + data: RawData::U16(data), + cfa: n.cfa.clone(), + bits: 14, + black: BlackLevel::uniform(black), + white: vec![white], + active_area: Rect::new(0, 0, n.width, n.height), + crop: Rect::new(0, 0, n.width, n.height), + orientation: Orientation::Normal, + color: ColorData::default(), + wb_multipliers: None, + linearized: false, + opcodes: OpcodeLists::default(), + metadata: Metadata::default(), + } + } + + #[test] + fn binning_matches_block_means_and_keeps_clipping() { + let img = Rgb32f::from_fn(24, 18, |x, y| [0.1 + x as f32 * 0.01, 0.2 + y as f32 * 0.01, 0.3]); + let mut raw = raw_from(&mosaic_from_rgb(&img, &Cfa::bayer("GRBG").unwrap()), 512.0, 16383.0); + assert!(raw.can_bin(2) && raw.can_bin(4) && !raw.can_bin(3)); + let b = raw.develop_binned(2, 0.99).unwrap().unwrap(); + assert_eq!((b.width, b.height), (12, 9)); + let p = b.get(3, 2); + // red sample of block (3,2) under GRBG is at (7,4); blue at (6,5); greens average (6,4)+(7,5) + assert!((p[0] - img.get(7, 4)[0]).abs() < 1e-3, "{p:?}"); + assert!((p[2] - 0.3).abs() < 1e-3); + assert!((p[1] - (img.get(6, 4)[1] + img.get(7, 5)[1]) / 2.0).abs() < 1e-3); + // a single clipped green sample keeps the block's green clipped + if let RawData::U16(d) = &mut raw.data { + d[4 * 24 + 6] = 16383; + } + assert!(raw.develop_binned(2, 0.99).unwrap().unwrap().get(3, 2)[1] >= 0.99); + // default crop is honoured at the binned scale + raw.crop = Rect::new(2, 2, 20, 12); + let c = raw.develop_binned(2, 0.99).unwrap().unwrap(); + assert_eq!((c.width, c.height), (10, 6)); + // X-Trans: 3×3 blocks hold all colours, 2×2 don't + let x = raw_from(&mosaic_from_rgb(&img, &Cfa::xtrans()), 0.0, 1000.0); + assert!(x.can_bin(3) && x.can_bin(6) && !x.can_bin(2)); + let xb = x.develop_binned(3, 0.99).unwrap().unwrap(); + assert_eq!((xb.width, xb.height), (8, 6)); + assert!((xb.get(4, 3)[2] - 0.3).abs() < 2e-3); + } + + #[test] + fn binning_matches_demosaic_downscaled_on_smooth_images() { + let img = Rgb32f::from_fn(64, 48, |x, y| [0.2 + 0.3 * (x as f32 / 64.0), 0.4 - 0.2 * (y as f32 / 48.0), 0.25]); + let raw = raw_from(&mosaic_from_rgb(&img, &Cfa::bayer("RGGB").unwrap()), 0.0, 4095.0); + let b = raw.develop_binned(2, 0.99).unwrap().unwrap(); + let full = lightcraft_raster::resample::half(&raw.develop(Method::Ahd).unwrap()); + for y in 2..b.height - 2 { + for x in 2..b.width - 2 { + for c in 0..3 { + assert!((b.get(x, y)[c] - full.get(x, y)[c]).abs() < 0.01, "({x},{y})"); + } + } + } + } +} diff --git a/crates/raw/src/demosaic/xtrans.rs b/crates/raw/src/demosaic/xtrans.rs index c26f768b..2b50b66c 100644 --- a/crates/raw/src/demosaic/xtrans.rs +++ b/crates/raw/src/demosaic/xtrans.rs @@ -13,24 +13,115 @@ use lightcraft_raster::par_rows; const EPS: f32 = 1e-3; -fn green(m: &Mosaic) -> Vec { +/// Same-colour neighbour offsets for every pattern phase, in scan order (so sums run in the same +/// order as a plain 2-D loop and interior pixels match border pixels bit for bit). +struct Offsets { + pw: usize, + ph: usize, + /// Per phase: green offsets at radius 1, then those at radius 2 (ring order as the loop below). + green: Vec<(Vec<(isize, isize)>, Vec<(isize, isize)>)>, + /// Per phase and colour (0 = red, 1 = blue): 5×5 offsets of that colour. + rb: Vec<[Vec<(isize, isize)>; 2]>, +} + +impl Offsets { + fn new(m: &Mosaic) -> Offsets { + let (pw, ph) = (m.cfa.width.max(1), m.cfa.height.max(1)); + // phases are evaluated far from the borders (pattern coordinates wrap) + let col = |x: isize, y: isize| m.cfa.color_at(x.rem_euclid(pw as isize) as usize, y.rem_euclid(ph as isize) as usize); + let mut green = Vec::new(); + let mut rb = Vec::new(); + for py in 0..ph as isize { + for px in 0..pw as isize { + let ring = |r: isize| { + let mut v = Vec::new(); + for dy in -r..=r { + for dx in -r..=r { + if (dx != 0 || dy != 0) && col(px + dx, py + dy) == 1 { + v.push((dx, dy)); + } + } + } + v + }; + green.push((ring(1), ring(2))); + let of = |c: u8| { + let mut v = Vec::new(); + for dy in -2isize..=2 { + for dx in -2isize..=2 { + if col(px + dx, py + dy) == c { + v.push((dx, dy)); + } + } + } + v + }; + rb.push([of(0), of(2)]); + } + } + Offsets { pw, ph, green, rb } + } + #[inline] + fn phase(&self, x: usize, y: usize) -> usize { + (y % self.ph) * self.pw + x % self.pw + } +} + +/// Robust weighted mean of `vals` around their median. +#[inline] +fn robust_mean(vals: &[f32]) -> f32 { + let mut sorted = [0f32; 24]; + let sorted = &mut sorted[..vals.len()]; + sorted.copy_from_slice(vals); + sorted.sort_unstable_by(|a, b| a.total_cmp(b)); + let med = sorted[sorted.len() / 2]; + let (mut s, mut ws) = (0.0f32, 0.0f32); + for &v in vals { + let wt = 1.0 / (EPS + (v - med).abs()); + s += wt * v; + ws += wt; + } + s / ws +} + +fn green(m: &Mosaic, off: &Offsets) -> Vec { let (w, h) = (m.w, m.h); let mut g = vec![0f32; w * h]; par_rows(&mut g, w, |y, row| { - let y = y as isize; + let interior_y = y >= 2 && y + 2 < h; + let yi = y as isize; let mut vals: Vec = Vec::with_capacity(24); for (x, o) in row.iter_mut().enumerate() { + if interior_y && x >= 2 && x + 2 < w { + // fast path: precomputed offsets, no reflection + let ph = off.phase(x, y); + let i = y * w + x; + if m.cfa.color_at(x, y) == 1 { + *o = m.data[i]; + continue; + } + let (r1, r2) = &off.green[ph]; + let ring = if r1.is_empty() { r2 } else { r1 }; + if ring.is_empty() { + *o = m.data[i]; + continue; + } + vals.clear(); + vals.extend(ring.iter().map(|&(dx, dy)| m.data[(i as isize + dy * w as isize + dx) as usize])); + *o = robust_mean(&vals); + continue; + } let x = x as isize; - if m.color(x, y) == 1 { - *o = m.get(x, y); + if m.color(x, yi) == 1 { + *o = m.get(x, yi); continue; } vals.clear(); for r in [1isize, 2] { for dy in -r..=r { for dx in -r..=r { - if (dx != 0 || dy != 0) && m.color(x + dx, y + dy) == 1 { - vals.push(m.get(x + dx, y + dy)); + if (dx != 0 || dy != 0) && m.color(x + dx, yi + dy) == 1 { + vals.push(m.get(x + dx, yi + dy)); } } } @@ -39,19 +130,10 @@ fn green(m: &Mosaic) -> Vec { } } if vals.is_empty() { - *o = m.get(x, y); + *o = m.get(x, yi); continue; } - let mut sorted = vals.clone(); - sorted.sort_by(|a, b| a.total_cmp(b)); - let med = sorted[sorted.len() / 2]; - let (mut s, mut ws) = (0.0f32, 0.0f32); - for &v in &vals { - let wt = 1.0 / (EPS + (v - med).abs()); - s += wt * v; - ws += wt; - } - *o = s / ws; + *o = robust_mean(&vals); } }); g @@ -59,31 +141,57 @@ fn green(m: &Mosaic) -> Vec { pub(crate) fn directional(m: &Mosaic) -> Rgb32f { let (w, h) = (m.w, m.h); - let g = green(m); + let off = Offsets::new(m); + let g = green(m, &off); let gg = |x: isize, y: isize| g[super::reflect(y, h) * w + super::reflect(x, w)]; let mut out = Rgb32f::new(w, h); par_rows(&mut out.data, w, |y, row| { - let y = y as isize; + let interior_y = y >= 2 && y + 2 < h; + let yi = y as isize; for (x, px) in row.iter_mut().enumerate() { + if interior_y && x >= 2 && x + 2 < w { + let i = y * w + x; + let own = m.cfa.color_at(x, y) as usize; + let gp = g[i]; + px[1] = gp; + let lists = &off.rb[off.phase(x, y)]; + for (li, c) in [0usize, 2].into_iter().enumerate() { + if c == own { + px[c] = m.data[i]; + continue; + } + let (mut s, mut ws) = (0.0f32, 0.0f32); + for &(dx, dy) in &lists[li] { + let j = (i as isize + dy * w as isize + dx) as usize; + let gq = g[j]; + let d2 = (dx * dx + dy * dy) as f32; + let wt = 1.0 / (d2 * (EPS + (gq - gp).abs())); + s += wt * (m.data[j] - gq); + ws += wt; + } + px[c] = if ws > 0.0 { gp + s / ws } else { gp }; + } + continue; + } let x = x as isize; - let own = m.color(x, y) as usize; - let gp = gg(x, y); + let own = m.color(x, yi) as usize; + let gp = gg(x, yi); px[1] = gp; for c in [0usize, 2] { if c == own { - px[c] = m.get(x, y); + px[c] = m.get(x, yi); continue; } let (mut s, mut ws) = (0.0f32, 0.0f32); for dy in -2isize..=2 { for dx in -2isize..=2 { - if m.color(x + dx, y + dy) as usize != c { + if m.color(x + dx, yi + dy) as usize != c { continue; } - let gq = gg(x + dx, y + dy); + let gq = gg(x + dx, yi + dy); let d2 = (dx * dx + dy * dy) as f32; let wt = 1.0 / (d2 * (EPS + (gq - gp).abs())); - s += wt * (m.get(x + dx, y + dy) - gq); + s += wt * (m.get(x + dx, yi + dy) - gq); ws += wt; } } @@ -93,3 +201,81 @@ pub(crate) fn directional(m: &Mosaic) -> Rgb32f { }); out } + +#[cfg(test)] +mod tests { + use super::*; + use crate::Cfa; + + /// The straightforward per-pixel formulation (reflection at every access): the precomputed + /// interior path must match it bit for bit. + fn reference(m: &Mosaic) -> Rgb32f { + let (w, h) = (m.w, m.h); + let mut g = vec![0f32; w * h]; + for y in 0..h as isize { + for x in 0..w as isize { + let o = &mut g[y as usize * w + x as usize]; + if m.color(x, y) == 1 { + *o = m.get(x, y); + continue; + } + let mut vals = Vec::new(); + for r in [1isize, 2] { + for dy in -r..=r { + for dx in -r..=r { + if (dx != 0 || dy != 0) && m.color(x + dx, y + dy) == 1 { + vals.push(m.get(x + dx, y + dy)); + } + } + } + if !vals.is_empty() { + break; + } + } + *o = if vals.is_empty() { m.get(x, y) } else { robust_mean(&vals) }; + } + } + let gg = |x: isize, y: isize| g[super::super::reflect(y, h) * w + super::super::reflect(x, w)]; + Rgb32f::from_fn(w, h, |x, y| { + let (x, y) = (x as isize, y as isize); + let own = m.color(x, y) as usize; + let gp = gg(x, y); + let mut px = [0f32, gp, 0.0]; + for c in [0usize, 2] { + if c == own { + px[c] = m.get(x, y); + continue; + } + let (mut s, mut ws) = (0.0f32, 0.0f32); + for dy in -2isize..=2 { + for dx in -2isize..=2 { + if m.color(x + dx, y + dy) as usize == c { + let gq = gg(x + dx, y + dy); + let wt = 1.0 / ((dx * dx + dy * dy) as f32 * (EPS + (gq - gp).abs())); + s += wt * (m.get(x + dx, y + dy) - gq); + ws += wt; + } + } + } + px[c] = if ws > 0.0 { gp + s / ws } else { gp }; + } + px + }) + } + + #[test] + fn fast_path_is_bit_exact() { + let cfa = Cfa::xtrans().shifted(1, 2); + let (w, h) = (37, 29); + let mut seed = 12345u32; + let data: Vec = (0..w * h) + .map(|_| { + seed = seed.wrapping_mul(1_664_525).wrapping_add(1_013_904_223); + (seed >> 8) as f32 / (1u32 << 24) as f32 + }) + .collect(); + let m = Mosaic { w, h, data: &data, cfa: &cfa }; + let (a, b) = (directional(&m), reference(&m)); + assert!(a.data.iter().zip(&b.data).all(|(p, q)| p.map(f32::to_bits) == q.map(f32::to_bits))); + } +} diff --git a/crates/raw/src/dng.rs b/crates/raw/src/dng.rs index f25eb604..4ffb0666 100644 --- a/crates/raw/src/dng.rs +++ b/crates/raw/src/dng.rs @@ -14,7 +14,7 @@ pub(crate) fn raw_ifd(tiff: &Tiff) -> Option<&Ifd> { .into_iter() .filter(|i| i.u32(t::NEW_SUBFILE_TYPE).unwrap_or(0) == 0) .filter(|i| matches!(i.u16(t::PHOTOMETRIC), Some(photometric::CFA) | Some(photometric::LINEAR_RAW))) - .max_by_key(|i| i.u64(t::IMAGE_WIDTH).unwrap_or(0) * i.u64(t::IMAGE_LENGTH).unwrap_or(0)) + .max_by_key(|i| i.u64(t::IMAGE_WIDTH).unwrap_or(0).saturating_mul(i.u64(t::IMAGE_LENGTH).unwrap_or(0))) } fn mat3(v: Option>) -> Option { diff --git a/crates/raw/src/highlight.rs b/crates/raw/src/highlight.rs index 388174ee..40d78751 100644 --- a/crates/raw/src/highlight.rs +++ b/crates/raw/src/highlight.rs @@ -7,6 +7,7 @@ //! normalised-convolution fill). Fully clipped pixels become neutral at the brightest plausible level. use lightcraft_raster::Rgb32f; +use rayon::prelude::*; /// Clip every channel of `wb ⊙ img` at `min_c(wb_c) · clip` — i.e. at the lowest channel's clip level after WB — /// then divide the multipliers back out. `clip` is the sensor clip in normalised units (≈ 1.0). @@ -28,25 +29,28 @@ fn fill_invalid(w: usize, h: usize, values: &mut [[f32; 3]], valid: &mut [bool]) let (cw, ch) = (w.div_ceil(2), h.div_ceil(2)); let mut cv = vec![[0f32; 3]; cw * ch]; let mut cval = vec![false; cw * ch]; - for y in 0..ch { - for x in 0..cw { - let (mut s, mut n) = ([0f32; 3], 0); - for dy in 0..2 { - for dx in 0..2 { - let (sx, sy) = (2 * x + dx, 2 * y + dy); - if sx < w && sy < h && valid[sy * w + sx] { - for c in 0..3 { - s[c] += values[sy * w + sx][c]; + { + let (values, valid) = (&*values, &*valid); + cv.par_chunks_mut(cw).zip(cval.par_chunks_mut(cw)).enumerate().for_each(|(y, (vrow, okrow))| { + for x in 0..cw { + let (mut s, mut n) = ([0f32; 3], 0); + for dy in 0..2 { + for dx in 0..2 { + let (sx, sy) = (2 * x + dx, 2 * y + dy); + if sx < w && sy < h && valid[sy * w + sx] { + for c in 0..3 { + s[c] += values[sy * w + sx][c]; + } + n += 1; } - n += 1; } } + if n > 0 { + vrow[x] = s.map(|v| v / n as f32); + okrow[x] = true; + } } - if n > 0 { - cv[y * cw + x] = s.map(|v| v / n as f32); - cval[y * cw + x] = true; - } - } + }); } if cw * ch < w * h { fill_invalid(cw, ch, &mut cv, &mut cval); @@ -54,10 +58,9 @@ fn fill_invalid(w: usize, h: usize, values: &mut [[f32; 3]], valid: &mut [bool]) // cannot shrink further (1×1): nothing valid anywhere handled above return; } - for y in 0..h { + values.par_chunks_mut(w).zip(valid.par_chunks_mut(w)).enumerate().for_each(|(y, (vrow, okrow))| { for x in 0..w { - let i = y * w + x; - if !valid[i] { + if !okrow[x] { // bilinear upsample of the coarse level let fx = ((x as f32 + 0.5) / 2.0 - 0.5).clamp(0.0, (cw - 1) as f32); let fy = ((y as f32 + 0.5) / 2.0 - 0.5).clamp(0.0, (ch - 1) as f32); @@ -71,11 +74,11 @@ fn fill_invalid(w: usize, h: usize, values: &mut [[f32; 3]], valid: &mut [bool]) let bot = g(x0, y1)[c] + (g(x1, y1)[c] - g(x0, y1)[c]) * tx; v[c] = top + (bot - top) * ty; } - values[i] = v; - valid[i] = true; + vrow[x] = v; + okrow[x] = true; } } - } + }); } /// Reconstruct partially clipped channels. `clip` is the sensor clip level in normalised units (use slightly @@ -84,52 +87,56 @@ fn fill_invalid(w: usize, h: usize, values: &mut [[f32; 3]], valid: &mut [bool]) pub fn reconstruct(img: &mut Rgb32f, wb: [f32; 3], clip: f32) -> usize { let (w, h) = (img.width, img.height); let n = w * h; - let clipped: Vec<[bool; 3]> = img.data.iter().map(|p| [p[0] >= clip, p[1] >= clip, p[2] >= clip]).collect(); - let count = clipped.iter().filter(|c| c.iter().any(|&b| b)).count(); + let is_clipped = |p: &[f32; 3]| p[0] >= clip || p[1] >= clip || p[2] >= clip; + let count = img.data.par_iter().filter(|p| is_clipped(p)).count(); if count == 0 { return 0; } // chromaticity (white-balanced ratios to the channel mean) of unclipped pixels let mut chroma = vec![[1f32 / 3.0; 3]; n]; let mut valid = vec![false; n]; - for i in 0..n { - if clipped[i].iter().any(|&b| b) { - continue; + chroma.par_iter_mut().zip(valid.par_iter_mut()).zip(img.data.par_iter()).for_each(|((r, v), p)| { + if is_clipped(p) { + return; } - let p = img.data[i]; let q = [p[0] * wb[0], p[1] * wb[1], p[2] * wb[2]]; let s = q[0] + q[1] + q[2]; if s > 1e-4 { - chroma[i] = q.map(|v| v.max(0.0) / s); - valid[i] = true; + *r = q.map(|v| v.max(0.0) / s); + *v = true; } - } + }); fill_invalid(w, h, &mut chroma, &mut valid); let max_level = wb.iter().cloned().fold(0.0f32, f32::max) * clip; - for i in 0..n { - let cl = clipped[i]; + img.data.par_iter_mut().zip(chroma.par_iter()).for_each(|(px, r)| { + let p = *px; + let cl = [p[0] >= clip, p[1] >= clip, p[2] >= clip]; if !cl.iter().any(|&b| b) { - continue; + return; } - let p = img.data[i]; let q = [p[0] * wb[0], p[1] * wb[1], p[2] * wb[2]]; - let r = chroma[i]; // estimate the white-balanced channel sum from the unclipped channels - let ests: Vec = (0..3).filter(|&c| !cl[c] && r[c] > 1e-3).map(|c| q[c] / r[c]).collect(); + let (mut sum, mut k) = (0f32, 0usize); + for c in 0..3 { + if !cl[c] && r[c] > 1e-3 { + sum += q[c] / r[c]; + k += 1; + } + } let mut out = q; - if ests.is_empty() { + if k == 0 { let v = q.iter().cloned().fold(0.0f32, f32::max).max(max_level); out = [v; 3]; } else { - let sum = ests.iter().sum::() / ests.len() as f32; + let sum = sum / k as f32; for c in 0..3 { if cl[c] { out[c] = q[c].max(sum * r[c]); } } } - img.data[i] = [out[0] / wb[0], out[1] / wb[1], out[2] / wb[2]]; - } + *px = [out[0] / wb[0], out[1] / wb[1], out[2] / wb[2]]; + }); count } diff --git a/crates/raw/src/lib.rs b/crates/raw/src/lib.rs index 2803d4b5..453c7523 100644 --- a/crates/raw/src/lib.rs +++ b/crates/raw/src/lib.rs @@ -5,7 +5,8 @@ //! opcode lists, [`Metadata`]); [`embedded_preview`] returns the largest embedded JPEG. //! - [`RawImage::normalized`] subtracts black, scales white to 1.0 and crops to the active area (applying DNG //! `OpcodeList1`/`OpcodeList2`); [`demosaic`] turns CFA data into camera-RGB [`Rgb32f`]; -//! [`RawImage::develop`] does all of it plus `OpcodeList3` and the default crop. +//! [`RawImage::develop`] does all of it plus `OpcodeList3` and the default crop; [`RawImage::develop_binned`] +//! produces the same at 1/k of the size straight from the mosaic (previews, thumbnails). //! - [`color`] implements the DNG colour model (dual-illuminant interpolation, forward matrices, white balance) //! and produces camera → linear Rec.2020 D65 matrices. //! @@ -18,6 +19,7 @@ //! colour matrix: [`color`] falls back to a documented neutral model. The decoders never panic on malformed input. #![forbid(unsafe_code)] +mod binned; pub mod color; pub mod demosaic; mod dng; diff --git a/crates/raw/src/vendor/arw.rs b/crates/raw/src/vendor/arw.rs index eb7c8e68..89020038 100644 --- a/crates/raw/src/vendor/arw.rs +++ b/crates/raw/src/vendor/arw.rs @@ -91,7 +91,7 @@ fn raw_ifd(tiff: &Tiff) -> Option<&Ifd> { tiff.all_ifds() .into_iter() .filter(|i| i.u16(t::PHOTOMETRIC) == Some(photometric::CFA) || i.contains(TONE_CURVE)) - .max_by_key(|i| i.u64(t::IMAGE_WIDTH).unwrap_or(0) * i.u64(t::IMAGE_LENGTH).unwrap_or(0)) + .max_by_key(|i| i.u64(t::IMAGE_WIDTH).unwrap_or(0).saturating_mul(i.u64(t::IMAGE_LENGTH).unwrap_or(0))) } pub(crate) fn decode(bytes: &[u8]) -> Result { diff --git a/crates/raw/src/vendor/nef.rs b/crates/raw/src/vendor/nef.rs index bcd8dd1a..11700bd1 100644 --- a/crates/raw/src/vendor/nef.rs +++ b/crates/raw/src/vendor/nef.rs @@ -22,7 +22,7 @@ fn raw_ifd(tiff: &Tiff) -> Option<&Ifd> { tiff.all_ifds() .into_iter() .filter(|i| i.u16(t::PHOTOMETRIC) == Some(photometric::CFA)) - .max_by_key(|i| i.u64(t::IMAGE_WIDTH).unwrap_or(0) * i.u64(t::IMAGE_LENGTH).unwrap_or(0)) + .max_by_key(|i| i.u64(t::IMAGE_WIDTH).unwrap_or(0).saturating_mul(i.u64(t::IMAGE_LENGTH).unwrap_or(0))) } /// Width without the optically masked columns some bodies append on the right: trailing columns (at most 64) diff --git a/crates/raw/tests/robust.proptest-regressions b/crates/raw/tests/robust.proptest-regressions new file mode 100644 index 00000000..21260e7d --- /dev/null +++ b/crates/raw/tests/robust.proptest-regressions @@ -0,0 +1,7 @@ +# Seeds for failure cases proptest has generated in the past. It is +# automatically read and these particular cases re-run before any +# novel cases are generated. +# +# It is recommended to check this file in to source control so that +# everyone who runs the test benefits from these saved cases. +cc 96fbbf63f3e7e0d186dc8aee69f5452640ec49e763149965b40707abbe37065c # shrinks to flips = [(16380191992317618766, 16), (5757657205609716137, 1)], cut = 11088808457807325322, lj = true