mirror of
https://github.com/storytold/lightcraft.git
synced 2026-10-24 10:09:01 -06:00
M13.1: lightcraft-merge crate — HDR merge, panorama stitching core
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 <noreply@anthropic.com>
This commit is contained in:
Brandon Thomas
co-authored by
Claude Opus 5.5
Claude Opus 5.5
parent
b13274f711
commit
543fdde41e
Generated
+16
@@ -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"
|
||||
|
||||
@@ -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" }
|
||||
|
||||
@@ -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 }
|
||||
@@ -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<f32> =
|
||||
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<Point> = 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<Point> = m.iter().map(|(i, _)| Point::new(ref_feat.points[*i].x as f64, ref_feat.points[*i].y as f64)).collect();
|
||||
let b: Vec<Point> = 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<bool>) {
|
||||
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));
|
||||
}
|
||||
}
|
||||
@@ -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<Keypoint>,
|
||||
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<Keypoint>, 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::<Vec<_>>(), 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<f32> = 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<Vec<(f32, f32, f32)>> = (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<usize> {
|
||||
if pts.len() <= max {
|
||||
return (0..pts.len()).collect();
|
||||
}
|
||||
let radii: Vec<f32> = (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<usize> = (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::<f32>() / DESC_LEN as f32;
|
||||
let var = d.iter().map(|v| (v - mean) * (v - mean)).sum::<f32>() / 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());
|
||||
}
|
||||
}
|
||||
@@ -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<ColorData>),
|
||||
/// 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<f64>,
|
||||
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<f64>) -> 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<f64> {
|
||||
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<usize>, orient: bool) -> Result<Frame> {
|
||||
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);
|
||||
}
|
||||
}
|
||||
@@ -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<f64>,
|
||||
/// The same from EXIF, when every frame had it.
|
||||
pub ev_exif: Option<Vec<f64>>,
|
||||
pub alignments: Vec<Alignment>,
|
||||
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<Frame>, opts: &HdrOptions, progress: &Progress) -> Result<HdrResult> {
|
||||
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<Vec<f64>> = frames.iter().map(|f| f.exposure).collect();
|
||||
let medians: Vec<f32> = frames.par_iter().map(|f| align::median_luma(&f.image, f.clip)).collect();
|
||||
let brightness: Vec<f64> = 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<usize> = (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<Option<Vec<bool>>> = vec![None; n];
|
||||
if opts.align {
|
||||
check(progress, 0.05, "Aligning")?;
|
||||
let g0 = 0.18 / medians[reference].max(1e-6);
|
||||
let rel: Vec<f32> = (0..n).map(|i| 2f64.powf(brightness[reference] - brightness[i]) as f32).collect();
|
||||
let rel: Vec<f32> = 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<bool>)> = (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<f32> = ln_k.iter().map(|v| v.exp() as f32).collect();
|
||||
let ev: Vec<f64> = 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<Plane>, 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<f32> = 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<f64> {
|
||||
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<f32> = (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<f32>) {
|
||||
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<Vec<bool>>],
|
||||
k: &[f32],
|
||||
order: &[usize],
|
||||
reference: usize,
|
||||
thr: f32,
|
||||
r: usize,
|
||||
) -> (Vec<Plane>, 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<f32>)> = (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<Plane> = (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
|
||||
}
|
||||
@@ -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<T> = std::result::Result<T, MergeError>;
|
||||
|
||||
/// 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) }
|
||||
}
|
||||
@@ -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<M3> {
|
||||
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<f64>, mut b: Vec<f64>, n: usize) -> Option<Vec<f64>> {
|
||||
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<f64>, Vec<Vec<f64>>) {
|
||||
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::<f64>().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<usize> = (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);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -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<Vec<u8>> {
|
||||
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());
|
||||
}
|
||||
}
|
||||
@@ -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<f32>,
|
||||
}
|
||||
|
||||
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<Buf>,
|
||||
wsum: Vec<Buf>,
|
||||
}
|
||||
|
||||
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<Buf> = 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<f32> = (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));
|
||||
}
|
||||
}
|
||||
@@ -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<u8>,
|
||||
}
|
||||
|
||||
/// Match every pair of images. `feats[k]` are in work coordinates scaled by `scales[k]`.
|
||||
pub fn match_pairs(feats: &[Features], scales: &[f64]) -> Vec<PairMatch> {
|
||||
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<PairMatch> = 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<Point> = m.iter().map(|(p, _)| Point::new(feats[i].points[*p].x as f64, feats[i].points[*p].y as f64)).collect();
|
||||
let b: Vec<Point> = 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<usize> {
|
||||
let mut parent: Vec<usize> = (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<usize, Vec<usize>> = 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<f64>, Option<f64>) {
|
||||
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<f64> {
|
||||
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<Option<Camera>> {
|
||||
let n = sizes.len();
|
||||
let mut cams: Vec<Option<Camera>> = vec![None; n];
|
||||
let deg = |k: usize| pairs.iter().filter(|p| p.i == k || p.j == k).map(|p| p.pts.len()).sum::<usize>();
|
||||
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<Option<Camera>>,
|
||||
/// 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<Camera>], pairs: &[PairMatch], anchor: usize, max_per_pair: usize) -> Adjusted {
|
||||
let n = cams.len();
|
||||
let idx: Vec<Option<usize>> = {
|
||||
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<usize> = (0..p.pts.len()).filter(|&k| p.levels.get(k).is_none_or(|l| *l <= 1)).collect();
|
||||
let use_idx: Vec<usize> = 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<Option<Camera>> = cams.to_vec();
|
||||
let sigma = 2.0f64;
|
||||
let cost = |cs: &[Option<Camera>]| -> (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<Camera>]) {
|
||||
let rs: Vec<M3> = 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::<f64>().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:?}");
|
||||
}
|
||||
}
|
||||
@@ -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<f64>, Vec<f64>)> {
|
||||
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<usize> = (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<f64> {
|
||||
(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::<f64>() / (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<bool> = (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<usize> = 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<f32> = 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));
|
||||
}
|
||||
}
|
||||
@@ -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<usize>,
|
||||
/// Cameras of the inputs (`None` = not used), in their full-resolution pixels, straightened.
|
||||
pub cameras: Vec<Option<Camera>>,
|
||||
/// Per-input exposure gains applied.
|
||||
pub gains: Vec<f64>,
|
||||
/// 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<Camera>]) -> (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<Camera>], 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<Camera>]) -> 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>) -> 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<f32>,
|
||||
}
|
||||
|
||||
/// Register the frames: cameras for the connected set.
|
||||
pub fn register(frames: &[Frame], progress: &Progress) -> Result<(Vec<Option<Camera>>, 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<f64> = 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<Frame>, opts: &PanoOptions, progress: &Progress) -> Result<PanoResult> {
|
||||
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<usize> = (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<f64> = (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<Option<Tile>> = (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<Tile>], 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::<f64>() / total;
|
||||
if mean <= 0.0 {
|
||||
return;
|
||||
}
|
||||
let (sn, sg) = (0.1f64, 1.0f64);
|
||||
let g0: Vec<f64> = 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);
|
||||
}
|
||||
}
|
||||
@@ -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());
|
||||
}
|
||||
}
|
||||
@@ -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<usize>,
|
||||
/// 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<Homography> {
|
||||
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::<f64>() / n, p.iter().map(|q| q.y).sum::<f64>() / 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::<f64>() / 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<Homography> {
|
||||
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<Fit> {
|
||||
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<usize> = 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<usize> = (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<Point> = best.iter().map(|&i| a[i]).collect();
|
||||
let ib: Vec<Point> = best.iter().map(|&i| b[i]).collect();
|
||||
let Some(hh) = fit_model(model, &ia, &ib) else { break };
|
||||
let inl: Vec<usize> = (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::<f64>() / 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::<Vec<_>>(), 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::<Vec<_>>()).unwrap();
|
||||
assert!((s.0[0]).abs() < 1e-9 && (s.0[1] + 2.0).abs() < 1e-9 && (s.0[3] - 2.0).abs() < 1e-9);
|
||||
}
|
||||
}
|
||||
@@ -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<f32> = 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<f32>, 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<Frame> = (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<Frame> {
|
||||
(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::<f32>() / 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::<f32>() / high.ghost.data.len() as f32;
|
||||
assert!(mean < 0.08, "the overlay is local (mean {mean})");
|
||||
}
|
||||
@@ -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<f32> = 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<M3>, Rgb32f) {
|
||||
let e = env();
|
||||
let rs: Vec<M3> = shots.iter().map(|s| rot(s.yaw, s.pitch, s.roll)).collect();
|
||||
let frames: Vec<Frame> = 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<f32> =
|
||||
(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<f32>], x: usize| -> Option<f32> { (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<Option<f32>> {
|
||||
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<Frame> = 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());
|
||||
}
|
||||
@@ -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<u8> {
|
||||
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<Vec<u8>> {
|
||||
}
|
||||
|
||||
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<Vec<u8>> = (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<Vec<u8>> {
|
||||
.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<Vec<u8>> = (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 {
|
||||
|
||||
@@ -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<f32> = (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);
|
||||
|
||||
Reference in New Issue
Block a user