From 9abf0d42103d1d9e1102f7771269f102ec8c54f3 Mon Sep 17 00:00:00 2001 From: Alex Chambers Date: Sat, 10 Oct 2026 00:39:09 +0200 Subject: [PATCH] LR-EDIT-LIGHT-HIGHLIGHTS: clipped highlights stay neutral when darkened (#523) Sensor-clipped pixels are rebuilt at raw decode. The tone map rolls them off to white at Neutral, so any colour they carry shows only when Highlights, a negative Exposure or Whites darkens them. They took colour from the nearest unclipped pixel (dark foliage, lens fringes along branches, pixels demosaiced from clipped samples), left two clipped channels at their white-balance-scaled clip levels (red or pink), and were neutral only in camera RGB, not after the colour model (a look fitted to the camera's JPEG, a camera profile). - highlight::reconstruct_with: chromaticity only from bright unclipped content away from edges and from the clip, through a normalised push-pull pyramid; a neutral estimate blended in where reliable content is scarce at every scale up to a quarter of the long edge, applied only as far as the pixel lies past the nearest near-clip unclipped pixels (so a surface still runs into clipping without a step); clipped pairs keep the estimate's ratio at their clip floors; partly clipped pixels fade into the clipped neutral on their way from the first channel clipping to the last; fully clipped pixels become the colour model's neutral. Pixels without a clipped channel are left exactly as they were. - engine: the decode colour step as CameraColourModel, with clipped_white, which keeps the camera's neutral for nearly singular or implausible models. - Tests in raw (continuity along a ramp through every clip, also as lines 1 to 16 rows thick with no reliable colour, on canvases 32 to 1024 rows), engine (pipeline renders at Highlights -100, #523's edit and Exposure -2 through the first clip, over the same sweep; a CC0 ILCE-7RM4 corpus file added to RAW_SAMPLES) and the GPU equivalence suite. - Render cache version 27. Co-Authored-By: Claude Opus 5.5 --- Cargo.lock | 1 + crates/engine/src/files.rs | 289 ++++++++++- crates/engine/src/media.rs | 2 +- crates/gpu/Cargo.toml | 1 + crates/gpu/tests/equivalence.rs | 53 ++ crates/raw/src/highlight.rs | 832 +++++++++++++++++++++----------- docs/parity.md | 2 +- xtask/src/main.rs | 6 + 8 files changed, 892 insertions(+), 294 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index 45d8ff15..e80e1c26 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -3181,6 +3181,7 @@ dependencies = [ "lightcraft-geom", "lightcraft-pipeline", "lightcraft-raster", + "lightcraft-raw", "lightcraft-scenes", "log", "pollster 1.0.1", diff --git a/crates/engine/src/files.rs b/crates/engine/src/files.rs index db762f00..52177ae5 100644 --- a/crates/engine/src/files.rs +++ b/crates/engine/src/files.rs @@ -200,6 +200,86 @@ fn preview_reason(bytes: &[u8], e: lightcraft_raw::RawError) -> Result { + /// White-balance multipliers (as shot). + pub wb: [f32; 3], + /// White-balanced camera RGB → linear Rec.2020: the camera matrix, or a look fitted to the camera's JPEG on top. + pub matrix: lightcraft_color::Mat3, + /// Baseline exposure, as a gain. + pub gain: f32, + /// A DNG profile's hue/saturation map and look table. + pub tables: Option<&'a lightcraft_raw::profile::ProfileTables>, + /// A fitted camera look's hue/saturation table. + pub hue_sat: Option<&'a crate::camera_preview::HueSat<'a>>, +} + +impl CameraColourModel<'_> { + /// Rebuild sensor-clipped channels; fully clipped highlights become [`clipped_white`] of the matrix, so they + /// come out of [`Self::apply`] neutral. + pub fn rebuild_highlights(&self, img: &mut Rgb32f) { + lightcraft_raw::highlight::reconstruct_with(img, self.wb, HIGHLIGHT_CLIP, clipped_white(&self.matrix)); + } + + /// 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)) + }); + } +} + +/// The white-balanced camera RGB that `model` (white-balanced camera RGB → linear Rec.2020) renders +/// neutral, as bright as it renders `[1, 1, 1]`: what fully clipped highlights are rebuilt to, so they come out +/// neutral whatever the colour model does to the camera's own neutral (a look fitted to the camera's JPEG or a +/// camera profile often tints it by a few per cent, which any darkening of the highlights would show). +/// `[1, 1, 1]` (the camera's own neutral, as before this existed) when the model can't be trusted with it: not +/// invertible or nearly so, an answer far from the camera's neutral (a sane colour model moves white by a few per +/// cent), or one that isn't neutral through the `f32` matrix the decode applies. +pub(crate) fn clipped_white(model: &lightcraft_color::Mat3) -> [f32; 3] { + const CAMERA: [f32; 3] = [1.0; 3]; + /// Largest condition number (max row sum norm) of a model trusted for this. + const CONDITION: f64 = 1e3; + /// Largest factor between the answer and the camera's neutral, per channel. + const SPREAD: f64 = 4.0; + let norm = |m: &lightcraft_color::Mat3| m.0.iter().map(|row| row.iter().map(|v| v.abs()).sum::()).fold(0.0, f64::max); + let Some(inverse) = model.inverse() else { return CAMERA }; + let condition = norm(model) * norm(&inverse); + if !condition.is_finite() || condition > CONDITION { + return CAMERA; + } + let luma = |v: [f64; 3]| { + let c = model.apply(v); + (0..3).map(|i| c[i] * f64::from(lightcraft_color::LUMA_2020[i])).sum::() + }; + let n = inverse.apply([1.0; 3]); + let (target, own) = (luma([1.0; 3]), luma(n)); + if !(target > 0.0 && own > 0.0 && target.is_finite() && own.is_finite()) { + return CAMERA; + } + let white = n.map(|v| v * target / own); + if !white.iter().all(|v| v.is_finite() && (1.0 / SPREAD..=SPREAD).contains(v)) { + return CAMERA; + } + let white = white.map(|v| v as f32); + // neutral through the matrix as the decode applies it + let m = model.to_f32(); + let out: [f32; 3] = std::array::from_fn(|i| m[i][0] * white[0] + m[i][1] * white[1] + m[i][2] * white[2]); + let (lo, hi) = (out[0].min(out[1]).min(out[2]), out[0].max(out[1]).max(out[2])); + if lo > 0.0 && hi.is_finite() && hi - lo <= 1e-3 * hi { white } else { CAMERA } +} + /// 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). X-Trans can only bin 3× (its 6×6 pattern), and its full demosaic is ~5× the @@ -302,32 +382,24 @@ fn load_bytes_now( let twin = denoise.and_then(|spec| denoised_camera_rgb(&raw, spec, factor)).filter(|d| d.width == img.width && d.height == img.height); // the samples aren't needed any more (the colour model below reads only the tags) raw.data = lightcraft_raw::RawData::U16(Vec::new()); - let wb = t.wb; - let m = camera_look.as_ref().map(|p| p.matrix.mul(&t.matrix)).unwrap_or(t.matrix).to_f32(); let hue_sat = camera_look.as_ref().and_then(|p| p.hue_sat.as_ref()).and_then(crate::camera_preview::HueSat::new); - let gain = 2f32.powf(t.baseline_exposure as f32); // A DNG's own profile look (hue/saturation map, look table), DNG spec chapter 6. let tables = lightcraft_raw::profile::ProfileTables::new(&raw.color.profile, lightcraft_raw::color::illuminant_weight(&raw.color, xy)); + let colour = CameraColourModel { + wb: t.wb, + matrix: camera_look.as_ref().map(|p| p.matrix.mul(&t.matrix)).unwrap_or(t.matrix), + gain: 2f32.powf(t.baseline_exposure as f32), + tables: tables.as_ref(), + hue_sat: hue_sat.as_ref(), + }; // camera RGB → what the pipeline takes: clipped highlights rebuilt, white balance and the colour model applied, // fitted to `max_edge` and upright (the same for the plain and the denoised picture) let finish = |mut img: Rgb32f| { let mut stages = vec![("develop", t0.elapsed())]; stages.push(("transform", t0.elapsed())); - lightcraft_raw::highlight::reconstruct(&mut img, wb, HIGHLIGHT_CLIP); + colour.rebuild_highlights(&mut img); stages.push(("highlights", t0.elapsed())); - 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 &tables { - Some(tables) => tables.apply(rgb, gain), - None => rgb.map(|v| v * gain), - }; - hue_sat.as_ref().map_or(rgb, |h| h.apply(rgb)).map(|v| v.max(0.0)) - }); + colour.apply(&mut img); stages.push(("colour", t0.elapsed())); let img = fit(&img, max_edge, max_edge, Filter::Box); stages.push(("fit", t0.elapsed())); @@ -868,6 +940,189 @@ mod tests { assert!((mean(&dim) / mean(&before) - 0.5).abs() < 0.02, "{} vs {}", mean(&dim), mean(&before)); } + /// A look fitted to a camera's JPEG, as `camera_preview` makes them: its rows don't sum to one, so the camera's + /// own white-balanced neutral comes out yellow (blue 12 % low). + fn yellowing_look() -> lightcraft_color::Mat3 { + lightcraft_color::Mat3([[0.42, 0.548, 0.04], [0.246, 0.728, 0.033], [-0.247, 0.679, 0.447]]) + } + + #[test] + fn clipped_white_renders_neutral_through_the_colour_model() { + use lightcraft_color::{LUMA_2020, Mat3}; + let look = yellowing_look(); + let white = clipped_white(&look); + let out = look.apply(white.map(f64::from)); + assert!(out.iter().all(|v| (v - out[1]).abs() < 1e-5), "{white:?} renders {out:?}"); + let luma = |c: [f64; 3]| (0..3).map(|i| c[i] * f64::from(LUMA_2020[i])).sum::(); + assert!((luma(out) - luma(look.apply([1.0; 3]))).abs() < 1e-5, "as bright as the camera's neutral"); + assert_eq!(clipped_white(&Mat3::IDENTITY), [1.0; 3]); + // a model that can't be inverted, or only just, keeps the camera's neutral + assert_eq!(clipped_white(&Mat3([[1.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, 1.0]])), [1.0; 3]); + assert_eq!(clipped_white(&Mat3([[1.0, 1.0, 0.0], [1.0, 1.0 + 1e-8, 0.0], [0.0, 0.0, 1.0]])), [1.0; 3]); + // so does one whose white would overflow the rebuilt highlights (2.5e38 in red here) + assert_eq!(clipped_white(&Mat3([[3e-39, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]])), [1.0; 3]); + // or move the camera's neutral implausibly far + assert_eq!(clipped_white(&Mat3([[0.1, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]])), [1.0; 3]); + // what it answers renders neutral through the f32 matrix, and fully clipped highlights stay finite + for model in [look, Mat3([[1.9, -0.7, -0.2], [-0.3, 1.4, -0.1], [0.0, -0.4, 1.4]])] { + let w = clipped_white(&model); + let m = model.to_f32(); + let out: [f32; 3] = std::array::from_fn(|i| m[i][0] * w[0] + m[i][1] * w[1] + m[i][2] * w[2]); + assert!(out.iter().all(|v| (v - out[1]).abs() <= 1e-3 * out[1]), "{w:?} renders {out:?}"); + let colour = CameraColourModel { wb: [2.5, 1.0, 1.6], matrix: model, gain: 1.0, tables: None, hue_sat: None }; + let mut img = Rgb32f::filled(4, 4, [1.0; 3]); + colour.rebuild_highlights(&mut img); + colour.apply(&mut img); + assert!(img.data.iter().flatten().all(|v| v.is_finite()), "{:?}", img.get(0, 0)); + } + } + + /// Issue #523: fully clipped highlights are neutral in the white-balanced camera RGB, but a colour model that + /// doesn't map that neutral to neutral tinted them. The tone map's roll-off to white hid the tint at Neutral; + /// Highlights -100, or the issue's edit (Exposure -0.88, Highlights -75…), showed it. + #[test] + fn clipped_highlights_stay_neutral_when_darkened() { + use lightcraft_pipeline::{RenderRequest, render}; + let wb = [2.4f32, 1.0, 1.6]; + let colour = CameraColourModel { wb, matrix: yellowing_look(), gain: 1.0, tables: None, hue_sat: None }; + // a grey wall (camera neutral) with a fully clipped light and a light clipped in green only + let mut img = Rgb32f::from_fn(96, 64, |x, y| { + let d = |cx: f32| ((x as f32 - cx).powi(2) + (y as f32 - 32.0).powi(2)).sqrt(); + if d(28.0) < 14.0 { + [1.0; 3] + } else if d(70.0) < 12.0 { + [0.6 / wb[0], 1.0, 0.7 / wb[2]] + } else { + [0.3 / wb[0], 0.3, 0.3 / wb[2]] + } + }); + colour.rebuild_highlights(&mut img); + colour.apply(&mut img); + let info = SourceInfo { raw: true, ..Default::default() }; + let mut issue = lightcraft_develop::DevelopSettings::default(); + let l = &mut issue.light; + (l.exposure, l.contrast, l.highlights, l.shadows, l.whites, l.blacks) = (-0.88, 5.0, -75.0, 52.0, 6.0, -37.0); + let mut h100 = lightcraft_develop::DevelopSettings::default(); + h100.light.highlights = -100.0; + for (name, s) in [("Highlights -100", h100), ("#523's edit", issue)] { + let out = render(&img, &info, &s, &RenderRequest::fit(96, 64)).image; + let p = out.get(28, 32); + assert!(p[1] < 250, "{name}: the clipped light is darkened ({p:?})"); + let spread = p[0].max(p[1]).max(p[2]) - p[0].min(p[1]).min(p[2]); + assert!(spread <= 1, "{name}: clipped light {p:?} is not neutral"); + } + } + + /// A smooth pink gradient (`[0.9k, k, 0.8k]` in the sensor, k from 0.8 to 1.35, white balance `[2.5, 1, 1.6]`) + /// through the point where green clips, 1600 px long, in `rows` of a 1600 × `height` picture that is otherwise + /// near black, rebuilt and rendered at Exposure -2. Returns the rendered colour of the middle row's last + /// unclipped pixel and the next two, and the largest step between neighbouring pixels along that row. + fn darkened_ramp(rows: std::ops::Range, height: usize) -> ([[u8; 4]; 3], u8) { + use lightcraft_pipeline::{RenderRequest, render}; + let wb = [2.5f32, 1.0, 1.6]; + let colour = CameraColourModel { wb, matrix: lightcraft_color::Mat3::IDENTITY, gain: 1.0, tables: None, hue_sat: None }; + let w = 1600usize; + let k = |x: usize| 0.8 + 0.55 * x as f32 / (w - 1) as f32; + let sensor = + Rgb32f::from_fn(w, height, |x, y| if rows.contains(&y) { [0.9 * k(x), k(x), 0.8 * k(x)].map(|v| v.min(1.0)) } else { [0.01; 3] }); + let y = (rows.start + rows.end) / 2; + let first = (0..w).find(|&x| sensor.get(x, y)[1] >= HIGHLIGHT_CLIP).unwrap(); + let mut img = sensor.clone(); + colour.rebuild_highlights(&mut img); + colour.apply(&mut img); + let mut s = lightcraft_develop::DevelopSettings::default(); + s.light.exposure = -2.0; + let out = render(&img, &SourceInfo { raw: true, ..Default::default() }, &s, &RenderRequest::fit(w, height)).image; + let step = |x: usize| { + let (a, b) = (out.get(x - 1, y), out.get(x, y)); + (0..3).map(|c| a[c].abs_diff(b[c])).max().unwrap_or(0) + }; + ([out.get(first - 1, y), out.get(first, y), out.get(first + 1, y)], (1..w).map(step).max().unwrap_or(0)) + } + + /// The darkened pink gradient carries on across the point where green clips: before the fix the fade towards + /// the clipped neutral, driven by how close the other channels were to their clip, made a step of 37 levels + /// there. + #[test] + fn rendered_colour_is_continuous_where_the_first_channel_clips() { + let ([a, b, c], worst) = darkened_ramp(0..16, 16); + let step = |p: [u8; 4], q: [u8; 4]| (0..3).map(|i| p[i].abs_diff(q[i])).max().unwrap_or(0); + assert!(step(a, b) <= 1 && step(b, c) <= 1, "at the first clip: {a:?} {b:?} {c:?}"); + assert!(worst <= 3, "largest step between neighbouring pixels: {worst} levels"); + } + + /// The same gradient as a thin bright line on near black: every pixel of it is next to a strong edge, so there + /// is no reliable colour around it. It still carries on across the first clip, however thin the line and however + /// much dark background there is (a 4-row line on 32 rows used to jump from [225, 126, 153] to [198, 177, 163]; + /// on 512 rows, from [225, 126, 153] to [218, 145, 151]). + #[test] + fn rendered_thin_line_is_continuous_where_the_first_channel_clips() { + let step = |p: [u8; 4], q: [u8; 4]| (0..3).map(|i| p[i].abs_diff(q[i])).max().unwrap_or(0); + for height in [32usize, 128, 512, 1024] { + for width in [1usize, 2, 4, 16] { + let top = (height - width) / 2; + let ([a, b, c], worst) = darkened_ramp(top..top + width, height); + let case = format!("{width} rows of {height}"); + assert!(step(a, b) <= 1 && step(b, c) <= 1, "{case}, at the first clip: {a:?} {b:?} {c:?}"); + assert!(worst <= 3, "{case}: largest step between neighbouring pixels {worst} levels"); + } + } + } + + /// Issue #523 on a public ILCE-7RM4 raw (raw.pixls.us, CC0; skipped without the corpus): a dusk sky clipped in + /// green behind a poplar. Its rebuilt colour used to come from the tree's foliage and lens fringes, a green halo + /// on one side of the tree and a magenta one on the other (Highlights -100 or the issue's edit show them). + #[test] + fn corpus_clipped_sky_beside_a_tree_has_no_coloured_halo() { + let dir = std::env::var_os("LIGHTCRAFT_CORPUS") + .map(std::path::PathBuf::from) + .unwrap_or_else(|| std::path::Path::new(env!("CARGO_MANIFEST_DIR")).join("../../corpus")) + .join("raw"); + let Ok(bytes) = std::fs::read(dir.join("arw-sony-a7rm4-14bit-compressed.arw")) else { + eprintln!("skip: {} absent", dir.display()); + return; + }; + let (img, _) = load_bytes(&bytes, 1600).unwrap(); + let (w, h) = (img.width, img.height); + let lum = |p: [f32; 3]| lightcraft_color::luminance_2020(p); + let chroma = |p: [f32; 3]| { + let s = (p[0] + p[1] + p[2]).max(1e-6); + [p[0] / s, p[2] / s] + }; + // the poplar stands at about 30 % of the width, its top 25 to 34 % down the frame + let px = |f: f32, n: usize| (f * n as f32) as usize; + let (mut near, mut far) = ([[0f64; 2]; 2], [[0f64; 2]; 2]); + let mut rows = 0; + for y in px(0.24, h)..px(0.34, h) { + let row = &img.data[y * w..(y + 1) * w]; + let sky = lum(row[px(0.2, w)]); + let dark: Vec = (px(0.25, w)..px(0.35, w)).filter(|&x| lum(row[x]) < 0.5 * sky).collect(); + let (Some(&left), Some(&right)) = (dark.first(), dark.last()) else { continue }; + let k = w as f32 / 6000.0; + let mean = |a: usize, b: usize| { + let n = (b - a) as f64; + (a..b).map(|x| chroma(row[x])).fold([0f64; 2], |s, c| [s[0] + c[0] as f64 / n, s[1] + c[1] as f64 / n]) + }; + let (n0, n1, f0, f1) = (px(8.0 / 6000.0, w).max(2), px(40.0 / 6000.0, w), (200.0 * k) as usize, (320.0 * k) as usize); + for (side, (n, f)) in + [(mean(left - n1, left - n0), mean(left - f1, left - f0)), (mean(right + n0, right + n1), mean(right + f0, right + f1))] + .into_iter() + .enumerate() + { + for c in 0..2 { + near[side][c] += n[c]; + far[side][c] += f[c]; + } + } + rows += 1; + } + assert!(rows > 50, "the tree was not found ({rows} rows)"); + for side in 0..2 { + let d = (0..2).map(|c| ((near[side][c] - far[side][c]) / rows as f64).abs()).fold(0.0, f64::max); + assert!(d < 0.004, "{} of the tree the clipped sky's chromaticity is off by {d:.4}", ["left", "right"][side]); + } + } + /// iPhone ProRAW-style semantic masks drive the Sky mask: Apple's sky matte (here the right /// half of the sensor, i.e. the bottom of the photo once turned 90° clockwise — where the sky /// heuristic would never look) is used; a matte with an unknown name is ignored. diff --git a/crates/engine/src/media.rs b/crates/engine/src/media.rs index 1231221b..b73d3f1d 100644 --- a/crates/engine/src/media.rs +++ b/crates/engine/src/media.rs @@ -62,7 +62,7 @@ impl SettingsHashes { } /// Bump when the pipeline's output changes, to invalidate cached thumbnails. -pub const RENDER_CACHE_VERSION: u64 = 22; +pub const RENDER_CACHE_VERSION: u64 = 27; /// Thumbnails render at one of these long edges (so window/cell size changes reuse the cache). pub const THUMB_SIZES: [usize; 4] = [128, 256, 384, 512]; diff --git a/crates/gpu/Cargo.toml b/crates/gpu/Cargo.toml index b19acf54..dc348ee8 100644 --- a/crates/gpu/Cargo.toml +++ b/crates/gpu/Cargo.toml @@ -30,3 +30,4 @@ denoise = ["dep:lightcraft-denoise", "lightcraft-denoise/runtime"] [dev-dependencies] lightcraft-scenes = { workspace = true } +lightcraft-raw = { workspace = true } diff --git a/crates/gpu/tests/equivalence.rs b/crates/gpu/tests/equivalence.rs index cb6e7c02..b10836c3 100644 --- a/crates/gpu/tests/equivalence.rs +++ b/crates/gpu/tests/equivalence.rs @@ -406,6 +406,59 @@ fn gpu_matches_cpu() { eprintln!("worst: mean {:.4} LSB, max {} LSB", worst.0, worst.1); } +/// Issue #523: a raw's clipped highlights rebuilt as the raw decode does it (fully clipped ones to the neutral of +/// the colour model applied after them, here a look whose rows don't sum to one) stay neutral when darkened, on +/// the GPU as on the CPU, and the two renders agree. +#[test] +fn clipped_highlights_stay_neutral_on_both() { + if !gpu() { + return; + } + let wb = [2.4f32, 1.0, 1.6]; + let look = lightcraft_color::Mat3([[0.42, 0.548, 0.04], [0.246, 0.728, 0.033], [-0.247, 0.679, 0.447]]); + let luma = |c: [f64; 3]| (0..3).map(|i| c[i] * f64::from(lightcraft_color::LUMA_2020[i])).sum::(); + let n = look.inverse().unwrap().apply([1.0; 3]); + let white = n.map(|v| (v * luma(look.apply([1.0; 3])) / luma(look.apply(n))) as f32); + // a grey wall with a fully clipped light, a light clipped in green only, and a dark bar with a purple fringe + // in front of the second + let mut cam = Rgb32f::from_fn(320, 240, |x, y| { + let d = |cx: f32| ((x as f32 - cx).powi(2) + (y as f32 - 120.0).powi(2)).sqrt(); + if d(90.0) < 50.0 { + [1.0; 3] + } else if (226..232).contains(&x) && (70..170).contains(&y) { + [0.02, 0.03, 0.02] + } else if (x == 225 || x == 232) && (70..170).contains(&y) { + [0.3, 0.25, 0.6] + } else if d(228.0) < 50.0 { + [0.6 / wb[0], 1.0, 0.7 / wb[2]] + } else { + [0.3 / wb[0], 0.3, 0.3 / wb[2]] + } + }); + lightcraft_raw::highlight::reconstruct_with(&mut cam, wb, 0.99, white); + let m = look.to_f32(); + let src = Arc::new(cam.map(|p| { + let c = [p[0] * wb[0], p[1] * wb[1], p[2] * wb[2]]; + std::array::from_fn(|i| (m[i][0] * c[0] + m[i][1] * c[1] + m[i][2] * c[2]).max(0.0)) + })); + let info = SourceInfo { raw: true, ..Default::default() }; + let req = RenderRequest::fit(320, 240); + let mut h100 = DevelopSettings::default(); + h100.light.highlights = -100.0; + let mut issue = DevelopSettings::default(); + let l = &mut issue.light; + (l.exposure, l.contrast, l.highlights, l.shadows, l.whites, l.blacks) = (-0.88, 5.0, -75.0, 52.0, 6.0, -37.0); + for (name, s) in [("clipped highlights -100", h100), ("clipped highlights #523", issue)] { + check(name, &src, &info, &s, &req); + let cpu = render(&src, &info, &s, &req).image; + let gpu = lightcraft_gpu::render(&src, &info, &s, &req, None).expect("gpu render").image; + for (who, img) in [("CPU", &cpu), ("GPU", &gpu)] { + let p = img.get(90, 120); + assert!(p[1] < 250 && p[0].max(p[1]).max(p[2]) - p[0].min(p[1]).min(p[2]) <= 1, "{name} {who}: clipped light {p:?}"); + } + } +} + #[test] fn rendered_sources_and_other_scenes() { if !gpu() { diff --git a/crates/raw/src/highlight.rs b/crates/raw/src/highlight.rs index 1731cb25..950be3af 100644 --- a/crates/raw/src/highlight.rs +++ b/crates/raw/src/highlight.rs @@ -2,9 +2,18 @@ //! //! - [`clip_neutral`]: white-balance-aware clipping so that sensor-clipped areas render neutral instead of //! magenta/cyan (each channel limited to the smallest white-balanced clip level). -//! - [`reconstruct`]: where only some channels are clipped, rebuild them from the unclipped channels using the -//! local chromaticity of nearby unclipped pixels (diffused into the clipped region with a coarse-to-fine -//! normalised-convolution fill). Fully clipped pixels become neutral at the brightest plausible level. +//! - [`reconstruct`] / [`reconstruct_with`]: where only some channels are clipped, rebuild them from the +//! unclipped channels using the chromaticity of reliable unclipped content nearby (bright, away from edges and +//! from the clipped pixels themselves; diffused into the clipped region coarse to fine), and neutral where there +//! is none within reach. Fully clipped pixels become neutral at the brightest plausible level; partly clipped +//! ones fade into it on their way from the first channel clipping to the last. +//! +//! Why the care: the white-balance gains push clipped channels apart (green clips first, red and blue keep +//! rising), so whatever colour the rebuilt pixels carry is hidden at Neutral only by the tone map's roll-off to +//! white. Any darkening (Highlights, negative Exposure or Whites) brings it out. Two things used to leak colour +//! in: the chromaticity of whatever unclipped pixel happened to be nearest (dark foliage, purple and green +//! fringes along branches, pixels demosaiced from clipped neighbours), and a colour model that does not render +//! the camera's white-balanced neutral as neutral (a look fitted to the camera's JPEG, a camera profile). use lightcraft_raster::Rgb32f; use rayon::prelude::*; @@ -20,119 +29,335 @@ pub fn clip_neutral(img: &mut Rgb32f, wb: [f32; 3], clip: f32) { } } -/// Fill `values` (per-pixel vectors) where `valid` is false from valid neighbours, coarse to fine. -fn fill_invalid(w: usize, h: usize, values: &mut [[f32; 3]], valid: &mut [bool]) { - if w == 0 || h == 0 || valid.iter().all(|&v| v) || !valid.iter().any(|&v| v) { - return; - } - // Build a pyramid of weighted means, fill each level from the next coarser one. - 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]; - { - let (values, valid) = (&*values, &*valid); - downsample(w, h, &mut cv, &mut cval, |i| valid[i].then(|| values[i])); - } - if cw * ch < w * h { - fill_invalid(cw, ch, &mut cv, &mut cval); - } else { - // cannot shrink further (1×1): nothing valid anywhere handled above - return; - } - values.par_chunks_mut(w).zip(valid.par_chunks_mut(w)).enumerate().for_each(|(y, (vrow, okrow))| { - for x in 0..w { - if !okrow[x] { - vrow[x] = upsample(&cv, cw, ch, x, y); - okrow[x] = true; - } - } - }); +/// Unclipped pixels lend their chromaticity in proportion to how close they come to clipping (their largest +/// channel, relative to the clip level): from nothing below `SOURCE_LO` to fully above `SOURCE_HI`. Dark content +/// (foliage, branches, buildings against a clipped sky) says nothing about the sky's colour. +const SOURCE_LO: f32 = 0.25; +const SOURCE_HI: f32 = 0.6; +/// Sources on a strong edge (a 3 × 3 half-resolution neighbour brighter or darker by this many stops) are +/// left out, fading from `EDGE_LO` to `EDGE_HI`: edge pixels mix two surfaces and carry lens fringes. +const EDGE_LO: f32 = 0.15; +const EDGE_HI: f32 = 0.4; +/// Source cells within this many half-resolution cells of a clipped pixel are left out: their channels were +/// demosaiced partly from clipped samples. (At least 1: the same neighbourhood tells which cells clipped +/// pixels sample.) +const MARGIN: usize = 1; +/// Weight of the coarser estimate against a cell's own sources when the pyramid is pushed back down (in units +/// of source density: a cell half covered by full-weight sources counts 0.5). +const PRIOR: f32 = 0.01; +/// The rebuilt colour fades to neutral where reliable content covers less than `SUPPORT` of the surroundings at +/// every scale up to `REACH` of the long edge (so previews and exports agree), gradually, as the coverage of the +/// widest one falls: a clipped sky between bare branches is better white than the colour of a river far away. +const REACH: f32 = 1.0 / 4.0; +const SUPPORT: f32 = 0.05; +/// An unclipped channel above `FLOOR_LO` of its clip level counts more and more as the clip floor it is about to +/// become (see [`reconstruct_with`]). +const FLOOR_LO: f32 = 0.75; +/// Without reliable colour around, a clipped pixel is rebuilt as neutral only as far as it is brighter (or +/// otherwise unlike) the unclipped pixels nearest its clip: the "rim", unclipped pixels above `RIM_LO` of the +/// clip level, weighted towards the clip. Up to `NEAR_RIM` stops past the rim (the rim is an average, a little +/// below the clip) the clipped channels stay at their clip, and over the next `DEPTH` stops they go over to the +/// neutral guess, so a surface runs into clipping without a step whatever its colour; well past it, or where the +/// rim is another surface (its colour differs), the neutral guess applies. +const RIM_LO: f32 = 0.85; +const NEAR_RIM: f32 = 0.03; +const DEPTH: f32 = 0.2; + +#[inline] +fn smoothstep(e0: f32, e1: f32, x: f32) -> f32 { + let t = ((x - e0) / (e1 - e0)).clamp(0.0, 1.0); + t * t * (3.0 - 2.0 * t) } -/// One pyramid step: each `2 × 2` block of a `w × h` level becomes the mean of its valid -/// entries (`at(i)`: the value at index `i` when valid) in `cv`, `cval` (`⌈w/2⌉ × ⌈h/2⌉`). -fn downsample(w: usize, h: usize, cv: &mut [[f32; 3]], cval: &mut [bool], at: impl Fn(usize) -> Option<[f32; 3]> + Sync) { - let cw = w.div_ceil(2); - cv.par_chunks_mut(cw).zip(cval.par_chunks_mut(cw)).enumerate().for_each(|(y, (vrow, okrow))| { +/// One level of a pyramid, `w × h` cells of `N - 1` values and a weight. Pulled: per cell the weighted sums of +/// the values and the weight, as densities (per pixel of the cell). Pushed: the estimated values, and still the +/// weight. The chromaticity pyramid holds red and blue chromaticity, the rim pyramid white-balanced RGB. +struct Level { + w: usize, + h: usize, + cells: Vec<[f32; N]>, +} + +/// Bilinear sample of a level at the centre of pixel `(x, y)` of a grid `scale` times finer. +#[inline] +fn sample(l: &Level, scale: f32, x: usize, y: usize) -> [f32; N] { + let (cv, cw, ch) = (&l.cells, l.w, l.h); + let fx = ((x as f32 + 0.5) / scale - 0.5).clamp(0.0, cw.saturating_sub(1) as f32); + let fy = ((y as f32 + 0.5) / scale - 0.5).clamp(0.0, ch.saturating_sub(1) as f32); + let (x0, y0) = (fx.floor() as usize, fy.floor() as usize); + let (x1, y1) = ((x0 + 1).min(cw.saturating_sub(1)), (y0 + 1).min(ch.saturating_sub(1))); + let (tx, ty) = (fx - x0 as f32, fy - y0 as f32); + let g = |xx: usize, yy: usize| cv.get(yy * cw + xx).copied().unwrap_or([0.0; N]); + let (a, b, c, d) = (g(x0, y0), g(x1, y0), g(x0, y1), g(x1, y1)); + std::array::from_fn(|k| { + let top = a[k] + (b[k] - a[k]) * tx; + let bot = c[k] + (d[k] - c[k]) * tx; + top + (bot - top) * ty + }) +} + +/// The half-resolution source level of `img`: every `2 × 2` block's source-weighted chromaticity sums, with +/// blocks near clipped pixels or on strong edges left out. Also returns, per cell, whether a clipped pixel lies +/// in its row within [`MARGIN`] cells (the vertical half of that test is left to the readers). +fn sources(img: &Rgb32f, wb: [f32; 3], clip: f32) -> (Level<3>, Vec) { + let (w, h) = (img.width, img.height); + let (cw, ch) = (w.div_ceil(2), h.div_ceil(2)); + // per cell: [red sum, blue sum, weight] densities, luminance, whether it holds a clipped pixel + let mut cells = vec![[0f32; 3]; cw * ch]; + let mut lum = vec![0f32; cw * ch]; + let mut near = vec![false; cw * ch]; + cells.par_chunks_mut(cw).zip(lum.par_chunks_mut(cw)).zip(near.par_chunks_mut(cw)).enumerate().for_each(|(y, ((crow, lrow), nrow))| { + let mut touched = vec![false; cw]; + let row = |r: usize| img.data.get(r * w..(r + 1) * w).unwrap_or(&[]); + let (r0, r1) = (row(2 * y), row(2 * y + 1)); + let lo = SOURCE_LO * clip; 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 - && let Some(v) = at(sy * w + sx) - { - for c in 0..3 { - s[c] += v[c]; - } - n += 1; - } + let (mut a, mut sum, mut n, mut hit) = ([0f32; 3], 0f32, 0usize, false); + let span = 2 * x..(2 * x + 2).min(w); + for p in r0.get(span.clone()).unwrap_or(&[]).iter().chain(r1.get(span).unwrap_or(&[])) { + let q = [p[0] * wb[0], p[1] * wb[1], p[2] * wb[2]]; + let s = q[0].max(0.0) + q[1].max(0.0) + q[2].max(0.0); + sum += s; + n += 1; + let top = p[0].max(p[1]).max(p[2]); + if top >= clip { + hit = true; + continue; + } + if top > lo && s > 1e-4 { + let b = smoothstep(SOURCE_LO, SOURCE_HI, top / clip) / s; + a[0] += b * q[0].max(0.0); + a[1] += b * q[2].max(0.0); + a[2] += b * s; } } - if n > 0 { - vrow[x] = s.map(|v| v / n as f32); - okrow[x] = true; - } + let k = 1.0 / n.max(1) as f32; + crow[x] = a.map(|v| v * k); + lrow[x] = (sum * k).max(1e-6); + touched[x] = hit; + } + for (x, nr) in nrow.iter_mut().enumerate() { + *nr = touched.get(x.saturating_sub(MARGIN)..=(x + MARGIN).min(cw - 1)).is_some_and(|t| t.iter().any(|&t| t)); } }); + // leave out cells near clipped pixels and on strong edges (luminance ratios to the neighbours) + let (lum, near_ref) = (&lum, &near); + let (edge_lo, edge_hi) = (EDGE_LO.exp2(), EDGE_HI.exp2()); + cells.par_chunks_mut(cw).enumerate().for_each(|(y, row)| { + let (ya, yb) = (y.saturating_sub(MARGIN), (y + MARGIN).min(ch - 1)); + let (ea, eb) = (y.saturating_sub(1), (y + 1).min(ch - 1)); + for (x, cell) in row.iter_mut().enumerate() { + if cell[2] <= 0.0 { + continue; + } + if (ya..=yb).any(|yy| near_ref.get(yy * cw + x).copied().unwrap_or(false)) { + *cell = [0.0; 3]; + continue; + } + let l = lum.get(y * cw + x).copied().unwrap_or(1.0); + let (xa, xb) = (x.saturating_sub(1), (x + 1).min(cw - 1)); + let (mut lo, mut hi) = (l, l); + for yy in ea..=eb { + for &v in lum.get(yy * cw + xa..=yy * cw + xb).unwrap_or(&[]) { + lo = lo.min(v); + hi = hi.max(v); + } + } + let keep = 1.0 - smoothstep(edge_lo, edge_hi, (hi / l).max(l / lo)); + *cell = cell.map(|v| v * keep); + } + }); + (Level { w: cw, h: ch, cells }, near) } -/// Bilinear sample of the `cw × ch` coarse level at fine pixel `(x, y)` (twice the resolution). -fn upsample(cv: &[[f32; 3]], cw: usize, ch: usize, x: usize, y: usize) -> [f32; 3] { - 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); - let (x0, y0) = (fx.floor() as usize, fy.floor() as usize); - let (x1, y1) = ((x0 + 1).min(cw - 1), (y0 + 1).min(ch - 1)); - let (tx, ty) = (fx - x0 as f32, fy - y0 as f32); - let g = |xx: usize, yy: usize| cv[yy * cw + xx]; - let mut v = [0f32; 3]; - for c in 0..3 { - let top = g(x0, y0)[c] + (g(x1, y0)[c] - g(x0, y0)[c]) * tx; - let bot = g(x0, y1)[c] + (g(x1, y1)[c] - g(x0, y1)[c]) * tx; - v[c] = top + (bot - top) * ty; +/// The quarter-resolution rim level of `img`: every `4 × 4` block's white-balanced RGB of its unclipped pixels +/// above [`RIM_LO`] of the clip level, weighted (steeply) towards the clip, as densities. Pixels next to clipped +/// ones (`margin`, by half-resolution cell) are left out, like the colour sources: they were demosaiced partly from +/// clipped samples, and lens fringes along clipped content would pass for its colour. +fn rim(img: &Rgb32f, wb: [f32; 3], clip: f32, margin: &(dyn Fn(usize, usize) -> bool + Sync)) -> Level<4> { + let (w, h) = (img.width, img.height); + let (cw, ch) = (w.div_ceil(4), h.div_ceil(4)); + let mut cells = vec![[0f32; 4]; cw * ch]; + cells.par_chunks_mut(cw).enumerate().for_each(|(y, row)| { + for (x, cell) in row.iter_mut().enumerate() { + let (mut a, mut n) = ([0f32; 4], 0usize); + for sy in 4 * y..(4 * y + 4).min(h) { + for (i, p) in img.data.get(sy * w + 4 * x..sy * w + (4 * x + 4).min(w)).unwrap_or(&[]).iter().enumerate() { + n += 1; + let top = p[0].max(p[1]).max(p[2]) / clip; + if !(RIM_LO..1.0).contains(&top) || margin((4 * x + i) / 2, sy / 2) { + continue; + } + let k = ((top - RIM_LO) / (1.0 - RIM_LO)).powi(4); + for c in 0..3 { + a[c] += k * p[c] * wb[c]; + } + a[3] += k; + } + } + *cell = a.map(|v| v / n.max(1) as f32); + } + }); + Level { w: cw, h: ch, cells } +} + +/// The next coarser level: each `2 × 2` block's mean densities. +fn pull(l: &Level) -> Level { + let (w, h) = (l.w.div_ceil(2), l.h.div_ceil(2)); + let mut cells = vec![[0f32; N]; w * h]; + cells.par_chunks_mut(w).enumerate().for_each(|(y, row)| { + for (x, cell) in row.iter_mut().enumerate() { + let (mut s, mut n) = ([0f32; N], 0usize); + for (sx, sy) in [(2 * x, 2 * y), (2 * x + 1, 2 * y), (2 * x, 2 * y + 1), (2 * x + 1, 2 * y + 1)] { + if let Some(c) = (sx < l.w && sy < l.h).then(|| l.cells.get(sy * l.w + sx)).flatten() { + for k in 0..N { + s[k] += c[k]; + } + n += 1; + } + } + *cell = s.map(|v| v / n.max(1) as f32); + } + }); + Level { w, h, cells } +} + +/// The estimated values of the cells of `levels[0]` that `needed` marks (all of them when `None`), and of every +/// cell of the coarser levels: each level's own sources blended with the next coarser estimate. The top level is +/// its own sources' weighted mean (normalised: however sparse they are in a large empty picture, the empty part +/// never pulls the estimate anywhere), and `empty` (at zero weight) when it has none at all. +fn push(levels: &mut [Level], empty: [f32; N], needed: Option<&(dyn Fn(usize, usize) -> bool + Sync)>) { + let estimate = |cell: [f32; N], prior: [f32; N]| -> [f32; N] { + let weight = cell[N - 1]; + let d = weight + PRIOR; + std::array::from_fn(|k| if k + 1 == N { weight } else { (cell[k] + PRIOR * prior[k]) / d }) + }; + let Some(last) = levels.last_mut() else { return }; + last.cells.iter_mut().for_each(|c| { + let weight = c[N - 1]; + *c = if weight > 0.0 { std::array::from_fn(|k| if k + 1 == N { weight } else { c[k] / weight }) } else { empty }; + }); + for i in (0..levels.len().saturating_sub(1)).rev() { + let (fine, coarse) = levels.split_at_mut(i + 1); + let (Some(fine), Some(coarse)) = (fine.last_mut(), coarse.first()) else { continue }; + let w = fine.w; + let needed = needed.filter(|_| i == 0); + fine.cells.par_chunks_mut(w).enumerate().for_each(|(y, row)| { + for (x, cell) in row.iter_mut().enumerate() { + if needed.is_some_and(|f| !f(x, y)) { + continue; + } + *cell = estimate(*cell, sample(coarse, 2.0, x, y)); + } + }); } - v } -/// White-balanced chromaticity of an unclipped pixel (`None`: clipped or too dark to tell). -#[inline] -fn chroma_of(p: &[f32; 3], wb: [f32; 3], clip: f32) -> Option<[f32; 3]> { - if p[0] >= clip || p[1] >= clip || p[2] >= clip { +/// The white-balanced colour of a partly clipped pixel (sensor values `p`, white-balanced `q`, clipped channels +/// `cl`) whose colour is estimated at chromaticity `r`: clipped channels rebuilt, unclipped ones kept. `None` +/// when no unclipped channel can tell the level. +fn rebuild(p: [f32; 3], q: [f32; 3], cl: [bool; 3], r: [f32; 3], clip: f32, max_level: f32) -> Option<[f32; 3]> { + // The white-balanced channel sum the unclipped channels imply: their mean, except that a channel about to clip + // counts more and more as the floor it is about to become (continuous with the neighbours where it has). + let (mut sum, mut near_clip, mut k) = (0f32, 0f32, 0usize); + for c in 0..3 { + if !cl[c] && r[c] > 1e-3 { + let s = q[c] / r[c]; + sum += s; + near_clip = near_clip.max(s * smoothstep(FLOOR_LO, 1.0, p[c] / clip)); + k += 1; + } + } + if k == 0 { return None; } - let q = [p[0] * wb[0], p[1] * wb[1], p[2] * wb[2]]; - let s = q[0] + q[1] + q[2]; - (s > 1e-4).then(|| q.map(|v| v.max(0.0) / s)) + // Two clipped channels must not be left at their clip levels: those differ by the white-balance gains (green + // clips first), and the pair would show as magenta or red. They keep the estimate's ratio instead, at the level + // that lifts each to its clip (no channel above the clipped neutral). + let floor = (0..3).filter(|&c| cl[c]).map(|c| q[c] / r[c].max(1e-3)).fold(0.0f32, f32::max); + let top = r[0].max(r[1]).max(r[2]).max(1e-3); + let sum = (sum / k as f32).max(near_clip).max(floor.min(max_level / top)); + Some(std::array::from_fn(|c| if cl[c] { q[c].max(sum * r[c]) } else { q[c] })) +} + +/// Reconstruct partially clipped channels, with a neutral that is neutral in camera RGB. See +/// [`reconstruct_with`]. +pub fn reconstruct(img: &mut Rgb32f, wb: [f32; 3], clip: f32) -> usize { + reconstruct_with(img, wb, clip, [1.0; 3]) } /// Reconstruct partially clipped channels. `clip` is the sensor clip level in normalised units (use slightly -/// below 1.0, e.g. 0.99), `wb` the white-balance multipliers that will be applied afterwards. -/// Returns the number of pixels that had at least one clipped channel. -/// -/// The chromaticity field is diffused coarse to fine from the unclipped pixels. A clipped pixel is -/// never valid at full resolution, so its chromaticity is always the bilinear sample of the -/// half-resolution level: we start the pyramid there, straight from the image, and never build -/// the full-resolution field (same result, a quarter of the memory and far less work). -pub fn reconstruct(img: &mut Rgb32f, wb: [f32; 3], clip: f32) -> usize { +/// below 1.0, e.g. 0.99), `wb` the white-balance multipliers that will be applied afterwards, `white` the +/// white-balanced camera RGB that the colour model applied afterwards renders neutral (`[1, 1, 1]` for a model +/// that maps the white-balanced neutral to neutral), scaled to the brightness that `[1, 1, 1]` renders at. +/// Fully clipped pixels become `white` times the brightest plausible level; partly clipped ones take the +/// chromaticity of reliable unclipped content nearby, and that of `white` where there is none. +/// Returns the number of pixels that had at least one clipped channel; the others are left exactly as they were. +pub fn reconstruct_with(img: &mut Rgb32f, wb: [f32; 3], clip: f32, white: [f32; 3]) -> usize { let (w, h) = (img.width, img.height); 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 { + if count == 0 || w == 0 || h == 0 { return 0; } - // chromaticity (white-balanced ratios to the channel sum) of unclipped pixels, at half resolution + // a white far from the camera's neutral (or not a number) is a broken colour model: keep the camera's + let plausible = white.iter().all(|v| v.is_finite() && (1.0 / 16.0..=16.0).contains(v)); + let white = if plausible { white } else { [1.0; 3] }; + let total = white[0] + white[1] + white[2]; + let neutral = [white[0] / total, white[1] / total, white[2] / total]; + // A clipped pixel is never a source, so its chromaticity is always the bilinear sample of the + // half-resolution estimate: the pyramid starts there, straight from the image. + let (first, near) = sources(img, wb, clip); + let mut levels = vec![first]; + while let Some(l) = levels.last().filter(|l| l.w > 1 || l.h > 1) { + let next = pull(l); + levels.push(next); + } + // a clipped pixel samples the cells within one of its own, which are within one of a clipped cell let (cw, ch) = (w.div_ceil(2), h.div_ceil(2)); - let mut cv = vec![[1f32 / 3.0; 3]; cw * ch]; - let mut cval = vec![false; cw * ch]; - { - let data = &img.data; - downsample(w, h, &mut cv, &mut cval, |i| chroma_of(&data[i], wb, clip)); + let needed = |x: usize, y: usize| (y.saturating_sub(1)..=(y + 1).min(ch - 1)).any(|yy| near.get(yy * cw + x).copied().unwrap_or(false)); + push(&mut levels, [neutral[0], neutral[2], 0.0], Some(&needed)); + // How much reliable content lies around each needed cell: its density at any scale up to [`REACH`] (level + // i's cells span 2^i cells of level 0). Near reliable content some fine level is dense with it; far from it + // only the widest one may still hold a little, and the trust fades with that. Kept in the third entry of the + // level 0 estimates, which clipped pixels sample with their colour. + let reach = REACH * w.max(h) as f32 / 2.0; + let widest = (0..levels.len()).find(|&i| (1u64 << i.min(62)) as f32 >= reach).unwrap_or(levels.len().saturating_sub(1)); + if let Some((first, coarser)) = levels.split_first_mut() { + let support: Vec<(&Level<3>, f32)> = coarser.iter().take(widest).enumerate().map(|(i, l)| (l, (1u64 << (i + 1).min(62)) as f32)).collect(); + let fw = first.w; + first.cells.par_chunks_mut(fw).enumerate().for_each(|(y, row)| { + for (x, cell) in row.iter_mut().enumerate() { + if !needed(x, y) { + continue; + } + let mut t = 0f32; + for &(l, scale) in &support { + t = t.max(smoothstep(0.0, SUPPORT, sample(l, scale, x, y)[2])); + if t >= 1.0 { + break; + } + } + cell[2] = t; + } + }); } - if cw * ch < w * h { - fill_invalid(cw, ch, &mut cv, &mut cval); - } - let any_valid = cw * ch < w * h && cval.iter().any(|&v| v); + let Some(field) = levels.first() else { return count }; + // the rim, only where some clipped pixel lacks reliable colour (level 0 of its pyramid holds the estimate) + let untrusted = field.cells.par_iter().enumerate().any(|(i, c)| c[2] < 1.0 && needed(i % field.w, i / field.w)); + let rims = untrusted + .then(|| { + let mut levels = vec![rim(img, wb, clip, &needed)]; + while let Some(l) = levels.last().filter(|l| l.w > 1 || l.h > 1) { + let next = pull(l); + levels.push(next); + } + push(&mut levels, [0.0; 4], None); + levels + }) + // a rim anywhere at all (else the neutral guess applies throughout) + .filter(|levels| levels.last().is_some_and(|top| top.cells.iter().any(|c| c[3] > 0.0))); + let rim_field = rims.as_ref().and_then(|l| l.first()); let max_level = wb.iter().cloned().fold(0.0f32, f32::max) * clip; img.data.par_chunks_mut(w).enumerate().for_each(|(y, row)| { for (x, px) in row.iter_mut().enumerate() { @@ -141,27 +366,45 @@ pub fn reconstruct(img: &mut Rgb32f, wb: [f32; 3], clip: f32) -> usize { if !cl.iter().any(|&b| b) { continue; } - let r = if any_valid { upsample(&cv, cw, ch, x, y) } else { [1f32 / 3.0; 3] }; let q = [p[0] * wb[0], p[1] * wb[1], p[2] * wb[2]]; - // estimate the white-balanced channel sum from the unclipped channels - 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 k == 0 { - let v = q.iter().cloned().fold(0.0f32, f32::max).max(max_level); - out = [v; 3]; - } else { - let sum = sum / k as f32; - for c in 0..3 { - if cl[c] { - out[c] = q[c].max(sum * r[c]); - } + // rebuilt with the colour of the reliable content around, and with neutral, blended by how much of it + // there is (blending the results, not the colours, keeps the estimate consistent with the pixel's own + // unclipped channels all the way) + let e = sample(field, 2.0, x, y); + let trust = e[2].clamp(0.0, 1.0); + let local = (trust > 0.0).then(|| rebuild(p, q, cl, [e[0], 1.0 - e[0] - e[1], e[1]], clip, max_level)).flatten(); + let plain = (trust < 1.0).then(|| rebuild(p, q, cl, neutral, clip, max_level)).flatten().map(|b| { + // only as far as the pixel lies past the rim nearest its clip (in every unclipped channel) + let Some(rim) = rim_field else { return b }; + let edge = sample(rim, 4.0, x, y); + let past = (0..3) + .filter(|&c| !cl[c]) + .map(|c| if edge[c] > 1e-6 && q[c] > 1e-6 { (q[c] / edge[c]).log2().abs() } else { f32::INFINITY }) + .fold(0.0f32, f32::max); + let k = smoothstep(NEAR_RIM, NEAR_RIM + DEPTH, past); + std::array::from_fn(|c| q[c] + (b[c] - q[c]) * k) + }); + let mut out = match (local, plain) { + (Some(a), Some(b)) => std::array::from_fn(|c| b[c] + (a[c] - b[c]) * trust), + (Some(a), None) | (None, Some(a)) => a, + (None, None) => { + let v = q.iter().cloned().fold(0.0f32, f32::max).max(max_level); + white.map(|c| c * v) } + }; + // Fade into the fully clipped neutral on the way from the first channel clipping (no change: continuous + // with the unclipped neighbours) to the last one clipping (all of it: continuous with the fully clipped + // ones). The way is measured in stops: how far the rebuilt clipped channels lie above their clip, + // against how far the lowest unclipped channel still lies below its own. + let above = (0..3).filter(|&c| cl[c]).map(|c| (out[c] / q[c].max(1e-9)).max(1.0).log2()).fold(0.0f32, f32::max); + let lowest = (0..3).filter(|&c| !cl[c]).map(|c| p[c]).fold(f32::MAX, f32::min); + let below = (clip / lowest.max(1e-6)).log2().max(0.0); + let t = if above > 0.0 && lowest < clip { smoothstep(0.5, 1.0, above / (above + below)) } else { 0.0 }; + if t > 0.0 { + let level = (out[0] + out[1] + out[2]) / total; + let full = q.iter().cloned().fold(0.0f32, f32::max).max(max_level); + let v = level + (full - level) * t; + out = std::array::from_fn(|c| out[c] + (white[c] * v - out[c]) * t); } *px = [out[0] / wb[0], out[1] / wb[1], out[2] / wb[2]]; } @@ -173,154 +416,170 @@ pub fn reconstruct(img: &mut Rgb32f, wb: [f32; 3], clip: f32) -> usize { mod tests { use super::*; - // The previous full-resolution implementation, kept as the reference for `reconstruct`. - /// Fill `values` (per-pixel vectors) where `valid` is false from valid neighbours, coarse to fine. - fn fill_invalid_ref(w: usize, h: usize, values: &mut [[f32; 3]], valid: &mut [bool]) { - if w == 0 || h == 0 || valid.iter().all(|&v| v) || !valid.iter().any(|&v| v) { - return; - } - // Build a pyramid of weighted means, fill each level from the next coarser one. - 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]; - { - 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; - } - } - } - if n > 0 { - vrow[x] = s.map(|v| v / n as f32); - okrow[x] = true; - } - } - }); - } - if cw * ch < w * h { - fill_invalid_ref(cw, ch, &mut cv, &mut cval); - } else { - // cannot shrink further (1×1): nothing valid anywhere handled above - return; - } - values.par_chunks_mut(w).zip(valid.par_chunks_mut(w)).enumerate().for_each(|(y, (vrow, okrow))| { - for x in 0..w { - 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); - let (x0, y0) = (fx.floor() as usize, fy.floor() as usize); - let (x1, y1) = ((x0 + 1).min(cw - 1), (y0 + 1).min(ch - 1)); - let (tx, ty) = (fx - x0 as f32, fy - y0 as f32); - let g = |xx: usize, yy: usize| cv[yy * cw + xx]; - let mut v = [0f32; 3]; - for c in 0..3 { - let top = g(x0, y0)[c] + (g(x1, y0)[c] - g(x0, y0)[c]) * tx; - let bot = g(x0, y1)[c] + (g(x1, y1)[c] - g(x0, y1)[c]) * tx; - v[c] = top + (bot - top) * ty; - } - vrow[x] = v; - okrow[x] = true; - } - } - }); + /// Daylight-like multipliers: green clips first, red and blue keep rising. + const WB: [f32; 3] = [2.5, 1.0, 1.6]; + + fn balanced(p: [f32; 3]) -> [f32; 3] { + [p[0] * WB[0], p[1] * WB[1], p[2] * WB[2]] } - /// Reconstruct partially clipped channels. `clip` is the sensor clip level in normalised units (use slightly - /// below 1.0, e.g. 0.99), `wb` the white-balance multipliers that will be applied afterwards. - /// Returns the number of pixels that had at least one clipped channel. - fn reconstruct_ref(img: &mut Rgb32f, wb: [f32; 3], clip: f32) -> usize { - let (w, h) = (img.width, img.height); - let n = w * h; - 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]; - 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 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 { - *r = q.map(|v| v.max(0.0) / s); - *v = true; - } - }); - fill_invalid_ref(w, h, &mut chroma, &mut valid); - let max_level = wb.iter().cloned().fold(0.0f32, f32::max) * clip; - 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) { - return; - } - let q = [p[0] * wb[0], p[1] * wb[1], p[2] * wb[2]]; - // estimate the white-balanced channel sum from the unclipped channels - 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 k == 0 { - let v = q.iter().cloned().fold(0.0f32, f32::max).max(max_level); - out = [v; 3]; + fn sensor(q: [f32; 3]) -> [f32; 3] { + [q[0] / WB[0], q[1] / WB[1], q[2] / WB[2]].map(|v| v.min(1.0)) + } + + fn chroma(q: [f32; 3]) -> [f32; 3] { + let s = q[0] + q[1] + q[2]; + q.map(|v| v / s) + } + + /// A pale blue sky whose green channel clipped, an unclipped strip of the same sky along the top, and a dark + /// branch with purple fringes on both sides standing in the clipped part (white-balanced sky colour returned). + fn branch_scene() -> (Rgb32f, [f32; 3]) { + let sky = [1.10, 1.25, 1.45]; + let img = Rgb32f::from_fn(192, 128, |x, y| { + sensor(if y < 24 { + sky.map(|v| v * 0.7) + } else if (90..100).contains(&x) && y >= 40 { + [0.06, 0.05, 0.03] + } else if (x == 89 || x == 100) && y >= 40 { + [0.75, 0.25, 0.96] } else { - let sum = sum / k as f32; - for c in 0..3 { - if cl[c] { - out[c] = q[c].max(sum * r[c]); - } - } + sky + }) + }); + (img, sky) + } + + /// Issue #523: the clipped sky beside a branch took the colour of the branch's fringes and of the branch + /// itself (magenta, green), which any darkening of the highlights then showed. + #[test] + fn fringes_and_dark_branches_do_not_tint_a_clipped_sky() { + let (mut img, sky) = branch_scene(); + assert!(reconstruct(&mut img, WB, 0.99) > 0); + let want = chroma(sky); + let mut worst = 0f32; + for y in 40..128 { + for x in [84, 86, 88, 101, 103, 105] { + let got = chroma(balanced(img.get(x, y))); + worst = worst.max((0..3).map(|c| (got[c] - want[c]).abs()).fold(0.0, f32::max)); } - *px = [out[0] / wb[0], out[1] / wb[1], out[2] / wb[2]]; - }); - count + } + assert!(worst < 0.01, "clipped sky beside the branch is off the sky's colour by {worst}"); } - /// `cargo test --release -p lightcraft-raw --lib -- --ignored bench_reconstruct --nocapture` + /// A warm light with red and green clipped: left at their clip levels the pair reads red-magenta (the white + /// balance multiplies red's clip level by 2.5, green's by 1). #[test] - #[ignore] - fn bench_reconstruct() { - let img = Rgb32f::from_fn(6000, 4000, |x, y| { - let v = 0.4 + 0.8 * ((x as f32 * 0.003).sin() * (y as f32 * 0.002).cos()).abs(); - [v.min(1.0), (v * 0.7).min(1.0), (v * 0.5).min(1.0)] + fn two_clipped_channels_are_not_left_at_their_clip_levels() { + let light = [3.5, 2.6, 0.8]; + let mut img = Rgb32f::from_fn(64, 64, |x, y| { + let d = ((x as f32 - 32.0).powi(2) + (y as f32 - 32.0).powi(2)).sqrt(); + if d < 20.0 { sensor(light) } else { [0.02; 3] } }); - let wb = [2.0, 1.0, 1.5]; - let time = |f: &dyn Fn(&mut Rgb32f) -> usize| { - (0..3) - .map(|_| { - let mut i = img.clone(); - let t = std::time::Instant::now(); - f(&mut i); - t.elapsed().as_secs_f64() * 1e3 - }) - .fold(f64::MAX, f64::min) - }; - let new = time(&|i| reconstruct(i, wb, 0.99)); - let old = time(&|i| reconstruct_ref(i, wb, 0.99)); - eprintln!("reconstruct 24 MP: {new:.0} ms (full-resolution reference: {old:.0} ms)"); + reconstruct(&mut img, WB, 0.99); + let q = balanced(img.get(32, 32)); + assert!((q[2] - 0.8).abs() < 0.1, "the unclipped channel is about kept: {q:?}"); + let ratio = q[1] / q[0]; + assert!((0.6..=1.05).contains(&ratio), "green / red {ratio} ({q:?}): not a warm white"); } - /// The half-resolution start gives exactly the full-resolution algorithm's result. #[test] - fn matches_the_full_resolution_reference() { + fn fully_clipped_pixels_follow_the_colour_model_white() { + let white = [1.03, 0.98, 1.06]; + let mut img = Rgb32f::filled(8, 8, [1.0; 3]); + img.data[0] = [0.3, 0.6, 0.4]; + reconstruct_with(&mut img, WB, 0.99, white); + let q = balanced(img.get(4, 4)); + let k = q[1] / white[1]; + assert!((0..3).all(|c| (q[c] - white[c] * k).abs() < 1e-4 * k), "{q:?} is not along {white:?}"); + assert!(k >= 2.5 * 0.99 - 1e-4, "clipped neutral at least as bright as the highest clip level: {k}"); + assert_eq!(img.data[0], [0.3, 0.6, 0.4]); + // a white that cannot be right falls back to camera neutral, and nothing overflows + for bad in [[f32::NAN, 1.0, 1.0], [2.5e38, 0.74, 0.74], [1.0, 0.0, 1.0], [40.0, 1.0, 1.0]] { + let mut img = Rgb32f::filled(4, 4, [1.0; 3]); + img.data[0] = [1.0, 1.0, 0.5]; + reconstruct_with(&mut img, WB, 0.99, bad); + let q = balanced(img.get(1, 1)); + assert!((q[0] - q[1]).abs() < 1e-4 && (q[2] - q[1]).abs() < 1e-4, "{bad:?}: {q:?}"); + assert!(img.data.iter().flatten().all(|v| v.is_finite()), "{bad:?}"); + } + } + + /// Red rises through its clip level across a sky whose green and blue clipped: the rebuilt colour changes + /// smoothly (no edge where red clips too) and arrives at the clipped neutral. + #[test] + fn rebuilt_colour_is_continuous_as_more_channels_clip() { + let blue = [1.0, 1.3, 1.9]; + let img0 = Rgb32f::from_fn(240, 40, |x, y| { + if y < 8 { + return sensor(blue.map(|v| v * 0.5)); // unclipped sky: the source of colour + } + let level = 0.6 + 2.2 * x as f32 / 239.0; // red at 0.24 to 1.12 of its clip level + sensor([level * WB[0] * 0.4, 1.6, 2.0]) + }); + let mut img = img0.clone(); + reconstruct(&mut img, WB, 0.99); + let row: Vec<[f32; 3]> = (0..240).map(|x| chroma(balanced(img.get(x, 30)))).collect(); + let step = row.windows(2).map(|w| (0..3).map(|c| (w[1][c] - w[0][c]).abs()).fold(0.0, f32::max)).fold(0.0, f32::max); + assert!(step < 0.01, "largest chromaticity step between neighbours {step}"); + let end = row[239]; + assert!(end.iter().all(|v| (v - 1.0 / 3.0).abs() < 1e-3), "fully clipped end {end:?}"); + } + + /// A smooth pink ramp (`[0.9k, k, 0.8k]`, k from 0.8 to 1.35) whose green clips first, then red, then blue, + /// in `rows` of a `h`-row picture that is otherwise near black: the rebuilt middle row, and where green first + /// clips on it. + fn rebuilt_ramp(rows: std::ops::Range, h: usize) -> (Vec<[f32; 3]>, usize) { + let (w, k0, k1) = (1600usize, 0.8f32, 1.35f32); + let k = |x: usize| k0 + (k1 - k0) * x as f32 / (w - 1) as f32; + let img0 = Rgb32f::from_fn(w, h, |x, y| if rows.contains(&y) { [0.9 * k(x), k(x), 0.8 * k(x)].map(|v| v.min(1.0)) } else { [0.01; 3] }); + let mut img = img0.clone(); + reconstruct(&mut img, WB, 0.99); + let y = (rows.start + rows.end) / 2; + let row: Vec<[f32; 3]> = (0..w).map(|x| balanced(img.get(x, y))).collect(); + (row, (0..w).find(|&x| img0.get(x, y)[1] >= 0.99).unwrap_or(0)) + } + + /// No channel changes by more than `at_clip` (relative) where green first clips, nor by 3 % anywhere. + fn assert_continuous(row: &[[f32; 3]], first: usize, at_clip: f32, case: &str) { + // across the first clip the colour carries on as it was (the sensor's own colour, not a fade) + let (a, b) = (row[first - 1], row[first]); + assert!((0..3).all(|c| (b[c] / a[c] - 1.0).abs() < at_clip), "{case}, first clip: {a:?} then {b:?}"); + let worst = row.windows(2).map(|p| (0..3).map(|c| (p[1][c] / p[0][c] - 1.0).abs()).fold(0.0, f32::max)).fold(0.0, f32::max); + // (a smooth ramp moves 0.04 % per pixel here; the rebuilt colour turns faster where it leaves the colour of + // the unclipped part for neutral, but never by a step) + assert!(worst < 0.03, "{case}: largest relative step between neighbours {worst}"); + } + + /// The rebuilt colour has no step anywhere along the ramp: not where green first clips (the neutral fade starts + /// from nothing there), nor where red does, nor where the last channel clips (the fade ends at the fully + /// clipped neutral). + #[test] + fn rebuilt_ramp_is_continuous_from_the_first_clip_to_the_last() { + let (row, first) = rebuilt_ramp(0..32, 32); + assert_continuous(&row, first, 0.002, "full ramp"); + } + + /// The same ramp as a thin line on near black has no reliable colour around (every pixel of it is on a strong + /// edge), so its clipped part falls back to the neutral guess. That guess only applies as far as the line runs + /// past its own unclipped part: no step where it first clips, however thin the line and however much dark + /// background surrounds it. + #[test] + fn thin_line_without_reliable_colour_is_continuous_where_it_clips() { + for height in [32usize, 512, 1024] { + for width in [1usize, 4, 16] { + let top = (height - width) / 2; + let (row, first) = rebuilt_ramp(top..top + width, height); + assert_continuous(&row, first, 0.002, &format!("{width} rows of {height}")); + // and well past the first clip the neutral guess has taken over (green as high as red and blue imply) + let far = row[1500]; + assert!(far[1] > 0.95 * far[2], "{width} rows of {height}: {far:?}"); + } + } + } + + #[test] + fn unclipped_pixels_are_left_exactly() { let mut seed = 0x2545_f491_4f6c_dd1du64; let mut rnd = move || { seed ^= seed << 13; @@ -328,28 +587,28 @@ mod tests { seed ^= seed << 17; (seed >> 40) as f32 / (1u64 << 24) as f32 }; - for (w, h) in [(1, 1), (1, 2), (3, 1), (7, 5), (64, 48), (129, 77)] { - for density in [0.02f32, 0.3, 0.97] { - // smooth colour with clipped blobs, some dark pixels, some fully clipped ones - let noise: Vec = (0..w * h).map(|_| rnd()).collect(); - let img = Rgb32f::from_fn(w, h, |x, y| { - let t = ((x as f32 * 0.21).sin() + (y as f32 * 0.17).cos()) * 0.5; - let base = [0.5 + 0.4 * t, 0.45 - 0.2 * t, 0.3 + 0.1 * t]; - let r = noise[y * w + x]; - if r < density * 0.3 { - [1.0, 1.0, 1.0] - } else if r < density { - [1.0, base[1] * 1.6, base[2]] - } else if r > 0.995 { - [0.0, 0.0, 0.00001] - } else { - base - } - }); - let wb = [2.1, 1.0, 1.6]; - let (mut a, mut b) = (img.clone(), img.clone()); - assert_eq!(reconstruct(&mut a, wb, 0.99), reconstruct_ref(&mut b, wb, 0.99)); - assert!(a.data == b.data, "{w}×{h} at {density}"); + for (w, h) in [(1, 1), (3, 1), (7, 5), (64, 48), (129, 77)] { + let noise: Vec = (0..w * h).map(|_| rnd()).collect(); + let img = Rgb32f::from_fn(w, h, |x, y| { + let t = ((x as f32 * 0.21).sin() + (y as f32 * 0.17).cos()) * 0.5; + let base = [0.5 + 0.4 * t, 0.45 - 0.2 * t, 0.3 + 0.1 * t]; + match noise[y * w + x] { + r if r < 0.1 => [1.0; 3], + r if r < 0.4 => [1.0, base[1] * 1.6, base[2]], + r if r > 0.995 => [0.0, 0.0, 0.00001], + _ => base, + } + }); + let mut out = img.clone(); + let n = reconstruct_with(&mut out, [2.1, 1.0, 1.6], 0.99, [1.02, 0.99, 1.05]); + let clipped = img.data.iter().filter(|p| p.iter().any(|&v| v >= 0.99)).count(); + assert_eq!(n, clipped); + for (a, b) in img.data.iter().zip(&out.data) { + if a.iter().all(|&v| v < 0.99) { + assert_eq!(a, b, "{w}×{h}"); + } else { + assert!(b.iter().all(|v| v.is_finite() && *v >= 0.0), "{w}×{h}: {b:?}"); + } } } } @@ -368,19 +627,23 @@ mod tests { #[test] fn reconstruct_restores_clipped_channel() { - // a warm gradient whose red channel clips in the bright half - let truth = Rgb32f::from_fn(64, 16, |x, _| { - let v = 0.3 + 1.2 * x as f32 / 63.0; + // a warm gradient whose red channel clips at its bright end + let truth = Rgb32f::from_fn(256, 16, |x, _| { + let v = 0.3 + 1.0 * x as f32 / 255.0; [v * 1.0, v * 0.6, v * 0.35] }); let mut img = truth.map(|p| p.map(|v| v.min(1.0))); - let before: f32 = img.data.iter().zip(&truth.data).map(|(a, b)| (a[0] - b[0]).abs()).sum(); + let first = (0..256).find(|&x| img.get(x, 8)[0] >= 0.999).unwrap(); + // near the unclipped gradient (within the reach of its colour) red is rebuilt from it + let err = |img: &Rgb32f| -> f32 { (first..first + 12).map(|x| (img.get(x, 8)[0] - truth.get(x, 8)[0]).abs()).sum() }; + let before = err(&img); let n = reconstruct(&mut img, [1.0, 1.0, 1.0], 0.999); assert!(n > 0); - let after: f32 = img.data.iter().zip(&truth.data).map(|(a, b)| (a[0] - b[0]).abs()).sum(); + let after = err(&img); assert!(after < before * 0.2, "before {before} after {after}"); // unclipped pixels untouched assert_eq!(img.get(0, 0), truth.get(0, 0)); + assert_eq!(img.get(first - 1, 8), truth.get(first - 1, 8)); } #[test] @@ -394,4 +657,23 @@ mod tests { let mut empty = Rgb32f::new(0, 0); assert_eq!(reconstruct(&mut empty, [1.0; 3], 0.99), 0); } + + /// `cargo test --release -p lightcraft-raw --lib -- --ignored bench_reconstruct --nocapture` + #[test] + #[ignore] + fn bench_reconstruct() { + let img = Rgb32f::from_fn(6000, 4000, |x, y| { + let v = 0.4 + 0.8 * ((x as f32 * 0.003).sin() * (y as f32 * 0.002).cos()).abs(); + [v.min(1.0), (v * 0.7).min(1.0), (v * 0.5).min(1.0)] + }); + let ms = (0..3) + .map(|_| { + let mut i = img.clone(); + let t = std::time::Instant::now(); + reconstruct(&mut i, [2.0, 1.0, 1.5], 0.99); + t.elapsed().as_secs_f64() * 1e3 + }) + .fold(f64::MAX, f64::min); + eprintln!("reconstruct 24 MP: {ms:.0} ms"); + } } diff --git a/docs/parity.md b/docs/parity.md index 4c27755b..8c243fd5 100644 --- a/docs/parity.md +++ b/docs/parity.md @@ -241,7 +241,7 @@ under the UI command that wraps it (a no-op error); the UI command now wins. | LR-EDIT-HDR-MODE | HDR editing | P2 | ⬜ | | | | LR-EDIT-LIGHT-EXPOSURE | Exposure | P0 | ✅ | `ctl:light.exposure` | | | LR-EDIT-LIGHT-CONTRAST | Contrast | P0 | ✅ | `ctl:light.contrast` | | -| LR-EDIT-LIGHT-HIGHLIGHTS | Highlights | P0 | ✅ | `ctl:light.highlights` | | +| LR-EDIT-LIGHT-HIGHLIGHTS | Highlights | P0 | ✅ | `ctl:light.highlights`, `crates/raw/src/highlight.rs` | sensor-clipped highlights stay neutral when Highlights (or Exposure, Whites) darkens them: fully clipped content is rebuilt to the neutral of the colour model applied after it, partly clipped content takes its colour from reliable unclipped content nearby (not from dark objects or lens fringes beside it) and white where there is none (#523) | | LR-EDIT-LIGHT-SHADOWS | Shadows | P0 | ✅ | `ctl:light.shadows` | | | LR-EDIT-LIGHT-WHITES | Whites | P0 | ✅ | `ctl:light.whites` | | | LR-EDIT-LIGHT-BLACKS | Blacks | P0 | ✅ | `ctl:light.blacks` | | diff --git a/xtask/src/main.rs b/xtask/src/main.rs index cfdec92a..4e633293 100644 --- a/xtask/src/main.rs +++ b/xtask/src/main.rs @@ -340,6 +340,12 @@ const RAW_SAMPLES: &[(&str, &str, &str)] = &[ "https://raw.pixls.us/data/Sony/ILCE-7M4/ILCE-7M4_DSC06676_FullFrame-LossLess-Compressed-Small.ARW", "cbbd0930c7d8706dff84c68a2004454266e6fd0d8354f5f76a106b5d776e0223", ), + // a dusk sky clipped in green behind a poplar: the clipped-highlight colour of issue #523 + ( + "arw-sony-a7rm4-14bit-compressed.arw", + "https://raw.pixls.us/getfile.php/3480/nice/Sony%20-%20ILCE-7RM4%20-%2014bit%2014bit%20compressed%20%283:2%29.ARW", + "e6dafe42643f69ab9d1fd00414b7a1f5104df354bb589201600ac934b29b5e4a", + ), // pre-2017 bodies: white balance only in the enciphered maker note, black level only in the SR2SubIFD (#148) ( "arw-sony-rx100m3.arw",