M5.4b: faster raw path — bin previews/thumbnails straight from the mosaic

- 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 <noreply@anthropic.com>
This commit is contained in:
Brandon Thomas 2026-10-01 09:01:53 -04:00
parent d1e00bc3cc
commit ec6637f238
13 changed files with 561 additions and 75 deletions

1
Cargo.lock generated
View File

@ -2093,6 +2093,7 @@ dependencies = [
"lightcraft-raw",
"lightcraft-scenes",
"log",
"rayon",
"serde",
"serde_json",
"thiserror 2.0.21",

View File

@ -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]

View File

@ -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");

View File

@ -167,8 +167,29 @@ pub fn probe_bytes(name: &str, bytes: &[u8]) -> Result<ProbeInfo, String> {
})
}
/// 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<usize> {
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<String> = 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 }));
}

View File

@ -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<bool> = std::sync::OnceLock::new();
*ON.get_or_init(|| std::env::var_os("LIGHTCRAFT_PROFILE").is_some())
}

194
crates/raw/src/binned.rs Normal file
View File

@ -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<Option<Rgb32f>> {
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<Normalized> = 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<u8>, [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<u8>, [u32; 3]) {
let mut n = [0u32; 3];
let l: Vec<u8> = (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})");
}
}
}
}
}

View File

@ -13,24 +13,115 @@ use lightcraft_raster::par_rows;
const EPS: f32 = 1e-3;
fn green(m: &Mosaic) -> Vec<f32> {
/// 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<f32> {
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<f32> = 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<f32> {
}
}
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<f32> {
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<f32> = (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)));
}
}

View File

@ -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<Vec<f64>>) -> Option<Mat3> {

View File

@ -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<f32> = (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::<f32>() / 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
}

View File

@ -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;

View File

@ -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<RawImage> {

View File

@ -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)

View File

@ -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