From 543fdde41e06f536058d318e98ca58d7817abd91 Mon Sep 17 00:00:00 2001 From: Brandon Thomas Date: Thu, 1 Oct 2026 10:06:00 -0400 Subject: [PATCH 1/3] =?UTF-8?q?M13.1:=20lightcraft-merge=20crate=20?= =?UTF-8?q?=E2=80=94=20HDR=20merge,=20panorama=20stitching=20core?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit New L3 crate `crates/merge` (no UI deps), clean-room from published papers: - features: Harris/MOPS-style oriented patch descriptors, ANMS, ratio-test matching - ransac: normalised DLT homographies, similarity/translation, adaptive RANSAC - hdr: brightness ordering, auto-align (model selection), robust exposure-ratio estimation chained to the reference, saturation-aware weighted radiance merge, deghost None/Low/Medium/High with feathered masks + overlay - pano: pairwise matching, Brown–Lowe verification, focal from homographies, rotation+focal bundle adjustment (LM), straightening, spherical/cylindrical/ perspective + auto, gain compensation, max-weight seams, multi-band blending in log space, boundary warp, auto crop, fill edges (push–pull) - output: LinearRaw DNG carrying source colour tags; raw crate gains Deflate tiles (16/32-bit float with FP predictor, 16-bit int with predictor 2) Tests: synthetic brackets (alignment < 0.3 px, EV < 0.02, radiance error noise-limited, ghost removal) and synthetic panoramas (relative rotations < 0.1°, focal < 1 %, error vs ground truth, seam test). Co-Authored-By: Claude Opus 5.5 --- Cargo.lock | 16 ++ Cargo.toml | 1 + crates/merge/Cargo.toml | 23 ++ crates/merge/src/align.rs | 142 ++++++++++ crates/merge/src/features.rs | 312 +++++++++++++++++++++ crates/merge/src/frame.rs | 200 ++++++++++++++ crates/merge/src/hdr.rs | 449 ++++++++++++++++++++++++++++++ crates/merge/src/lib.rs | 57 ++++ crates/merge/src/linalg.rs | 247 +++++++++++++++++ crates/merge/src/output.rs | 142 ++++++++++ crates/merge/src/pano/blend.rs | 300 ++++++++++++++++++++ crates/merge/src/pano/camera.rs | 446 ++++++++++++++++++++++++++++++ crates/merge/src/pano/finish.rs | 210 ++++++++++++++ crates/merge/src/pano/mod.rs | 461 +++++++++++++++++++++++++++++++ crates/merge/src/pano/project.rs | 83 ++++++ crates/merge/src/ransac.rs | 291 +++++++++++++++++++ crates/merge/tests/hdr.rs | 163 +++++++++++ crates/merge/tests/pano.rs | 247 +++++++++++++++++ crates/raw/src/dngwrite.rs | 121 ++++++++ crates/raw/tests/dng.rs | 34 +++ 20 files changed, 3945 insertions(+) create mode 100644 crates/merge/Cargo.toml create mode 100644 crates/merge/src/align.rs create mode 100644 crates/merge/src/features.rs create mode 100644 crates/merge/src/frame.rs create mode 100644 crates/merge/src/hdr.rs create mode 100644 crates/merge/src/lib.rs create mode 100644 crates/merge/src/linalg.rs create mode 100644 crates/merge/src/output.rs create mode 100644 crates/merge/src/pano/blend.rs create mode 100644 crates/merge/src/pano/camera.rs create mode 100644 crates/merge/src/pano/finish.rs create mode 100644 crates/merge/src/pano/mod.rs create mode 100644 crates/merge/src/pano/project.rs create mode 100644 crates/merge/src/ransac.rs create mode 100644 crates/merge/tests/hdr.rs create mode 100644 crates/merge/tests/pano.rs diff --git a/Cargo.lock b/Cargo.lock index 01e13890..9a384cca 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -2144,6 +2144,22 @@ dependencies = [ "serde_json", ] +[[package]] +name = "lightcraft-merge" +version = "0.1.0" +dependencies = [ + "lightcraft-codecs", + "lightcraft-color", + "lightcraft-geom", + "lightcraft-meta", + "lightcraft-raster", + "lightcraft-raw", + "lightcraft-scenes", + "rayon", + "serde", + "thiserror 2.0.21", +] + [[package]] name = "lightcraft-meta" version = "0.1.0" diff --git a/Cargo.toml b/Cargo.toml index 76abfa4c..836fd417 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -24,6 +24,7 @@ lightcraft-pipeline = { path = "crates/pipeline" } lightcraft-gpu = { path = "crates/gpu" } lightcraft-catalog = { path = "crates/catalog" } lightcraft-preview = { path = "crates/preview" } +lightcraft-merge = { path = "crates/merge" } lightcraft-engine = { path = "crates/engine" } lightcraft-ui-egui = { path = "crates/ui-egui" } lightcraft-mcp = { path = "crates/mcp" } diff --git a/crates/merge/Cargo.toml b/crates/merge/Cargo.toml new file mode 100644 index 00000000..f0f5fb1f --- /dev/null +++ b/crates/merge/Cargo.toml @@ -0,0 +1,23 @@ +[package] +name = "lightcraft-merge" +version.workspace = true +edition.workspace = true +license.workspace = true +description = "LightCraft Photo Merge: HDR merge (align, deghost) and panorama stitching (L3)" + +[lints] +workspace = true + +[dependencies] +lightcraft-geom = { workspace = true } +lightcraft-color = { workspace = true } +lightcraft-raster = { workspace = true } +lightcraft-raw = { workspace = true } +lightcraft-codecs = { workspace = true } +lightcraft-meta = { workspace = true } +serde = { workspace = true } +thiserror = { workspace = true } +rayon = { workspace = true } + +[dev-dependencies] +lightcraft-scenes = { workspace = true } diff --git a/crates/merge/src/align.rs b/crates/merge/src/align.rs new file mode 100644 index 00000000..4298f810 --- /dev/null +++ b/crates/merge/src/align.rs @@ -0,0 +1,142 @@ +//! Frame alignment for HDR brackets (handheld, small motion): features on exposure-normalised, +//! gamma-encoded luminance at ≤ 1200 px, matched to the reference frame, then a RANSAC homography +//! (falling back to similarity, then translation, then identity when the scene gives too little). + +use lightcraft_geom::{Homography, Point}; +use lightcraft_raster::resample::{Filter, fit}; +use lightcraft_raster::{Plane, Rgb32f}; +use rayon::prelude::*; + +use crate::features::{self, Features}; +use crate::ransac::{self, Model}; + +/// Long edge of the images features are detected on. +pub const WORK_EDGE: usize = 1200; + +/// A grey, perceptually encoded work image: `(luminance · gain)^(1/2.2)`, clamped to 0..1, and +/// the scale from full-resolution pixels to work pixels. +pub fn work_image(img: &Rgb32f, gain: f32, max_edge: usize) -> (Plane, f64) { + let small = if img.width.max(img.height) > max_edge { fit(img, max_edge, max_edge, Filter::Box) } else { img.clone() }; + let scale = small.width as f64 / img.width as f64; + let p = small.map(|c| { + let y = (0.25 * c[0] + 0.6 * c[1] + 0.15 * c[2]) * gain; + y.clamp(0.0, 1.0).powf(1.0 / 2.2) + }); + (p, scale) +} + +/// Median luminance of the (unsaturated, non-black) pixels: a brightness estimate for frames +/// without EXIF. +pub fn median_luma(img: &Rgb32f, clip: f32) -> f32 { + let step = (img.data.len() / 100_000).max(1); + let mut v: Vec = + img.data.iter().step_by(step).filter(|p| p[0].max(p[1]).max(p[2]) < clip * 0.98).map(|c| 0.25 * c[0] + 0.6 * c[1] + 0.15 * c[2]).collect(); + if v.is_empty() { + return clip; + } + v.sort_by(|a, b| a.total_cmp(b)); + v[v.len() / 2].max(1e-6) +} + +#[derive(Clone, Debug)] +pub struct Alignment { + /// Maps reference pixel coordinates to this frame's coordinates (full resolution). + pub h: Homography, + pub model: &'static str, + pub inliers: usize, + pub rms: f64, +} + +impl Alignment { + pub fn identity() -> Alignment { + Alignment { h: Homography::IDENTITY, model: "identity", inliers: 0, rms: 0.0 } + } +} + +/// Conjugate a work-resolution homography to full resolution: S⁻¹ · H · S with S = diag(s, s, 1). +pub fn scale_homography(h: &Homography, s: f64) -> Homography { + let m = h.0; + Homography([m[0], m[1], m[2] / s, m[3], m[4], m[5] / s, m[6] * s, m[7] * s, m[8]]) +} + +/// Whether `h` is a plausible handheld motion for a `w × h` frame: corners move less than +/// `max_frac` of the diagonal and the area changes by < 20 %. +pub fn plausible(h: &Homography, w: f64, ht: f64, max_frac: f64) -> bool { + let diag = (w * w + ht * ht).sqrt(); + let corners = [Point::new(0.0, 0.0), Point::new(w, 0.0), Point::new(w, ht), Point::new(0.0, ht)]; + let moved: Vec = corners.iter().map(|&c| h.apply(c)).collect(); + if corners.iter().zip(&moved).any(|(a, b)| a.dist(*b) > max_frac * diag || !b.x.is_finite() || !b.y.is_finite()) { + return false; + } + let area = |p: &[Point]| { + let mut s = 0.0; + for i in 0..p.len() { + let (a, b) = (p[i], p[(i + 1) % p.len()]); + s += a.x * b.y - b.x * a.y; + } + s.abs() / 2.0 + }; + let r = area(&moved) / (w * ht); + (0.8..1.25).contains(&r) +} + +/// Align `frame` to `reference` (both as work images with features already detected). +pub fn align_pair(ref_feat: &Features, feat: &Features, scale: f64, full_w: usize, full_h: usize, seed: u64) -> Alignment { + let m = features::match_features(ref_feat, feat, 0.8); + let a: Vec = m.iter().map(|(i, _)| Point::new(ref_feat.points[*i].x as f64, ref_feat.points[*i].y as f64)).collect(); + let b: Vec = m.iter().map(|(_, j)| Point::new(feat.points[*j].x as f64, feat.points[*j].y as f64)).collect(); + let (ww, wh) = (full_w as f64 * scale, full_h as f64 * scale); + // Fit every model; keep the simplest one that explains (almost) as many matches as the most + // general one — extra degrees of freedom fitted to clustered points extrapolate badly. + let mut fits: Vec<(&'static str, ransac::Fit)> = Vec::new(); + for (model, name, min) in [(Model::Translation, "translation", 5), (Model::Similarity, "similarity", 8), (Model::Homography, "homography", 16)] { + if let Some(f) = ransac::ransac(model, &a, &b, 1.5, min, seed) + && plausible(&f.h, ww, wh, 0.12) + { + fits.push((name, f)); + } + } + let Some(most) = fits.iter().map(|(_, f)| f.inliers.len()).max() else { return Alignment::identity() }; + let (name, f) = fits.into_iter().find(|(_, f)| f.inliers.len() as f64 >= most as f64 * 0.97).expect("the best fit qualifies"); + Alignment { h: scale_homography(&f.h, scale), model: name, inliers: f.inliers.len(), rms: f.rms / scale } +} + +/// Resample `img` into the reference geometry: `out(p) = img(h(p))`. Returns the warped image and +/// a validity mask (false where `h(p)` falls outside `img`). +pub fn warp(img: &Rgb32f, h: &Homography) -> (Rgb32f, Vec) { + let (w, ht) = (img.width, img.height); + if *h == Homography::IDENTITY { + return (img.clone(), vec![true; w * ht]); + } + let mut out = Rgb32f::new(w, ht); + let mut valid = vec![false; w * ht]; + out.data.par_chunks_mut(w).zip(valid.par_chunks_mut(w)).enumerate().for_each(|(y, (row, vrow))| { + for x in 0..w { + let p = h.apply(Point::new(x as f64 + 0.5, y as f64 + 0.5)); + let (sx, sy) = (p.x as f32, p.y as f32); + if sx >= 0.0 && sy >= 0.0 && sx <= w as f32 && sy <= ht as f32 { + row[x] = img.sample_bilinear(sx, sy); + vrow[x] = true; + } + } + }); + (out, valid) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn homography_scaling_conjugates() { + let h = Homography([1.01, 0.02, 5.0, -0.01, 0.99, -3.0, 1e-5, 2e-5, 1.0]); + let s = 0.25; + let full = scale_homography(&h, s); + let p = Point::new(400.0, 300.0); + let w = h.apply(Point::new(p.x * s, p.y * s)); + let f = full.apply(p); + assert!((f.x * s - w.x).abs() < 1e-9 && (f.y * s - w.y).abs() < 1e-9); + assert!(plausible(&Homography::IDENTITY, 100.0, 100.0, 0.1)); + assert!(!plausible(&Homography([1.0, 0.0, 50.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]), 100.0, 100.0, 0.1)); + } +} diff --git a/crates/merge/src/features.rs b/crates/merge/src/features.rs new file mode 100644 index 00000000..353bed8f --- /dev/null +++ b/crates/merge/src/features.rs @@ -0,0 +1,312 @@ +//! Feature points and descriptors (our own implementation, after the published MOPS method: +//! M. Brown, R. Szeliski, S. Winder, "Multi-Image Matching using Multi-Scale Oriented Patches", +//! CVPR 2005). +//! +//! - **Detector:** Harris corners on a Gaussian pyramid (octaves of ×½). The corner strength is the +//! harmonic mean of the structure tensor's eigenvalues, `det / trace` (MOPS); integration scale +//! σᵢ = 1.5, derivative scale σd = 1.0. Local maxima over 3×3 above a threshold, refined to +//! sub-pixel by a parabola fit in x and y. +//! - **Selection:** adaptive non-maximal suppression (ANMS) keeps the `max` points with the largest +//! suppression radius (distance to the nearest point that is ≥ 10 % stronger), which spreads +//! points evenly over the image. +//! - **Orientation:** the direction of the gradient of the level blurred with σ = 4.5. +//! - **Descriptor:** 8 × 8 samples with 5-pixel spacing (a 40 × 40 window) taken from the level +//! blurred with σ = 2.5, on a grid rotated by the orientation; normalised to zero mean and unit +//! variance (bias/gain invariant) — 64 floats. +//! - **Matching:** brute-force nearest neighbour in L2 with Lowe's ratio test, kept only when +//! mutually consistent. + +use lightcraft_raster::Plane; +use lightcraft_raster::blur::gaussian; +use lightcraft_raster::resample::half; +use rayon::prelude::*; + +pub const DESC_LEN: usize = 64; + +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct Keypoint { + /// Position in pixels of the input image (pixel centres at +0.5). + pub x: f32, + pub y: f32, + /// Pyramid level (0 = full resolution). + pub level: u8, + /// Orientation (radians). + pub angle: f32, + pub strength: f32, +} + +#[derive(Clone, Debug, Default)] +pub struct Features { + pub points: Vec, + pub desc: Vec<[f32; DESC_LEN]>, +} + +impl Features { + pub fn len(&self) -> usize { + self.points.len() + } + pub fn is_empty(&self) -> bool { + self.points.is_empty() + } +} + +/// Detect up to `max` features in a grey image (values roughly perceptual, e.g. gamma-encoded 0..1). +pub fn detect(img: &Plane, max: usize) -> Features { + let mut levels = vec![img.clone()]; + while levels.len() < 4 { + let l = levels.last().expect("non-empty"); + if l.width.min(l.height) < 96 { + break; + } + levels.push(half(l)); + } + let per_level: Vec<(Vec, Vec<[f32; DESC_LEN]>)> = levels + .par_iter() + .enumerate() + .map(|(li, lvl)| { + let cands = harris(lvl); + let orient_img = gaussian(lvl, 4.5); + let desc_img = gaussian(lvl, 2.5); + let (gx, gy) = gradients(&orient_img); + let scale = (1u32 << li) as f32; + let mut pts = Vec::new(); + let mut descs = Vec::new(); + for (x, y, s) in cands { + let (xi, yi) = (x as usize, y as usize); + let (dx, dy) = (gx.get(xi.min(gx.width - 1), yi.min(gx.height - 1)), gy.get(xi.min(gx.width - 1), yi.min(gx.height - 1))); + let angle = dy.atan2(dx); + let Some(d) = descriptor(&desc_img, x, y, angle) else { continue }; + pts.push(Keypoint { x: x * scale, y: y * scale, level: li as u8, angle, strength: s }); + descs.push(d); + } + (pts, descs) + }) + .collect(); + let mut all: Vec<(Keypoint, [f32; DESC_LEN])> = per_level.into_iter().flat_map(|(p, d)| p.into_iter().zip(d)).collect(); + all.sort_by(|a, b| b.0.strength.total_cmp(&a.0.strength)); + all.truncate(max.saturating_mul(8).max(64)); + let keep = anms(&all.iter().map(|(k, _)| *k).collect::>(), max); + let mut f = Features::default(); + for i in keep { + f.points.push(all[i].0); + f.desc.push(all[i].1); + } + f +} + +fn gradients(p: &Plane) -> (Plane, Plane) { + let (w, h) = (p.width, p.height); + let gx = Plane::from_fn(w, h, |x, y| (p.get_clamped(x as isize + 1, y as isize) - p.get_clamped(x as isize - 1, y as isize)) * 0.5); + let gy = Plane::from_fn(w, h, |x, y| (p.get_clamped(x as isize, y as isize + 1) - p.get_clamped(x as isize, y as isize - 1)) * 0.5); + (gx, gy) +} + +/// Harris corners of one level: (x, y in level pixel coordinates (centre +0.5), strength). +fn harris(lvl: &Plane) -> Vec<(f32, f32, f32)> { + let (w, h) = (lvl.width, lvl.height); + if w < 24 || h < 24 { + return vec![]; + } + let s = gaussian(lvl, 1.0); + let (gx, gy) = gradients(&s); + let xx = gaussian(&gx.zip_map(&gx, |a, b| a * b), 1.5); + let yy = gaussian(&gy.zip_map(&gy, |a, b| a * b), 1.5); + let xy = gaussian(&gx.zip_map(&gy, |a, b| a * b), 1.5); + let mut r = Plane::new(w, h); + for i in 0..r.data.len() { + let (a, b, c) = (xx.data[i], yy.data[i], xy.data[i]); + let tr = a + b; + r.data[i] = if tr > 1e-12 { (a * b - c * c) / tr } else { 0.0 }; + } + // threshold relative to the level's strongest response (robust to image contrast) + let mut sorted: Vec = r.data.iter().copied().step_by(7).collect(); + sorted.sort_by(|a, b| b.total_cmp(a)); + let top = sorted.get(sorted.len() / 200).copied().unwrap_or(0.0); + let thr = (top * 0.01).max(1e-7); + let border = 22usize; + let rows: Vec> = (border..h.saturating_sub(border)) + .into_par_iter() + .map(|y| { + let mut out = Vec::new(); + for x in border..w.saturating_sub(border) { + let v = r.get(x, y); + if v <= thr { + continue; + } + let mut is_max = true; + 'n: for dy in -1i32..=1 { + for dx in -1i32..=1 { + if (dx, dy) != (0, 0) { + let o = r.get((x as i32 + dx) as usize, (y as i32 + dy) as usize); + if o > v || (o == v && (dy < 0 || (dy == 0 && dx < 0))) { + is_max = false; + break 'n; + } + } + } + } + if !is_max { + continue; + } + let sub = |a: f32, b: f32, c: f32| { + let d = a - 2.0 * b + c; + if d.abs() > 1e-12 { (0.5 * (a - c) / d).clamp(-0.5, 0.5) } else { 0.0 } + }; + let ox = sub(r.get(x - 1, y), v, r.get(x + 1, y)); + let oy = sub(r.get(x, y - 1), v, r.get(x, y + 1)); + out.push((x as f32 + 0.5 + ox, y as f32 + 0.5 + oy, v)); + } + out + }) + .collect(); + rows.into_iter().flatten().collect() +} + +/// Adaptive non-maximal suppression: indices (into `pts`, sorted by decreasing strength) of the `max` +/// points with the largest suppression radii. +fn anms(pts: &[Keypoint], max: usize) -> Vec { + if pts.len() <= max { + return (0..pts.len()).collect(); + } + let radii: Vec = (0..pts.len()) + .into_par_iter() + .map(|i| { + let p = pts[i]; + let mut best = f32::INFINITY; + for q in &pts[..i] { + if q.strength * 0.9 > p.strength { + let d = (q.x - p.x).powi(2) + (q.y - p.y).powi(2); + best = best.min(d); + } + } + best + }) + .collect(); + let mut idx: Vec = (0..pts.len()).collect(); + idx.sort_by(|&a, &b| radii[b].total_cmp(&radii[a]).then(a.cmp(&b))); + idx.truncate(max); + idx +} + +fn descriptor(img: &Plane, x: f32, y: f32, angle: f32) -> Option<[f32; DESC_LEN]> { + let (s, c) = angle.sin_cos(); + let mut d = [0f32; DESC_LEN]; + let (w, h) = (img.width as f32, img.height as f32); + for j in 0..8 { + for i in 0..8 { + let u = (i as f32 - 3.5) * 5.0; + let v = (j as f32 - 3.5) * 5.0; + let sx = x + c * u - s * v; + let sy = y + s * u + c * v; + if sx < 0.5 || sy < 0.5 || sx > w - 0.5 || sy > h - 0.5 { + return None; + } + d[j * 8 + i] = bilinear(img, sx, sy); + } + } + let mean = d.iter().sum::() / DESC_LEN as f32; + let var = d.iter().map(|v| (v - mean) * (v - mean)).sum::() / DESC_LEN as f32; + if var < 1e-10 { + return None; + } + let inv = 1.0 / var.sqrt(); + for v in d.iter_mut() { + *v = (*v - mean) * inv; + } + Some(d) +} + +#[inline] +pub fn bilinear(img: &Plane, x: f32, y: f32) -> f32 { + let fx = x - 0.5; + let fy = y - 0.5; + let x0 = fx.floor(); + let y0 = fy.floor(); + let (tx, ty) = (fx - x0, fy - y0); + let (x0, y0) = (x0 as isize, y0 as isize); + let a = img.get_clamped(x0, y0); + let b = img.get_clamped(x0 + 1, y0); + let c = img.get_clamped(x0, y0 + 1); + let d = img.get_clamped(x0 + 1, y0 + 1); + let top = a + (b - a) * tx; + let bot = c + (d - c) * tx; + top + (bot - top) * ty +} + +fn dist2(a: &[f32; DESC_LEN], b: &[f32; DESC_LEN]) -> f32 { + let mut s = 0.0; + for i in 0..DESC_LEN { + let d = a[i] - b[i]; + s += d * d; + } + s +} + +/// Nearest and second-nearest neighbour of each descriptor of `a` in `b`: (index, d1², d2²). +fn nearest(a: &Features, b: &Features) -> Vec<(usize, f32, f32)> { + a.desc + .par_iter() + .map(|da| { + let (mut bi, mut d1, mut d2) = (usize::MAX, f32::INFINITY, f32::INFINITY); + for (j, db) in b.desc.iter().enumerate() { + let d = dist2(da, db); + if d < d1 { + d2 = d1; + d1 = d; + bi = j; + } else if d < d2 { + d2 = d; + } + } + (bi, d1, d2) + }) + .collect() +} + +/// Putative matches `(index in a, index in b)`: ratio test (`ratio` on distances) and mutual +/// consistency. +pub fn match_features(a: &Features, b: &Features, ratio: f32) -> Vec<(usize, usize)> { + if a.is_empty() || b.is_empty() { + return vec![]; + } + let ab = nearest(a, b); + let ba = nearest(b, a); + let r2 = ratio * ratio; + ab.iter().enumerate().filter(|(i, (j, d1, d2))| *j != usize::MAX && *d1 < r2 * *d2 && ba[*j].0 == *i).map(|(i, (j, _, _))| (i, *j)).collect() +} + +#[cfg(test)] +mod tests { + use super::*; + + fn pattern(w: usize, h: usize, dx: f32, dy: f32) -> Plane { + Plane::from_fn(w, h, |x, y| { + let (u, v) = (x as f32 + 0.5 - dx, y as f32 + 0.5 - dy); + let a = ((u * 0.071).sin() * (v * 0.053).cos() + + (u * 0.023 + v * 0.031).sin() * 0.7 + + ((u * 0.011).floor() as i32 + (v * 0.013).floor() as i32).rem_euclid(2) as f32 * 0.6) + * 0.25 + + 0.5; + a.clamp(0.0, 1.0) + }) + } + + #[test] + fn shifted_images_match_with_the_shift() { + let a = pattern(400, 300, 0.0, 0.0); + let b = pattern(400, 300, 7.0, -4.0); + let fa = detect(&a, 300); + let fb = detect(&b, 300); + assert!(fa.len() > 50, "{} features", fa.len()); + let m = match_features(&fa, &fb, 0.8); + assert!(m.len() > 20, "{} matches", m.len()); + let good = m + .iter() + .filter(|(i, j)| { + let (p, q) = (fa.points[*i], fb.points[*j]); + ((q.x - p.x - 7.0).powi(2) + (q.y - p.y + 4.0).powi(2)).sqrt() < 1.5 + }) + .count(); + assert!(good as f32 > m.len() as f32 * 0.8, "{good}/{}", m.len()); + } +} diff --git a/crates/merge/src/frame.rs b/crates/merge/src/frame.rs new file mode 100644 index 00000000..df016b5d --- /dev/null +++ b/crates/merge/src/frame.rs @@ -0,0 +1,200 @@ +//! Merge inputs: scene-linear frames decoded from raw files (camera RGB, before white balance) or +//! standard images (linear light in the file's primaries), with their exposure from EXIF. + +use lightcraft_color::Mat3; +use lightcraft_geom::Orientation; +use lightcraft_meta::Metadata; +use lightcraft_raster::Rgb32f; +use lightcraft_raster::resample::{Filter, fit}; +use lightcraft_raw::ColorData; + +use crate::{MergeError, Result}; + +/// What the RGB values of a frame mean. +#[derive(Clone, Debug, PartialEq)] +pub enum FrameColor { + /// Camera RGB of a raw file (not white balanced): the DNG colour tags of the source. + Camera(Box), + /// Linear RGB with these primaries: linear RGB → XYZ (D50-adapted, ICC PCS). + Linear { to_xyz_d50: Mat3 }, +} + +/// One input photo. +#[derive(Clone, Debug)] +pub struct Frame { + /// Scene-linear pixels; 1.0 = the sensor/file white level. + pub image: Rgb32f, + /// Values at or above this (any channel) are saturated. + pub clip: f32, + /// Capture brightness from EXIF, in EV: log2(t · ISO / N²) (higher = brighter frame); `None` + /// when the file doesn't say. + pub exposure: Option, + pub raw: bool, + pub color: FrameColor, + /// How `image` must be rotated/flipped for display (`Normal` once applied). + pub orientation: Orientation, + pub metadata: Metadata, + /// Raw files: the source's `BaselineExposure` (EV), carried into the output. + pub baseline_exposure: f64, +} + +impl Frame { + /// A frame from linear pixels in the develop working space (linear Rec.2020), e.g. a rendered + /// procedural scene. `exposure` as in [`Frame::exposure`]. + pub fn from_linear_rec2020(image: Rgb32f, exposure: Option) -> Frame { + let to_xyz_d50 = lightcraft_color::bradford(lightcraft_color::D65, lightcraft_color::D50).mul(&lightcraft_color::REC2020.to_xyz()); + let clip = detect_clip(&image); + Frame { + image, + clip, + exposure, + raw: false, + color: FrameColor::Linear { to_xyz_d50 }, + orientation: Orientation::Normal, + metadata: Metadata::default(), + baseline_exposure: 0.0, + } + } + + pub fn width(&self) -> usize { + self.image.width + } + pub fn height(&self) -> usize { + self.image.height + } + + /// Apply the display orientation to the pixels. + pub fn orient(&mut self) { + if self.orientation != Orientation::Normal { + self.image = self.image.oriented(self.orientation); + self.orientation = Orientation::Normal; + } + } +} + +/// EV of a capture from its EXIF: log2(t · ISO / N²). ISO defaults to 100 and N to 1 when absent +/// (only ratios between frames matter); `None` without an exposure time. +pub fn exif_exposure(m: &Metadata) -> Option { + let t = m.exposure_time.filter(|t| *t > 0.0 && t.is_finite())?; + let iso = m.iso.filter(|i| *i > 0).unwrap_or(100) as f64; + let n = m.f_number.filter(|n| *n > 0.0 && n.is_finite()).unwrap_or(1.0); + Some((t * iso / (n * n)).log2()) +} + +/// The saturation level of a frame: 1.0, or the (lower) level at which values pile up when the +/// file's white level is optimistic. Returned slightly below the level. +pub fn detect_clip(img: &Rgb32f) -> f32 { + let maxv = img.data.iter().map(|p| p[0].max(p[1]).max(p[2])).fold(0.0f32, f32::max); + if maxv.is_nan() || maxv <= 0.0 { + return 1.0; + } + let level = maxv.min(1.0); + // a pile-up: the top 0.5 % of the range is ≥ 4× denser than the 2 % just below it + let (mut near, mut below) = (0usize, 0usize); + for p in &img.data { + let m = p[0].max(p[1]).max(p[2]); + if m >= level * 0.995 { + near += 1; + } else if m >= level * 0.97 && m < level * 0.99 { + below += 1; + } + } + let piled = near as f64 > img.data.len() as f64 * 1e-4 && near >= 4 && near as f64 / 0.005 > 4.0 * below as f64 / 0.02; + if maxv >= 1.0 || piled { level * 0.995 } else { 1.0 } +} + +/// Decode a photo (raw or standard format) into a frame no larger than `max_edge` (when given). +/// `orient` applies the EXIF orientation to the pixels. +pub fn load_frame(bytes: &[u8], max_edge: Option, orient: bool) -> Result { + let metadata = lightcraft_meta::extract(bytes); + let exposure = exif_exposure(&metadata); + if lightcraft_raw::probe(bytes).is_some() { + let mut raw = lightcraft_raw::decode(bytes).map_err(|e| MergeError::Decode(e.to_string()))?; + // geometric lens corrections stay with the develop settings ("profile corrections") + raw.opcodes + .list3 + .retain(|op| !matches!(op, lightcraft_raw::Opcode::WarpRectilinear { .. } | lightcraft_raw::Opcode::FixVignetteRadial { .. })); + let small = max_edge.is_some_and(|m| m <= 1200); + let method = if small { lightcraft_raw::Method::Bilinear } else { lightcraft_raw::Method::Ahd }; + let mut image = raw.develop(method).map_err(|e| MergeError::Decode(e.to_string()))?; + for p in image.data.iter_mut() { + *p = p.map(|v| v.max(0.0)); + } + let clip = detect_clip(&image); + if let Some(m) = max_edge + && image.width.max(image.height) > m + { + image = fit(&image, m, m, Filter::Box); + } + let mut f = Frame { + image, + clip, + exposure, + raw: true, + color: FrameColor::Camera(Box::new(raw.color.clone())), + orientation: raw.orientation, + metadata: if raw.metadata == Metadata::default() { metadata } else { raw.metadata.clone() }, + baseline_exposure: raw.color.baseline_exposure, + }; + if f.exposure.is_none() { + f.exposure = exif_exposure(&f.metadata); + } + if orient { + f.orient(); + } + return Ok(f); + } + let opts = match max_edge { + Some(m) => lightcraft_codecs::DecodeOptions::fit(m as u32, m as u32), + None => lightcraft_codecs::DecodeOptions::default(), + }; + let d = lightcraft_codecs::decode(bytes, opts).map_err(|e| MergeError::Decode(e.to_string()))?; + let mut image = d.image; + if let Some(m) = max_edge + && image.width.max(image.height) > m + { + image = fit(&image, m, m, Filter::Box); + } + let clip = detect_clip(&image); + let mut f = Frame { + image, + clip, + exposure, + raw: false, + color: FrameColor::Linear { to_xyz_d50: d.space.to_xyz_d50 }, + orientation: Orientation::from_exif(d.orientation), + metadata, + baseline_exposure: 0.0, + }; + if orient { + f.orient(); + } + Ok(f) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn exposure_from_exif() { + let m = Metadata { exposure_time: Some(1.0 / 125.0), f_number: Some(8.0), iso: Some(100), ..Default::default() }; + let e = exif_exposure(&m).unwrap(); + let m2 = Metadata { exposure_time: Some(4.0 / 125.0), ..m.clone() }; + assert!((exif_exposure(&m2).unwrap() - e - 2.0).abs() < 1e-9); + assert!(exif_exposure(&Metadata::default()).is_none()); + } + + #[test] + fn clip_detection() { + let mut img = Rgb32f::from_fn(100, 100, |x, _| [x as f32 / 100.0; 3]); + assert_eq!(detect_clip(&img), 1.0); + for p in img.data.iter_mut().take(500) { + *p = [0.9, 0.9, 0.9]; + } + for p in img.data.iter_mut() { + *p = p.map(|v| v.min(0.9)); + } + assert!((detect_clip(&img) - 0.9 * 0.995).abs() < 1e-6); + } +} diff --git a/crates/merge/src/hdr.rs b/crates/merge/src/hdr.rs new file mode 100644 index 00000000..7c596c4a --- /dev/null +++ b/crates/merge/src/hdr.rs @@ -0,0 +1,449 @@ +//! HDR merge of exposure brackets. +//! +//! 1. **Order and reference.** Frames are ordered by brightness (EXIF exposure, or median +//! luminance without EXIF); the middle one is the reference (geometry, colour and the output's +//! default brightness come from it). +//! 2. **Auto align** ([`crate::align`]): each frame is registered to the reference and resampled +//! into its geometry. +//! 3. **Exposure normalisation.** The exposure ratio between brightness-adjacent frames is +//! measured robustly — the median of log ratios over pixels well exposed in both — and chained +//! to the reference; EXIF is the fallback when too few pixels overlap. Radiance +//! `Rᵢ = vᵢ / kᵢ` is expressed in reference units. +//! 4. **Weights.** `wᵢ = kᵢ · (1 − smoothstep(0.80·clip, 0.96·clip, max channel))`: longer +//! exposures get more weight (shot-noise optimal for a linear sensor), saturated pixels none +//! (all channels together, so a clipped channel never shifts colour). The darkest frame keeps a +//! tiny weight so fully saturated pixels still get a value. +//! 5. **Deghost** (None/Low/Medium/High): on a ≤ 1024 px version, a per-pixel consistency test +//! against the reference: the reference radiance (or, where the reference is clipped/too dark, +//! the nearest well-exposed frame's) predicts what each frame should have recorded; frames that +//! disagree by more than a threshold (in log units, with a noise floor) are masked there. Masks +//! are dilated and feathered (Gaussian), upsampled, and remove the frame's weight. Their union +//! is the deghost overlay. + +use lightcraft_raster::blur::gaussian; +use lightcraft_raster::{Plane, Rgb32f}; +use rayon::prelude::*; +use serde::{Deserialize, Serialize}; + +use crate::align::{self, Alignment}; +use crate::frame::{Frame, FrameColor}; +use crate::{MergeError, Progress, Result, check}; + +#[derive(Clone, Copy, Debug, Default, PartialEq, Eq, Serialize, Deserialize)] +#[serde(rename_all = "lowercase")] +pub enum Deghost { + #[default] + None, + Low, + Medium, + High, +} + +impl Deghost { + /// (log-ratio threshold, dilation/feather radius in deghost-map pixels). + fn params(self) -> Option<(f32, usize)> { + match self { + Deghost::None => None, + Deghost::Low => Some((1.0, 2)), + Deghost::Medium => Some((0.6, 3)), + Deghost::High => Some((0.35, 5)), + } + } +} + +#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)] +#[serde(default, rename_all = "camelCase")] +pub struct HdrOptions { + pub align: bool, + pub deghost: Deghost, +} + +impl Default for HdrOptions { + fn default() -> Self { + HdrOptions { align: true, deghost: Deghost::None } + } +} + +/// The merged result, in the reference frame's geometry and units. +#[derive(Clone, Debug)] +pub struct HdrResult { + /// Scene-linear radiance; 1.0 = the reference frame's white level. + pub radiance: Rgb32f, + /// Deghost overlay (0..1, the union of the masks) at the radiance's size; empty without deghosting. + pub ghost: Plane, + /// Index (into the inputs) of the reference frame. + pub reference: usize, + /// Exposure of each input relative to the reference, in EV (as used for the merge). + pub ev: Vec, + /// The same from EXIF, when every frame had it. + pub ev_exif: Option>, + pub alignments: Vec, + pub color: FrameColor, + pub raw: bool, + pub orientation: lightcraft_geom::Orientation, + pub metadata: lightcraft_meta::Metadata, + pub baseline_exposure: f64, +} + +#[inline] +fn smoothstep(a: f32, b: f32, x: f32) -> f32 { + let t = ((x - a) / (b - a)).clamp(0.0, 1.0); + t * t * (3.0 - 2.0 * t) +} + +#[inline] +fn maxc(p: [f32; 3]) -> f32 { + p[0].max(p[1]).max(p[2]) +} + +/// Merge exposure brackets. Frames must share their pixel size (and orientation). +pub fn merge_hdr(mut frames: Vec, opts: &HdrOptions, progress: &Progress) -> Result { + let n = frames.len(); + if n < 2 { + return Err(MergeError::TooFew(n, 2)); + } + let (w, h) = (frames[0].width(), frames[0].height()); + if frames.iter().any(|f| f.width() != w || f.height() != h) { + return Err(MergeError::Mismatch("HDR merge needs photos of the same size".into())); + } + if frames.iter().any(|f| f.orientation != frames[0].orientation) { + for f in frames.iter_mut() { + f.orient(); + } + let (w0, h0) = (frames[0].width(), frames[0].height()); + if frames.iter().any(|f| f.width() != w0 || f.height() != h0) { + return Err(MergeError::Mismatch("HDR merge needs photos of the same size".into())); + } + } + let (w, h) = (frames[0].width(), frames[0].height()); + check(progress, 0.02, "Analysing exposures")?; + + // 1. brightness order and reference + let exif: Option> = frames.iter().map(|f| f.exposure).collect(); + let medians: Vec = frames.par_iter().map(|f| align::median_luma(&f.image, f.clip)).collect(); + let brightness: Vec = match &exif { + Some(e) if e.iter().any(|v| (v - e[0]).abs() > 0.1) => e.clone(), + _ => medians.iter().map(|m| (*m as f64).log2()).collect(), + }; + let mut order: Vec = (0..n).collect(); + order.sort_by(|&a, &b| brightness[a].total_cmp(&brightness[b]).then(a.cmp(&b))); + let reference = order[n / 2]; + + // 2. alignment + let mut alignments = vec![Alignment::identity(); n]; + let mut valid: Vec>> = vec![None; n]; + if opts.align { + check(progress, 0.05, "Aligning")?; + let g0 = 0.18 / medians[reference].max(1e-6); + let rel: Vec = (0..n).map(|i| 2f64.powf(brightness[reference] - brightness[i]) as f32).collect(); + let rel: Vec = if exif.is_some() && brightness.iter().any(|b| (b - brightness[0]).abs() > 0.1) { + rel + } else { + (0..n).map(|i| medians[reference] / medians[i].max(1e-6)).collect() + }; + let feats: Vec<(crate::features::Features, f64)> = frames + .par_iter() + .enumerate() + .map(|(i, f)| { + let (p, s) = align::work_image(&f.image, g0 * rel[i], align::WORK_EDGE); + (crate::features::detect(&p, 1500), s) + }) + .collect(); + check(progress, 0.25, "Aligning")?; + let al: Vec<(usize, Alignment)> = (0..n) + .into_par_iter() + .filter(|&i| i != reference) + .map(|i| (i, align::align_pair(&feats[reference].0, &feats[i].0, feats[i].1, w, h, i as u64 + 1))) + .collect(); + for (i, a) in al { + alignments[i] = a; + } + check(progress, 0.35, "Resampling")?; + let warped: Vec<(usize, Rgb32f, Vec)> = (0..n) + .into_par_iter() + .filter(|&i| i != reference && alignments[i].model != "identity") + .map(|i| { + let (img, v) = align::warp(&frames[i].image, &alignments[i].h); + (i, img, v) + }) + .collect(); + for (i, img, v) in warped { + frames[i].image = img; + valid[i] = Some(v); + } + } + check(progress, 0.45, "Normalising exposures")?; + + // 3. exposure ratios between brightness neighbours, chained to the reference (ln units) + let mut step_ln = vec![0f64; n.saturating_sub(1)]; + for p in 0..n - 1 { + let (a, b) = (order[p], order[p + 1]); + let est = log_ratio(&frames[a], valid[a].as_deref(), &frames[b], valid[b].as_deref()); + step_ln[p] = match (est, &exif) { + (Some(r), _) => r, + (None, Some(e)) => (e[b] - e[a]) * std::f64::consts::LN_2, + (None, None) => ((medians[b] / medians[a].max(1e-6)) as f64).max(1e-6).ln(), + }; + } + let mut ln_k = vec![0f64; n]; + let rpos = n / 2; + for p in (0..rpos).rev() { + ln_k[order[p]] = ln_k[order[p + 1]] - step_ln[p]; + } + for p in rpos + 1..n { + ln_k[order[p]] = ln_k[order[p - 1]] + step_ln[p - 1]; + } + let k: Vec = ln_k.iter().map(|v| v.exp() as f32).collect(); + let ev: Vec = ln_k.iter().map(|v| v / std::f64::consts::LN_2).collect(); + let ev_exif = exif.as_ref().map(|e| e.iter().map(|v| v - e[reference]).collect()); + + // 4./5. deghost maps (small scale) + check(progress, 0.55, "Deghosting")?; + let ghost_small: Option<(Vec, usize)> = opts.deghost.params().map(|(thr, r)| deghost_maps(&frames, &valid, &k, &order, reference, thr, r)); + check(progress, 0.7, "Merging")?; + + let darkest = order[0]; + let clips: Vec = frames.iter().map(|f| f.clip).collect(); + let mut out = Rgb32f::new(w, h); + let fr = &frames; + let vd = &valid; + let gs = ghost_small.as_ref(); + out.data.par_chunks_mut(w).enumerate().for_each(|(y, row)| { + for x in 0..w { + let i_px = y * w + x; + let mut acc = [0f64; 3]; + let mut ws = 0f64; + for i in 0..n { + if vd[i].as_ref().is_some_and(|v| !v[i_px]) { + continue; + } + let v = fr[i].image.data[i_px]; + let m = maxc(v); + let mut wt = k[i] * (1.0 - smoothstep(0.80 * clips[i], 0.96 * clips[i], m)); + if i == darkest { + wt += 1e-6; + } + if i == reference { + wt += 1e-9; + } + if i != reference + && let Some((maps, fac)) = gs + { + let g = sample_plane(&maps[i], (x as f32 + 0.5) / *fac as f32, (y as f32 + 0.5) / *fac as f32); + wt *= 1.0 - g; + } + if wt <= 0.0 { + continue; + } + let inv = 1.0 / k[i]; + for c in 0..3 { + acc[c] += (wt * v[c] * inv) as f64; + } + ws += wt as f64; + } + if ws > 0.0 { + row[x] = acc.map(|a| (a / ws) as f32); + } else { + row[x] = fr[reference].image.data[i_px]; + } + } + }); + check(progress, 0.92, "Merging")?; + + let ghost = match &ghost_small { + Some((maps, fac)) => { + let fac = *fac as f32; + let mut g = Plane::new(w, h); + g.data.par_chunks_mut(w).enumerate().for_each(|(y, row)| { + for (x, o) in row.iter_mut().enumerate() { + let (sx, sy) = ((x as f32 + 0.5) / fac, (y as f32 + 0.5) / fac); + *o = (0..n).filter(|&i| i != reference).map(|i| sample_plane(&maps[i], sx, sy)).fold(0.0, f32::max); + } + }); + g + } + None => Plane::new(0, 0), + }; + let rf = &frames[reference]; + Ok(HdrResult { + radiance: out, + ghost, + reference, + ev, + ev_exif, + alignments, + color: rf.color.clone(), + raw: rf.raw, + orientation: rf.orientation, + metadata: rf.metadata.clone(), + baseline_exposure: rf.baseline_exposure, + }) +} + +/// Bilinear sample of a plane at continuous coordinates (pixel centres at +0.5). +pub(crate) fn sample_plane(p: &Plane, x: f32, y: f32) -> f32 { + if p.width == 0 || p.height == 0 { + return 0.0; + } + crate::features::bilinear(p, x, y) +} + +/// Median of ln(b / a) on the green channel over pixels well exposed in both (≥ 500 samples). +fn log_ratio(a: &Frame, va: Option<&[bool]>, b: &Frame, vb: Option<&[bool]>) -> Option { + let len = a.image.data.len(); + let step = (len / 400_000).max(1); + let (ca, cb) = (a.clip, b.clip); + let lo = |f: &Frame| if f.raw { 0.01 } else { 0.04 }; + let (la, lb) = (lo(a) * ca, lo(b) * cb); + let mut v: Vec = (0..len) + .into_par_iter() + .step_by(step) + .filter_map(|i| { + if va.is_some_and(|v| !v[i]) || vb.is_some_and(|v| !v[i]) { + return None; + } + let (pa, pb) = (a.image.data[i], b.image.data[i]); + if maxc(pa) >= 0.8 * ca || maxc(pb) >= 0.8 * cb || pa[1] <= la || pb[1] <= lb { + return None; + } + Some((pb[1] / pa[1]).ln()) + }) + .collect(); + if v.len() < 500 { + return None; + } + v.sort_by(|x, y| x.total_cmp(y)); + Some(v[v.len() / 2] as f64) +} + +/// Box-downsample by `fac` (with the fraction of valid source pixels). +fn shrink(img: &Rgb32f, valid: Option<&[bool]>, fac: usize) -> (Rgb32f, Vec) { + let (w, h) = (img.width.div_ceil(fac), img.height.div_ceil(fac)); + let mut out = Rgb32f::new(w, h); + let mut vf = vec![0f32; w * h]; + out.data.par_chunks_mut(w).zip(vf.par_chunks_mut(w)).enumerate().for_each(|(y, (row, vrow))| { + for x in 0..w { + let mut s = [0f32; 3]; + let (mut nv, mut nt) = (0usize, 0usize); + for sy in y * fac..((y + 1) * fac).min(img.height) { + for sx in x * fac..((x + 1) * fac).min(img.width) { + let i = sy * img.width + sx; + nt += 1; + if valid.is_none_or(|v| v[i]) { + let p = img.data[i]; + for c in 0..3 { + s[c] += p[c]; + } + nv += 1; + } + } + } + if nv > 0 { + row[x] = s.map(|v| v / nv as f32); + } + vrow[x] = nv as f32 / nt.max(1) as f32; + } + }); + (out, vf) +} + +/// Per-frame ghost masks (0..1) at 1/`fac` scale; the reference's is all zero. +fn deghost_maps( + frames: &[Frame], + valid: &[Option>], + k: &[f32], + order: &[usize], + reference: usize, + thr: f32, + r: usize, +) -> (Vec, usize) { + let n = frames.len(); + let (w, h) = (frames[0].width(), frames[0].height()); + let fac = w.max(h).div_ceil(1024).max(1); + let small: Vec<(Rgb32f, Vec)> = (0..n).into_par_iter().map(|i| shrink(&frames[i].image, valid[i].as_deref(), fac)).collect(); + let (sw, sh) = (small[0].0.width, small[0].0.height); + let rpos = order.iter().position(|&i| i == reference).unwrap_or(0); + let good = |i: usize, p: usize| -> bool { + let f = &frames[i]; + let m = maxc(small[i].0.data[p]); + let dark = if f.raw { 0.01 } else { 0.04 }; + small[i].1[p] > 0.99 && m < 0.8 * f.clip && m > dark * f.clip + }; + // predicted reference radiance per pixel + let pred: Vec<[f32; 3]> = (0..sw * sh) + .into_par_iter() + .map(|p| { + let rad = |i: usize| small[i].0.data[p].map(|v| v / k[i]); + if good(reference, p) { + return rad(reference); + } + let bright = maxc(small[reference].0.data[p]) >= 0.8 * frames[reference].clip; + if bright { + for q in (0..rpos).rev() { + if good(order[q], p) { + return rad(order[q]); + } + } + rad(order[0]) + } else { + for &i in &order[rpos + 1..] { + if good(i, p) { + return rad(i); + } + } + rad(reference) + } + }) + .collect(); + let maps: Vec = (0..n) + .into_par_iter() + .map(|i| { + if i == reference { + return Plane::new(sw, sh); + } + let f = &frames[i]; + let e = if f.raw { 0.01 } else { 0.03 } * f.clip; + let mut m = Plane::new(sw, sh); + for p in 0..sw * sh { + if small[i].1[p] < 0.5 { + continue; + } + let obs = small[i].0.data[p]; + let mut d = 0f32; + for c in 0..3 { + let o = obs[c].min(f.clip); + let pr = (pred[p][c] * k[i]).min(f.clip); + d = d.max(((o + e) / (pr + e)).ln().abs()); + } + m.data[p] = smoothstep(thr * 0.7, thr, d); + } + // dilate (max filter, square r) then feather + let dil = dilate(&m, r); + let mut g = gaussian(&dil, r as f32); + for v in g.data.iter_mut() { + *v = v.clamp(0.0, 1.0); + } + g + }) + .collect(); + (maps, fac) +} + +fn dilate(p: &Plane, r: usize) -> Plane { + let (w, h) = (p.width, p.height); + let mut tmp = Plane::new(w, h); + for y in 0..h { + for x in 0..w { + let (a, b) = (x.saturating_sub(r), (x + r).min(w - 1)); + tmp.data[y * w + x] = p.data[y * w + a..=y * w + b].iter().copied().fold(0.0, f32::max); + } + } + let mut out = Plane::new(w, h); + for y in 0..h { + let (a, b) = (y.saturating_sub(r), (y + r).min(h - 1)); + for x in 0..w { + out.data[y * w + x] = (a..=b).map(|yy| tmp.data[yy * w + x]).fold(0.0, f32::max); + } + } + out +} diff --git a/crates/merge/src/lib.rs b/crates/merge/src/lib.rs new file mode 100644 index 00000000..93fafa6f --- /dev/null +++ b/crates/merge/src/lib.rs @@ -0,0 +1,57 @@ +//! LightCraft Photo Merge: HDR merge and panorama stitching (layer L3, no UI dependencies). +//! +//! - [`frame`]: inputs — scene-linear frames from raw files (camera RGB) or standard images. +//! - [`hdr`]: exposure brackets → one scene-linear radiance image (auto-align, exposure +//! normalisation, saturation-aware weighting, deghosting with an overlay). +//! - [`pano`]: overlapping photos → one panorama (features, matching, RANSAC, rotation bundle +//! adjustment, spherical/cylindrical/perspective projection, exposure compensation, seams, +//! multi-band blending, boundary warp, auto crop, fill edges). +//! - [`output`]: results as linear DNGs (16-bit float for HDR) with the source's colour tags. +//! +//! Every algorithm here is our own implementation from published papers (cited in each module); +//! no code from GPL projects was consulted. +#![forbid(unsafe_code)] + +pub mod align; +pub mod features; +pub mod frame; +pub mod hdr; +pub mod linalg; +pub mod output; +pub mod pano; +pub mod ransac; + +pub use frame::{Frame, FrameColor, load_frame}; +pub use hdr::{Deghost, HdrOptions, HdrResult, merge_hdr}; +pub use output::{DngSamples, write_linear_dng}; +pub use pano::{PanoOptions, PanoResult, Projection, stitch}; + +#[derive(Debug, Clone, PartialEq, thiserror::Error)] +pub enum MergeError { + #[error("{0} photo(s) selected; at least {1} are needed")] + TooFew(usize, usize), + #[error("{0}")] + Mismatch(String), + #[error("could not read a photo: {0}")] + Decode(String), + #[error("{0}")] + NoOverlap(String), + #[error("cancelled")] + Cancelled, + #[error("could not write the result: {0}")] + Output(String), +} + +pub type Result = std::result::Result; + +/// Progress callback: `(fraction 0..1, stage)`; return `false` to cancel. +pub type Progress<'a> = dyn Fn(f32, &str) -> bool + Sync + 'a; + +/// A progress callback that ignores progress and never cancels. +pub fn no_progress(_: f32, _: &str) -> bool { + true +} + +pub(crate) fn check(p: &Progress, f: f32, stage: &str) -> Result<()> { + if p(f.clamp(0.0, 1.0), stage) { Ok(()) } else { Err(MergeError::Cancelled) } +} diff --git a/crates/merge/src/linalg.rs b/crates/merge/src/linalg.rs new file mode 100644 index 00000000..5c8edf7c --- /dev/null +++ b/crates/merge/src/linalg.rs @@ -0,0 +1,247 @@ +//! Small dense linear algebra: Gaussian elimination, symmetric eigen-decomposition (cyclic Jacobi), +//! 3×3 rotations (Rodrigues). Textbook methods (Golub & Van Loan, "Matrix Computations", §8.5 for +//! Jacobi); sizes here are tiny (≤ a few hundred unknowns), so clarity beats speed. + +/// 3×3 matrix, row-major. +pub type M3 = [[f64; 3]; 3]; + +pub const I3: M3 = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]; + +pub fn mul3(a: &M3, b: &M3) -> M3 { + let mut r = [[0.0; 3]; 3]; + for i in 0..3 { + for j in 0..3 { + r[i][j] = (0..3).map(|k| a[i][k] * b[k][j]).sum(); + } + } + r +} + +pub fn apply3(a: &M3, v: [f64; 3]) -> [f64; 3] { + [0, 1, 2].map(|i| a[i][0] * v[0] + a[i][1] * v[1] + a[i][2] * v[2]) +} + +pub fn transpose3(a: &M3) -> M3 { + [0, 1, 2].map(|i| [a[0][i], a[1][i], a[2][i]]) +} + +pub fn det3(a: &M3) -> f64 { + a[0][0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1]) - a[0][1] * (a[1][0] * a[2][2] - a[1][2] * a[2][0]) + + a[0][2] * (a[1][0] * a[2][1] - a[1][1] * a[2][0]) +} + +pub fn inv3(a: &M3) -> Option { + let d = det3(a); + if d.abs() < 1e-300 || !d.is_finite() { + return None; + } + let c = |r0: usize, c0: usize, r1: usize, c1: usize| a[r0][c0] * a[r1][c1] - a[r0][c1] * a[r1][c0]; + Some([ + [c(1, 1, 2, 2) / d, -c(0, 1, 2, 2) / d, c(0, 1, 1, 2) / d], + [-c(1, 0, 2, 2) / d, c(0, 0, 2, 2) / d, -c(0, 0, 1, 2) / d], + [c(1, 0, 2, 1) / d, -c(0, 0, 2, 1) / d, c(0, 0, 1, 1) / d], + ]) +} + +pub fn dot(a: [f64; 3], b: [f64; 3]) -> f64 { + a[0] * b[0] + a[1] * b[1] + a[2] * b[2] +} + +pub fn cross(a: [f64; 3], b: [f64; 3]) -> [f64; 3] { + [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]] +} + +pub fn norm(a: [f64; 3]) -> f64 { + dot(a, a).sqrt() +} + +pub fn normalize(a: [f64; 3]) -> [f64; 3] { + let n = norm(a); + if n > 0.0 { a.map(|v| v / n) } else { a } +} + +/// Rotation matrix from a rotation vector (axis × angle), Rodrigues' formula. +pub fn rodrigues(w: [f64; 3]) -> M3 { + let th = norm(w); + if th < 1e-12 { + // first order: I + [w]x + return [[1.0, -w[2], w[1]], [w[2], 1.0, -w[0]], [-w[1], w[0], 1.0]]; + } + let k = w.map(|v| v / th); + let (s, c) = th.sin_cos(); + let v = 1.0 - c; + [ + [c + k[0] * k[0] * v, k[0] * k[1] * v - k[2] * s, k[0] * k[2] * v + k[1] * s], + [k[1] * k[0] * v + k[2] * s, c + k[1] * k[1] * v, k[1] * k[2] * v - k[0] * s], + [k[2] * k[0] * v - k[1] * s, k[2] * k[1] * v + k[0] * s, c + k[2] * k[2] * v], + ] +} + +/// Rotation vector of a rotation matrix (inverse of [`rodrigues`]). +pub fn rotation_vector(r: &M3) -> [f64; 3] { + let tr = (r[0][0] + r[1][1] + r[2][2] - 1.0) / 2.0; + let th = tr.clamp(-1.0, 1.0).acos(); + let a = [r[2][1] - r[1][2], r[0][2] - r[2][0], r[1][0] - r[0][1]]; + if th < 1e-9 { + return a.map(|v| v / 2.0); + } + if (std::f64::consts::PI - th).abs() < 1e-6 { + // axis from the largest diagonal element of (R + I) / 2 + let b = [0, 1, 2].map(|i| [0, 1, 2].map(|j| (r[i][j] + if i == j { 1.0 } else { 0.0 }) / 2.0)); + let i = (0..3).max_by(|&x, &y| b[x][x].total_cmp(&b[y][y])).unwrap_or(0); + let axis = normalize(b[i]); + return axis.map(|v| v * th); + } + let s = th / (2.0 * th.sin()); + a.map(|v| v * s) +} + +/// Nearest rotation to `m` (polar decomposition via a few Newton iterations X ← (X + X⁻ᵀ)/2). +pub fn orthonormalize(m: &M3) -> M3 { + let mut x = *m; + for _ in 0..20 { + let Some(inv) = inv3(&x) else { return I3 }; + let it = transpose3(&inv); + let mut n = [[0.0; 3]; 3]; + for i in 0..3 { + for j in 0..3 { + n[i][j] = 0.5 * (x[i][j] + it[i][j]); + } + } + x = n; + } + if det3(&x) < 0.0 { + // a reflection: flip the least significant axis + for row in x.iter_mut() { + row[2] = -row[2]; + } + } + x +} + +/// Solve `a · x = b` (`a` is `n × n`, row-major) by Gaussian elimination with partial pivoting. +pub fn solve(mut a: Vec, mut b: Vec, n: usize) -> Option> { + for col in 0..n { + let piv = (col..n).max_by(|&i, &j| a[i * n + col].abs().total_cmp(&a[j * n + col].abs()))?; + let pv = a[piv * n + col]; + if pv.abs() < 1e-300 || !pv.is_finite() { + return None; + } + if piv != col { + for k in 0..n { + a.swap(piv * n + k, col * n + k); + } + b.swap(piv, col); + } + for r in col + 1..n { + let f = a[r * n + col] / a[col * n + col]; + if f != 0.0 { + for k in col..n { + a[r * n + k] -= f * a[col * n + k]; + } + b[r] -= f * b[col]; + } + } + } + let mut x = vec![0.0; n]; + for r in (0..n).rev() { + let s: f64 = (r + 1..n).map(|k| a[r * n + k] * x[k]).sum(); + x[r] = (b[r] - s) / a[r * n + r]; + } + x.iter().all(|v| v.is_finite()).then_some(x) +} + +/// Eigen-decomposition of a symmetric `n × n` matrix (cyclic Jacobi). Returns eigenvalues ascending +/// and the matching unit eigenvectors. +pub fn sym_eigen(a: &[f64], n: usize) -> (Vec, Vec>) { + let mut m = a.to_vec(); + let mut v = vec![0.0; n * n]; + for i in 0..n { + v[i * n + i] = 1.0; + } + for _sweep in 0..100 { + let off: f64 = (0..n).flat_map(|i| (0..n).filter(move |&j| j != i).map(move |j| (i, j))).map(|(i, j)| m[i * n + j] * m[i * n + j]).sum(); + let scale: f64 = (0..n).map(|i| m[i * n + i] * m[i * n + i]).sum::().max(1e-300); + if off <= 1e-24 * scale { + break; + } + for p in 0..n { + for q in p + 1..n { + let apq = m[p * n + q]; + if apq.abs() < 1e-300 { + continue; + } + let theta = (m[q * n + q] - m[p * n + p]) / (2.0 * apq); + let t = theta.signum() / (theta.abs() + (theta * theta + 1.0).sqrt()); + let t = if theta == 0.0 { 1.0 } else { t }; + let c = 1.0 / (t * t + 1.0).sqrt(); + let s = t * c; + for k in 0..n { + let (mkp, mkq) = (m[k * n + p], m[k * n + q]); + m[k * n + p] = c * mkp - s * mkq; + m[k * n + q] = s * mkp + c * mkq; + } + for k in 0..n { + let (mpk, mqk) = (m[p * n + k], m[q * n + k]); + m[p * n + k] = c * mpk - s * mqk; + m[q * n + k] = s * mpk + c * mqk; + } + for k in 0..n { + let (vkp, vkq) = (v[k * n + p], v[k * n + q]); + v[k * n + p] = c * vkp - s * vkq; + v[k * n + q] = s * vkp + c * vkq; + } + } + } + } + let mut idx: Vec = (0..n).collect(); + idx.sort_by(|&a, &b| m[a * n + a].total_cmp(&m[b * n + b])); + let vals = idx.iter().map(|&i| m[i * n + i]).collect(); + let vecs = idx.iter().map(|&i| (0..n).map(|k| v[k * n + i]).collect()).collect(); + (vals, vecs) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn rodrigues_round_trip() { + for w in [[0.1, -0.2, 0.3], [0.0, 0.0, 0.0], [1.5, 0.2, -0.7], [0.0, 3.0, 0.0]] { + let r = rodrigues(w); + assert!((det3(&r) - 1.0).abs() < 1e-9); + let back = rotation_vector(&r); + let r2 = rodrigues(back); + for i in 0..3 { + for j in 0..3 { + assert!((r[i][j] - r2[i][j]).abs() < 1e-9, "{w:?}"); + } + } + } + } + + #[test] + fn solve_and_eigen() { + let a = vec![4.0, 1.0, 0.0, 1.0, 3.0, 1.0, 0.0, 1.0, 2.0]; + let x = solve(a.clone(), vec![1.0, 2.0, 3.0], 3).unwrap(); + for r in 0..3 { + let s: f64 = (0..3).map(|k| a[r * 3 + k] * x[k]).sum(); + assert!((s - [1.0, 2.0, 3.0][r]).abs() < 1e-12); + } + let (vals, vecs) = sym_eigen(&a, 3); + assert!(vals[0] <= vals[1] && vals[1] <= vals[2]); + for (l, v) in vals.iter().zip(&vecs) { + for r in 0..3 { + let s: f64 = (0..3).map(|k| a[r * 3 + k] * v[k]).sum(); + assert!((s - l * v[r]).abs() < 1e-9); + } + } + let r = orthonormalize(&[[1.02, 0.01, 0.0], [-0.01, 0.98, 0.02], [0.0, -0.02, 1.01]]); + let rt = mul3(&r, &transpose3(&r)); + for i in 0..3 { + for j in 0..3 { + assert!((rt[i][j] - I3[i][j]).abs() < 1e-9); + } + } + } +} diff --git a/crates/merge/src/output.rs b/crates/merge/src/output.rs new file mode 100644 index 00000000..3d0bbfe5 --- /dev/null +++ b/crates/merge/src/output.rs @@ -0,0 +1,142 @@ +//! Merge results as linear DNGs (`PhotometricInterpretation = LinearRaw`, three samples per +//! pixel), written with LightCraft's own TIFF/DNG writer. +//! +//! - Pixel values are scaled so the brightest pixel sits just below 1.0 (the white level); the +//! scale is carried in `BaselineExposure`, so a default rendering looks like the reference frame. +//! - Raw sources: the reference frame's DNG colour tags (`ColorMatrix*`, `ForwardMatrix*`, +//! `CameraCalibration*`, `AsShotNeutral`, illuminants) are carried over unchanged — the data is +//! still camera RGB. +//! - Other sources: the "camera" is the file's linear RGB space; `ColorMatrix1` = XYZ(D65) → RGB, +//! `ForwardMatrix1` = RGB → XYZ(D50), `AsShotNeutral` = (1, 1, 1). + +use lightcraft_color::{D50, D65, Mat3, bradford}; +use lightcraft_geom::Orientation; +use lightcraft_meta::Metadata; +use lightcraft_raster::Rgb32f; +use lightcraft_raw::{BlackLevel, ColorData, DngCompression, DngWriteOptions, OpcodeLists, RawData, RawFormat, RawImage, Rect}; + +use crate::frame::FrameColor; +use crate::{MergeError, Result}; + +/// Sample encoding of the output DNG. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum DngSamples { + /// 16-bit floating point, Deflate + floating-point predictor (HDR). + Half, + /// 32-bit floating point, Deflate + floating-point predictor. + Float, + /// 16-bit integer, Deflate + horizontal predictor (panoramas of normal-range photos). + U16, +} + +/// The DNG colour tags for frames of `color`. +pub fn color_data(color: &FrameColor) -> ColorData { + match color { + FrameColor::Camera(c) => (**c).clone(), + FrameColor::Linear { to_xyz_d50 } => { + let to_xyz_d65 = bradford(D50, D65).mul(to_xyz_d50); + let cm = to_xyz_d65.inverse().unwrap_or(Mat3::IDENTITY); + ColorData { + illuminant: [21, 0], + color_matrix: [Some(cm), None], + forward_matrix: [Some(*to_xyz_d50), None], + camera_calibration: [None, None], + analog_balance: None, + as_shot_neutral: Some([1.0, 1.0, 1.0]), + as_shot_white_xy: None, + baseline_exposure: 0.0, + } + } + } +} + +/// Encode `img` (scene-linear, 1.0 = source white) as a LinearRaw DNG. +pub fn write_linear_dng( + img: &Rgb32f, + color: &FrameColor, + orientation: Orientation, + metadata: &Metadata, + baseline_exposure: f64, + samples: DngSamples, +) -> Result> { + let (w, h) = (img.width, img.height); + if w == 0 || h == 0 { + return Err(MergeError::Output("empty image".into())); + } + let peak = img.data.iter().flat_map(|p| p.iter()).copied().filter(|v| v.is_finite()).fold(0.0f32, f32::max); + let target = 0.95f32; + // never brighten: values ≤ 1 keep their scale (so a plain panorama keeps its white level) + let scale = if peak > target { target / peak } else { 1.0 }; + let mut cd = color_data(color); + cd.baseline_exposure = baseline_exposure + (1.0 / scale as f64).log2(); + let clean = |v: f32| if v.is_finite() { (v * scale).max(0.0) } else { 0.0 }; + let (data, bits, white) = match samples { + DngSamples::Half | DngSamples::Float => { + (RawData::F32(img.data.iter().flat_map(|p| p.map(clean)).collect()), if samples == DngSamples::Half { 16 } else { 32 }, 1.0) + } + DngSamples::U16 => { + (RawData::U16(img.data.iter().flat_map(|p| p.map(|v| (clean(v) * 65535.0).round().min(65535.0) as u16)).collect()), 16, 65535.0) + } + }; + let mut meta = metadata.clone(); + meta.width = Some(w as u32); + meta.height = Some(h as u32); + meta.orientation = Some(orientation); + let raw = RawImage { + format: RawFormat::Dng, + width: w, + height: h, + cpp: 3, + data, + cfa: None, + bits, + black: BlackLevel::uniform(0.0), + white: vec![white], + active_area: Rect::new(0, 0, w, h), + crop: Rect::new(0, 0, w, h), + orientation, + color: cd, + wb_multipliers: None, + linearized: false, + opcodes: OpcodeLists::default(), + metadata: meta, + }; + let compression = DngCompression::Deflate { tile: 256, half: samples == DngSamples::Half }; + lightcraft_raw::write_dng(&raw, &DngWriteOptions { compression, ..Default::default() }).map_err(|e| MergeError::Output(e.to_string())) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn linear_dng_round_trips_through_the_raw_decoder() { + let img = Rgb32f::from_fn(80, 60, |x, y| [x as f32 / 20.0, y as f32 / 30.0, 0.3]); + let color = FrameColor::Linear { to_xyz_d50: bradford(D65, D50).mul(&lightcraft_color::SRGB.to_xyz()) }; + let bytes = write_linear_dng(&img, &color, Orientation::Normal, &Metadata::default(), 0.0, DngSamples::Half).unwrap(); + let raw = lightcraft_raw::decode(&bytes).unwrap(); + assert_eq!((raw.width, raw.height, raw.cpp), (80, 60, 3)); + let gain = 2f32.powf(raw.color.baseline_exposure as f32); + let back = raw.develop(lightcraft_raw::Method::Bilinear).unwrap(); + for (a, b) in img.data.iter().zip(&back.data) { + for c in 0..3 { + assert!((a[c] - b[c] * gain).abs() <= a[c] * 2e-3 + 1e-5, "{a:?} {b:?}"); + } + } + // the colour model maps the source white to neutral (no white-balance shift) + let xy = lightcraft_raw::color::as_shot_white_xy(&raw); + assert!((xy.x - D65.x).abs() < 2e-3 && (xy.y - D65.y).abs() < 2e-3, "{xy:?}"); + let t = lightcraft_raw::color::camera_transform(&raw, xy); + let white = t.matrix.apply([t.wb[0] as f64, t.wb[1] as f64, t.wb[2] as f64]); + assert!(white.iter().all(|v| (v - white[0]).abs() < 1e-3), "{white:?}"); + // sRGB red maps to the Rec.2020 coordinates of sRGB red + let red = t.matrix.apply([t.wb[0] as f64, 0.0, 0.0]); + let want = lightcraft_color::SRGB.to_space(&lightcraft_color::REC2020).apply([1.0, 0.0, 0.0]); + let k = white[0]; + for c in 0..3 { + assert!((red[c] / k - want[c]).abs() < 5e-3, "{red:?} vs {want:?}"); + } + let u16 = write_linear_dng(&img, &color, Orientation::Normal, &Metadata::default(), 0.0, DngSamples::U16).unwrap(); + assert!(lightcraft_raw::decode(&u16).is_ok()); + } +} diff --git a/crates/merge/src/pano/blend.rs b/crates/merge/src/pano/blend.rs new file mode 100644 index 00000000..e3147da1 --- /dev/null +++ b/crates/merge/src/pano/blend.rs @@ -0,0 +1,300 @@ +//! Multi-band blending (P. J. Burt & E. H. Adelson, "A Multiresolution Spline With Application to +//! Image Mosaics", ACM ToG 1983), as used by Brown & Lowe: each image's Laplacian pyramid is +//! blended with the Gaussian pyramid of its (binary) seam mask, so low frequencies blend over wide +//! regions and high frequencies over narrow ones. Blending happens in log space, which turns +//! residual exposure differences into additive offsets the coarse bands smooth away. +//! +//! Images arrive as tiles placed on a canvas padded to a multiple of 2^levels; tile origins are +//! multiples of 2^levels so every pyramid level of a tile lines up with the canvas pyramid. + +use rayon::prelude::*; + +/// A planar-interleaved float buffer with `c` channels. +#[derive(Clone, Debug, PartialEq)] +pub struct Buf { + pub w: usize, + pub h: usize, + pub c: usize, + pub data: Vec, +} + +impl Buf { + pub fn new(w: usize, h: usize, c: usize) -> Buf { + Buf { w, h, c, data: vec![0.0; w * h * c] } + } + #[inline] + fn at(&self, x: isize, y: isize, ch: usize) -> f32 { + let x = x.clamp(0, self.w as isize - 1) as usize; + let y = y.clamp(0, self.h as isize - 1) as usize; + self.data[(y * self.w + x) * self.c + ch] + } +} + +const K: [f32; 5] = [1.0 / 16.0, 4.0 / 16.0, 6.0 / 16.0, 4.0 / 16.0, 1.0 / 16.0]; + +/// Binomial 5-tap blur then decimate by 2 (output `⌈w/2⌉ × ⌈h/2⌉`). +pub fn down(b: &Buf) -> Buf { + let (w2, h2) = (b.w.div_ceil(2), b.h.div_ceil(2)); + // horizontal pass at decimated x + let mut tmp = Buf::new(w2, b.h, b.c); + tmp.data.par_chunks_mut(w2 * b.c).enumerate().for_each(|(y, row)| { + for x in 0..w2 { + for ch in 0..b.c { + let mut s = 0.0; + for (k, wk) in K.iter().enumerate() { + s += wk * b.at(2 * x as isize + k as isize - 2, y as isize, ch); + } + row[x * b.c + ch] = s; + } + } + }); + let mut out = Buf::new(w2, h2, b.c); + out.data.par_chunks_mut(w2 * b.c).enumerate().for_each(|(y, row)| { + for x in 0..w2 { + for ch in 0..b.c { + let mut s = 0.0; + for (k, wk) in K.iter().enumerate() { + s += wk * tmp.at(x as isize, 2 * y as isize + k as isize - 2, ch); + } + row[x * b.c + ch] = s; + } + } + }); + out +} + +/// Expand to `w × h` (bilinear, pixel-centre aligned with [`down`]'s grid). +pub fn up(b: &Buf, w: usize, h: usize) -> Buf { + let mut out = Buf::new(w, h, b.c); + out.data.par_chunks_mut(w * b.c).enumerate().for_each(|(y, row)| { + let fy = (y as f32 + 0.5) / 2.0 - 0.5; + let y0 = fy.floor(); + let ty = fy - y0; + let y0 = y0 as isize; + for x in 0..w { + let fx = (x as f32 + 0.5) / 2.0 - 0.5; + let x0 = fx.floor(); + let tx = fx - x0; + let x0 = x0 as isize; + for ch in 0..b.c { + let a = b.at(x0, y0, ch); + let bb = b.at(x0 + 1, y0, ch); + let c = b.at(x0, y0 + 1, ch); + let d = b.at(x0 + 1, y0 + 1, ch); + let top = a + (bb - a) * tx; + let bot = c + (d - c) * tx; + row[x * b.c + ch] = top + (bot - top) * ty; + } + } + }); + out +} + +/// Fill pixels with weight 0 from their neighbourhood (push–pull: weighted pyramid down, fill +/// from coarser levels up). `vals` has `c` channels; `wgt` one. +pub fn push_pull(vals: &mut Buf, wgt: &[f32]) { + if wgt.iter().all(|&w| w > 0.0) || !wgt.iter().any(|&w| w > 0.0) { + return; + } + let (w, h, c) = (vals.w, vals.h, vals.c); + // premultiplied values + weight as the last channel + let mut pm = Buf::new(w, h, c + 1); + for i in 0..w * h { + let a = wgt[i].clamp(0.0, 1.0); + for ch in 0..c { + pm.data[i * (c + 1) + ch] = vals.data[i * c + ch] * a; + } + pm.data[i * (c + 1) + c] = a; + } + let mut levels = vec![pm]; + while levels.last().is_some_and(|l| l.w > 1 || l.h > 1) { + let d = down(levels.last().expect("level")); + levels.push(d); + } + // normalise the coarsest, then fill each finer level's holes from the coarser one + let mut filled = { + let l = levels.last().expect("level"); + let mut f = Buf::new(l.w, l.h, c); + for i in 0..l.w * l.h { + let a = l.data[i * (c + 1) + c]; + for ch in 0..c { + f.data[i * c + ch] = if a > 1e-12 { l.data[i * (c + 1) + ch] / a } else { 0.0 }; + } + } + f + }; + for li in (0..levels.len() - 1).rev() { + let l = &levels[li]; + let coarse = up(&filled, l.w, l.h); + let mut f = Buf::new(l.w, l.h, c); + for i in 0..l.w * l.h { + let a = l.data[i * (c + 1) + c].min(1.0); + for ch in 0..c { + let v = if a > 1e-12 { l.data[i * (c + 1) + ch] / l.data[i * (c + 1) + c] } else { 0.0 }; + f.data[i * c + ch] = v * a + coarse.data[i * c + ch] * (1.0 - a); + } + } + filled = f; + } + for i in 0..w * h { + if wgt[i] <= 0.0 { + for ch in 0..c { + vals.data[i * c + ch] = filled.data[i * c + ch]; + } + } + } +} + +/// Accumulates tiles into the canvas Laplacian pyramid. +pub struct Blender { + pub levels: usize, + pub w: usize, + pub h: usize, + acc: Vec, + wsum: Vec, +} + +impl Blender { + /// A blender for a `w × h` canvas (both multiples of 2^levels). + pub fn new(w: usize, h: usize, levels: usize) -> Blender { + let mut acc = Vec::new(); + let mut wsum = Vec::new(); + for l in 0..=levels { + acc.push(Buf::new(w >> l, h >> l, 3)); + wsum.push(Buf::new(w >> l, h >> l, 1)); + } + Blender { levels, w, h, acc, wsum } + } + + /// Add a tile at canvas offset (`x0`, `y0`) (multiples of 2^levels; tile size too): `vals` + /// (3 channels, holes already filled) and its seam mask (0..1). + pub fn add(&mut self, x0: usize, y0: usize, vals: Buf, mask: Buf) { + let mut g = vals; + let mut m = mask; + for l in 0..=self.levels { + let (ox, oy) = (x0 >> l, y0 >> l); + let (gn, mn) = if l < self.levels { (Some(down(&g)), Some(down(&m))) } else { (None, None) }; + let lap = match &gn { + Some(gn) => { + let u = up(gn, g.w, g.h); + let mut lap = g.clone(); + for (a, b) in lap.data.iter_mut().zip(&u.data) { + *a -= b; + } + lap + } + None => g.clone(), + }; + let (acc, ws) = (&mut self.acc[l], &mut self.wsum[l]); + let aw = acc.w; + let tw = lap.w; + acc.data.par_chunks_mut(aw * 3).zip(ws.data.par_chunks_mut(aw)).enumerate().for_each(|(y, (arow, wrow))| { + if y < oy || y >= oy + lap.h { + return; + } + let ty = y - oy; + for tx in 0..tw { + let x = ox + tx; + if x >= aw { + break; + } + let mk = m.data[ty * tw + tx]; + if mk <= 0.0 { + continue; + } + for ch in 0..3 { + arow[x * 3 + ch] += lap.data[(ty * tw + tx) * 3 + ch] * mk; + } + wrow[x] += mk; + } + }); + if let (Some(gn), Some(mn)) = (gn, mn) { + g = gn; + m = mn; + } + } + } + + /// Normalise and collapse the pyramid into the canvas. + pub fn finish(self) -> Buf { + let mut res: Option = None; + for l in (0..=self.levels).rev() { + let mut b = self.acc[l].clone(); + let ws = &self.wsum[l]; + for i in 0..b.w * b.h { + let s = ws.data[i]; + for ch in 0..3 { + b.data[i * 3 + ch] = if s > 1e-8 { b.data[i * 3 + ch] / s } else { 0.0 }; + } + } + if let Some(r) = res { + let u = up(&r, b.w, b.h); + for (a, v) in b.data.iter_mut().zip(&u.data) { + *a += v; + } + } + res = Some(b); + } + res.expect("level 0") + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn single_tile_reconstructs_exactly() { + let (w, h) = (64, 48); + let mut v = Buf::new(w, h, 3); + for (i, x) in v.data.iter_mut().enumerate() { + *x = ((i * 37) % 101) as f32 / 101.0; + } + let mut m = Buf::new(w, h, 1); + m.data.fill(1.0); + let mut b = Blender::new(w, h, 3); + b.add(0, 0, v.clone(), m); + let out = b.finish(); + for (a, b) in v.data.iter().zip(&out.data) { + assert!((a - b).abs() < 1e-4); + } + } + + #[test] + fn two_tiles_blend_without_a_step() { + // left tile value 0, right tile value 1, seam in the middle: the result is a smooth ramp + let (w, h) = (128, 32); + let mut b = Blender::new(w, h, 4); + for (val, side) in [(0.0f32, 0), (1.0, 1)] { + let mut v = Buf::new(w, h, 3); + v.data.fill(val); + let mut m = Buf::new(w, h, 1); + for y in 0..h { + for x in 0..w { + m.data[y * w + x] = if (x >= w / 2) == (side == 1) { 1.0 } else { 0.0 }; + } + } + b.add(0, 0, v, m); + } + let out = b.finish(); + let row: Vec = (0..w).map(|x| out.data[(16 * w + x) * 3]).collect(); + let max_step = row.windows(2).map(|p| (p[1] - p[0]).abs()).fold(0.0, f32::max); + assert!(max_step < 0.15, "max step {max_step}"); + assert!(row[2].abs() < 0.02 && (row[w - 3] - 1.0).abs() < 0.02); + } + + #[test] + fn push_pull_fills_holes() { + let (w, h) = (20, 10); + let mut v = Buf::new(w, h, 1); + let mut wg = vec![0.0; w * h]; + for y in 0..h { + for x in 0..w / 2 { + v.data[y * w + x] = 0.7; + wg[y * w + x] = 1.0; + } + } + push_pull(&mut v, &wg); + assert!(v.data.iter().all(|x| (x - 0.7).abs() < 1e-4)); + } +} diff --git a/crates/merge/src/pano/camera.rs b/crates/merge/src/pano/camera.rs new file mode 100644 index 00000000..1257b952 --- /dev/null +++ b/crates/merge/src/pano/camera.rs @@ -0,0 +1,446 @@ +//! Panorama registration: pairwise matching → camera rotations and focal lengths. +//! +//! After M. Brown & D. G. Lowe, "Automatic Panoramic Image Stitching using Invariant Features" +//! (IJCV 2007): +//! - pairwise feature matches verified with a RANSAC homography; a pair is accepted when the +//! inliers explain enough of the matches (`nᵢ > 8 + 0.3·n_f`); +//! - the largest connected set of images is kept; +//! - focal length initialised from the homographies (Szeliski & Shum 1997 relations, derived here +//! from the orthonormality of `K₁⁻¹ H K₀`), rotations by chaining the pairwise rotations along a +//! maximum spanning tree (most inliers first); +//! - bundle adjustment of every camera's rotation and focal length (Levenberg–Marquardt on the +//! symmetric reprojection error, Huber-weighted, with a numeric Jacobian per match); +//! - automatic straightening: the world "up" is the direction most perpendicular to every +//! camera's x axis. + +use lightcraft_geom::{Homography, Point}; +use rayon::prelude::*; + +use crate::features::{self, Features}; +use crate::linalg::{self, M3}; +use crate::ransac::{self, Model}; + +/// A pinhole camera: `u = f · (R d)ₓᵧ / (R d)_z + c` for a world direction `d`. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct Camera { + /// Focal length (pixels of the image it belongs to). + pub f: f64, + /// World → camera rotation. + pub r: M3, + pub cx: f64, + pub cy: f64, +} + +impl Camera { + /// World direction of image pixel `(x, y)`. + #[inline] + pub fn ray(&self, x: f64, y: f64) -> [f64; 3] { + let d = [(x - self.cx) / self.f, (y - self.cy) / self.f, 1.0]; + linalg::apply3(&linalg::transpose3(&self.r), d) + } + /// Image pixel of world direction `d` (`None` behind the camera). + #[inline] + pub fn project(&self, d: [f64; 3]) -> Option<(f64, f64)> { + let p = linalg::apply3(&self.r, d); + (p[2] > 1e-9).then(|| (self.f * p[0] / p[2] + self.cx, self.f * p[1] / p[2] + self.cy)) + } +} + +/// Verified matches between two images (full-resolution pixel coordinates). +#[derive(Clone, Debug)] +pub struct PairMatch { + pub i: usize, + pub j: usize, + /// Maps image `i` pixels to image `j` pixels. + pub h: Homography, + pub pts: Vec<(Point, Point)>, + /// Coarsest pyramid level of the two features of each match (finer = more precise). + pub levels: Vec, +} + +/// Match every pair of images. `feats[k]` are in work coordinates scaled by `scales[k]`. +pub fn match_pairs(feats: &[Features], scales: &[f64]) -> Vec { + let n = feats.len(); + let pairs: Vec<(usize, usize)> = (0..n).flat_map(|i| (i + 1..n).map(move |j| (i, j))).collect(); + let mut out: Vec = pairs + .par_iter() + .filter_map(|&(i, j)| { + let m = features::match_features(&feats[i], &feats[j], 0.8); + if m.len() < 12 { + return None; + } + let a: Vec = m.iter().map(|(p, _)| Point::new(feats[i].points[*p].x as f64, feats[i].points[*p].y as f64)).collect(); + let b: Vec = m.iter().map(|(_, q)| Point::new(feats[j].points[*q].x as f64, feats[j].points[*q].y as f64)).collect(); + let fit = ransac::ransac(Model::Homography, &a, &b, 3.0, 12, (i * 1000 + j) as u64)?; + // Brown & Lowe's probabilistic verification + if (fit.inliers.len() as f64) <= 8.0 + 0.3 * m.len() as f64 { + return None; + } + let (si, sj) = (scales[i], scales[j]); + let pts = fit.inliers.iter().map(|&k| (Point::new(a[k].x / si, a[k].y / si), Point::new(b[k].x / sj, b[k].y / sj))).collect(); + let levels = fit.inliers.iter().map(|&k| feats[i].points[m[k].0].level.max(feats[j].points[m[k].1].level)).collect(); + // conjugate to full resolution: H_full = S_j⁻¹ · H · S_i + let s_i = Homography([si, 0.0, 0.0, 0.0, si, 0.0, 0.0, 0.0, 1.0]); + let s_j_inv = Homography([1.0 / sj, 0.0, 0.0, 0.0, 1.0 / sj, 0.0, 0.0, 0.0, 1.0]); + Some(PairMatch { i, j, h: s_j_inv.mul(&fit.h).mul(&s_i), pts, levels }) + }) + .collect(); + out.sort_by_key(|p| (p.i, p.j)); + out +} + +/// The largest connected component of the match graph (sorted image indices). +pub fn largest_component(n: usize, pairs: &[PairMatch]) -> Vec { + let mut parent: Vec = (0..n).collect(); + fn find(p: &mut [usize], x: usize) -> usize { + let mut r = x; + while p[r] != r { + r = p[r]; + } + let mut c = x; + while p[c] != r { + let nx = p[c]; + p[c] = r; + c = nx; + } + r + } + for m in pairs { + let (a, b) = (find(&mut parent, m.i), find(&mut parent, m.j)); + if a != b { + parent[a] = b; + } + } + let mut groups: std::collections::BTreeMap> = Default::default(); + for i in 0..n { + let r = find(&mut parent, i); + groups.entry(r).or_default().push(i); + } + groups.into_values().max_by_key(|g| g.len()).unwrap_or_default() +} + +/// Focal length candidates (f₀, f₁) from a homography between centred image coordinates. +pub fn focals_from_homography(h: &Homography) -> (Option, Option) { + let m = h.0; + let (h00, h01, h02, h10, h11, h12, h20, h21) = (m[0], m[1], m[2], m[3], m[4], m[5], m[6], m[7]); + let pick = |n1: f64, d1: f64, n2: f64, d2: f64| -> Option { + let a = (d1.abs() > 1e-12).then(|| n1 / d1).filter(|v| *v > 0.0); + let b = (d2.abs() > 1e-12).then(|| n2 / d2).filter(|v| *v > 0.0); + let v = match (a, b) { + (Some(a), Some(b)) => { + if d1.abs() > d2.abs() { + a + } else { + b + } + } + (Some(a), None) => a, + (None, Some(b)) => b, + _ => return None, + }; + Some(v.sqrt()) + }; + // rows of K₁⁻¹HK₀ orthogonal / equal length → f₀; columns → f₁ + let f0 = pick(-h02 * h12, h00 * h10 + h01 * h11, h12 * h12 - h02 * h02, h00 * h00 + h01 * h01 - h10 * h10 - h11 * h11); + let f1 = pick(-(h00 * h01 + h10 * h11), h20 * h21, h01 * h01 + h11 * h11 - h00 * h00 - h10 * h10, h20 * h20 - h21 * h21); + (f0, f1) +} + +fn centred(h: &Homography, ci: (f64, f64), cj: (f64, f64)) -> Homography { + let ti = Homography([1.0, 0.0, ci.0, 0.0, 1.0, ci.1, 0.0, 0.0, 1.0]); + let tj_inv = Homography([1.0, 0.0, -cj.0, 0.0, 1.0, -cj.1, 0.0, 0.0, 1.0]); + tj_inv.mul(h).mul(&ti) +} + +/// Initial cameras for `ids` (a connected set): shared focal `f0`, rotations chained along a +/// maximum spanning tree rooted at the best-connected image. +pub fn initial_cameras(ids: &[usize], sizes: &[(usize, usize)], pairs: &[PairMatch], f0: f64) -> Vec> { + let n = sizes.len(); + let mut cams: Vec> = vec![None; n]; + let deg = |k: usize| pairs.iter().filter(|p| p.i == k || p.j == k).map(|p| p.pts.len()).sum::(); + let Some(&root) = ids.iter().max_by_key(|&&k| (deg(k), usize::MAX - k)) else { return cams }; + let c = |k: usize| (sizes[k].0 as f64 / 2.0, sizes[k].1 as f64 / 2.0); + cams[root] = Some(Camera { f: f0, r: linalg::I3, cx: c(root).0, cy: c(root).1 }); + let mut edges: Vec<&PairMatch> = pairs.iter().filter(|p| ids.contains(&p.i) && ids.contains(&p.j)).collect(); + edges.sort_by_key(|p| std::cmp::Reverse(p.pts.len())); + loop { + let mut grew = false; + for e in &edges { + let (known, new, h) = match (cams[e.i].is_some(), cams[e.j].is_some()) { + (true, false) => (e.i, e.j, e.h), + (false, true) => ( + e.j, + e.i, + match e.h.inverse() { + Some(h) => h, + None => continue, + }, + ), + _ => continue, + }; + // h maps `known` → `new`; R_new R_knownᵀ ≈ K_new⁻¹ H K_known + let hc = centred(&h, c(known), c(new)); + let k = [[f0, 0.0, 0.0], [0.0, f0, 0.0], [0.0, 0.0, 1.0]]; + let kinv = [[1.0 / f0, 0.0, 0.0], [0.0, 1.0 / f0, 0.0], [0.0, 0.0, 1.0]]; + let hm = [[hc.0[0], hc.0[1], hc.0[2]], [hc.0[3], hc.0[4], hc.0[5]], [hc.0[6], hc.0[7], hc.0[8]]]; + let rel = linalg::mul3(&linalg::mul3(&kinv, &hm), &k); + let d = linalg::det3(&rel); + if !d.is_finite() || d.abs() < 1e-12 { + continue; + } + let s = d.signum() * d.abs().cbrt(); + let rel = rel.map(|r| r.map(|v| v / s)); + let rel = linalg::orthonormalize(&rel); + let rk = cams[known].expect("known").r; + cams[new] = Some(Camera { f: f0, r: linalg::mul3(&rel, &rk), cx: c(new).0, cy: c(new).1 }); + grew = true; + } + if !grew { + break; + } + } + cams +} + +/// Bundle adjustment result. +#[derive(Clone, Debug)] +pub struct Adjusted { + pub cameras: Vec>, + /// RMS reprojection error (pixels) before and after. + pub rms_before: f64, + pub rms: f64, + pub iterations: usize, +} + +struct Obs { + i: usize, + j: usize, + a: Point, + b: Point, +} + +/// Residuals of one match: image-j point projected into i, and image-i point into j. +fn residual(ci: &Camera, cj: &Camera, a: Point, b: Point) -> [f64; 4] { + let big = 1e3; + let p = ci.project(cj.ray(b.x, b.y)); + let q = cj.project(ci.ray(a.x, a.y)); + let (r0, r1) = p.map_or((big, big), |(x, y)| (x - a.x, y - a.y)); + let (r2, r3) = q.map_or((big, big), |(x, y)| (x - b.x, y - b.y)); + [r0, r1, r2, r3] +} + +fn perturb(c: &Camera, k: usize, eps: f64) -> Camera { + let mut c = *c; + if k < 3 { + let mut w = [0.0; 3]; + w[k] = eps; + c.r = linalg::mul3(&linalg::rodrigues(w), &c.r); + } else { + c.f += eps; + } + c +} + +/// Levenberg–Marquardt over (rotation increment, focal) of every camera except the anchor's +/// rotation (gauge freedom). At most `max_per_pair` matches per pair are used. +pub fn bundle_adjust(cams: &[Option], pairs: &[PairMatch], anchor: usize, max_per_pair: usize) -> Adjusted { + let n = cams.len(); + let idx: Vec> = { + let mut k = 0; + cams.iter() + .map(|c| { + c.map(|_| { + k += 1; + k - 1 + }) + }) + .collect() + }; + let m = idx.iter().flatten().count(); + let np = 4 * m; + let mut obs = Vec::new(); + for p in pairs { + if cams[p.i].is_none() || cams[p.j].is_none() { + continue; + } + // prefer the precisely located (fine-level) matches when there are enough of them + let fine: Vec = (0..p.pts.len()).filter(|&k| p.levels.get(k).is_none_or(|l| *l <= 1)).collect(); + let use_idx: Vec = if fine.len() >= 16 { fine } else { (0..p.pts.len()).collect() }; + let step = (use_idx.len() / max_per_pair.max(1)).max(1); + for &k in use_idx.iter().step_by(step) { + let (a, b) = p.pts[k]; + obs.push(Obs { i: p.i, j: p.j, a, b }); + } + } + let mut cur: Vec> = cams.to_vec(); + let sigma = 2.0f64; + let cost = |cs: &[Option]| -> (f64, f64) { + let (mut c, mut sq) = (0.0, 0.0); + for o in &obs { + let r = residual(cs[o.i].as_ref().expect("cam"), cs[o.j].as_ref().expect("cam"), o.a, o.b); + for v in r { + let a = v.abs(); + sq += v * v; + c += if a <= sigma { 0.5 * v * v } else { sigma * (a - 0.5 * sigma) }; + } + } + (c, (sq / (obs.len().max(1) * 4) as f64).sqrt()) + }; + let (mut c0, rms_before) = cost(&cur); + let mut lambda = 1e-3; + let mut iterations = 0; + let _ = n; + for _ in 0..100 { + iterations += 1; + // normal equations (Huber IRLS weights) + let (jtj, jtr) = obs + .par_iter() + .fold( + || (vec![0.0f64; np * np], vec![0.0f64; np]), + |(mut jtj, mut jtr), o| { + let (ci, cj) = (cur[o.i].as_ref().expect("cam"), cur[o.j].as_ref().expect("cam")); + let r = residual(ci, cj, o.a, o.b); + let mut cols: Vec<(usize, [f64; 4])> = Vec::with_capacity(8); + for (cam, which) in [(o.i, 0), (o.j, 1)] { + let base = 4 * idx[cam].expect("idx"); + for k in 0..4 { + if k < 3 && cam == anchor { + continue; + } + let c = if which == 0 { ci } else { cj }; + let eps = if k < 3 { 1e-6 } else { 1e-4 * c.f }; + let pc = perturb(c, k, eps); + let rr = if which == 0 { residual(&pc, cj, o.a, o.b) } else { residual(ci, &pc, o.a, o.b) }; + let mut d = [0.0; 4]; + for t in 0..4 { + d[t] = (rr[t] - r[t]) / eps; + } + cols.push((base + k, d)); + } + } + for t in 0..4 { + let a = r[t].abs(); + let w = if a <= sigma { 1.0 } else { sigma / a }; + for &(p, dp) in &cols { + jtr[p] += w * dp[t] * r[t]; + for &(q, dq) in &cols { + jtj[p * np + q] += w * dp[t] * dq[t]; + } + } + } + (jtj, jtr) + }, + ) + .reduce( + || (vec![0.0f64; np * np], vec![0.0f64; np]), + |(mut a, mut b), (c, d)| { + for (x, y) in a.iter_mut().zip(&c) { + *x += y; + } + for (x, y) in b.iter_mut().zip(&d) { + *x += y; + } + (a, b) + }, + ); + let mut improved = false; + for _ in 0..10 { + let mut a = jtj.clone(); + for p in 0..np { + let d = a[p * np + p]; + a[p * np + p] = d + lambda * d.max(1e-9) + 1e-12; + } + let Some(delta) = linalg::solve(a, jtr.iter().map(|v| -v).collect(), np) else { + lambda *= 10.0; + continue; + }; + let mut trial = cur.clone(); + for (cam, c) in trial.iter_mut().enumerate() { + let (Some(c), Some(b)) = (c.as_mut(), idx[cam]) else { continue }; + let b = 4 * b; + if cam != anchor { + c.r = linalg::orthonormalize(&linalg::mul3(&linalg::rodrigues([delta[b], delta[b + 1], delta[b + 2]]), &c.r)); + } + c.f = (c.f + delta[b + 3]).max(1.0); + } + let (c1, _) = cost(&trial); + if c1 < c0 { + let rel = (c0 - c1) / c0.max(1e-300); + cur = trial; + c0 = c1; + lambda = (lambda * 0.3).max(1e-9); + improved = rel > 1e-9; + break; + } + lambda *= 10.0; + } + if !improved { + break; + } + } + let (_, rms) = cost(&cur); + Adjusted { cameras: cur, rms_before, rms, iterations } +} + +/// Rotate the world so that it is level: "down" (+y) is the direction most perpendicular to every +/// camera's x axis, and "forward" (+z) the mean viewing direction. +pub fn straighten(cams: &mut [Option]) { + let rs: Vec = cams.iter().flatten().map(|c| c.r).collect(); + if rs.is_empty() { + return; + } + let mut cov = [0.0f64; 9]; + let (mut ysum, mut zsum, mut xsum) = ([0.0; 3], [0.0; 3], [0.0; 3]); + for r in &rs { + let (x, y, z) = (r[0], r[1], r[2]); + for a in 0..3 { + for b in 0..3 { + cov[a * 3 + b] += x[a] * x[b]; + } + ysum[a] += y[a]; + zsum[a] += z[a]; + xsum[a] += x[a]; + } + } + let (vals, vecs) = linalg::sym_eigen(&cov, 3); + let total = vals.iter().sum::().max(1e-12); + let mut down = [vecs[0][0], vecs[0][1], vecs[0][2]]; + if vals[1] < 1e-3 * total { + // all x axes (nearly) parallel (a vertical sweep): use the mean camera "down" made + // perpendicular to the common x axis + let xa = linalg::normalize(xsum); + let y = linalg::normalize(ysum); + let k = linalg::dot(y, xa); + down = linalg::normalize([y[0] - k * xa[0], y[1] - k * xa[1], y[2] - k * xa[2]]); + } + if linalg::dot(down, ysum) < 0.0 { + down = down.map(|v| -v); + } + let z = linalg::normalize(zsum); + let k = linalg::dot(z, down); + let fwd = linalg::normalize([z[0] - k * down[0], z[1] - k * down[1], z[2] - k * down[2]]); + let right = linalg::cross(down, fwd); + let w = [right, down, fwd]; + let wt = linalg::transpose3(&w); + for c in cams.iter_mut().flatten() { + c.r = linalg::mul3(&c.r, &wt); + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn focal_from_a_rotation_homography() { + let f = 900.0; + let r = linalg::rodrigues([0.05, 0.3, 0.02]); + let k = [[f, 0.0, 0.0], [0.0, f, 0.0], [0.0, 0.0, 1.0]]; + let kinv = [[1.0 / f, 0.0, 0.0], [0.0, 1.0 / f, 0.0], [0.0, 0.0, 1.0]]; + let h = linalg::mul3(&linalg::mul3(&k, &r), &kinv); + let hh = Homography([h[0][0], h[0][1], h[0][2], h[1][0], h[1][1], h[1][2], h[2][0], h[2][1], h[2][2]]); + let (a, b) = focals_from_homography(&hh); + assert!((a.unwrap() - f).abs() < 1.0, "{a:?}"); + assert!((b.unwrap() - f).abs() < 1.0, "{b:?}"); + } +} diff --git a/crates/merge/src/pano/finish.rs b/crates/merge/src/pano/finish.rs new file mode 100644 index 00000000..ac57a18d --- /dev/null +++ b/crates/merge/src/pano/finish.rs @@ -0,0 +1,210 @@ +//! Panorama finishing: boundary warp, auto crop, fill edges. +//! +//! - **Boundary warp** (0–100): a separable stretch that maps, per column, the smoothed top and +//! bottom edges of the valid area to the canvas edges, then per row the left and right edges; +//! `amount` blends between the identity and that stretch. (A simple, documented alternative to +//! content-preserving mesh warps.) +//! - **Auto crop**: the largest axis-aligned rectangle of fully valid pixels (maximal rectangle in a +//! binary matrix via a per-row histogram and a monotone stack, O(w·h)). +//! - **Fill edges**: transparent areas are filled by push–pull diffusion of the surrounding pixels +//! in log space (no patch synthesis). + +use lightcraft_raster::{Plane, Rgb32f}; +use rayon::prelude::*; + +use super::blend::{Buf, push_pull}; + +fn edges(alpha: &Plane, vertical: bool) -> Option<(Vec, Vec)> { + let (w, h) = (alpha.width, alpha.height); + let (n, len) = if vertical { (w, h) } else { (h, w) }; + let at = |i: usize, k: usize| if vertical { alpha.get(i, k) } else { alpha.get(k, i) }; + let mut lo = vec![f64::NAN; n]; + let mut hi = vec![f64::NAN; n]; + for i in 0..n { + if let Some(a) = (0..len).find(|&k| at(i, k) >= 0.5) { + let b = (0..len).rev().find(|&k| at(i, k) >= 0.5).unwrap_or(a); + lo[i] = a as f64; + hi[i] = b as f64 + 1.0; + } + } + let known: Vec = (0..n).filter(|&i| lo[i].is_finite()).collect(); + if known.len() < 2 { + return None; + } + // fill gaps by linear interpolation / nearest at the ends + for i in 0..n { + if lo[i].is_finite() { + continue; + } + let p = known.iter().rev().find(|&&k| k < i).copied(); + let q = known.iter().find(|&&k| k > i).copied(); + let (a, b) = match (p, q) { + (Some(p), Some(q)) => { + let t = (i - p) as f64 / (q - p) as f64; + (lo[p] + (lo[q] - lo[p]) * t, hi[p] + (hi[q] - hi[p]) * t) + } + (Some(p), None) => (lo[p], hi[p]), + (None, Some(q)) => (lo[q], hi[q]), + _ => (0.0, len as f64), + }; + lo[i] = a; + hi[i] = b; + } + // smooth (box, radius n/40) so the stretch follows the overall shape, not every jag + let r = (n / 40).max(1) as isize; + let smooth = |v: &[f64]| -> Vec { + (0..n as isize) + .map(|i| { + let (a, b) = ((i - r).max(0) as usize, ((i + r) as usize).min(n - 1)); + v[a..=b].iter().sum::() / (b - a + 1) as f64 + }) + .collect() + }; + Some((smooth(&lo), smooth(&hi))) +} + +fn stretch(img: &Rgb32f, alpha: &Plane, amount: f64, vertical: bool) -> (Rgb32f, Plane) { + let Some((lo, hi)) = edges(alpha, vertical) else { return (img.clone(), alpha.clone()) }; + let (w, h) = (img.width, img.height); + let len = if vertical { h } else { w } as f64; + let mut out = Rgb32f::new(w, h); + let mut oa = Plane::new(w, h); + out.data.par_chunks_mut(w).zip(oa.data.par_chunks_mut(w)).enumerate().for_each(|(y, (row, arow))| { + for x in 0..w { + let (i, k) = if vertical { (x, y) } else { (y, x) }; + let t = (k as f64 + 0.5) / len; + let src = lo[i] + t * (hi[i] - lo[i]); + let s = (1.0 - amount) * (k as f64 + 0.5) + amount * src; + let (sx, sy) = if vertical { (x as f32 + 0.5, s as f32) } else { (s as f32, y as f32 + 0.5) }; + row[x] = img.sample_bilinear(sx, sy); + arow[x] = crate::features::bilinear(alpha, sx, sy); + } + }); + (out, oa) +} + +/// Boundary warp by `amount` (0..1). +pub fn boundary_warp(img: &Rgb32f, alpha: &Plane, amount: f64) -> (Rgb32f, Plane) { + if amount <= 0.0 { + return (img.clone(), alpha.clone()); + } + let amount = amount.min(1.0); + let (i1, a1) = stretch(img, alpha, amount, true); + stretch(&i1, &a1, amount, false) +} + +/// Largest axis-aligned rectangle where `alpha ≥ 0.999`: (x, y, w, h) in pixels. +pub fn auto_crop(alpha: &Plane) -> Option<(usize, usize, usize, usize)> { + let (w, h) = (alpha.width, alpha.height); + if w == 0 || h == 0 { + return None; + } + // work on a ≤ 1500 px grid with conservative (min) pooling + let f = w.max(h).div_ceil(1500).max(1); + let (gw, gh) = (w.div_ceil(f), h.div_ceil(f)); + let ok = |gx: usize, gy: usize| -> bool { + for y in gy * f..((gy + 1) * f).min(h) { + for x in gx * f..((gx + 1) * f).min(w) { + if alpha.get(x, y) < 0.999 { + return false; + } + } + } + true + }; + let grid: Vec = (0..gh).into_par_iter().flat_map_iter(|gy| (0..gw).map(move |gx| (gx, gy))).map(|(gx, gy)| ok(gx, gy)).collect(); + let mut heights = vec![0usize; gw]; + let mut best = (0usize, 0, 0, 0, 0); // area, x, y, w, h + for gy in 0..gh { + for gx in 0..gw { + heights[gx] = if grid[gy * gw + gx] { heights[gx] + 1 } else { 0 }; + } + let mut stack: Vec = Vec::new(); + for i in 0..=gw { + let cur = if i < gw { heights[i] } else { 0 }; + while let Some(&top) = stack.last() { + if heights[top] <= cur { + break; + } + stack.pop(); + let hgt = heights[top]; + let left = stack.last().map_or(0, |&s| s + 1); + let width = i - left; + let area = hgt * width; + if area > best.0 { + best = (area, left, gy + 1 - hgt, width, hgt); + } + } + stack.push(i); + } + } + if best.0 == 0 { + return None; + } + let (x, y) = (best.1 * f, best.2 * f); + Some((x, y, (best.3 * f).min(w - x), (best.4 * f).min(h - y))) +} + +/// Fill transparent pixels (alpha < 1) from their surroundings; returns the image with alpha 1. +pub fn fill_edges(img: &Rgb32f, alpha: &Plane) -> Rgb32f { + let (w, h) = (img.width, img.height); + let mut b = Buf::new(w, h, 3); + for (i, p) in img.data.iter().enumerate() { + for c in 0..3 { + b.data[i * 3 + c] = (p[c].max(1e-6)).ln(); + } + } + let wgt: Vec = alpha.data.iter().map(|a| if *a >= 0.999 { 1.0 } else { 0.0 }).collect(); + push_pull(&mut b, &wgt); + let mut out = img.clone(); + for (i, p) in out.data.iter_mut().enumerate() { + let a = alpha.data[i].clamp(0.0, 1.0); + if a < 0.999 { + let f = [0, 1, 2].map(|c| b.data[i * 3 + c].exp()); + // the warped/blended edge pixels are premultiplied-like: un-premultiply then mix + let v = if a > 1e-3 { p.map(|v| v / a) } else { [0.0; 3] }; + *p = [0, 1, 2].map(|c| v[c] * a + f[c] * (1.0 - a)); + } + } + out +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn crop_finds_the_largest_rectangle() { + let mut a = Plane::new(100, 60); + for y in 10..50 { + for x in 5..90 { + a.set(x, y, 1.0); + } + } + for y in 0..60 { + a.set(40, y, 1.0); // a thin spike doesn't matter + } + assert_eq!(auto_crop(&a), Some((5, 10, 85, 40))); + } + + #[test] + fn warp_fills_a_curved_boundary() { + let (w, h) = (200, 100); + let mut a = Plane::new(w, h); + let img = Rgb32f::from_fn(w, h, |x, _| [(x + 1) as f32 / w as f32; 3]); + for x in 0..w { + let dip = (10.0 * (1.0 - ((x as f32 - 100.0) / 100.0).powi(2))) as usize; + for y in dip..h - dip { + a.set(x, y, 1.0); + } + } + let before = a.data.iter().filter(|v| **v >= 0.5).count(); + let (_, a2) = boundary_warp(&img, &a, 1.0); + let after = a2.data.iter().filter(|v| **v >= 0.5).count(); + assert!(after as f64 > 0.97 * (w * h) as f64 && after > before, "{before} → {after}"); + let (_, a0) = boundary_warp(&img, &a, 0.0); + assert_eq!(a0, a); + let filled = fill_edges(&img, &a); + assert!(filled.data.iter().all(|p| p[0] > 0.0)); + } +} diff --git a/crates/merge/src/pano/mod.rs b/crates/merge/src/pano/mod.rs new file mode 100644 index 00000000..d81abda9 --- /dev/null +++ b/crates/merge/src/pano/mod.rs @@ -0,0 +1,461 @@ +//! Panorama stitching. +//! +//! 1. Features ([`crate::features`]) on each frame at ≤ 1000 px, pairwise matching and RANSAC +//! homographies, largest connected set ([`camera`]). +//! 2. Camera rotations + focal lengths: initialised from the homographies, bundle-adjusted, +//! straightened ([`camera`]). +//! 3. Projection: spherical, cylindrical or perspective (rectilinear), or chosen automatically +//! from the field of view ([`choose_projection`]). +//! 4. Exposure compensation: EXIF exposure, then per-image gains solved from the mean intensities +//! of every overlap (Brown & Lowe's gain compensation, a small linear least-squares problem +//! with a prior pulling gains towards 1). +//! 5. Seams: each canvas pixel belongs to the image whose centre-weighted feather weight is +//! highest there (Voronoi-like seams away from image borders); multi-band blending in log space +//! ([`blend`]). +//! 6. Boundary warp, auto crop, fill edges ([`finish`]). + +pub mod blend; +pub mod camera; +pub mod finish; +pub mod project; + +use lightcraft_raster::{Plane, Rgb32f}; +use rayon::prelude::*; +use serde::{Deserialize, Serialize}; + +pub use project::Projection; + +use crate::frame::{Frame, FrameColor}; +use crate::linalg; +use crate::{MergeError, Progress, Result, check}; +use blend::{Blender, Buf, push_pull}; +use camera::Camera; + +#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)] +#[serde(default, rename_all = "camelCase")] +pub struct PanoOptions { + pub projection: Projection, + /// 0..100. + pub boundary_warp: f64, + pub auto_crop: bool, + pub fill_edges: bool, + /// Output size limit (megapixels); the panorama is scaled down to fit. + pub max_megapixels: f64, +} + +impl Default for PanoOptions { + fn default() -> Self { + PanoOptions { projection: Projection::Auto, boundary_warp: 0.0, auto_crop: false, fill_edges: false, max_megapixels: 40.0 } + } +} + +/// The stitched panorama. +#[derive(Clone, Debug)] +pub struct PanoResult { + /// Scene-linear pixels (same units as the frames, after exposure compensation). + pub image: Rgb32f, + /// Coverage (1 = image data, 0 = empty canvas). + pub alpha: Plane, + /// Auto-crop rectangle (normalised x, y, w, h), when requested. + pub crop: Option<[f64; 4]>, + /// The projection used (never `Auto`). + pub projection: Projection, + /// Indices (into the inputs) of the frames in the panorama; the others didn't connect. + pub used: Vec, + /// Cameras of the inputs (`None` = not used), in their full-resolution pixels, straightened. + pub cameras: Vec>, + /// Per-input exposure gains applied. + pub gains: Vec, + /// Bundle adjustment RMS reprojection error (pixels) before/after. + pub rms_before: f64, + pub rms: f64, + /// Horizontal and vertical field of view (degrees). + pub fov: (f64, f64), + /// Output focal length and the surface coordinates of the canvas's top-left corner: canvas + /// pixel (x, y) (before boundary warp) shows surface point (x + 0.5 + origin.0, y + 0.5 + origin.1). + pub focal: f64, + pub origin: (f64, f64), + pub color: FrameColor, + pub metadata: lightcraft_meta::Metadata, + pub raw: bool, + pub baseline_exposure: f64, +} + +/// Field of view (degrees) of the cameras' union, measured in the spherical projection. +fn fov(frames: &[Frame], cams: &[Option]) -> (f64, f64) { + let (mut tmin, mut tmax, mut pmin, mut pmax) = (f64::MAX, f64::MIN, f64::MAX, f64::MIN); + for (f, c) in frames.iter().zip(cams) { + let Some(c) = c else { continue }; + for (x, y) in border(f.width(), f.height(), 16) { + if let Some((t, p)) = Projection::Spherical.forward(c.ray(x, y), 1.0) { + tmin = tmin.min(t); + tmax = tmax.max(t); + pmin = pmin.min(p); + pmax = pmax.max(p); + } + } + } + ((tmax - tmin).to_degrees(), (pmax - pmin).to_degrees()) +} + +/// Points along an image's border. +fn border(w: usize, h: usize, per_edge: usize) -> Vec<(f64, f64)> { + let (w, h) = (w as f64, h as f64); + let mut v = Vec::new(); + for k in 0..=per_edge { + let t = k as f64 / per_edge as f64; + v.extend([(t * w, 0.0), (t * w, h), (0.0, t * h), (w, t * h)]); + } + v +} + +/// Surface bounding box of every used image in `proj` (output focal `f`). +fn bounds(frames: &[Frame], cams: &[Option], proj: Projection, f: f64) -> Option<(f64, f64, f64, f64)> { + let (mut x0, mut y0, mut x1, mut y1) = (f64::MAX, f64::MAX, f64::MIN, f64::MIN); + for (fr, c) in frames.iter().zip(cams) { + let Some(c) = c else { continue }; + for (x, y) in border(fr.width(), fr.height(), 24) { + let (u, v) = proj.forward(c.ray(x, y), f)?; + x0 = x0.min(u); + y0 = y0.min(v); + x1 = x1.max(u); + y1 = y1.max(v); + } + } + (x1 > x0 && y1 > y0).then_some((x0, y0, x1, y1)) +} + +/// Auto projection: perspective when the field of view is moderate and the rectilinear canvas +/// stays compact; cylindrical for wide but not tall panoramas; spherical otherwise. +pub fn choose_projection(frames: &[Frame], cams: &[Option]) -> Projection { + let (h, v) = fov(frames, cams); + let area: f64 = frames.iter().zip(cams).filter(|(_, c)| c.is_some()).map(|(f, _)| (f.width() * f.height()) as f64).sum(); + let fmed = median(cams.iter().flatten().map(|c| c.f).collect()); + if h < 110.0 + && v < 110.0 + && let Some((x0, y0, x1, y1)) = bounds(frames, cams, Projection::Perspective, fmed) + && (x1 - x0) * (y1 - y0) < 2.5 * area + { + return Projection::Perspective; + } + if v < 90.0 { Projection::Cylindrical } else { Projection::Spherical } +} + +fn median(mut v: Vec) -> f64 { + if v.is_empty() { + return 1.0; + } + v.sort_by(|a, b| a.total_cmp(b)); + v[v.len() / 2] +} + +/// One image resampled onto its canvas tile. +struct Tile { + x0: usize, + y0: usize, + w: usize, + h: usize, + rgb: Vec<[f32; 3]>, + /// Feather weight (0 = outside the image). + wgt: Vec, +} + +/// Register the frames: cameras for the connected set. +pub fn register(frames: &[Frame], progress: &Progress) -> Result<(Vec>, f64, f64)> { + let n = frames.len(); + let work: Vec<(crate::features::Features, f64)> = frames + .par_iter() + .map(|f| { + let med = crate::align::median_luma(&f.image, f.clip); + let (p, s) = crate::align::work_image(&f.image, 0.18 / med.max(1e-6), 1000); + (crate::features::detect(&p, 1200), s) + }) + .collect(); + check(progress, 0.2, "Matching")?; + let feats: Vec<_> = work.iter().map(|(f, _)| f.clone()).collect(); + let scales: Vec = work.iter().map(|(_, s)| *s).collect(); + let pairs = camera::match_pairs(&feats, &scales); + let ids = camera::largest_component(n, &pairs); + if ids.len() < 2 { + return Err(MergeError::NoOverlap("the photos don't overlap enough to be stitched".into())); + } + check(progress, 0.3, "Aligning")?; + // focal length: EXIF 35 mm equivalent, else the homographies, else ~50° horizontal + let sizes: Vec<(usize, usize)> = frames.iter().map(|f| (f.width(), f.height())).collect(); + let mut fs = Vec::new(); + for p in &pairs { + if !(ids.contains(&p.i) && ids.contains(&p.j)) { + continue; + } + let ci = (sizes[p.i].0 as f64 / 2.0, sizes[p.i].1 as f64 / 2.0); + let cj = (sizes[p.j].0 as f64 / 2.0, sizes[p.j].1 as f64 / 2.0); + let ti = lightcraft_geom::Homography([1.0, 0.0, ci.0, 0.0, 1.0, ci.1, 0.0, 0.0, 1.0]); + let tj = lightcraft_geom::Homography([1.0, 0.0, -cj.0, 0.0, 1.0, -cj.1, 0.0, 0.0, 1.0]); + let hc = tj.mul(&p.h).mul(&ti); + if let (Some(a), Some(b)) = camera::focals_from_homography(&hc) { + fs.push((a * b).sqrt()); + } + } + let long = sizes[ids[0]].0.max(sizes[ids[0]].1) as f64; + let exif_f = frames[ids[0]].metadata.focal_length_35mm.filter(|v| *v > 1.0).map(|f35| f35 / 36.0 * long); + let hom_f = (!fs.is_empty()).then(|| median(fs.clone())).filter(|f| *f > 0.1 * long && *f < 20.0 * long); + let f0 = exif_f.or(hom_f).unwrap_or(long * 1.07); + let cams = camera::initial_cameras(&ids, &sizes, &pairs, f0); + let anchor = ids.iter().copied().find(|&i| cams[i].is_some()).unwrap_or(ids[0]); + check(progress, 0.35, "Aligning")?; + let adj = camera::bundle_adjust(&cams, &pairs, anchor, 60); + let mut cams = adj.cameras; + camera::straighten(&mut cams); + Ok((cams, adj.rms_before, adj.rms)) +} + +/// Stitch `frames` (oriented, linear) into a panorama. +pub fn stitch(frames: Vec, opts: &PanoOptions, progress: &Progress) -> Result { + let n = frames.len(); + if n < 2 { + return Err(MergeError::TooFew(n, 2)); + } + check(progress, 0.01, "Finding features")?; + let (cams, rms_before, rms) = register(&frames, progress)?; + let used: Vec = (0..n).filter(|&i| cams[i].is_some()).collect(); + let fov = fov(&frames, &cams); + let projection = match opts.projection { + Projection::Auto => choose_projection(&frames, &cams), + p => p, + }; + let fmed = median(used.iter().map(|&i| cams[i].expect("used").f).collect()); + let Some((bx0, by0, bx1, by1)) = bounds(&frames, &cams, projection, fmed) else { + return Err(MergeError::NoOverlap(format!("the field of view ({:.0}°) is too wide for a perspective projection", fov.0))); + }; + let area: f64 = used.iter().map(|&i| (frames[i].width() * frames[i].height()) as f64).sum(); + if projection == Projection::Perspective && (bx1 - bx0) * (by1 - by0) > 12.0 * area { + return Err(MergeError::NoOverlap(format!("the field of view ({:.0}°) is too wide for a perspective projection", fov.0))); + } + // output scale: native resolution, limited to max_megapixels + let max_px = (opts.max_megapixels.max(0.1) * 1e6).min(400e6); + let s = (max_px / ((bx1 - bx0) * (by1 - by0))).sqrt().min(1.0); + let f_out = fmed * s; + let (ox, oy) = (bx0 * s, by0 * s); + let (cw, ch) = (((bx1 - bx0) * s).ceil() as usize + 1, ((by1 - by0) * s).ceil() as usize + 1); + let levels = ((cw.min(ch) as f64).log2() - 3.0).floor().clamp(1.0, 7.0) as usize; + let unit = 1usize << levels; + let (pw, ph) = (cw.div_ceil(unit) * unit, ch.div_ceil(unit) * unit); + check(progress, 0.45, "Projecting")?; + + // exposure from EXIF relative to the first used frame + let ev0 = frames[used[0]].exposure; + let mut gains: Vec = (0..n) + .map(|i| match (frames[i].exposure, ev0) { + (Some(e), Some(e0)) => 2f64.powf(e0 - e), + _ => 1.0, + }) + .collect(); + + // resample every used frame onto its tile + let tiles: Vec> = (0..n) + .into_par_iter() + .map(|i| { + let c = cams[i]?; + let fr = &frames[i]; + let (mut x0, mut y0, mut x1, mut y1) = (f64::MAX, f64::MAX, f64::MIN, f64::MIN); + for (x, y) in border(fr.width(), fr.height(), 32) { + let (u, v) = projection.forward(c.ray(x, y), f_out)?; + x0 = x0.min(u - ox); + y0 = y0.min(v - oy); + x1 = x1.max(u - ox); + y1 = y1.max(v - oy); + } + let pad = 2.0 * unit as f64; + let tx0 = (((x0 - pad).max(0.0) as usize) / unit) * unit; + let ty0 = (((y0 - pad).max(0.0) as usize) / unit) * unit; + let tx1 = (((x1 + pad).max(0.0) as usize).min(pw)).div_ceil(unit) * unit; + let ty1 = (((y1 + pad).max(0.0) as usize).min(ph)).div_ceil(unit) * unit; + let (tw, th) = (tx1.min(pw).saturating_sub(tx0), ty1.min(ph).saturating_sub(ty0)); + if tw == 0 || th == 0 { + return None; + } + let (fw, fh) = (fr.width() as f64, fr.height() as f64); + let mut rgb = vec![[0f32; 3]; tw * th]; + let mut wgt = vec![0f32; tw * th]; + rgb.par_chunks_mut(tw).zip(wgt.par_chunks_mut(tw)).enumerate().for_each(|(ty, (row, wrow))| { + let v = (ty0 + ty) as f64 + 0.5 + oy; + for tx in 0..tw { + let u = (tx0 + tx) as f64 + 0.5 + ox; + let d = projection.inverse(u, v, f_out); + let Some((px, py)) = c.project(d) else { continue }; + if px < 0.0 || py < 0.0 || px > fw || py > fh { + continue; + } + let wx = 1.0 - (2.0 * px / fw - 1.0).abs(); + let wy = 1.0 - (2.0 * py / fh - 1.0).abs(); + let w = (wx * wy).max(1e-6) as f32; + row[tx] = fr.image.sample_bilinear(px as f32, py as f32); + wrow[tx] = w; + } + }); + Some(Tile { x0: tx0, y0: ty0, w: tw, h: th, rgb, wgt }) + }) + .collect(); + check(progress, 0.6, "Exposure compensation")?; + compensate_gains(&tiles, &mut gains, &frames); + + // seams: label = max feather weight + let mut best = vec![0f32; pw * ph]; + let mut label = vec![u16::MAX; pw * ph]; + for (i, t) in tiles.iter().enumerate() { + let Some(t) = t else { continue }; + for ty in 0..t.h { + let row = (t.y0 + ty) * pw + t.x0; + for tx in 0..t.w { + let w = t.wgt[ty * t.w + tx]; + if w > best[row + tx] { + best[row + tx] = w; + label[row + tx] = i as u16; + } + } + } + } + check(progress, 0.7, "Blending")?; + let mut blender = Blender::new(pw, ph, levels); + for (i, t) in tiles.into_iter().enumerate() { + let Some(t) = t else { continue }; + let g = gains[i] as f32; + let mut vals = Buf::new(t.w, t.h, 3); + let mut mask = Buf::new(t.w, t.h, 1); + for ty in 0..t.h { + for tx in 0..t.w { + let k = ty * t.w + tx; + if t.wgt[k] > 0.0 { + for c in 0..3 { + vals.data[k * 3 + c] = (t.rgb[k][c] * g).max(1e-6).ln(); + } + } + if label[(t.y0 + ty) * pw + t.x0 + tx] == i as u16 { + mask.data[k] = 1.0; + } + } + } + push_pull(&mut vals, &t.wgt); + blender.add(t.x0, t.y0, vals, mask); + check(progress, 0.7 + 0.2 * (i + 1) as f32 / n as f32, "Blending")?; + } + let out = blender.finish(); + let mut image = Rgb32f::new(cw, ch); + let mut alpha = Plane::new(cw, ch); + for y in 0..ch { + for x in 0..cw { + let k = y * pw + x; + if label[k] != u16::MAX { + image.data[y * cw + x] = [0, 1, 2].map(|c| out.data[k * 3 + c].exp()); + alpha.data[y * cw + x] = 1.0; + } + } + } + check(progress, 0.92, "Finishing")?; + let (mut image, alpha) = finish::boundary_warp(&image, &alpha, opts.boundary_warp / 100.0); + let crop = if opts.auto_crop { + finish::auto_crop(&alpha).map(|(x, y, w, h)| [x as f64 / cw as f64, y as f64 / ch as f64, w as f64 / cw as f64, h as f64 / ch as f64]) + } else { + None + }; + if opts.fill_edges { + image = finish::fill_edges(&image, &alpha); + } + let r = &frames[used[0]]; + Ok(PanoResult { + image, + alpha, + crop, + projection, + used, + cameras: cams, + gains, + rms_before, + rms, + fov, + focal: f_out, + origin: (ox, oy), + color: r.color.clone(), + metadata: r.metadata.clone(), + raw: r.raw, + baseline_exposure: r.baseline_exposure, + }) +} + +/// Gain compensation (Brown & Lowe §6): minimise +/// Σᵢⱼ Nᵢⱼ [ (gᵢ Īᵢⱼ − gⱼ Īⱼᵢ)² / σN² + (1 − gᵢ/g⁰ᵢ)² / σg² ] over the overlaps, with intensities +/// relative to the mean overlap level (σN = 0.1; σg = 1, a weak prior that only fixes the overall +/// level) and `g⁰` the EXIF-based gains. +fn compensate_gains(tiles: &[Option], gains: &mut [f64], frames: &[Frame]) { + let n = tiles.len(); + let lum = |p: [f32; 3]| 0.25 * p[0] as f64 + 0.6 * p[1] as f64 + 0.15 * p[2] as f64; + let mut nij = vec![0f64; n * n]; + let mut iij = vec![0f64; n * n]; + for i in 0..n { + let Some(a) = &tiles[i] else { continue }; + for j in 0..n { + let Some(b) = tiles[j].as_ref().filter(|_| i != j) else { continue }; + let (x0, y0) = (a.x0.max(b.x0), a.y0.max(b.y0)); + let (x1, y1) = ((a.x0 + a.w).min(b.x0 + b.w), (a.y0 + a.h).min(b.y0 + b.h)); + let (mut cnt, mut sum) = (0usize, 0f64); + for y in (y0..y1).step_by(2) { + for x in (x0..x1).step_by(2) { + let ka = (y - a.y0) * a.w + x - a.x0; + let kb = (y - b.y0) * b.w + x - b.x0; + if a.wgt[ka] > 0.0 && b.wgt[kb] > 0.0 { + let pa = a.rgb[ka]; + if pa.iter().any(|v| *v >= frames[i].clip * 0.95) || b.rgb[kb].iter().any(|v| *v >= frames[j].clip * 0.95) { + continue; + } + cnt += 1; + sum += lum(pa); + } + } + } + if cnt > 50 { + nij[i * n + j] = cnt as f64; + iij[i * n + j] = sum / cnt as f64; + } + } + } + let total: f64 = nij.iter().sum(); + if total <= 0.0 { + return; + } + let mean = (0..n * n).map(|k| nij[k] * iij[k] * gains[k / n]).sum::() / total; + if mean <= 0.0 { + return; + } + let (sn, sg) = (0.1f64, 1.0f64); + let g0: Vec = gains.to_vec(); + let mut a = vec![0f64; n * n]; + let mut b = vec![0f64; n]; + for i in 0..n { + if tiles[i].is_none() { + a[i * n + i] = 1.0; + b[i] = g0[i]; + continue; + } + for j in 0..n { + let nn = nij[i * n + j]; + if nn <= 0.0 || nij[j * n + i] <= 0.0 { + continue; + } + let (ii, ij) = (iij[i * n + j] / mean, iij[j * n + i] / mean); + // d/dgᵢ of Nᵢⱼ (gᵢ Iᵢⱼ − gⱼ Iⱼᵢ)²/σN² (both orderings of the pair are in the sum) + a[i * n + i] += 2.0 * nn * ii * ii / (sn * sn); + a[i * n + j] -= 2.0 * nn * ii * ij / (sn * sn); + a[i * n + i] += 2.0 * nn / (sg * sg * g0[i] * g0[i]); + b[i] += 2.0 * nn / (sg * sg * g0[i]); + } + if a[i * n + i] == 0.0 { + a[i * n + i] = 1.0; + b[i] = g0[i]; + } + } + if let Some(g) = linalg::solve(a, b, n) + && g.iter().all(|v| v.is_finite() && *v > 0.0) + { + gains.copy_from_slice(&g); + } +} diff --git a/crates/merge/src/pano/project.rs b/crates/merge/src/pano/project.rs new file mode 100644 index 00000000..0407f472 --- /dev/null +++ b/crates/merge/src/pano/project.rs @@ -0,0 +1,83 @@ +//! Panorama projections. World directions use camera-style axes: x right, y down, z forward. +//! Surface coordinates are in units of the output focal length `f` (pixels). + +use serde::{Deserialize, Serialize}; + +#[derive(Clone, Copy, Debug, Default, PartialEq, Eq, Serialize, Deserialize)] +#[serde(rename_all = "lowercase")] +pub enum Projection { + /// Pick from the field of view (see [`crate::pano::choose_projection`]). + #[default] + Auto, + /// Longitude/latitude (equirectangular): any field of view; straight lines bend. + Spherical, + /// Longitude/height on a cylinder: wide horizontal panoramas; verticals stay straight. + Cylindrical, + /// Rectilinear (gnomonic): straight lines stay straight; fields of view < ~120°. + Perspective, +} + +impl Projection { + /// Project a world direction to surface coordinates (`None` when not representable). + #[inline] + pub fn forward(self, d: [f64; 3], f: f64) -> Option<(f64, f64)> { + match self { + Projection::Perspective => { + if d[2] <= 1e-6 { + return None; + } + Some((f * d[0] / d[2], f * d[1] / d[2])) + } + Projection::Cylindrical => { + let r = (d[0] * d[0] + d[2] * d[2]).sqrt(); + if r < 1e-12 { + return None; + } + Some((f * d[0].atan2(d[2]), f * d[1] / r)) + } + Projection::Spherical | Projection::Auto => { + let n = (d[0] * d[0] + d[1] * d[1] + d[2] * d[2]).sqrt(); + if n < 1e-12 { + return None; + } + Some((f * d[0].atan2(d[2]), f * (d[1] / n).clamp(-1.0, 1.0).asin())) + } + } + } + + /// The world direction (not normalised) of surface coordinates. + #[inline] + pub fn inverse(self, x: f64, y: f64, f: f64) -> [f64; 3] { + match self { + Projection::Perspective => [x / f, y / f, 1.0], + Projection::Cylindrical => { + let t = x / f; + [t.sin(), y / f, t.cos()] + } + Projection::Spherical | Projection::Auto => { + let (t, p) = (x / f, y / f); + [t.sin() * p.cos(), p.sin(), t.cos() * p.cos()] + } + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn projections_round_trip() { + for p in [Projection::Perspective, Projection::Cylindrical, Projection::Spherical] { + for d in [[0.1, -0.2, 1.0], [-0.5, 0.3, 0.8], [0.0, 0.0, 1.0]] { + let (x, y) = p.forward(d, 1000.0).unwrap(); + let back = p.inverse(x, y, 1000.0); + let (n1, n2) = (crate::linalg::normalize(d), crate::linalg::normalize(back)); + for c in 0..3 { + assert!((n1[c] - n2[c]).abs() < 1e-9, "{p:?} {d:?} {back:?}"); + } + } + } + assert!(Projection::Perspective.forward([1.0, 0.0, -0.1], 1.0).is_none()); + } +} diff --git a/crates/merge/src/ransac.rs b/crates/merge/src/ransac.rs new file mode 100644 index 00000000..b0937643 --- /dev/null +++ b/crates/merge/src/ransac.rs @@ -0,0 +1,291 @@ +//! Robust 2-D motion models: normalised DLT homographies (Hartley & Zisserman, "Multiple View +//! Geometry", §4.4) and RANSAC (Fischler & Bolles 1981) with the adaptive iteration count. +//! Deterministic: the sampler is a fixed-seed xorshift, so results are reproducible. + +use lightcraft_geom::{Homography, Point}; + +use crate::linalg::sym_eigen; + +/// The kind of motion to fit. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum Model { + Translation, + /// Rotation + uniform scale + translation (4 dof). + Similarity, + Homography, +} + +impl Model { + fn min_samples(self) -> usize { + match self { + Model::Translation => 1, + Model::Similarity => 2, + Model::Homography => 4, + } + } +} + +#[derive(Clone, Debug)] +pub struct Fit { + /// Maps points of the first set to the second. + pub h: Homography, + pub inliers: Vec, + /// RMS transfer error of the inliers (pixels). + pub rms: f64, +} + +/// Small deterministic PRNG (xorshift64*). +pub struct Rng(u64); + +impl Rng { + pub fn new(seed: u64) -> Rng { + Rng(seed.max(1)) + } + pub fn next_u64(&mut self) -> u64 { + let mut x = self.0; + x ^= x >> 12; + x ^= x << 25; + x ^= x >> 27; + self.0 = x; + x.wrapping_mul(0x2545_f491_4f6c_dd1d) + } + pub fn below(&mut self, n: usize) -> usize { + (self.next_u64() % n.max(1) as u64) as usize + } + /// Uniform in [0, 1). + pub fn unit(&mut self) -> f64 { + (self.next_u64() >> 11) as f64 / (1u64 << 53) as f64 + } +} + +/// Least-squares fit of `model` to all pairs (`a[i]` → `b[i]`). +pub fn fit_model(model: Model, a: &[Point], b: &[Point]) -> Option { + let n = a.len().min(b.len()); + if n < model.min_samples() { + return None; + } + match model { + Model::Translation => { + let (mut dx, mut dy) = (0.0, 0.0); + for i in 0..n { + dx += b[i].x - a[i].x; + dy += b[i].y - a[i].y; + } + Some(Homography([1.0, 0.0, dx / n as f64, 0.0, 1.0, dy / n as f64, 0.0, 0.0, 1.0])) + } + Model::Similarity => { + // closed form (Umeyama without reflection): b ≈ s R a + t + let ca = centroid(&a[..n]); + let cb = centroid(&b[..n]); + let (mut sxx, mut sxy, mut var) = (0.0, 0.0, 0.0); + for i in 0..n { + let (ax, ay) = (a[i].x - ca.x, a[i].y - ca.y); + let (bx, by) = (b[i].x - cb.x, b[i].y - cb.y); + sxx += ax * bx + ay * by; + sxy += ax * by - ay * bx; + var += ax * ax + ay * ay; + } + if var < 1e-12 { + return None; + } + let (c, s) = (sxx / var, sxy / var); + let tx = cb.x - (c * ca.x - s * ca.y); + let ty = cb.y - (s * ca.x + c * ca.y); + Some(Homography([c, -s, tx, s, c, ty, 0.0, 0.0, 1.0])) + } + Model::Homography => dlt(&a[..n], &b[..n]), + } +} + +fn centroid(p: &[Point]) -> Point { + let n = p.len().max(1) as f64; + Point::new(p.iter().map(|q| q.x).sum::() / n, p.iter().map(|q| q.y).sum::() / n) +} + +/// Hartley normalisation: translate to the centroid, scale to mean distance √2. +fn normaliser(p: &[Point]) -> [f64; 3] { + let c = centroid(p); + let md = p.iter().map(|q| ((q.x - c.x).powi(2) + (q.y - c.y).powi(2)).sqrt()).sum::() / p.len().max(1) as f64; + let s = if md > 1e-12 { std::f64::consts::SQRT_2 / md } else { 1.0 }; + [s, -s * c.x, -s * c.y] +} + +/// Normalised direct linear transform (≥ 4 points). +pub fn dlt(a: &[Point], b: &[Point]) -> Option { + let n = a.len(); + if n < 4 { + return None; + } + let ta = normaliser(a); + let tb = normaliser(b); + let mut ata = [0.0f64; 81]; + for i in 0..n { + let (x, y) = (a[i].x * ta[0] + ta[1], a[i].y * ta[0] + ta[2]); + let (u, v) = (b[i].x * tb[0] + tb[1], b[i].y * tb[0] + tb[2]); + let r1 = [-x, -y, -1.0, 0.0, 0.0, 0.0, u * x, u * y, u]; + let r2 = [0.0, 0.0, 0.0, -x, -y, -1.0, v * x, v * y, v]; + for r in [r1, r2] { + for p in 0..9 { + if r[p] == 0.0 { + continue; + } + for q in 0..9 { + ata[p * 9 + q] += r[p] * r[q]; + } + } + } + } + let (_, vecs) = sym_eigen(&ata, 9); + let h = &vecs[0]; + let hn = Homography([h[0], h[1], h[2], h[3], h[4], h[5], h[6], h[7], h[8]]); + // denormalise: H = Tb⁻¹ · Hn · Ta + let ta_m = Homography([ta[0], 0.0, ta[1], 0.0, ta[0], ta[2], 0.0, 0.0, 1.0]); + let tb_inv = Homography([1.0 / tb[0], 0.0, -tb[1] / tb[0], 0.0, 1.0 / tb[0], -tb[2] / tb[0], 0.0, 0.0, 1.0]); + let h = tb_inv.mul(&hn).mul(&ta_m); + let s = h.0[8]; + if s.abs() < 1e-14 || h.0.iter().any(|v| !v.is_finite()) { + return None; + } + Some(Homography(h.0.map(|v| v / s))) +} + +pub fn transfer_error2(h: &Homography, a: Point, b: Point) -> f64 { + let p = h.apply(a); + (p.x - b.x).powi(2) + (p.y - b.y).powi(2) +} + +/// RANSAC fit of `model` mapping `a` → `b`, inlier threshold `thresh` pixels. Refits on the +/// inliers (twice) and returns `None` when fewer than `min_inliers` agree. +pub fn ransac(model: Model, a: &[Point], b: &[Point], thresh: f64, min_inliers: usize, seed: u64) -> Option { + let n = a.len().min(b.len()); + let k = model.min_samples(); + if n < k.max(min_inliers) { + return None; + } + let t2 = thresh * thresh; + let mut rng = Rng::new(seed ^ 0x9e37_79b9_7f4a_7c15); + let mut best: Vec = Vec::new(); + let mut iters = 2000usize; + let mut it = 0; + let mut sa = Vec::with_capacity(k); + let mut sb = Vec::with_capacity(k); + while it < iters { + it += 1; + sa.clear(); + sb.clear(); + let mut idx = [usize::MAX; 4]; + let mut ok = true; + for s in 0..k { + let mut tries = 0; + let i = loop { + let i = rng.below(n); + if !idx[..s].contains(&i) { + break i; + } + tries += 1; + if tries > 20 { + ok = false; + break i; + } + }; + idx[s] = i; + sa.push(a[i]); + sb.push(b[i]); + } + if !ok { + continue; + } + let h = match model { + Model::Homography => { + let src = [sa[0], sa[1], sa[2], sa[3]]; + let dst = [sb[0], sb[1], sb[2], sb[3]]; + if degenerate(&src) || degenerate(&dst) { + continue; + } + match Homography::from_quads(&src, &dst) { + Some(h) => h, + None => continue, + } + } + _ => match fit_model(model, &sa, &sb) { + Some(h) => h, + None => continue, + }, + }; + let inl: Vec = (0..n).filter(|&i| transfer_error2(&h, a[i], b[i]) < t2).collect(); + if inl.len() > best.len() { + best = inl; + let w = best.len() as f64 / n as f64; + let p_fail = 1.0 - w.powi(k as i32); + if p_fail <= 1e-12 { + iters = it; + } else { + let need = ((1e-3f64).ln() / p_fail.ln()).ceil(); + iters = iters.min(need.max(16.0) as usize); + } + } + } + if best.len() < min_inliers.max(k) { + return None; + } + let mut h = None; + for _ in 0..3 { + let ia: Vec = best.iter().map(|&i| a[i]).collect(); + let ib: Vec = best.iter().map(|&i| b[i]).collect(); + let Some(hh) = fit_model(model, &ia, &ib) else { break }; + let inl: Vec = (0..n).filter(|&i| transfer_error2(&hh, a[i], b[i]) < t2).collect(); + h = Some(hh); + if inl.len() < best.len() || inl == best { + break; + } + best = inl; + } + let h = h?; + if best.len() < min_inliers.max(k) { + return None; + } + let rms = (best.iter().map(|&i| transfer_error2(&h, a[i], b[i])).sum::() / best.len() as f64).sqrt(); + Some(Fit { h, inliers: best, rms }) +} + +/// Three (nearly) collinear points among the four. +fn degenerate(p: &[Point; 4]) -> bool { + for (i, j, k) in [(0, 1, 2), (0, 1, 3), (0, 2, 3), (1, 2, 3)] { + let area = ((p[j].x - p[i].x) * (p[k].y - p[i].y) - (p[j].y - p[i].y) * (p[k].x - p[i].x)).abs(); + let scale = p[i].dist(p[j]).max(p[i].dist(p[k])).max(1e-9); + if area < 1e-3 * scale * scale { + return true; + } + } + false +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn ransac_recovers_a_homography_with_outliers() { + let truth = Homography([1.02, 0.03, 12.0, -0.02, 0.99, -7.0, 1e-5, -2e-5, 1.0]); + let mut rng = Rng::new(5); + let mut a = Vec::new(); + let mut b = Vec::new(); + for i in 0..200 { + let p = Point::new(rng.unit() * 800.0, rng.unit() * 600.0); + let q = if i % 4 == 0 { Point::new(rng.unit() * 800.0, rng.unit() * 600.0) } else { truth.apply(p) }; + let q = Point::new(q.x + (rng.unit() - 0.5) * 0.5, q.y + (rng.unit() - 0.5) * 0.5); + a.push(p); + b.push(q); + } + let f = ransac(Model::Homography, &a, &b, 2.0, 20, 1).unwrap(); + assert!(f.inliers.len() >= 145, "{}", f.inliers.len()); + assert!(f.rms < 0.5); + for p in [Point::new(0.0, 0.0), Point::new(800.0, 600.0), Point::new(400.0, 100.0)] { + assert!(f.h.apply(p).dist(truth.apply(p)) < 0.5); + } + let t = ransac(Model::Translation, &a, &a.iter().map(|p| Point::new(p.x + 3.0, p.y - 1.0)).collect::>(), 1.0, 10, 2).unwrap(); + assert!((t.h.0[2] - 3.0).abs() < 1e-9 && (t.h.0[5] + 1.0).abs() < 1e-9); + let s = fit_model(Model::Similarity, &a[..10], &a[..10].iter().map(|p| Point::new(-p.y * 2.0 + 5.0, p.x * 2.0)).collect::>()).unwrap(); + assert!((s.0[0]).abs() < 1e-9 && (s.0[1] + 2.0).abs() < 1e-9 && (s.0[3] - 2.0).abs() < 1e-9); + } +} diff --git a/crates/merge/tests/hdr.rs b/crates/merge/tests/hdr.rs new file mode 100644 index 00000000..68b9ba4b --- /dev/null +++ b/crates/merge/tests/hdr.rs @@ -0,0 +1,163 @@ +//! HDR merge on synthetic brackets: a procedural scene (scene-linear truth with real highlight +//! headroom) is "photographed" at several exposures with clipping, sensor noise, 14-bit +//! quantisation and small handheld shifts; the merge must recover the radiance. + +use lightcraft_merge::frame::Frame; +use lightcraft_merge::hdr::{Deghost, HdrOptions, merge_hdr}; +use lightcraft_merge::no_progress; +use lightcraft_raster::Rgb32f; + +const W: usize = 640; +const H: usize = 420; +const M: usize = 16; + +struct Rng(u64); +impl Rng { + fn unit(&mut self) -> f64 { + self.0 ^= self.0 << 13; + self.0 ^= self.0 >> 7; + self.0 ^= self.0 << 17; + (self.0 >> 11) as f64 / (1u64 << 53) as f64 + } + fn gauss(&mut self) -> f32 { + let (a, b) = (self.unit().max(1e-12), self.unit()); + ((-2.0 * a.ln()).sqrt() * (std::f64::consts::TAU * b).cos()) as f32 + } +} + +/// Scene radiance (with margin), scaled so the median is ≈ 0.18 at exposure 1. +fn truth() -> Rgb32f { + let scene = lightcraft_scenes::demo_library().into_iter().find(|s| s.kind == lightcraft_scenes::Kind::OceanSunset).unwrap(); + let img = scene.render(W + 2 * M, H + 2 * M); + let mut l: Vec = img.data.iter().map(|p| p[1]).collect(); + l.sort_by(|a, b| a.total_cmp(b)); + let s = 0.18 / l[l.len() / 2].max(1e-6); + img.map(|p| p.map(|v| v * s)) +} + +/// A disk of radiance `val` centred at (cx, cy), radius r (truth coordinates). +fn with_object(t: &Rgb32f, cx: f32, cy: f32, r: f32, val: [f32; 3]) -> Rgb32f { + Rgb32f::from_fn(t.width, t.height, |x, y| { + let d = ((x as f32 + 0.5 - cx).powi(2) + (y as f32 + 0.5 - cy).powi(2)).sqrt(); + if d < r { val } else { t.get(x, y) } + }) +} + +/// Photograph `t` at exposure `k`, content shifted by (dx, dy). +fn shoot(t: &Rgb32f, k: f32, dx: f32, dy: f32, seed: u64) -> Frame { + let mut rng = Rng(seed * 0x9e37_79b9 + 7); + let mut img = Rgb32f::from_fn(W, H, |x, y| t.sample_bilinear(x as f32 + 0.5 + M as f32 - dx, y as f32 + 0.5 + M as f32 - dy).map(|v| v * k)); + // shot + read noise, clipping at the white level, 14-bit quantisation + for p in img.data.iter_mut() { + for c in 0..3 { + let s = p[c].max(0.0); + let sigma = (1e-4 * s + 4e-8).sqrt(); + p[c] = ((s + sigma * rng.gauss()).clamp(0.0, 1.0) * 16383.0).round() / 16383.0; + } + } + let mut f = Frame::from_linear_rec2020(img, Some((k as f64).log2())); + f.raw = true; + f.clip = 1.0; + f +} + +fn log_err(a: [f32; 3], b: [f32; 3]) -> f32 { + (0..3).map(|c| ((a[c] + 0.003).ln() - (b[c] + 0.003).ln()).abs()).fold(0.0, f32::max) +} + +fn pct(mut v: Vec, p: f32) -> f32 { + v.sort_by(|a, b| a.total_cmp(b)); + v[((v.len() - 1) as f32 * p) as usize] +} + +#[test] +fn bracket_merge_recovers_radiance_and_alignment() { + let t = truth(); + let shifts = [(3.25f32, -2.0f32), (0.0, 0.0), (-4.0, 1.5)]; + let ks = [0.25f32, 1.0, 4.0]; + let frames: Vec = (0..3).map(|i| shoot(&t, ks[i], shifts[i].0, shifts[i].1, i as u64 + 1)).collect(); + let t0 = std::time::Instant::now(); + let r = merge_hdr(frames, &HdrOptions { align: true, deghost: Deghost::None }, &no_progress).unwrap(); + eprintln!("merge {W}x{H}: {:?}", t0.elapsed()); + assert_eq!(r.reference, 1); + // exposure estimates within 0.02 EV of the truth + for (i, k) in ks.iter().enumerate() { + assert!((r.ev[i] - (*k as f64).log2()).abs() < 0.02, "ev {:?}", r.ev); + } + // alignment: the reference → frame mapping is the content shift + for i in [0, 2] { + let a = &r.alignments[i]; + for (x, y) in [(50.0, 50.0), (W as f64 - 50.0, H as f64 - 50.0), (320.0, 200.0)] { + let p = a.h.apply(lightcraft_geom::Point::new(x, y)); + let err = ((p.x - x - shifts[i].0 as f64).powi(2) + (p.y - y - shifts[i].1 as f64).powi(2)).sqrt(); + assert!(err < 0.3, "frame {i} ({}, {} inliers, rms {:.3}): alignment error {err:.3} px", a.model, a.inliers, a.rms); + } + } + // radiance vs truth (reference geometry), away from the borders, inside the bracket's range + let mut errs = Vec::new(); + let mut recovered = Vec::new(); + let mut bias = 0f64; + for y in 12..H - 12 { + for x in 12..W - 12 { + let want = t.get(x + M, y + M); + let m = want[0].max(want[1]).max(want[2]); + if !(0.002..3.5).contains(&m) { + continue; + } + let e = log_err(r.radiance.get(x, y), want); + let g = r.radiance.get(x, y)[1]; + bias += ((g + 0.003).ln() - (want[1] + 0.003).ln()) as f64; + errs.push(e); + if m > 1.2 { + recovered.push(e); + } + } + } + let n_err = errs.len(); + eprintln!("mean signed log error (green): {:.4}", bias / n_err as f64); + let (med, p95) = (pct(errs.clone(), 0.5), pct(errs, 0.95)); + eprintln!("log error: median {med:.4}, p95 {p95:.4}; highlights recovered: {} px", recovered.len()); + // noise-limited (the error is the worst of three channels); no bias + assert!(med < 0.03 && p95 < 0.1, "median {med}, p95 {p95}"); + assert!((bias / n_err as f64).abs() < 0.01); + assert!(recovered.len() > 200, "the scene should have highlights above the reference's clip"); + let rmed = pct(recovered, 0.5); + assert!(rmed < 0.05, "recovered highlights: median log error {rmed}"); +} + +#[test] +fn deghost_removes_a_moving_object() { + let t = truth(); + let ks = [0.25f32, 1.0, 4.0]; + // a dark object at a different place in each frame (truth coordinates) + let pos = [(200.0f32, 260.0f32), (330.0, 260.0), (460.0, 260.0)]; + let frames = |_: ()| -> Vec { + (0..3).map(|i| shoot(&with_object(&t, pos[i].0, pos[i].1, 28.0, [0.01, 0.012, 0.02]), ks[i], 0.0, 0.0, 10 + i as u64)).collect() + }; + let ref_truth = with_object(&t, pos[1].0, pos[1].1, 28.0, [0.01, 0.012, 0.02]); + let err_at = |img: &Rgb32f, cx: f32, cy: f32| -> f32 { + let mut v = Vec::new(); + for y in (cy - 20.0) as usize..(cy + 20.0) as usize { + for x in (cx - 20.0) as usize..(cx + 20.0) as usize { + let (tx, ty) = (x - M, y - M); + v.push(log_err(img.get(tx, ty), ref_truth.get(x, y))); + } + } + v.iter().sum::() / v.len() as f32 + }; + let none = merge_hdr(frames(()), &HdrOptions { align: false, deghost: Deghost::None }, &no_progress).unwrap(); + let high = merge_hdr(frames(()), &HdrOptions { align: false, deghost: Deghost::High }, &no_progress).unwrap(); + let med = merge_hdr(frames(()), &HdrOptions { align: false, deghost: Deghost::Medium }, &no_progress).unwrap(); + for (cx, cy) in [pos[0], pos[2]] { + let (en, eh, em) = (err_at(&none.radiance, cx, cy), err_at(&high.radiance, cx, cy), err_at(&med.radiance, cx, cy)); + eprintln!("ghost at ({cx},{cy}): none {en:.3}, medium {em:.3}, high {eh:.3}"); + assert!(en > 0.15, "without deghosting the object leaves a ghost ({en})"); + assert!(eh < 0.08 && em < 0.12, "deghosting removes it (high {eh}, medium {em})"); + // the overlay marks the ghost + let g = high.ghost.get((cx - M as f32) as usize, (cy - M as f32) as usize); + assert!(g > 0.5, "overlay at the ghost: {g}"); + } + assert!(none.ghost.data.is_empty()); + let mean = high.ghost.data.iter().sum::() / high.ghost.data.len() as f32; + assert!(mean < 0.08, "the overlay is local (mean {mean})"); +} diff --git a/crates/merge/tests/pano.rs b/crates/merge/tests/pano.rs new file mode 100644 index 00000000..402037ca --- /dev/null +++ b/crates/merge/tests/pano.rs @@ -0,0 +1,247 @@ +//! Panorama stitching on synthetic input: a procedural scene is used as a spherical environment +//! (longitude/latitude map) and photographed by pinhole cameras with known rotations and focal +//! length. The stitch must recover the geometry and reproduce the environment without seams. + +use lightcraft_merge::frame::Frame; +use lightcraft_merge::linalg::{self, M3}; +use lightcraft_merge::no_progress; +use lightcraft_merge::pano::{PanoOptions, Projection, stitch}; +use lightcraft_raster::Rgb32f; + +const THETA: f64 = 200.0; // degrees of longitude covered by the environment map +const PHI: f64 = 70.0; + +fn env() -> Rgb32f { + let scene = lightcraft_scenes::demo_library().into_iter().find(|s| s.kind == lightcraft_scenes::Kind::Canyon).unwrap(); + let img = scene.render(2400, 840); + let mut l: Vec = img.data.iter().map(|p| p[1]).collect(); + l.sort_by(|a, b| a.total_cmp(b)); + let s = 0.25 / l[l.len() / 2].max(1e-6); + img.map(|p| p.map(|v| (v * s).min(0.98))) +} + +fn sample_env(e: &Rgb32f, d: [f64; 3]) -> Option<[f32; 3]> { + let n = linalg::norm(d); + let t = d[0].atan2(d[2]).to_degrees(); + let p = (d[1] / n).asin().to_degrees(); + if t.abs() > THETA / 2.0 || p.abs() > PHI / 2.0 { + return None; + } + let x = (t / THETA + 0.5) * e.width as f64; + let y = (p / PHI + 0.5) * e.height as f64; + Some(e.sample_bilinear(x as f32, y as f32)) +} + +fn rot(yaw: f64, pitch: f64, roll: f64) -> M3 { + // camera-to-world = Ry(yaw) · Rx(pitch) · Rz(roll); world→camera is its transpose + let ry = linalg::rodrigues([0.0, yaw.to_radians(), 0.0]); + let rx = linalg::rodrigues([pitch.to_radians(), 0.0, 0.0]); + let rz = linalg::rodrigues([0.0, 0.0, roll.to_radians()]); + linalg::transpose3(&linalg::mul3(&linalg::mul3(&ry, &rx), &rz)) +} + +fn shoot(e: &Rgb32f, r: &M3, f: f64, w: usize, h: usize, gain: f32) -> Frame { + let rt = linalg::transpose3(r); + let img = Rgb32f::from_fn(w, h, |x, y| { + let d = [(x as f64 + 0.5 - w as f64 / 2.0) / f, (y as f64 + 0.5 - h as f64 / 2.0) / f, 1.0]; + sample_env(e, linalg::apply3(&rt, d)).map(|p| p.map(|v| v * gain)).unwrap_or([0.0; 3]) + }); + Frame::from_linear_rec2020(img, None) +} + +fn angle_between(a: &M3, b: &M3) -> f64 { + let d = linalg::mul3(a, &linalg::transpose3(b)); + linalg::norm(linalg::rotation_vector(&d)).to_degrees() +} + +struct Shot { + yaw: f64, + pitch: f64, + roll: f64, + gain: f32, +} + +fn run(shots: &[Shot], f: f64, opts: &PanoOptions) -> (lightcraft_merge::pano::PanoResult, Vec, Rgb32f) { + let e = env(); + let rs: Vec = shots.iter().map(|s| rot(s.yaw, s.pitch, s.roll)).collect(); + let frames: Vec = shots.iter().zip(&rs).map(|(s, r)| shoot(&e, r, f, 800, 600, s.gain)).collect(); + let t0 = std::time::Instant::now(); + let res = stitch(frames, opts, &no_progress).unwrap(); + eprintln!( + "stitch {} frames → {}×{} ({:?}) in {:?}; BA rms {:.3} → {:.3} px; fov {:.0}°×{:.0}°", + shots.len(), + res.image.width, + res.image.height, + res.projection, + t0.elapsed(), + res.rms_before, + res.rms, + res.fov.0, + res.fov.1 + ); + (res, rs, e) +} + +/// Check the recovered geometry and compare the panorama with the environment rendered through the +/// recovered projection; returns (median, p99) of the log error over covered pixels. +fn check(res: &lightcraft_merge::pano::PanoResult, rs: &[M3], e: &Rgb32f, f: f64) -> (f32, f32, f32) { + let cams: Vec<_> = res.cameras.iter().map(|c| c.expect("all frames used")).collect(); + for c in &cams { + assert!((c.f - f).abs() / f < 0.01, "focal {} vs {f}", c.f); + } + for i in 0..cams.len() { + for j in i + 1..cams.len() { + let truth = linalg::mul3(&rs[i], &linalg::transpose3(&rs[j])); + let got = linalg::mul3(&cams[i].r, &linalg::transpose3(&cams[j].r)); + let err = angle_between(&truth, &got); + assert!(err < 0.1, "relative rotation {i}-{j} off by {err:.3}°"); + } + } + // world alignment: recovered world → environment world + let mut acc = [[0.0; 3]; 3]; + for (c, r) in cams.iter().zip(rs) { + let m = linalg::mul3(&linalg::transpose3(r), &c.r); + for a in 0..3 { + for b in 0..3 { + acc[a][b] += m[a][b]; + } + } + } + let wmap = linalg::orthonormalize(&acc); + let (w, h) = (res.image.width, res.image.height); + let scale = res.focal; + let (bx, by) = res.origin; + // the environment rendered through the recovered geometry (alpha 0 where it has no data) + let mut truth = Rgb32f::new(w, h); + let mut talpha = lightcraft_raster::Plane::new(w, h); + for y in 0..h { + for x in 0..w { + let d = res.projection.inverse(x as f64 + 0.5 + bx, y as f64 + 0.5 + by, scale); + if let Some(t) = sample_env(e, linalg::apply3(&wmap, d)) + && res.alpha.get(x, y) >= 1.0 + { + truth.set(x, y, t); + talpha.set(x, y, 1.0); + } + } + } + // global level (exposure compensation fixes relative, not absolute, exposure) + let mut lr: Vec = + (0..w * h).filter(|&k| talpha.data[k] > 0.0).map(|k| (res.image.data[k][1] + 0.01).ln() - (truth.data[k][1] + 0.01).ln()).collect(); + lr.sort_by(|a, b| a.total_cmp(b)); + let offset = lr[lr.len() / 2]; + let mut errs = Vec::new(); + let mut emap = Rgb32f::new(w, h); + for k in 0..w * h { + if talpha.data[k] == 0.0 { + continue; + } + let (g, t) = (res.image.data[k], truth.data[k]); + let er = (0..3).map(|c| ((g[c] + 0.01).ln() - offset - (t[c] + 0.01).ln()).abs()).fold(0.0, f32::max); + errs.push(er); + emap.data[k] = [er * 5.0, 1.0 - er * 5.0, 0.0]; + } + if let Ok(dir) = std::env::var("LC_MERGE_DUMP") { + let tag = format!("{:?}-{}", res.projection, w); + ppm(&format!("{dir}/{tag}.ppm"), &res.image); + ppm(&format!("{dir}/{tag}-err.ppm"), &emap); + } + // seams: the column-to-column change of mean log luminance must follow the truth's + let steps = column_steps(&res.image, &talpha); + let tsteps = column_steps(&truth, &talpha); + // summed over 5 columns: resampling legitimately spreads a sharp scene edge over neighbours + let win = |v: &[Option], x: usize| -> Option { (x..x + 5).map(|k| v.get(k).copied().flatten()).sum() }; + let seam = (0..steps.len().saturating_sub(5)).filter_map(|x| Some((win(&steps, x)? - win(&tsteps, x)?).abs())).fold(0.0, f32::max); + errs.sort_by(|a, b| a.total_cmp(b)); + (errs[errs.len() / 2], errs[errs.len() * 99 / 100], seam) +} + +fn ppm(path: &str, img: &Rgb32f) { + let mut b = format!("P6 {} {} 255\n", img.width, img.height).into_bytes(); + for p in &img.data { + for c in 0..3 { + b.push((p[c].clamp(0.0, 1.0).powf(1.0 / 2.2) * 255.0) as u8); + } + } + std::fs::write(path, b).unwrap(); +} + +/// Per column: mean change of log luminance to the next column (where both are covered). +fn column_steps(img: &Rgb32f, alpha: &lightcraft_raster::Plane) -> Vec> { + let (w, h) = (img.width, img.height); + (0..w - 1) + .map(|x| { + let mut s = 0.0; + let mut n = 0; + for y in 0..h { + if alpha.get(x, y) >= 1.0 && alpha.get(x + 1, y) >= 1.0 { + s += (img.get(x + 1, y)[1] + 0.01).ln() - (img.get(x, y)[1] + 0.01).ln(); + n += 1; + } + } + (n > h / 3).then(|| s / n as f32) + }) + .collect() +} + +#[test] +fn spherical_panorama_recovers_geometry_without_seams() { + let shots = [ + Shot { yaw: -60.0, pitch: 1.5, roll: 0.8, gain: 1.0 }, + Shot { yaw: -30.0, pitch: -1.0, roll: -0.5, gain: 1.0 }, + Shot { yaw: 0.0, pitch: 0.5, roll: 1.0, gain: 1.0 }, + Shot { yaw: 31.0, pitch: -1.5, roll: 0.0, gain: 1.3 }, // auto exposure drifted + Shot { yaw: 61.0, pitch: 1.0, roll: -1.0, gain: 1.0 }, + ]; + let f = 700.0; + let opts = PanoOptions { projection: Projection::Spherical, max_megapixels: 100.0, ..Default::default() }; + let (res, rs, e) = run(&shots, f, &opts); + assert_eq!(res.used.len(), 5); + assert!(res.rms < 0.5, "BA rms {}", res.rms); + assert!((res.fov.0 - 182.0).abs() < 8.0, "fov {:?}", res.fov); + // exposure compensation undid the 1.3× drift (relative to its neighbours) + let g = &res.gains; + let rel = g[3] / ((g[2] + g[4]) / 2.0); + assert!((rel - 1.0 / 1.3).abs() < 0.03, "gains {g:?}"); + let (med, p99, seam) = check(&res, &rs, &e, f); + eprintln!("log error vs environment: median {med:.4}, p99 {p99:.4}, seam step {seam:.4}"); + assert!(med < 0.02 && p99 < 0.15, "median {med}, p99 {p99}"); + // (an unblended 1.3× exposure seam would be a step of ln 1.3 ≈ 0.26) + assert!(seam < 0.05, "a column step of {seam} beyond the scene's suggests a seam"); +} + +#[test] +fn auto_projection_and_finishing() { + let shots = [ + Shot { yaw: -22.0, pitch: 0.5, roll: 0.3, gain: 1.0 }, + Shot { yaw: 0.0, pitch: -0.4, roll: 0.0, gain: 1.0 }, + Shot { yaw: 21.0, pitch: 0.8, roll: -0.4, gain: 1.0 }, + ]; + let f = 800.0; + let (res, rs, e) = run(&shots, f, &PanoOptions { projection: Projection::Auto, auto_crop: true, max_megapixels: 100.0, ..Default::default() }); + assert_eq!(res.projection, Projection::Perspective, "a ~90° set is shown rectilinear"); + let (med, p99, seam) = check(&res, &rs, &e, f); + eprintln!("perspective: median {med:.4}, p99 {p99:.4}, seam step {seam:.4}"); + assert!(med < 0.02 && p99 < 0.1 && seam < 0.05); + let c = res.crop.expect("auto crop"); + assert!(c[2] * c[3] > 0.6, "crop {c:?}"); + let (w, h) = (res.image.width as f64, res.image.height as f64); + for y in (c[1] * h) as usize..((c[1] + c[3]) * h) as usize { + for x in (c[0] * w) as usize..((c[0] + c[2]) * w) as usize { + assert!(res.alpha.get(x, y) >= 0.999, "crop must be inside the image data"); + } + } + // boundary warp fills the canvas; cylindrical also works + let (warped, _, _) = + run(&shots, f, &PanoOptions { projection: Projection::Cylindrical, boundary_warp: 100.0, fill_edges: true, ..Default::default() }); + let cover = warped.alpha.data.iter().filter(|a| **a >= 0.5).count() as f64 / warped.alpha.data.len() as f64; + assert!(cover > 0.97, "boundary warp coverage {cover}"); + assert!(warped.image.data.iter().all(|p| p.iter().all(|v| v.is_finite()))); +} + +#[test] +fn unrelated_photos_do_not_stitch() { + let a = lightcraft_scenes::demo_library(); + let frames: Vec = a.iter().take(2).map(|s| Frame::from_linear_rec2020(s.render(400, 300), None)).collect(); + assert!(stitch(frames, &PanoOptions::default(), &no_progress).is_err()); +} diff --git a/crates/raw/src/dngwrite.rs b/crates/raw/src/dngwrite.rs index 6bdeeb27..10bb8ded 100644 --- a/crates/raw/src/dngwrite.rs +++ b/crates/raw/src/dngwrite.rs @@ -16,6 +16,61 @@ pub enum DngCompression { Lj92 { tile: u32, }, + /// Deflate (zlib) tiles of `tile × tile` pixels. Floating-point data is stored as 16-bit (`half`) or + /// 32-bit IEEE floats with the floating-point predictor (TIFF predictor 3, Adobe TIFF Technical Note 3); + /// integer data with horizontal differencing (predictor 2). 16-bit float is what HDR merges usually use. + Deflate { + tile: u32, + half: bool, + }, +} + +/// IEEE 754 binary32 → binary16 (round to nearest even; overflow → ±inf, NaN kept). +pub fn f32_to_f16(v: f32) -> u16 { + let b = v.to_bits(); + let sign = ((b >> 16) & 0x8000) as u16; + let exp = ((b >> 23) & 0xff) as i32; + let man = b & 0x7f_ffff; + if exp == 0xff { + return sign | 0x7c00 | if man != 0 { 0x200 } else { 0 }; + } + let e = exp - 127 + 15; + if e >= 0x1f { + return sign | 0x7c00; + } + if e <= 0 { + if e < -10 { + return sign; + } + // subnormal half + let m = man | 0x80_0000; + let shift = (14 - e) as u32; + let half = m >> shift; + let rem = m & ((1 << shift) - 1); + let mid = 1 << (shift - 1); + let r = if rem > mid || (rem == mid && half & 1 == 1) { half + 1 } else { half }; + return sign | r as u16; + } + let half = ((e as u32) << 10) | (man >> 13); + let rem = man & 0x1fff; + let r = if rem > 0x1000 || (rem == 0x1000 && half & 1 == 1) { half + 1 } else { half }; + sign | r as u16 +} + +/// Apply the floating-point predictor (inverse of `unpack::undo_float_predictor`) to one row of +/// `n` big-endian samples of `bytes_per` bytes: split into byte planes (most significant first), then +/// byte-difference with `stride` samples. +fn float_predict(row: &[u8], n: usize, bytes_per: usize, stride: usize) -> Vec { + let mut planes = vec![0u8; row.len()]; + for i in 0..n { + for b in 0..bytes_per { + planes[b * n + i] = row[i * bytes_per + b]; + } + } + for i in (stride..planes.len()).rev() { + planes[i] = planes[i].wrapping_sub(planes[i - stride]); + } + planes } #[derive(Clone, Debug, PartialEq)] @@ -131,6 +186,42 @@ pub fn write_dng(raw: &RawImage, opts: &DngWriteOptions) -> Result> { } match &raw.data { + RawData::F32(v) if matches!(opts.compression, DngCompression::Deflate { .. }) => { + let DngCompression::Deflate { tile, half } = opts.compression else { unreachable!() }; + let tile = (tile.max(16) & !15) as usize; + let bp = if half { 2 } else { 4 }; + ifd.set(t::BITS_PER_SAMPLE, Value::Short(vec![(bp * 8) as u16; cpp])); + ifd.set(t::SAMPLE_FORMAT, Value::Short(vec![3; cpp])); + ifd.set(t::COMPRESSION, Value::Short(vec![compression::ADOBE_DEFLATE])); + ifd.set(t::PREDICTOR, Value::Short(vec![3])); + let (ta, td) = (w.div_ceil(tile), h.div_ceil(tile)); + let tiles: Vec> = (0..ta * td) + .into_par_iter() + .map(|i| { + let (tx, ty) = (i % ta * tile, i / ta * tile); + let mut buf = Vec::with_capacity(tile * tile * cpp * bp); + let mut row = Vec::with_capacity(tile * cpp * bp); + for y in 0..tile { + let sy = (ty + y).min(h - 1); + row.clear(); + for x in 0..tile { + let sx = (tx + x).min(w - 1); + for s in 0..cpp { + let f = v[(sy * w + sx) * cpp + s]; + if half { + row.extend_from_slice(&f32_to_f16(f).to_be_bytes()); + } else { + row.extend_from_slice(&f.to_bits().to_be_bytes()); + } + } + } + buf.extend_from_slice(&float_predict(&row, tile * cpp, bp, cpp)); + } + miniz_oxide::deflate::compress_to_vec_zlib(&buf, 6) + }) + .collect(); + ifd.set_image(ImageData::Tiles { tile_width: tile as u32, tile_height: tile as u32, tiles }); + } RawData::F32(v) => { ifd.set(t::BITS_PER_SAMPLE, Value::Short(vec![32; cpp])); ifd.set(t::SAMPLE_FORMAT, Value::Short(vec![3; cpp])); @@ -164,6 +255,36 @@ pub fn write_dng(raw: &RawImage, opts: &DngWriteOptions) -> Result> { .collect(); ifd.set_image(ImageData::Strips { rows_per_strip: rows as u32, strips }); } + DngCompression::Deflate { tile, .. } => { + let tile = (tile.max(16) & !15) as usize; + ifd.set(t::BITS_PER_SAMPLE, Value::Short(vec![16; cpp])); + ifd.set(t::COMPRESSION, Value::Short(vec![compression::ADOBE_DEFLATE])); + ifd.set(t::PREDICTOR, Value::Short(vec![2])); + let (ta, td) = (w.div_ceil(tile), h.div_ceil(tile)); + let tiles: Vec> = (0..ta * td) + .into_par_iter() + .map(|i| { + let (tx, ty) = (i % ta * tile, i / ta * tile); + let mut buf = Vec::with_capacity(tile * tile * cpp * 2); + let mut row = vec![0u16; tile * cpp]; + for y in 0..tile { + let sy = (ty + y).min(h - 1); + for x in 0..tile { + let sx = (tx + x).min(w - 1); + for s in 0..cpp { + row[x * cpp + s] = v[(sy * w + sx) * cpp + s]; + } + } + for i in (cpp..row.len()).rev() { + row[i] = row[i].wrapping_sub(row[i - cpp]); + } + row.iter().for_each(|&x| opts.order.put_u16(&mut buf, x)); + } + miniz_oxide::deflate::compress_to_vec_zlib(&buf, 6) + }) + .collect(); + ifd.set_image(ImageData::Tiles { tile_width: tile as u32, tile_height: tile as u32, tiles }); + } DngCompression::Lj92 { tile } => { let tile = (tile.max(2) & !1) as usize; if tile > 65534 { diff --git a/crates/raw/tests/dng.rs b/crates/raw/tests/dng.rs index 6fa7b670..04f59507 100644 --- a/crates/raw/tests/dng.rs +++ b/crates/raw/tests/dng.rs @@ -106,6 +106,40 @@ fn writer_roundtrip_all_layouts() { } } +#[test] +fn deflate_roundtrips_integer_and_float() { + for (w, h) in [(64, 48), (70, 45)] { + for cfa in [Some(Cfa::bayer("RGGB").unwrap()), None] { + let cpp = if cfa.is_some() { 1 } else { 3 }; + let raw = synthetic(w, h, cfa, cpp); + for order in [ByteOrder::Little, ByteOrder::Big] { + let bytes = + write_dng(&raw, &DngWriteOptions { compression: DngCompression::Deflate { tile: 32, half: false }, order, xmp: None }).unwrap(); + assert_same(&raw, &decode(&bytes).unwrap()); + } + } + } + // linear float data: 32-bit exact, 16-bit within half precision + let mut raw = synthetic(70, 45, None, 3); + let vals: Vec = (0..70 * 45 * 3).map(|i| ((i % 97) as f32 / 97.0).powi(3) * 0.9 + 1e-4).collect(); + raw.data = RawData::F32(vals.clone()); + raw.black = BlackLevel::uniform(0.0); + raw.white = vec![1.0]; + for half in [false, true] { + let bytes = write_dng(&raw, &DngWriteOptions { compression: DngCompression::Deflate { tile: 32, half }, ..Default::default() }).unwrap(); + let back = decode(&bytes).unwrap(); + let RawData::F32(b) = &back.data else { panic!("float data expected") }; + for (a, b) in vals.iter().zip(b) { + let tol = if half { a * 1e-3 + 1e-7 } else { 0.0 }; + assert!((a - b).abs() <= tol, "{a} vs {b} (half: {half})"); + } + } + assert_eq!(lightcraft_raw::dngwrite::f32_to_f16(1.0), 0x3c00); + assert_eq!(lightcraft_raw::dngwrite::f32_to_f16(-2.0), 0xc000); + assert_eq!(lightcraft_raw::dngwrite::f32_to_f16(1e9), 0x7c00); + assert_eq!(lightcraft_raw::dngwrite::f32_to_f16(2f32.powi(-24)), 0x0001); +} + #[test] fn normalized_and_develop() { let mut raw = synthetic(40, 30, Some(Cfa::bayer("RGGB").unwrap()), 1); From bc2cda97cfd5b07adc7bcf1ec3f13f0231a5efbf Mon Sep 17 00:00:00 2001 From: Brandon Thomas Date: Thu, 1 Oct 2026 10:30:03 -0400 Subject: [PATCH 2/3] M13.2: Photo Merge commands, dialog with live preview, CLI - engine: `merge.hdr`, `merge.panorama`, `merge.hdrPanorama` (bracket size from EXIF); MergeJob (Send) with progress/cancel, preview renders through the develop pipeline (deghost overlay, auto-crop frame); finish writes -HDR/-Pano.dng next to the first source, imports, selects, applies Auto Settings / auto crop - ui: Photo Merge dialog (HDR / Panorama / HDR Panorama) with live preview on a worker thread, background final merge with progress toast; context menu and Ctrl+H / Ctrl+M (Ctrl now means the physical Control key on macOS) - cli: `merge hdr|panorama|hdr-panorama`, `synth-merge` for synthetic inputs - merge::synth: procedural brackets (DNG) and panorama views; docs/merge.md Co-Authored-By: Claude Opus 5.5 --- Cargo.lock | 3 + apps/lightcraft-cli/Cargo.toml | 1 + apps/lightcraft-cli/src/main.rs | 119 ++++++++ crates/engine/Cargo.toml | 2 + crates/engine/src/cmd/merge.rs | 72 +++++ crates/engine/src/cmd/mod.rs | 2 + crates/engine/src/lib.rs | 3 + crates/engine/src/media.rs | 3 + crates/engine/src/merge.rs | 408 +++++++++++++++++++++++++++ crates/engine/src/tests_merge.rs | 82 ++++++ crates/merge/Cargo.toml | 4 +- crates/merge/src/lib.rs | 1 + crates/merge/src/synth.rs | 81 ++++++ crates/ui-egui/src/control.rs | 2 +- crates/ui-egui/src/headless.rs | 2 +- crates/ui-egui/src/lib.rs | 15 +- crates/ui-egui/src/menus.rs | 7 + crates/ui-egui/src/merge.rs | 389 +++++++++++++++++++++++++ crates/ui-egui/src/panels/dialogs.rs | 10 +- crates/ui-egui/src/panels/grid.rs | 8 + crates/ui-egui/src/shortcuts.rs | 5 +- crates/ui-egui/src/state.rs | 4 + docs/merge.md | 63 +++++ 23 files changed, 1277 insertions(+), 9 deletions(-) create mode 100644 crates/engine/src/cmd/merge.rs create mode 100644 crates/engine/src/merge.rs create mode 100644 crates/engine/src/tests_merge.rs create mode 100644 crates/merge/src/synth.rs create mode 100644 crates/ui-egui/src/merge.rs create mode 100644 docs/merge.md diff --git a/Cargo.lock b/Cargo.lock index 9a384cca..e2679dc2 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -2035,6 +2035,7 @@ dependencies = [ "lightcraft-codecs", "lightcraft-engine", "lightcraft-mcp", + "lightcraft-merge", "lightcraft-raster", "lightcraft-ui-egui", "serde_json", @@ -2097,6 +2098,7 @@ dependencies = [ "lightcraft-develop", "lightcraft-geom", "lightcraft-gpu", + "lightcraft-merge", "lightcraft-meta", "lightcraft-pipeline", "lightcraft-preview", @@ -2104,6 +2106,7 @@ dependencies = [ "lightcraft-raw", "lightcraft-scenes", "log", + "rayon", "serde", "serde_json", "thiserror 2.0.21", diff --git a/apps/lightcraft-cli/Cargo.toml b/apps/lightcraft-cli/Cargo.toml index 939bdefe..8bb6e114 100644 --- a/apps/lightcraft-cli/Cargo.toml +++ b/apps/lightcraft-cli/Cargo.toml @@ -19,4 +19,5 @@ lightcraft-mcp = { workspace = true } lightcraft-ui-egui = { workspace = true } lightcraft-codecs = { workspace = true } lightcraft-raster = { workspace = true } +lightcraft-merge = { workspace = true } serde_json = { workspace = true } diff --git a/apps/lightcraft-cli/src/main.rs b/apps/lightcraft-cli/src/main.rs index 46186fb4..0df0a917 100644 --- a/apps/lightcraft-cli/src/main.rs +++ b/apps/lightcraft-cli/src/main.rs @@ -4,6 +4,8 @@ //! lightcraft-cli mcp [--connect [ADDR]] [--demo] [--compact] [FILES/FOLDERS…] //! lightcraft-cli render -o [--set control=value]… [--settings FILE.json] [--preset ID] [--size N] [--quality Q] //! lightcraft-cli snapshot [--library DIR | --demo] [--script FILE.jsonl] [-o OUT.png] [--size WxH] [--scale S] [FILES…] +//! lightcraft-cli merge hdr|panorama|hdr-panorama [OPTIONS] FILES… +//! lightcraft-cli synth-merge hdr|panorama -o DIR //! lightcraft-cli commands [--json] //! lightcraft-cli controls [--json] //! ``` @@ -46,6 +48,15 @@ USAGE: -o, --output OUT PNG path (a final screenshot is written here if the script took none) --size WxH window size in points (default 1600x1000) --scale S pixels per point (default 1) + lightcraft-cli merge hdr|panorama|hdr-panorama [OPTIONS] FILES… + Photo Merge: writes -HDR.dng / -Pano.dng / -HDR-Pano.dng next to the first file and + prints the result as JSON. Options: + --deghost none|low|medium|high --no-align --bracket N (HDR panorama) + --projection auto|spherical|cylindrical|perspective --boundary-warp 0..100 + --auto-crop --fill-edges --no-auto-settings + --preview OUT.png only render a ≤ 1024 px preview (nothing written next to the files) + lightcraft-cli synth-merge hdr|panorama -o DIR + Write synthetic merge inputs (procedural scene; bracketed DNGs or overlapping PNG views). lightcraft-cli commands [--json] list every command id with its parameters lightcraft-cli controls [--json] list every develop control id with its range lightcraft-cli --version | --help @@ -58,6 +69,8 @@ fn main() -> ExitCode { Some("render") => render(&args[1..]), Some("snapshot") => snapshot(&args[1..]), Some("commands") => commands(&args[1..]), + Some("merge") => merge(&args[1..]), + Some("synth-merge") => synth_merge(&args[1..]), Some("controls") => controls(&args[1..]), Some("--version" | "-V" | "version") => { println!("lightcraft-cli {}", env!("CARGO_PKG_VERSION")); @@ -148,6 +161,112 @@ fn mcp(args: &[String]) -> Result<(), String> { server.serve(BufReader::new(stdin.lock()), std::io::stdout().lock()).map_err(|e| e.to_string()) } +fn merge(args: &[String]) -> Result<(), String> { + let cmd = match args.first().map(String::as_str) { + Some("hdr") => "merge.hdr", + Some("panorama" | "pano") => "merge.panorama", + Some("hdr-panorama") => "merge.hdrPanorama", + _ => return Err("merge: expected hdr, panorama or hdr-panorama".into()), + }; + let mut p = serde_json::Map::new(); + let mut files = Vec::new(); + let mut i = 1; + while i < args.len() { + match args[i].as_str() { + "--deghost" => { + p.insert("deghost".into(), json!(take_value(args, &mut i, "--deghost")?)); + } + "--no-align" => { + p.insert("align".into(), json!(false)); + } + "--bracket" => { + p.insert("bracket".into(), json!(take_value(args, &mut i, "--bracket")?.parse::().map_err(|e| e.to_string())?)); + } + "--projection" => { + p.insert("projection".into(), json!(take_value(args, &mut i, "--projection")?)); + } + "--boundary-warp" => { + p.insert("boundaryWarp".into(), json!(take_value(args, &mut i, "--boundary-warp")?.parse::().map_err(|e| e.to_string())?)); + } + "--auto-crop" => { + p.insert("autoCrop".into(), json!(true)); + } + "--fill-edges" => { + p.insert("fillEdges".into(), json!(true)); + } + "--no-auto-settings" => { + p.insert("autoSettings".into(), json!(false)); + } + "--preview" => { + p.insert("preview".into(), json!(true)); + p.insert("previewPath".into(), json!(take_value(args, &mut i, "--preview")?)); + p.insert("showOverlay".into(), json!(true)); + } + a if a.starts_with("--") => return Err(format!("unknown option `{a}`")), + f => files.push(f.to_string()), + } + i += 1; + } + let mut s = Session::new().with_fs(); + let paths = expand_paths(&files); + let r = s.execute("library.import", &json!({"paths": paths})).map_err(|e| e.to_string())?; + let mut ids: Vec = r["imported"].as_array().map(|a| a.iter().filter_map(Value::as_u64).collect()).unwrap_or_default(); + // keep the command line's order + ids.sort_by_key(|id| { + paths.iter().position(|f| { + s.catalog + .photo(lightcraft_engine::catalog::PhotoId(*id)) + .is_some_and(|ph| matches!(&ph.source, lightcraft_engine::catalog::Source::File { path } if path == f)) + }) + }); + p.insert("ids".into(), json!(ids)); + let t0 = std::time::Instant::now(); + let mut out = s.execute(cmd, &Value::Object(p)).map_err(|e| e.to_string())?; + out["seconds"] = json!(t0.elapsed().as_secs_f64()); + println!("{}", serde_json::to_string_pretty(&out).unwrap_or_default()); + Ok(()) +} + +fn synth_merge(args: &[String]) -> Result<(), String> { + let kind = args.first().map(String::as_str).unwrap_or(""); + let mut dir = None; + let mut i = 1; + while i < args.len() { + match args[i].as_str() { + "-o" | "--output" => dir = Some(take_value(args, &mut i, "-o")?.to_string()), + a => return Err(format!("unknown option `{a}`")), + } + i += 1; + } + let dir = dir.ok_or("synth-merge needs -o DIR")?; + std::fs::create_dir_all(&dir).map_err(|e| e.to_string())?; + let mut written = Vec::new(); + match kind { + "hdr" => { + for (k, b) in lightcraft_merge::synth::bracket_dngs(1800, 1200, &[-2.0, 0.0, 2.0]).into_iter().enumerate() { + let p = Path::new(&dir).join(format!("bracket-{k}.dng")); + std::fs::write(&p, b).map_err(|e| e.to_string())?; + written.push(p); + } + } + "panorama" | "pano" => { + for (k, v) in lightcraft_merge::synth::pano_views(1200, 900, 1000.0, &[-50.0, -25.0, 0.0, 25.0, 50.0]).into_iter().enumerate() { + let p = Path::new(&dir).join(format!("view-{k}.png")); + let img = v.to_srgb8(); + let png = lightcraft_codecs::encode_png(&lightcraft_codecs::EncodeImage::rgba8(&img), &lightcraft_codecs::EncodeMeta::default()) + .map_err(|e| e.to_string())?; + std::fs::write(&p, png).map_err(|e| e.to_string())?; + written.push(p); + } + } + _ => return Err("synth-merge: expected hdr or panorama".into()), + } + for p in written { + println!("{}", p.display()); + } + Ok(()) +} + fn take_value<'a>(args: &'a [String], i: &mut usize, flag: &str) -> Result<&'a str, String> { *i += 1; args.get(*i).map(String::as_str).ok_or_else(|| format!("{flag} needs a value")) diff --git a/crates/engine/Cargo.toml b/crates/engine/Cargo.toml index 59922a08..a514ecf2 100644 --- a/crates/engine/Cargo.toml +++ b/crates/engine/Cargo.toml @@ -21,6 +21,8 @@ lightcraft-preview = { workspace = true } lightcraft-codecs = { workspace = true } lightcraft-raw = { workspace = true } lightcraft-meta = { workspace = true } +lightcraft-merge = { workspace = true } +rayon = { workspace = true } serde = { workspace = true } serde_json = { workspace = true } thiserror = { workspace = true } diff --git a/crates/engine/src/cmd/merge.rs b/crates/engine/src/cmd/merge.rs new file mode 100644 index 00000000..102ff0d3 --- /dev/null +++ b/crates/engine/src/cmd/merge.rs @@ -0,0 +1,72 @@ +//! Photo Merge commands: `merge.hdr`, `merge.panorama`, `merge.hdrPanorama`. +//! +//! Run synchronously here (CLI, MCP, control channel); the desktop UI plans the same job with +//! [`Session::plan_merge`] and runs it on a worker thread. `preview: true` merges ≤ 1024 px frames +//! and (with `previewPath`) writes the rendered preview as PNG instead of adding a photo. + +use lightcraft_catalog::PhotoId; +use serde_json::{Value, json}; + +use super::{CommandSpec, bad, bool_or, cmd, str_param}; +use crate::{Result, Session}; + +fn can_merge(s: &Session) -> std::result::Result<(), String> { + if cfg!(target_arch = "wasm32") && s.media.file_bytes.is_none() { + return Err("Photo Merge needs access to the photo files".into()); + } + // the photo count is checked when run (an explicit `ids` parameter may name them) + if s.catalog.photos().nth(1).is_none() { Err("the library has fewer than 2 photos".into()) } else { Ok(()) } +} + +fn run(s: &mut Session, id: &str, p: &Value) -> Result { + let (kind, finish) = crate::merge::parse(id, p)?; + let ids: Vec = s.targets(p); + if !bool_or(p, "preview", false) { + return s.merge_now(kind, finish, &ids); + } + let job = s.plan_merge(kind, finish, &ids, true)?; + let out = job.run(&|_, _| true).map_err(|e| bad(id, e))?; + let mut info = out.info; + if let (Some(path), Some(img)) = (str_param(p, "previewPath"), &out.preview) { + let png = lightcraft_codecs::encode_png(&lightcraft_codecs::EncodeImage::rgba8(img), &lightcraft_codecs::EncodeMeta::default()) + .map_err(|e| bad(id, e.to_string()))?; + std::fs::write(path, png).map_err(|e| bad(id, format!("{path}: {e}")))?; + info["previewPath"] = json!(path); + } + if let Some(img) = &out.preview { + info["previewSize"] = json!([img.width, img.height]); + } + Ok(info) +} + +pub fn specs() -> Vec { + vec![ + cmd!( + "merge.hdr", + "HDR Merge", + [], + None, + "{ids?, align=true, deghost=none|low|medium|high, autoSettings=true, stack=false, preview=false, showOverlay=false, previewPath?}", + can_merge, + |s, p| run(s, "merge.hdr", p) + ), + cmd!( + "merge.panorama", + "Panorama Merge", + [], + None, + "{ids?, projection=auto|spherical|cylindrical|perspective, boundaryWarp=0..100, autoCrop=false, fillEdges=false, autoSettings=true, stack=false, maxMegapixels=40, preview=false, previewPath?}", + can_merge, + |s, p| run(s, "merge.panorama", p) + ), + cmd!( + "merge.hdrPanorama", + "HDR Panorama Merge", + [], + None, + "{ids?, bracket=0 (auto from EXIF), align, deghost, projection, boundaryWarp, autoCrop, fillEdges, autoSettings, stack, preview, previewPath?}", + can_merge, + |s, p| run(s, "merge.hdrPanorama", p) + ), + ] +} diff --git a/crates/engine/src/cmd/mod.rs b/crates/engine/src/cmd/mod.rs index 3e1b8c54..62331349 100644 --- a/crates/engine/src/cmd/mod.rs +++ b/crates/engine/src/cmd/mod.rs @@ -8,6 +8,7 @@ mod develop; mod edit; mod library; mod masks; +mod merge; mod preset_files; mod query; mod xmp; @@ -102,6 +103,7 @@ pub fn command_specs() -> &'static [CommandSpec] { v.extend(library::specs()); v.extend(develop::specs()); v.extend(masks::specs()); + v.extend(merge::specs()); v.extend(query::specs()); v.extend(xmp::specs()); v.extend(preset_files::specs()); diff --git a/crates/engine/src/lib.rs b/crates/engine/src/lib.rs index 2f971c97..c259d8bf 100644 --- a/crates/engine/src/lib.rs +++ b/crates/engine/src/lib.rs @@ -18,6 +18,7 @@ pub mod files; pub mod import; pub mod library; pub mod media; +pub mod merge; pub mod presets; pub mod sidecar; mod view; @@ -321,4 +322,6 @@ mod tests_import; #[cfg(test)] mod tests_library; #[cfg(test)] +mod tests_merge; +#[cfg(test)] mod tests_xmp; diff --git a/crates/engine/src/media.rs b/crates/engine/src/media.rs index 1c0a4070..f8907827 100644 --- a/crates/engine/src/media.rs +++ b/crates/engine/src/media.rs @@ -92,6 +92,8 @@ pub struct MediaCache { pub preview_capacity: usize, pub file_loader: Option, pub file_probe: Option, + /// Reads a photo file's bytes (Photo Merge); `None` = the local file system. + pub file_bytes: Option, scenes: Vec, /// Rendered thumbnails (memory, plus disk once a library is attached). pub rendered: Arc, @@ -106,6 +108,7 @@ impl Default for MediaCache { preview_capacity: 0, file_loader: None, file_probe: None, + file_bytes: None, scenes: Vec::new(), rendered: Arc::new(PreviewCache::memory(RENDERED_MEM_BYTES)), } diff --git a/crates/engine/src/merge.rs b/crates/engine/src/merge.rs new file mode 100644 index 00000000..608928ed --- /dev/null +++ b/crates/engine/src/merge.rs @@ -0,0 +1,408 @@ +//! Photo Merge (HDR, Panorama, HDR Panorama) in the library. +//! +//! A merge is planned from the session ([`crate::Session::plan_merge`] → [`MergeJob`], which is +//! `Send` and self-contained), run anywhere (a worker thread in the app, inline for the CLI/MCP: +//! [`MergeJob::run`], with progress and cancellation), and finished on the session +//! ([`crate::Session::finish_merge`]): the DNG is written next to the first source photo +//! (`-HDR.dng`, `-Pano.dng`, `-HDR-Pano.dng`), imported, selected, and given Auto Settings / +//! the auto-crop rectangle when asked. A *preview* run works on ≤ 1024 px frames and returns a +//! rendered image (with the deghost overlay when asked) instead of a file. + +use std::sync::Arc; + +use lightcraft_catalog::{PhotoId, Source}; +use lightcraft_develop::DevelopSettings; +use lightcraft_merge::{Deghost, Frame, HdrOptions, PanoOptions, Projection}; +use lightcraft_pipeline::{Quality, RenderRequest, SourceInfo}; +use lightcraft_raster::{Rgb32f, Rgba8}; +use serde::{Deserialize, Serialize}; +use serde_json::{Value, json}; + +use crate::{EngineError, Result, Session}; + +/// Long edge of the frames of a preview merge. +pub const PREVIEW_EDGE: usize = 1024; + +#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)] +#[serde(tag = "kind", rename_all = "camelCase")] +pub enum MergeKind { + Hdr { + hdr: HdrOptions, + }, + Panorama { + pano: PanoOptions, + }, + /// Merge each bracket (`bracket` consecutive photos; 0 = detect from EXIF), then stitch. + HdrPanorama { + hdr: HdrOptions, + pano: PanoOptions, + bracket: usize, + }, +} + +impl MergeKind { + fn suffix(&self) -> &'static str { + match self { + MergeKind::Hdr { .. } => "HDR", + MergeKind::Panorama { .. } => "Pano", + MergeKind::HdrPanorama { .. } => "HDR-Pano", + } + } + fn min_photos(&self) -> usize { + match self { + MergeKind::HdrPanorama { .. } => 4, + _ => 2, + } + } +} + +/// Options that apply after the merge. +#[derive(Clone, Debug, Default, PartialEq, Serialize, Deserialize)] +#[serde(default, rename_all = "camelCase")] +pub struct MergeFinish { + /// Run Auto Settings on the result. + pub auto_settings: bool, + /// Group the result with its sources. LightCraft has no stacks yet: recorded in the result's + /// caption instead (see `docs/merge.md`). + pub stack: bool, + /// Preview only: tint deghosted areas red. + pub show_overlay: bool, +} + +/// Reads a source photo's bytes (by catalog path). +pub type ByteReader = Arc std::result::Result, String> + Send + Sync>; + +/// A planned merge, detached from the session. +#[derive(Clone)] +pub struct MergeJob { + pub kind: MergeKind, + pub finish: MergeFinish, + pub sources: Vec<(PhotoId, String)>, + /// `Some(edge)`: preview at that frame size; `None`: the full-resolution merge. + pub preview: Option, + pub read: ByteReader, +} + +/// What a merge produced. +pub struct MergeOutput { + /// The DNG (also for previews, at preview size). + pub dng: Vec, + /// Rendered preview (previews only). + pub preview: Option, + /// Auto-crop rectangle (normalised x, y, w, h) for panoramas. + pub crop: Option<[f64; 4]>, + /// Statistics for the caller (exposures, alignment, projection, …). + pub info: Value, +} + +pub type ProgressFn<'a> = dyn Fn(f32, &str) -> bool + Sync + 'a; + +fn err(e: impl std::fmt::Display) -> String { + e.to_string() +} + +impl MergeJob { + fn load(&self, orient: bool, progress: &ProgressFn, span: (f32, f32)) -> std::result::Result, String> { + use rayon::prelude::*; + let done = std::sync::atomic::AtomicUsize::new(0); + let n = self.sources.len(); + let frames: std::result::Result, String> = self + .sources + .par_iter() + .map(|(_, path)| { + let bytes = (self.read)(path)?; + let f = lightcraft_merge::load_frame(&bytes, self.preview, orient).map_err(|e| format!("{path}: {e}")); + let k = done.fetch_add(1, std::sync::atomic::Ordering::Relaxed) + 1; + progress(span.0 + (span.1 - span.0) * k as f32 / n as f32, "Reading photos"); + f + }) + .collect(); + frames + } + + /// Run the merge. `progress(fraction, stage)` returns false to cancel. + pub fn run(&self, progress: &ProgressFn) -> std::result::Result { + if self.sources.len() < self.kind.min_photos() { + return Err(format!("select at least {} photos", self.kind.min_photos())); + } + if !progress(0.0, "Reading photos") { + return Err("cancelled".into()); + } + let sub = |a: f32, b: f32| move |f: f32, s: &str| progress(a + (b - a) * f, s); + let (image, color, orientation, metadata, baseline, samples, crop, mut info, ghost) = match &self.kind { + MergeKind::Hdr { hdr } => { + let frames = self.load(false, progress, (0.0, 0.3))?; + let p = sub(0.3, 0.9); + let r = lightcraft_merge::merge_hdr(frames, hdr, &p).map_err(err)?; + let info = json!({ + "reference": self.sources[r.reference].0.0, + "ev": r.ev, "evExif": r.ev_exif, + "alignment": r.alignments.iter().map(|a| json!({"model": a.model, "inliers": a.inliers, "rmsPx": a.rms})).collect::>(), + }); + let ghost = (!r.ghost.data.is_empty()).then(|| r.ghost.oriented(r.orientation)); + (r.radiance, r.color, r.orientation, r.metadata, r.baseline_exposure, lightcraft_merge::DngSamples::Half, None, info, ghost) + } + MergeKind::Panorama { pano } => { + let frames = self.load(true, progress, (0.0, 0.25))?; + let p = sub(0.25, 0.9); + let r = lightcraft_merge::stitch(frames, pano, &p).map_err(err)?; + let info = pano_info(&r, &self.sources); + ( + r.image, + r.color, + lightcraft_geom::Orientation::Normal, + r.metadata, + r.baseline_exposure, + lightcraft_merge::DngSamples::U16, + r.crop, + info, + None, + ) + } + MergeKind::HdrPanorama { hdr, pano, bracket } => { + let frames = self.load(false, progress, (0.0, 0.2))?; + let b = + if *bracket >= 2 { *bracket } else { detect_bracket(&frames).ok_or("could not detect the bracket size from EXIF; set it")? }; + if !frames.len().is_multiple_of(b) || frames.len() / b < 2 { + return Err(format!("{} photos don't form brackets of {b} for at least two positions", frames.len())); + } + let groups = frames.len() / b; + let mut merged = Vec::new(); + let mut it = frames.into_iter(); + for g in 0..groups { + let set: Vec = it.by_ref().take(b).collect(); + let p = sub(0.2 + 0.4 * g as f32 / groups as f32, 0.2 + 0.4 * (g + 1) as f32 / groups as f32); + let r = lightcraft_merge::merge_hdr(set, hdr, &p).map_err(err)?; + let mut f = Frame { + image: r.radiance, + clip: f32::INFINITY, + exposure: None, + raw: r.raw, + color: r.color, + orientation: r.orientation, + metadata: r.metadata, + baseline_exposure: r.baseline_exposure, + }; + f.orient(); + merged.push(f); + } + let p = sub(0.6, 0.9); + let r = lightcraft_merge::stitch(merged, pano, &p).map_err(err)?; + let mut info = pano_info(&r, &self.sources); + info["bracket"] = json!(b); + ( + r.image, + r.color, + lightcraft_geom::Orientation::Normal, + r.metadata, + r.baseline_exposure, + lightcraft_merge::DngSamples::Half, + r.crop, + info, + None, + ) + } + }; + if !progress(0.9, "Writing DNG") { + return Err("cancelled".into()); + } + let mut metadata = metadata; + if self.finish.stack { + let names: Vec<&str> = + self.sources.iter().map(|(_, p)| std::path::Path::new(p).file_name().and_then(|n| n.to_str()).unwrap_or(p)).collect(); + metadata.caption = Some(format!("Merged from {}", names.join(", "))); + } + metadata.software = Some("LightCraft Photo Merge".into()); + let dng = lightcraft_merge::write_linear_dng(&image, &color, orientation, &metadata, baseline, samples).map_err(err)?; + info["width"] = json!(image.width); + info["height"] = json!(image.height); + let preview = match self.preview { + Some(edge) => Some(render_preview(&dng, edge, self.finish.auto_settings, crop, ghost.as_ref().filter(|_| self.finish.show_overlay))?), + None => None, + }; + progress(1.0, "Done"); + Ok(MergeOutput { dng, preview, crop, info }) + } +} + +fn pano_info(r: &lightcraft_merge::PanoResult, sources: &[(PhotoId, String)]) -> Value { + json!({ + "projection": r.projection, + "used": r.used.iter().filter_map(|&i| sources.get(i).map(|s| s.0.0)).collect::>(), + "fov": [r.fov.0, r.fov.1], + "rmsPx": r.rms, + "gains": r.gains, + "crop": r.crop, + }) +} + +/// Bracket size from EXIF: the smallest period of the exposure sequence (≥ 2). +pub fn detect_bracket(frames: &[Frame]) -> Option { + let ev: Vec = frames.iter().map(|f| f.exposure).collect::>()?; + let n = ev.len(); + (2..=n / 2).find(|&b| n.is_multiple_of(b) && (0..n).all(|k| (ev[k] - ev[k % b]).abs() < 0.2) && (1..b).any(|k| (ev[k] - ev[0]).abs() > 0.5)) +} + +/// Render a merged DNG like the library will show it (default raw settings, optionally Auto +/// Settings and the auto crop), with the deghost overlay tinted red. +fn render_preview( + dng: &[u8], + edge: usize, + auto: bool, + crop: Option<[f64; 4]>, + ghost: Option<&lightcraft_raster::Plane>, +) -> std::result::Result { + let (src, info) = crate::files::load_bytes(dng, edge)?; + let mut s = DevelopSettings::for_raw(info.as_shot_temp, info.as_shot_tint); + if auto { + apply_auto(&mut s, &src, &info); + } + if let Some(c) = crop { + s.crop.geometry.rect = lightcraft_geom::Rect::from_xywh(c[0], c[1], c[2], c[3]); + } + let req = RenderRequest { max_w: edge, max_h: edge, quality: Quality::Full, apply_crop: false }; + let mut img = lightcraft_pipeline::render(&src, &info, &s, &req).image; + if let Some(g) = ghost { + let (w, h) = (img.width, img.height); + for y in 0..h { + for x in 0..w { + let gx = ((x as f32 + 0.5) / w as f32 * g.width as f32) as usize; + let gy = ((y as f32 + 0.5) / h as f32 * g.height as f32) as usize; + let a = g.get(gx.min(g.width - 1), gy.min(g.height - 1)).clamp(0.0, 1.0) * 0.6; + if a > 0.01 { + let p = &mut img.data[y * w + x]; + p[0] = (p[0] as f32 * (1.0 - a) + 255.0 * a) as u8; + p[1] = (p[1] as f32 * (1.0 - a)) as u8; + p[2] = (p[2] as f32 * (1.0 - a)) as u8; + } + } + } + } + if let Some(c) = crop { + // draw the auto-crop frame + let (w, h) = (img.width as f64, img.height as f64); + let (x0, y0, x1, y1) = ((c[0] * w) as usize, (c[1] * h) as usize, ((c[0] + c[2]) * w) as usize, ((c[1] + c[3]) * h) as usize); + for y in y0..y1.min(img.height) { + for x in x0..x1.min(img.width) { + if y == y0 || y + 1 == y1 || x == x0 || x + 1 == x1 { + img.data[y * img.width + x] = [255, 255, 255, 255]; + } + } + } + } + Ok(img) +} + +fn apply_auto(s: &mut DevelopSettings, src: &Rgb32f, info: &SourceInfo) { + let a = lightcraft_pipeline::auto::auto_tone(src, info, s); + s.light.exposure = a.exposure; + s.light.contrast = a.contrast; + s.light.highlights = a.highlights; + s.light.shadows = a.shadows; + s.light.whites = a.whites; + s.light.blacks = a.blacks; + s.color.vibrance = a.vibrance; + s.color.saturation = a.saturation; +} + +/// A free file name `dir/-.dng` (`-2`, `-3`… when taken). +fn output_path(first: &str, suffix: &str) -> String { + let p = std::path::Path::new(first); + let dir = p.parent().map(|d| d.to_path_buf()).unwrap_or_default(); + let stem = p.file_stem().map(|s| s.to_string_lossy().to_string()).unwrap_or_else(|| "Merge".into()); + for k in 1.. { + let name = if k == 1 { format!("{stem}-{suffix}.dng") } else { format!("{stem}-{suffix}-{k}.dng") }; + let cand = dir.join(name); + if !cand.exists() { + return cand.to_string_lossy().to_string(); + } + } + unreachable!() +} + +/// Parse `merge.*` command parameters into a kind and finishing options. +pub fn parse(id: &str, p: &Value) -> Result<(MergeKind, MergeFinish)> { + let bad = |m: String| EngineError::BadParams { cmd: id.into(), msg: m }; + let deghost: Deghost = serde_json::from_value(p.get("deghost").cloned().unwrap_or(json!("none"))).map_err(|e| bad(format!("deghost: {e}")))?; + let hdr = HdrOptions { align: p.get("align").and_then(Value::as_bool).unwrap_or(true), deghost }; + let projection: Projection = + serde_json::from_value(p.get("projection").cloned().unwrap_or(json!("auto"))).map_err(|e| bad(format!("projection: {e}")))?; + let pano = PanoOptions { + projection, + boundary_warp: p.get("boundaryWarp").and_then(Value::as_f64).unwrap_or(0.0).clamp(0.0, 100.0), + auto_crop: p.get("autoCrop").and_then(Value::as_bool).unwrap_or(false), + fill_edges: p.get("fillEdges").and_then(Value::as_bool).unwrap_or(false), + max_megapixels: p.get("maxMegapixels").and_then(Value::as_f64).unwrap_or(40.0), + }; + let kind = match id { + "merge.hdr" => MergeKind::Hdr { hdr }, + "merge.panorama" => MergeKind::Panorama { pano }, + _ => MergeKind::HdrPanorama { hdr, pano, bracket: p.get("bracket").and_then(Value::as_u64).unwrap_or(0) as usize }, + }; + let finish = MergeFinish { + auto_settings: p.get("autoSettings").and_then(Value::as_bool).unwrap_or(true), + stack: p.get("stack").and_then(Value::as_bool).unwrap_or(false), + show_overlay: p.get("showOverlay").and_then(Value::as_bool).unwrap_or(false), + }; + Ok((kind, finish)) +} + +impl Session { + /// Plan a merge of `ids` (catalog order of the given ids is kept). + pub fn plan_merge(&self, kind: MergeKind, finish: MergeFinish, ids: &[PhotoId], preview: bool) -> Result { + let mut sources = Vec::new(); + for id in ids { + let p = self.catalog.photo(*id).ok_or(lightcraft_catalog::CatalogError::NoPhoto(*id))?; + match &p.source { + Source::File { path } => sources.push((*id, path.clone())), + Source::Demo { .. } => return Err(EngineError::Other(format!("{} is a demo photo; Photo Merge needs photo files", p.file_name))), + } + } + if sources.len() < kind.min_photos() { + return Err(EngineError::Other(format!("select at least {} photos to merge", kind.min_photos()))); + } + let read: ByteReader = match &self.media.file_bytes { + Some(r) => r.clone(), + None => Arc::new(|path: &str| std::fs::read(path).map_err(|e| format!("{path}: {e}"))), + }; + Ok(MergeJob { kind, finish, sources, preview: preview.then_some(PREVIEW_EDGE), read }) + } + + /// Write a finished (non-preview) merge next to its first source, import it, select it and + /// apply Auto Settings / the auto crop. Returns `{id, path, …info}`. + pub fn finish_merge(&mut self, job: &MergeJob, out: MergeOutput) -> Result { + let first = &job.sources.first().ok_or_else(|| EngineError::Other("nothing merged".into()))?.1; + let path = output_path(first, job.kind.suffix()); + std::fs::write(&path, &out.dng).map_err(|e| EngineError::Other(format!("{path}: {e}")))?; + let report = crate::import::import(self, std::slice::from_ref(&path), crate::import::ImportMode::Add)?; + let id = report + .imported + .first() + .copied() + .map(PhotoId) + .ok_or_else(|| EngineError::Other(format!("the merged file {path} could not be imported: {:?}", report.failed)))?; + self.selection = crate::Selection::single(id); + if job.finish.auto_settings { + self.execute("develop.auto", &json!({}))?; + } + if let Some(c) = out.crop + && let Some(d) = self.develop_of(id) + { + let mut d = (*d).clone(); + d.crop.geometry.rect = lightcraft_geom::Rect::from_xywh(c[0], c[1], c[2], c[3]); + self.set_develop(id, d, "Auto Crop")?; + } + let mut info = out.info; + info["id"] = json!(id.0); + info["path"] = json!(path); + info["bytes"] = json!(out.dng.len()); + Ok(info) + } + + /// Plan, run and finish a merge synchronously (CLI, MCP, control channel). + pub fn merge_now(&mut self, kind: MergeKind, finish: MergeFinish, ids: &[PhotoId]) -> Result { + let job = self.plan_merge(kind, finish, ids, false)?; + let out = job.run(&|_, _| true).map_err(EngineError::Other)?; + self.finish_merge(&job, out) + } +} diff --git a/crates/engine/src/tests_merge.rs b/crates/engine/src/tests_merge.rs new file mode 100644 index 00000000..507df3e1 --- /dev/null +++ b/crates/engine/src/tests_merge.rs @@ -0,0 +1,82 @@ +//! Photo Merge through the command layer: synthetic brackets / panorama views on disk → merged +//! DNG next to the sources → imported, selected, Auto Settings / auto crop applied. + +use serde_json::json; + +use crate::Session; + +fn temp_dir(tag: &str) -> std::path::PathBuf { + let d = std::env::temp_dir().join(format!("lc-merge-{tag}-{}", std::process::id())); + let _ = std::fs::remove_dir_all(&d); + std::fs::create_dir_all(&d).unwrap(); + d +} + +fn import(s: &mut Session, paths: Vec) -> Vec { + let r = s.execute("library.import", &json!({"paths": paths})).unwrap(); + r["imported"].as_array().unwrap().iter().map(|v| v.as_u64().unwrap()).collect() +} + +#[test] +fn hdr_merge_command_creates_and_imports_a_dng() { + let dir = temp_dir("hdr"); + let mut paths = Vec::new(); + for (i, b) in lightcraft_merge::synth::bracket_dngs(480, 320, &[-2.0, 0.0, 2.0]).into_iter().enumerate() { + let p = dir.join(format!("IMG_{i}.dng")); + std::fs::write(&p, b).unwrap(); + paths.push(p.to_string_lossy().to_string()); + } + let mut s = Session::new().with_fs(); + let ids = import(&mut s, paths); + assert_eq!(ids.len(), 3); + // preview: rendered PNG, nothing added + let png = dir.join("preview.png"); + let pv = s.execute("merge.hdr", &json!({"ids": ids, "deghost": "medium", "preview": true, "showOverlay": true, "previewPath": png})).unwrap(); + assert!(png.exists() && pv["previewSize"][0].as_u64().unwrap() > 100, "{pv}"); + assert_eq!(s.catalog.photos().count(), 3); + let ev: Vec = pv["ev"].as_array().unwrap().iter().map(|v| v.as_f64().unwrap()).collect(); + assert!((ev[0] + 2.0).abs() < 0.1 && (ev[2] - 2.0).abs() < 0.1, "{ev:?}"); + // full merge + let r = s.execute("merge.hdr", &json!({"ids": ids, "autoSettings": true, "stack": true})).unwrap(); + let path = r["path"].as_str().unwrap(); + assert!(path.ends_with("IMG_0-HDR.dng"), "{path}"); + let id = lightcraft_catalog::PhotoId(r["id"].as_u64().unwrap()); + let p = s.catalog.photo(id).unwrap().clone(); + assert_eq!((p.width, p.height), (480, 320)); + assert_eq!(p.kind, lightcraft_catalog::MediaKind::Raw); + assert_eq!(s.selection.active, Some(id)); + let d = s.develop_of(id).unwrap(); + assert!(d.light.exposure != 0.0 || d.light.highlights != 0.0, "auto settings applied"); + let img = s.render_now(id, 256, 256).unwrap(); + assert_eq!(img.image.width, 256); + // the second merge doesn't overwrite the first + let r2 = s.execute("merge.hdr", &json!({"ids": ids, "autoSettings": false})).unwrap(); + assert!(r2["path"].as_str().unwrap().ends_with("IMG_0-HDR-2.dng")); + let _ = std::fs::remove_dir_all(&dir); +} + +#[test] +fn panorama_command_with_auto_crop() { + let dir = temp_dir("pano"); + let mut paths = Vec::new(); + for (i, v) in lightcraft_merge::synth::pano_views(480, 360, 420.0, &[-25.0, 0.0, 25.0]).into_iter().enumerate() { + let img = v.to_srgb8(); + let png = lightcraft_codecs::encode_png(&lightcraft_codecs::EncodeImage::rgba8(&img), &lightcraft_codecs::EncodeMeta::default()).unwrap(); + let p = dir.join(format!("P_{i}.png")); + std::fs::write(&p, png).unwrap(); + paths.push(p.to_string_lossy().to_string()); + } + let mut s = Session::new().with_fs(); + let ids = import(&mut s, paths); + let r = s.execute("merge.panorama", &json!({"ids": ids, "projection": "cylindrical", "autoCrop": true, "autoSettings": false})).unwrap(); + assert_eq!(r["projection"], "cylindrical"); + assert_eq!(r["used"].as_array().unwrap().len(), 3); + let id = lightcraft_catalog::PhotoId(r["id"].as_u64().unwrap()); + let p = s.catalog.photo(id).unwrap().clone(); + assert!(p.width > 800 && p.height > 300, "{}×{}", p.width, p.height); + let crop = s.develop_of(id).unwrap().crop.geometry.rect; + assert!(crop.width() < 1.0 && crop.width() * crop.height() > 0.5, "{crop:?}"); + // too few photos + assert!(s.execute("merge.panorama", &json!({"ids": [ids[0]]})).is_err()); + let _ = std::fs::remove_dir_all(&dir); +} diff --git a/crates/merge/Cargo.toml b/crates/merge/Cargo.toml index f0f5fb1f..296f06be 100644 --- a/crates/merge/Cargo.toml +++ b/crates/merge/Cargo.toml @@ -15,9 +15,7 @@ lightcraft-raster = { workspace = true } lightcraft-raw = { workspace = true } lightcraft-codecs = { workspace = true } lightcraft-meta = { workspace = true } +lightcraft-scenes = { workspace = true } serde = { workspace = true } thiserror = { workspace = true } rayon = { workspace = true } - -[dev-dependencies] -lightcraft-scenes = { workspace = true } diff --git a/crates/merge/src/lib.rs b/crates/merge/src/lib.rs index 93fafa6f..65bd86c4 100644 --- a/crates/merge/src/lib.rs +++ b/crates/merge/src/lib.rs @@ -20,6 +20,7 @@ pub mod linalg; pub mod output; pub mod pano; pub mod ransac; +pub mod synth; pub use frame::{Frame, FrameColor, load_frame}; pub use hdr::{Deghost, HdrOptions, HdrResult, merge_hdr}; diff --git a/crates/merge/src/synth.rs b/crates/merge/src/synth.rs new file mode 100644 index 00000000..1d8fd8aa --- /dev/null +++ b/crates/merge/src/synth.rs @@ -0,0 +1,81 @@ +//! Synthetic merge inputs from procedural scenes (tests, demos, UI checks — no media files): +//! exposure brackets as linear DNGs and overlapping panorama views as PNGs. + +use lightcraft_geom::Orientation; +use lightcraft_meta::Metadata; +use lightcraft_raster::Rgb32f; + +use crate::frame::FrameColor; +use crate::linalg::{self, M3}; +use crate::output::{DngSamples, write_linear_dng}; + +fn scene(kind: lightcraft_scenes::Kind) -> lightcraft_scenes::Scene { + lightcraft_scenes::demo_library().into_iter().find(|s| s.kind == kind).expect("scene kind in the demo library") +} + +fn srgb_color() -> FrameColor { + FrameColor::Linear { to_xyz_d50: lightcraft_color::bradford(lightcraft_color::D65, lightcraft_color::D50).mul(&lightcraft_color::SRGB.to_xyz()) } +} + +/// Exposure brackets (`evs`, e.g. `[-2, 0, 2]`) of a sunset scene as 16-bit linear DNGs, each +/// shifted by a few pixels (handheld) with EXIF exposure times; `w × h` pixels. +pub fn bracket_dngs(w: usize, h: usize, evs: &[f64]) -> Vec> { + let m = 12usize; + let truth = scene(lightcraft_scenes::Kind::OceanSunset).render(w + 2 * m, h + 2 * m); + let mut l: Vec = truth.data.iter().map(|p| p[1]).collect(); + l.sort_by(|a, b| a.total_cmp(b)); + let s = 0.18 / l[l.len() / 2].max(1e-6); + let to_srgb = lightcraft_color::REC2020.to_space(&lightcraft_color::SRGB).to_f32(); + evs.iter() + .enumerate() + .map(|(i, ev)| { + let k = 2f32.powf(*ev as f32) * s; + let (dx, dy) = ((i as f32 * 2.3) % 5.0 - 2.0, (i as f32 * 1.7) % 4.0 - 1.5); + let img = Rgb32f::from_fn(w, h, |x, y| { + let p = truth.sample_bilinear(x as f32 + 0.5 + m as f32 - dx, y as f32 + 0.5 + m as f32 - dy); + let q = [0, 1, 2].map(|r| to_srgb[r][0] * p[0] + to_srgb[r][1] * p[1] + to_srgb[r][2] * p[2]); + q.map(|v| (v * k).clamp(0.0, 1.0)) + }); + let meta = Metadata { + make: Some("LightCraft".into()), + model: Some("Synthetic Bracket".into()), + exposure_time: Some(2f64.powf(*ev) / 125.0), + f_number: Some(8.0), + iso: Some(100), + ..Default::default() + }; + write_linear_dng(&img, &srgb_color(), Orientation::Normal, &meta, 0.0, DngSamples::U16).expect("encode") + }) + .collect() +} + +/// `yaws.len()` overlapping views (degrees of yaw, small pitch/roll wobble) of a canyon scene used +/// as a 200° × 70° spherical environment, `w × h` pixels with focal length `f`, as linear sRGB. +pub fn pano_views(w: usize, h: usize, f: f64, yaws: &[f64]) -> Vec { + let env = scene(lightcraft_scenes::Kind::Canyon).render(2400, 840); + let mut l: Vec = env.data.iter().map(|p| p[1]).collect(); + l.sort_by(|a, b| a.total_cmp(b)); + let s = 0.25 / l[l.len() / 2].max(1e-6); + let (theta, phi) = (200.0f64, 70.0f64); + yaws.iter() + .enumerate() + .map(|(i, yaw)| { + let wob = (i as f64 * 1.3).sin(); + let ry = linalg::rodrigues([0.0, yaw.to_radians(), 0.0]); + let rx = linalg::rodrigues([(wob * 1.5).to_radians(), 0.0, 0.0]); + let rz = linalg::rodrigues([0.0, 0.0, (wob * 0.8).to_radians()]); + let c2w: M3 = linalg::mul3(&linalg::mul3(&ry, &rx), &rz); + Rgb32f::from_fn(w, h, |x, y| { + let d = linalg::apply3(&c2w, [(x as f64 + 0.5 - w as f64 / 2.0) / f, (y as f64 + 0.5 - h as f64 / 2.0) / f, 1.0]); + let t = d[0].atan2(d[2]).to_degrees(); + let p = (d[1] / linalg::norm(d)).asin().to_degrees(); + if t.abs() > theta / 2.0 || p.abs() > phi / 2.0 { + return [0.0; 3]; + } + let px = (t / theta + 0.5) * env.width as f64; + let py = (p / phi + 0.5) * env.height as f64; + env.sample_bilinear(px as f32, py as f32).map(|v| (v * s).min(1.0)) + }) + }) + .collect() +} diff --git a/crates/ui-egui/src/control.rs b/crates/ui-egui/src/control.rs index 8d6cbfae..53912d37 100644 --- a/crates/ui-egui/src/control.rs +++ b/crates/ui-egui/src/control.rs @@ -77,7 +77,7 @@ pub fn inspect(app: &LightcraftApp, ctx: &egui::Context) -> Value { "selection": app.session.selection.ids.iter().map(|p| p.0).collect::>(), "activeMask": app.session.active_mask, "widgetCount": app.widgets.len(), - "perf": {"frameMs": app.perf.frame_ms, "fps": app.perf.fps, "lastRenderMs": app.renderer.last_main_ms, "renderQueue": app.renderer.queued(), "rendersInFlight": app.renderer.in_flight(), "rendersDone": app.renderer.completed, "thumbTextures": app.renderer.thumb_textures(), "gpu": lightcraft_engine::gpu::available().then(lightcraft_engine::gpu::adapter_name).flatten()}, + "perf": {"frameMs": app.perf.frame_ms, "fps": app.perf.fps, "lastRenderMs": app.renderer.last_main_ms, "renderQueue": app.renderer.queued(), "rendersInFlight": app.renderer.in_flight(), "mergeRunning": app.merge.busy(), "lastMerge": app.merge.last_result, "rendersDone": app.renderer.completed, "thumbTextures": app.renderer.thumb_textures(), "gpu": lightcraft_engine::gpu::available().then(lightcraft_engine::gpu::adapter_name).flatten()}, "status": app.ui.status, }) } diff --git a/crates/ui-egui/src/headless.rs b/crates/ui-egui/src/headless.rs index ccac06f2..ae6398b3 100644 --- a/crates/ui-egui/src/headless.rs +++ b/crates/ui-egui/src/headless.rs @@ -171,7 +171,7 @@ impl Headless { /// Is anything still in progress (renders, queued input)? pub fn busy(&self) -> bool { - self.app.renderer.in_flight() > 0 || !self.app.synthetic.is_empty() || !self.events.is_empty() + self.app.renderer.in_flight() > 0 || self.app.merge.busy() || !self.app.synthetic.is_empty() || !self.events.is_empty() } /// Run frames until nothing is pending (renders finished, input consumed) for a few frames in diff --git a/crates/ui-egui/src/lib.rs b/crates/ui-egui/src/lib.rs index 05013ce1..3de5c6d0 100644 --- a/crates/ui-egui/src/lib.rs +++ b/crates/ui-egui/src/lib.rs @@ -9,6 +9,7 @@ pub mod control; pub mod headless; pub mod icons; pub mod menus; +pub mod merge; pub mod panels; pub mod render; pub mod shortcuts; @@ -80,6 +81,8 @@ pub struct LightcraftApp { pub widgets: Vec<(String, egui::Rect)>, /// In-progress on-canvas gesture (brush stroke points, gradient drag…). pub gesture: Option, + /// Photo Merge dialog previews and background merges. + pub merge: merge::MergeState, } impl LightcraftApp { @@ -105,6 +108,7 @@ impl LightcraftApp { image_rect: None, widgets: vec![], gesture: None, + merge: merge::MergeState::default(), } } @@ -168,7 +172,7 @@ impl LightcraftApp { return; } let now = now_ms(); - let busy = self.renderer.in_flight() > 0; + let busy = self.renderer.in_flight() > 0 || self.merge.busy(); let mut shots = std::mem::take(&mut self.pending_screenshots); let mut shadow_ticked = false; shots.retain_mut(|s| { @@ -225,7 +229,13 @@ impl LightcraftApp { let memory = main.memory(|m| m.clone()); view.ctx.memory_mut(|m| *m = memory); view.run(headless::HeadlessView::raw_input(size, ppp, time, vec![]), |ui| self.ui(ui)); - let img = capture.then(|| view.paint(&self.renderer.cpu_textures())); + let img = capture.then(|| { + let mut tex = self.renderer.cpu_textures(); + if let (Some((t, ..)), Some(px)) = (&self.merge.preview, &self.merge.preview_pixels) { + tex.insert(t.id(), crate::softpaint::CpuTexture::linear(px.clone())); + } + view.paint(&tex) + }); self.shadow = Some(view); img } @@ -275,6 +285,7 @@ impl LightcraftApp { ctx.request_repaint(); } self.renderer.poll(ctx, &mut self.session); + merge::poll(self, ctx); self.session.persist_if_dirty(); self.collect_screenshots(ctx); self.issue_screenshots(ctx); diff --git a/crates/ui-egui/src/menus.rs b/crates/ui-egui/src/menus.rs index 1dd7101f..a052fe25 100644 --- a/crates/ui-egui/src/menus.rs +++ b/crates/ui-egui/src/menus.rs @@ -55,6 +55,9 @@ pub const UI_COMMANDS: &[UiCommand] = &[ ("dialog.createPreset", "Create Preset…", Some("Cmd+Shift+P"), "Photo"), ("dialog.copySettings", "Choose Edit Settings to Copy…", Some("Cmd+Shift+C"), "Photo"), ("dialog.export", "Export…", Some("Cmd+Shift+E"), "File"), + ("dialog.mergeHdr", "Photo Merge ▸ HDR…", Some("Ctrl+H"), "Photo"), + ("dialog.mergePanorama", "Photo Merge ▸ Panorama…", Some("Ctrl+M"), "Photo"), + ("dialog.mergeHdrPanorama", "Photo Merge ▸ HDR Panorama…", None, "Photo"), ("file.addPhotos", "Add Photos…", Some("Cmd+Shift+I"), "File"), ("file.importPresets", "Import Presets…", None, "File"), ("file.exportPresets", "Export Presets…", None, "File"), @@ -281,6 +284,9 @@ pub fn run_ui_command(app: &mut LightcraftApp, id: &str, p: &Value) -> Option crate::merge::open(app, "merge.hdr"), + "dialog.mergePanorama" => crate::merge::open(app, "merge.panorama"), + "dialog.mergeHdrPanorama" => crate::merge::open(app, "merge.hdrPanorama"), "app.about" => { app.ui.dialog = Some(Dialog::About); Ok(Value::Null) @@ -363,6 +369,7 @@ pub fn ui_enabled(app: &LightcraftApp, id: &str) -> bool { s if s.starts_with("panel.") || s.starts_with("tool.") || s.starts_with("section.") => app.session.active().is_some() || s == "panel.close", "app.export" | "dialog.export" | "dialog.createPreset" | "dialog.copySettings" => app.session.active().is_some(), "file.exportPresets" => app.session.presets.iter().any(|p| !p.builtin), + s if s.starts_with("dialog.merge") => app.session.targets(&serde_json::json!({})).len() >= 2 && app.merge.final_task.is_none(), _ => true, } } diff --git a/crates/ui-egui/src/merge.rs b/crates/ui-egui/src/merge.rs new file mode 100644 index 00000000..507eee2f --- /dev/null +++ b/crates/ui-egui/src/merge.rs @@ -0,0 +1,389 @@ +//! Photo Merge in the UI: the HDR / Panorama / HDR Panorama dialog (options + live preview) and +//! the background merge jobs. +//! +//! Previews and the final merge run [`lightcraft_engine::merge::MergeJob`]s on worker threads with +//! progress and cancellation; the window stays responsive. Changing an option cancels the running +//! preview and starts a new one. The final merge keeps running after the dialog closes (progress +//! in a toast); when it finishes, the result is written, imported and selected +//! ([`lightcraft_engine::Session::finish_merge`]). + +use std::sync::atomic::{AtomicBool, Ordering}; +use std::sync::mpsc::{Receiver, channel}; +use std::sync::{Arc, Mutex}; + +use lightcraft_catalog::PhotoId; +use lightcraft_engine::merge::{MergeJob, MergeOutput}; +use serde::{Deserialize, Serialize}; +use serde_json::{Value, json}; + +use crate::LightcraftApp; + +/// The dialog's options (also the parameters of the `merge.*` commands). +#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)] +#[serde(default, rename_all = "camelCase")] +pub struct MergeDialog { + /// "merge.hdr" | "merge.panorama" | "merge.hdrPanorama". + pub command: String, + pub align: bool, + /// none | low | medium | high + pub deghost: String, + pub show_overlay: bool, + pub auto_settings: bool, + pub stack: bool, + /// auto | spherical | cylindrical | perspective + pub projection: String, + pub boundary_warp: f64, + pub auto_crop: bool, + pub fill_edges: bool, + /// HDR panorama: photos per bracket (0 = from EXIF). + pub bracket: u32, +} + +impl Default for MergeDialog { + fn default() -> Self { + MergeDialog { + command: "merge.hdr".into(), + align: true, + deghost: "none".into(), + show_overlay: false, + auto_settings: true, + stack: false, + projection: "auto".into(), + boundary_warp: 0.0, + auto_crop: true, + fill_edges: false, + bracket: 0, + } + } +} + +impl MergeDialog { + pub fn for_command(command: &str) -> MergeDialog { + MergeDialog { command: command.into(), ..Default::default() } + } + pub fn title(&self) -> &'static str { + match self.command.as_str() { + "merge.panorama" => "Panorama Merge Preview", + "merge.hdrPanorama" => "HDR Panorama Merge Preview", + _ => "HDR Merge Preview", + } + } + pub fn is_hdr(&self) -> bool { + self.command != "merge.panorama" + } + pub fn is_pano(&self) -> bool { + self.command != "merge.hdr" + } + pub fn params(&self) -> Value { + serde_json::to_value(self).unwrap_or_default() + } +} + +/// A merge running on a worker thread. +pub struct MergeTask { + pub preview: bool, + pub job: MergeJob, + /// The dialog options it was started with. + pub options: MergeDialog, + pub progress: Arc>, + pub cancel: Arc, + rx: Receiver>, +} + +/// Merge state of the app. +#[derive(Default)] +pub struct MergeState { + pub preview_task: Option, + pub final_task: Option, + /// The latest preview (texture, size, its options, stats) and the latest error. + pub preview: Option<(egui::TextureHandle, [usize; 2], MergeDialog, Value)>, + pub preview_pixels: Option>, + pub error: Option, + /// Options whose preview failed (not retried until they change). + pub failed_options: Option, + /// The last finished merge (`{id, path, …}`). + pub last_result: Option, + /// Photos the dialog merges (the selection when it opened). + pub ids: Vec, +} + +impl MergeState { + /// A preview or a merge is running. + pub fn busy(&self) -> bool { + self.preview_task.is_some() || self.final_task.is_some() + } +} + +fn spawn(app: &LightcraftApp, options: &MergeDialog, ids: &[PhotoId], preview: bool) -> Result { + let (kind, finish) = lightcraft_engine::merge::parse(&options.command, &options.params()).map_err(|e| e.to_string())?; + let job = app.session.plan_merge(kind, finish, ids, preview).map_err(|e| e.to_string())?; + let progress = Arc::new(Mutex::new((0.0, "Starting".to_string()))); + let cancel = Arc::new(AtomicBool::new(false)); + let (tx, rx) = channel(); + let (j, p, c) = (job.clone(), progress.clone(), cancel.clone()); + let work = move || { + let r = j.run(&|f, stage| { + if let Ok(mut g) = p.lock() { + *g = (f, stage.to_string()); + } + !c.load(Ordering::Relaxed) + }); + let _ = tx.send(r); + }; + #[cfg(not(target_arch = "wasm32"))] + std::thread::Builder::new().name("photo-merge".into()).spawn(work).map_err(|e| e.to_string())?; + #[cfg(target_arch = "wasm32")] + work(); + Ok(MergeTask { preview, job, options: options.clone(), progress, cancel, rx }) +} + +/// Open the merge dialog for the selection. +pub fn open(app: &mut LightcraftApp, command: &str) -> Result { + let ids = app.session.targets(&json!({})); + if ids.len() < 2 { + return Err("select at least 2 photos to merge".into()); + } + app.merge.ids = ids; + app.merge.error = None; + app.merge.preview = None; + app.merge.preview_pixels = None; + if let Some(t) = app.merge.preview_task.take() { + t.cancel.store(true, Ordering::Relaxed); + } + app.ui.dialog = Some(crate::state::Dialog::Merge { opts: MergeDialog::for_command(command) }); + Ok(json!({"photos": app.merge.ids.len()})) +} + +/// Start the full-resolution merge (the dialog's Merge button / `ui.dialog.confirm`). +pub fn start_final(app: &mut LightcraftApp, opts: &MergeDialog) -> Result { + if app.merge.final_task.is_some() { + return Err("a merge is already running".into()); + } + if let Some(t) = app.merge.preview_task.take() { + t.cancel.store(true, Ordering::Relaxed); + } + let ids = app.merge.ids.clone(); + let task = spawn(app, opts, &ids, false)?; + app.merge.final_task = Some(task); + Ok(json!({"started": true, "photos": ids.len()})) +} + +/// Per frame: keep the preview in sync with the dialog's options, collect finished jobs. +pub fn poll(app: &mut LightcraftApp, ctx: &egui::Context) { + // preview for the open dialog + let dialog_opts = match &app.ui.dialog { + Some(crate::state::Dialog::Merge { opts }) => Some(opts.clone()), + _ => None, + }; + match &dialog_opts { + Some(opts) => { + let current = app.merge.preview_task.as_ref().map(|t| &t.options).or(app.merge.preview.as_ref().map(|p| &p.2)); + let stale = current != Some(opts) && !last_failed_matches(app, opts); + if stale { + if let Some(t) = app.merge.preview_task.take() { + t.cancel.store(true, Ordering::Relaxed); + } + let ids = app.merge.ids.clone(); + match spawn(app, opts, &ids, true) { + Ok(t) => { + app.merge.error = None; + app.merge.preview_task = Some(t); + } + Err(e) => { + app.merge.error = Some(e); + app.merge.failed_options = Some(opts.clone()); + } + } + } + } + None => { + if let Some(t) = app.merge.preview_task.take() { + t.cancel.store(true, Ordering::Relaxed); + } + } + } + if let Some(t) = &app.merge.preview_task { + match t.rx.try_recv() { + Ok(Ok(out)) => { + let t = app.merge.preview_task.take().expect("task"); + if let Some(img) = out.preview { + let color = Arc::new(egui::ColorImage::from_rgba_unmultiplied([img.width, img.height], &img.as_bytes())); + let tex = ctx.load_texture("merge-preview", color.clone(), egui::TextureOptions::LINEAR); + app.merge.preview_pixels = Some(color); + app.merge.preview = Some((tex, [img.width, img.height], t.options, out.info)); + } + } + Ok(Err(e)) => { + let t = app.merge.preview_task.take().expect("task"); + if e != "cancelled" { + app.merge.error = Some(e); + app.merge.failed_options = Some(t.options); + } + } + Err(_) => ctx.request_repaint_after(std::time::Duration::from_millis(50)), + } + } + // the final merge + if let Some(t) = &app.merge.final_task { + match t.rx.try_recv() { + Ok(r) => { + let t = app.merge.final_task.take().expect("task"); + let what = if t.options.command == "merge.panorama" { + "Panorama" + } else if t.options.command == "merge.hdrPanorama" { + "HDR panorama" + } else { + "HDR" + }; + match r.map_err(lightcraft_engine::EngineError::Other).and_then(|out| app.session.finish_merge(&t.job, out)) { + Ok(v) => { + app.merge.last_result = Some(v); + app.toast(ctx, format!("{what} merge added")); + } + Err(e) => { + app.ui.status = e.to_string(); + app.toast(ctx, format!("{what} merge failed: {e}")); + } + } + } + Err(_) => { + let (f, stage) = t.progress.lock().map(|g| g.clone()).unwrap_or_default(); + let now = ctx.input(|i| i.time); + app.ui.toast = Some((format!("Merging… {stage} {:.0}%", f * 100.0), now + 0.5)); + ctx.request_repaint_after(std::time::Duration::from_millis(100)); + } + } + } +} + +fn last_failed_matches(app: &LightcraftApp, opts: &MergeDialog) -> bool { + app.merge.failed_options.as_ref() == Some(opts) +} + +/// The dialog body (options on the left, preview on the right). +pub fn body(app: &mut LightcraftApp, ui: &mut egui::Ui, opts: &mut MergeDialog) { + let t = crate::theme::Tokens::get(ui.ctx()); + ui.horizontal_top(|ui| { + // preview + ui.vertical(|ui| { + let box_size = egui::vec2(520.0, 360.0); + let (rect, _) = ui.allocate_exact_size(box_size, egui::Sense::hover()); + ui.painter().rect_filled(rect, 4.0, t.canvas); + crate::widgets::register(ui.ctx(), "merge.preview", rect); + if let Some((tex, size, _, _)) = &app.merge.preview { + let s = (box_size.x / size[0] as f32).min(box_size.y / size[1] as f32); + let r = egui::Rect::from_center_size(rect.center(), egui::vec2(size[0] as f32 * s, size[1] as f32 * s)); + ui.painter().image(tex.id(), r, egui::Rect::from_min_max(egui::pos2(0.0, 0.0), egui::pos2(1.0, 1.0)), egui::Color32::WHITE); + } + if let Some(task) = &app.merge.preview_task { + let (f, stage) = task.progress.lock().map(|g| g.clone()).unwrap_or_default(); + ui.painter().text( + rect.center_bottom() - egui::vec2(0.0, 14.0), + egui::Align2::CENTER_CENTER, + format!("{stage}… {:.0}%", f * 100.0), + t.font(12.0), + t.text, + ); + } + if let Some(e) = &app.merge.error { + ui.painter().text(rect.center(), egui::Align2::CENTER_CENTER, e, t.font(13.0), t.text); + } + let info = app.merge.preview.as_ref().map(|p| p.3.clone()); + if let Some(info) = info { + let txt = if opts.is_pano() { + format!( + "{} photos · {} · {:.0}° × {:.0}°", + info["used"].as_array().map_or(0, Vec::len), + info["projection"].as_str().unwrap_or(""), + info["fov"][0].as_f64().unwrap_or(0.0), + info["fov"][1].as_f64().unwrap_or(0.0) + ) + } else { + let ev: Vec = + info["ev"].as_array().map(|a| a.iter().filter_map(Value::as_f64).map(|v| format!("{v:+.1}")).collect()).unwrap_or_default(); + format!("{} photos · exposures {} EV", ev.len(), ev.join(" / ")) + }; + ui.label(egui::RichText::new(txt).color(t.text_dim)); + } + }); + ui.add_space(12.0); + // options + ui.vertical(|ui| { + ui.set_width(220.0); + if opts.is_pano() { + ui.label(egui::RichText::new("Projection").color(t.text_label)); + for (v, l) in [("auto", "Auto Select"), ("spherical", "Spherical"), ("cylindrical", "Cylindrical"), ("perspective", "Perspective")] { + if ui.radio(opts.projection == v, l).clicked() { + opts.projection = v.to_string(); + } + } + ui.add_space(4.0); + let mut bw = opts.boundary_warp; + if num(ui, &BOUNDARY_WARP, &mut bw) { + opts.boundary_warp = bw.round(); + } + ui.checkbox(&mut opts.fill_edges, "Fill Edges"); + ui.checkbox(&mut opts.auto_crop, "Auto Crop"); + ui.add_space(4.0); + } + if opts.is_hdr() { + ui.checkbox(&mut opts.align, "Auto Align"); + ui.label(egui::RichText::new("Deghost Amount").color(t.text_label)); + ui.horizontal(|ui| { + for (i, (v, l)) in [("none", "None"), ("low", "Low"), ("medium", "Med"), ("high", "High")].iter().enumerate() { + if crate::widgets::text_button(ui, &format!("mergeDeghost-{i}"), l, opts.deghost == *v).clicked() { + opts.deghost = v.to_string(); + } + } + }); + if opts.command == "merge.hdr" { + ui.add_enabled_ui(opts.deghost != "none", |ui| ui.checkbox(&mut opts.show_overlay, "Show Deghost Overlay")); + } + if opts.command == "merge.hdrPanorama" { + let mut b = opts.bracket as f64; + if num(ui, &BRACKET, &mut b) { + opts.bracket = b.round() as u32; + } + } + ui.add_space(4.0); + } + ui.checkbox(&mut opts.auto_settings, "Auto Settings"); + ui.checkbox(&mut opts.stack, "Create Stack"); + ui.add_space(8.0); + ui.label(egui::RichText::new(format!("{} photos", app.merge.ids.len())).color(t.text_dim)); + }); + }); +} + +const BOUNDARY_WARP: lightcraft_develop::ControlSpec = lightcraft_develop::ControlSpec { + id: "merge.boundaryWarp", + label: "Boundary Warp", + section: lightcraft_develop::Section::Light, + min: 0.0, + max: 100.0, + default: 0.0, + step: 1.0, + decimals: 0, + track: lightcraft_develop::Track::Plain, +}; +const BRACKET: lightcraft_develop::ControlSpec = lightcraft_develop::ControlSpec { + id: "merge.bracket", + label: "Photos per bracket (0 = auto)", + section: lightcraft_develop::Section::Light, + min: 0.0, + max: 9.0, + default: 0.0, + step: 1.0, + decimals: 0, + track: lightcraft_develop::Track::Plain, +}; + +fn num(ui: &mut egui::Ui, spec: &lightcraft_develop::ControlSpec, v: &mut f64) -> bool { + match crate::widgets::slider(ui, spec, *v, true, None).value { + Some(n) => { + *v = n; + true + } + None => false, + } +} diff --git a/crates/ui-egui/src/panels/dialogs.rs b/crates/ui-egui/src/panels/dialogs.rs index 0269aad2..466e48ac 100644 --- a/crates/ui-egui/src/panels/dialogs.rs +++ b/crates/ui-egui/src/panels/dialogs.rs @@ -25,6 +25,7 @@ pub fn show(app: &mut LightcraftApp, ctx: &egui::Context) { Dialog::CreatePreset { .. } => "Create Preset", Dialog::CopySettings { .. } => "Choose Edit Settings to Copy", Dialog::Export { .. } => "Export", + Dialog::Merge { opts } => opts.title(), Dialog::About => "About LightCraft", Dialog::Shortcuts => "Keyboard Shortcuts", }; @@ -156,6 +157,7 @@ pub fn show(app: &mut LightcraftApp, ctx: &egui::Context) { }); field(ui, "Folder", |ui| ui.add(egui::TextEdit::singleline(dir).desired_width(f32::INFINITY))); } + Dialog::Merge { opts } => crate::merge::body(app, ui, opts), Dialog::About => { ui.label(egui::RichText::new("LightCraft").font(t.semibold(20.0)).color(t.text)); ui.label(format!("Version {} — a clean-room, pure-Rust photo library and raw developer.", env!("CARGO_PKG_VERSION"))); @@ -190,7 +192,12 @@ pub fn show(app: &mut LightcraftApp, ctx: &egui::Context) { if !informational && ui.button("Cancel").clicked() { close = true; } - if ui.button(if informational { "Close" } else { "OK" }).clicked() { + let ok = match dlg { + Dialog::Merge { .. } => "Merge", + _ if informational => "Close", + _ => "OK", + }; + if ui.button(ok).clicked() { if informational { close = true; } else { @@ -229,6 +236,7 @@ pub fn confirm_dialog(app: &mut LightcraftApp, dlg: &Dialog) -> Result crate::merge::start_final(app, opts), Dialog::About | Dialog::Shortcuts => Ok(serde_json::Value::Null), } } diff --git a/crates/ui-egui/src/panels/grid.rs b/crates/ui-egui/src/panels/grid.rs index 6312c0cb..c4b14db4 100644 --- a/crates/ui-egui/src/panels/grid.rs +++ b/crates/ui-egui/src/panels/grid.rs @@ -256,6 +256,14 @@ pub fn context_menu(app: &mut LightcraftApp, ui: &mut egui::Ui, id: PhotoId) { if ui.button("Reset Edits").clicked() { let _ = app.run("develop.reset", json!({})); } + ui.menu_button("Photo Merge", |ui| { + let n = app.session.targets(&json!({})).len(); + for (id, label) in [("dialog.mergeHdr", "HDR…"), ("dialog.mergePanorama", "Panorama…"), ("dialog.mergeHdrPanorama", "HDR Panorama…")] { + if ui.add_enabled(n >= 2, egui::Button::new(label)).clicked() { + let _ = app.run(id, json!({})); + } + } + }); ui.separator(); if ui.button("Rotate Left").clicked() { let _ = app.run("photo.rotateLeft", json!({})); diff --git a/crates/ui-egui/src/shortcuts.rs b/crates/ui-egui/src/shortcuts.rs index af53a4bd..a808c99f 100644 --- a/crates/ui-egui/src/shortcuts.rs +++ b/crates/ui-egui/src/shortcuts.rs @@ -38,7 +38,10 @@ pub fn parse(s: &str) -> Option<(Modifiers, Key)> { fn matches(i: &egui::InputState, m: Modifiers, k: Key) -> bool { i.events.iter().any(|e| match e { egui::Event::Key { key, pressed: true, modifiers, .. } => { - *key == k && modifiers.command == m.command && modifiers.shift == m.shift && modifiers.alt == m.alt + // `Ctrl` is the physical Control key (on macOS distinct from Cmd; elsewhere Cmd = Ctrl) + let ctrl_ok = if m.ctrl { modifiers.ctrl } else { !modifiers.ctrl || modifiers.command }; + let cmd_ok = m.ctrl || modifiers.command == m.command; + *key == k && ctrl_ok && cmd_ok && modifiers.shift == m.shift && modifiers.alt == m.alt } _ => false, }) diff --git a/crates/ui-egui/src/state.rs b/crates/ui-egui/src/state.rs index 34b3b15a..c5af03f8 100644 --- a/crates/ui-egui/src/state.rs +++ b/crates/ui-egui/src/state.rs @@ -141,6 +141,10 @@ pub enum Dialog { limit_kb: u32, dir: String, }, + /// Photo Merge (HDR / Panorama / HDR Panorama) options; the preview lives in the app. + Merge { + opts: crate::merge::MergeDialog, + }, About, Shortcuts, } diff --git a/docs/merge.md b/docs/merge.md new file mode 100644 index 00000000..b802bacf --- /dev/null +++ b/docs/merge.md @@ -0,0 +1,63 @@ +# Photo Merge (HDR, Panorama, HDR Panorama) + +`crates/merge` (layer L3, no UI dependencies) implements the algorithms; `crates/engine/src/merge.rs` +plans, runs and finishes merges in the library; the desktop UI shows a dialog with a live preview. +All algorithms are our own implementations from published papers (cited in the module docs). + +## Using it + +- **App:** select the photos, then *Photo Merge ▸ HDR…* (`Ctrl+H`), *Panorama…* (`Ctrl+M`) or + *HDR Panorama…* from the photo context menu (or the `dialog.merge*` commands). The dialog previews + the merge at ≤ 1024 px and re-runs the preview when an option changes; *Merge* runs the full merge + in the background (progress in a toast). The result is written next to the first photo as + `-HDR.dng`, `-Pano.dng` or `-HDR-Pano.dng`, imported and selected. +- **Commands** (CLI, MCP, control channel): + - `merge.hdr {ids?, align=true, deghost=none|low|medium|high, autoSettings=true, stack=false, preview=false, showOverlay=false, previewPath?}` + - `merge.panorama {ids?, projection=auto|spherical|cylindrical|perspective, boundaryWarp=0..100, autoCrop=false, fillEdges=false, autoSettings=true, maxMegapixels=40, preview, previewPath?}` + - `merge.hdrPanorama {ids?, bracket=0 (from EXIF), …both sets of options}` +- **CLI:** `lightcraft-cli merge hdr|panorama|hdr-panorama [options] FILES…`; + `lightcraft-cli synth-merge hdr|panorama -o DIR` writes synthetic inputs (procedural scene). + +## What it does + +**HDR** (`hdr.rs`): frames ordered by EXIF exposure (or median luminance); the middle one is the +reference. Auto Align: features on exposure-normalised luminance, RANSAC translation / similarity / +homography (the simplest model that explains ≥ 97 % of the matches). Exposure ratios are measured +(median log ratio of well-exposed pixels between brightness neighbours), EXIF is the fallback. +Radiance is a weighted mean in linear light, weight `k · (1 − smoothstep(0.8·clip, 0.96·clip, max +channel))`. Deghost: a per-pixel consistency test against the reference radiance (with a noise +floor), dilated and feathered masks; the union is the deghost overlay. Output: 16-bit float DNG. + +**Panorama** (`pano/`): MOPS-style features, pairwise RANSAC homographies with Brown & Lowe's +verification, focal from homographies (or EXIF 35 mm), rotations chained along a maximum spanning +tree, Levenberg–Marquardt bundle adjustment of rotations + focal lengths, straightening. +Projections: spherical, cylindrical, perspective; *Auto* picks perspective for compact fields of +view < 110°, cylindrical when the vertical field is < 90°, else spherical. Exposure: EXIF, then +gain compensation from the overlaps. Seams: each pixel goes to the image with the highest +centre-weighted feather weight; multi-band blending (Laplacian pyramids) in log space. Boundary +Warp: a separable stretch of the smoothed valid-area edges to the canvas edges. Auto Crop: the +largest rectangle of image data (set as the develop crop, so it stays adjustable). Fill Edges: +push–pull diffusion. Output: 16-bit integer DNG (16-bit float for HDR panoramas). + +**DNG output** (`output.rs`): LinearRaw, 3 samples, Deflate tiles; values scaled below 1.0 with the +scale in `BaselineExposure`; raw sources keep their DNG colour tags (camera RGB), other sources get +`ColorMatrix1`/`ForwardMatrix1` of their RGB space. + +## Measured (synthetic tests, `cargo test -p lightcraft-merge`) + +- HDR, 3 brackets ±2 EV with 14-bit noise and handheld shifts: alignment error < 0.3 px, exposure + error < 0.02 EV, radiance median log error < 0.03 (noise-limited, bias < 0.01), recovered + highlights < 0.05; moving object: ghost error 0.2–2.0 without deghosting, < 0.06 with High. +- Panorama, 5 views over 182° with a 1.3× exposure drift: relative rotations < 0.1°, focal < 1 %, + BA RMS ≈ 0.35 px; median log error vs ground truth < 0.02; no seam steps. + +## Limits (known gaps) + +- No lens-distortion model in the panorama bundle adjustment (rectilinear lenses assumed); no + 360° wrap-around; flat (translating-camera) scans are fitted with the rotation model. +- No camera response calibration for JPEG brackets (sRGB decoding is taken as linear light). +- Boundary Warp is a separable stretch, not a content-preserving mesh warp; Fill Edges is diffusion, + not patch synthesis. +- *Create Stack*: LightCraft has no stacks yet; the result's caption lists its sources. +- Panorama output is limited to `maxMegapixels` (default 40) — blending memory is ~25 bytes/pixel. +- The browser build has no file access for merges yet. From b5a873a6312258e143694d203a0c3202306d8b88 Mon Sep 17 00:00:00 2001 From: Brandon Thomas Date: Thu, 1 Oct 2026 10:43:17 -0400 Subject: [PATCH 3/3] docs: Photo Merge in README and parity Co-Authored-By: Claude Opus 5.5 --- README.md | 1 + docs/parity.md | 10 +++++----- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/README.md b/README.md index 19493f0b..d95367da 100644 --- a/README.md +++ b/README.md @@ -238,6 +238,7 @@ Lightroom parity tracker — every feature, menu item and shortcut with its stat | MCP server (headless or live app, persistent libraries), CLI, control channel | ✅ | | XMP sidecars (read/write, auto-write), reading `crs:` develop settings, preset files (`.lcpreset`, XMP presets) | ✅ | | Optics (distortion, vignetting, auto + manual CA, defringe, DNG-embedded lens corrections), Geometry (transforms, Constrain Crop), Upright (Auto/Level/Vertical/Full/Guided) | ✅ · camera lens profiles (our own) ⬜ | +| Photo Merge: HDR (auto-align, deghost), Panorama (spherical/cylindrical/perspective, boundary warp, auto crop), HDR Panorama → DNG | ✅ | | GPU pipeline (wgpu compute, CPU-exact within 1/255) with CPU fallback | ✅ · WebGPU in the browser 🚧 | | Web build (same UI in the browser via WASM): persistent library in OPFS/IndexedDB, Web Worker rendering, export downloads | ✅ · WebGPU, Safari/Firefox testing 🚧 | diff --git a/docs/parity.md b/docs/parity.md index aa57982b..cd80f25d 100644 --- a/docs/parity.md +++ b/docs/parity.md @@ -35,7 +35,7 @@ and prints the summary; `cargo xtask parity --write` refreshes the summary table | L. Presets (PRE) | 1 | 4 | 2 | 1 | 0/2 (0%) | 1/2 (50%) | | M. Versions & history (VER) | 3 | 1 | 1 | 0 | 1/1 (100%) | 2/3 (67%) | | N. Copy / paste / sync (SYNC) | 4 | 0 | 1 | 0 | 3/3 (100%) | 1/1 (100%) | -| O. Merge (MERGE) | 0 | 0 | 4 | 0 | — | — | +| O. Merge (MERGE) | 3 | 0 | 1 | 0 | — | — | | P. Enhance (ENH) | 0 | 0 | 2 | 0 | — | — | | Q. HDR (HDR) | 0 | 0 | 5 | 0 | — | — | | R. Video (VID) | 0 | 0 | 5 | 1 | — | 0/2 (0%) | @@ -48,7 +48,7 @@ and prints the summary; `cargo xtask parity --write` refreshes the summary table | Y. Menus | 42 | 19 | 24 | 8 | 37/47 (79%) | 3/23 (13%) | | Z. Keyboard shortcuts (desktop) | 45 | 16 | 19 | 1 | 39/52 (75%) | 6/23 (26%) | | Lightroom Classic extras | 1 | 17 | 61 | 9 | — | 1/21 (5%) | -| **Total** | 187 | 107 | 206 | 37 | 150/198 (76%) | 34/144 (24%) | +| **Total** | 190 | 107 | 203 | 37 | 150/198 (76%) | 34/144 (24%) | ## Top gaps @@ -407,9 +407,9 @@ under the UI command that wraps it (a no-op error); the UI command now wins. | Id | Feature | Tier | Status | Evidence | Missing / notes | |---|---|---|---|---|---| -| LR-MERGE-HDR | HDR merge | P2 | ⬜ | | | -| LR-MERGE-PANO | Panorama | P2 | ⬜ | | | -| LR-MERGE-HDRPANO | HDR panorama | P2 | ⬜ | | | +| LR-MERGE-HDR | HDR merge | P2 | ✅ | `cmd:merge.hdr`, `cmd:dialog.mergeHdr` | auto align, deghost None–High + overlay, auto settings, float DNG; no Create Stack yet (no stacks), JPEG brackets treated as linear | +| LR-MERGE-PANO | Panorama | P2 | ✅ | `cmd:merge.panorama`, `cmd:dialog.mergePanorama` | spherical/cylindrical/perspective + auto, boundary warp, auto crop, fill edges (diffusion), DNG; no lens model / 360° wrap | +| LR-MERGE-HDRPANO | HDR panorama | P2 | ✅ | `cmd:merge.hdrPanorama`, `cmd:dialog.mergeHdrPanorama` | brackets grouped by EXIF | | LR-MERGE-HEADLESS | Merge with last settings | P2 | ⬜ | | | ## P. Enhance (ENH)