Merge main (Photo Merge); M1.14b: stack merge results with their sources

Photo Merge's Create Stack option now puts the merged DNG on top of a
collapsed stack of its sources (Catalog::stack_with_ops), so the grid
shows the result; test in tests_merge (and no stack without the option).
Photo Merge dialogs live in a Photo > Photo Merge submenu.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
Brandon Thomas 2026-10-01 11:00:51 -04:00
commit ce3dd5e156
43 changed files with 5240 additions and 15 deletions

18
Cargo.lock generated
View File

@ -2055,6 +2055,7 @@ dependencies = [
"lightcraft-codecs",
"lightcraft-engine",
"lightcraft-mcp",
"lightcraft-merge",
"lightcraft-raster",
"lightcraft-ui-egui",
"serde_json",
@ -2117,6 +2118,7 @@ dependencies = [
"lightcraft-develop",
"lightcraft-geom",
"lightcraft-gpu",
"lightcraft-merge",
"lightcraft-meta",
"lightcraft-pipeline",
"lightcraft-preview",
@ -2165,6 +2167,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"

View File

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

View File

@ -238,6 +238,7 @@ Lightroom parity tracker — every feature, menu item and shortcut with its stat
| MCP server (headless or live app, persistent libraries), CLI, control channel | ✅ |
| XMP sidecars (read/write, auto-write), reading `crs:` develop settings, preset files (`.lcpreset`, XMP presets) | ✅ |
| Optics (distortion, vignetting, auto + manual CA, defringe, DNG-embedded lens corrections), Geometry (transforms, Constrain Crop), Upright (Auto/Level/Vertical/Full/Guided) | ✅ · camera lens profiles (our own) ⬜ |
| Photo Merge: HDR (auto-align, deghost), Panorama (spherical/cylindrical/perspective, boundary warp, auto crop), HDR Panorama → DNG | ✅ |
| GPU pipeline (wgpu compute, CPU-exact within 1/255) with CPU fallback | ✅ · WebGPU in the browser 🚧 |
| Web build (same UI in the browser via WASM): persistent library in OPFS/IndexedDB, Web Worker rendering, export downloads | ✅ · WebGPU, Safari/Firefox testing 🚧 |

View File

@ -19,4 +19,5 @@ lightcraft-mcp = { workspace = true }
lightcraft-ui-egui = { workspace = true }
lightcraft-codecs = { workspace = true }
lightcraft-raster = { workspace = true }
lightcraft-merge = { workspace = true }
serde_json = { workspace = true }

View File

@ -4,6 +4,8 @@
//! lightcraft-cli mcp [--connect [ADDR]] [--demo] [--compact] [FILES/FOLDERS…]
//! lightcraft-cli render <in> -o <out> [--set control=value]… [--settings FILE.json] [--preset ID] [--size N] [--quality Q]
//! lightcraft-cli snapshot [--library DIR | --demo] [--script FILE.jsonl] [-o OUT.png] [--size WxH] [--scale S] [FILES…]
//! lightcraft-cli merge hdr|panorama|hdr-panorama [OPTIONS] FILES…
//! lightcraft-cli synth-merge hdr|panorama -o DIR
//! lightcraft-cli commands [--json]
//! lightcraft-cli controls [--json]
//! ```
@ -46,6 +48,15 @@ USAGE:
-o, --output OUT PNG path (a final screenshot is written here if the script took none)
--size WxH window size in points (default 1600x1000)
--scale S pixels per point (default 1)
lightcraft-cli merge hdr|panorama|hdr-panorama [OPTIONS] FILES…
Photo Merge: writes <first>-HDR.dng / -Pano.dng / -HDR-Pano.dng next to the first file and
prints the result as JSON. Options:
--deghost none|low|medium|high --no-align --bracket N (HDR panorama)
--projection auto|spherical|cylindrical|perspective --boundary-warp 0..100
--auto-crop --fill-edges --no-auto-settings
--preview OUT.png only render a ≤ 1024 px preview (nothing written next to the files)
lightcraft-cli synth-merge hdr|panorama -o DIR
Write synthetic merge inputs (procedural scene; bracketed DNGs or overlapping PNG views).
lightcraft-cli commands [--json] list every command id with its parameters
lightcraft-cli controls [--json] list every develop control id with its range
lightcraft-cli --version | --help
@ -58,6 +69,8 @@ fn main() -> ExitCode {
Some("render") => render(&args[1..]),
Some("snapshot") => snapshot(&args[1..]),
Some("commands") => commands(&args[1..]),
Some("merge") => merge(&args[1..]),
Some("synth-merge") => synth_merge(&args[1..]),
Some("controls") => controls(&args[1..]),
Some("--version" | "-V" | "version") => {
println!("lightcraft-cli {}", env!("CARGO_PKG_VERSION"));
@ -148,6 +161,112 @@ fn mcp(args: &[String]) -> Result<(), String> {
server.serve(BufReader::new(stdin.lock()), std::io::stdout().lock()).map_err(|e| e.to_string())
}
fn merge(args: &[String]) -> Result<(), String> {
let cmd = match args.first().map(String::as_str) {
Some("hdr") => "merge.hdr",
Some("panorama" | "pano") => "merge.panorama",
Some("hdr-panorama") => "merge.hdrPanorama",
_ => return Err("merge: expected hdr, panorama or hdr-panorama".into()),
};
let mut p = serde_json::Map::new();
let mut files = Vec::new();
let mut i = 1;
while i < args.len() {
match args[i].as_str() {
"--deghost" => {
p.insert("deghost".into(), json!(take_value(args, &mut i, "--deghost")?));
}
"--no-align" => {
p.insert("align".into(), json!(false));
}
"--bracket" => {
p.insert("bracket".into(), json!(take_value(args, &mut i, "--bracket")?.parse::<u64>().map_err(|e| e.to_string())?));
}
"--projection" => {
p.insert("projection".into(), json!(take_value(args, &mut i, "--projection")?));
}
"--boundary-warp" => {
p.insert("boundaryWarp".into(), json!(take_value(args, &mut i, "--boundary-warp")?.parse::<f64>().map_err(|e| e.to_string())?));
}
"--auto-crop" => {
p.insert("autoCrop".into(), json!(true));
}
"--fill-edges" => {
p.insert("fillEdges".into(), json!(true));
}
"--no-auto-settings" => {
p.insert("autoSettings".into(), json!(false));
}
"--preview" => {
p.insert("preview".into(), json!(true));
p.insert("previewPath".into(), json!(take_value(args, &mut i, "--preview")?));
p.insert("showOverlay".into(), json!(true));
}
a if a.starts_with("--") => return Err(format!("unknown option `{a}`")),
f => files.push(f.to_string()),
}
i += 1;
}
let mut s = Session::new().with_fs();
let paths = expand_paths(&files);
let r = s.execute("library.import", &json!({"paths": paths})).map_err(|e| e.to_string())?;
let mut ids: Vec<u64> = r["imported"].as_array().map(|a| a.iter().filter_map(Value::as_u64).collect()).unwrap_or_default();
// keep the command line's order
ids.sort_by_key(|id| {
paths.iter().position(|f| {
s.catalog
.photo(lightcraft_engine::catalog::PhotoId(*id))
.is_some_and(|ph| matches!(&ph.source, lightcraft_engine::catalog::Source::File { path } if path == f))
})
});
p.insert("ids".into(), json!(ids));
let t0 = std::time::Instant::now();
let mut out = s.execute(cmd, &Value::Object(p)).map_err(|e| e.to_string())?;
out["seconds"] = json!(t0.elapsed().as_secs_f64());
println!("{}", serde_json::to_string_pretty(&out).unwrap_or_default());
Ok(())
}
fn synth_merge(args: &[String]) -> Result<(), String> {
let kind = args.first().map(String::as_str).unwrap_or("");
let mut dir = None;
let mut i = 1;
while i < args.len() {
match args[i].as_str() {
"-o" | "--output" => dir = Some(take_value(args, &mut i, "-o")?.to_string()),
a => return Err(format!("unknown option `{a}`")),
}
i += 1;
}
let dir = dir.ok_or("synth-merge needs -o DIR")?;
std::fs::create_dir_all(&dir).map_err(|e| e.to_string())?;
let mut written = Vec::new();
match kind {
"hdr" => {
for (k, b) in lightcraft_merge::synth::bracket_dngs(1800, 1200, &[-2.0, 0.0, 2.0]).into_iter().enumerate() {
let p = Path::new(&dir).join(format!("bracket-{k}.dng"));
std::fs::write(&p, b).map_err(|e| e.to_string())?;
written.push(p);
}
}
"panorama" | "pano" => {
for (k, v) in lightcraft_merge::synth::pano_views(1200, 900, 1000.0, &[-50.0, -25.0, 0.0, 25.0, 50.0]).into_iter().enumerate() {
let p = Path::new(&dir).join(format!("view-{k}.png"));
let img = v.to_srgb8();
let png = lightcraft_codecs::encode_png(&lightcraft_codecs::EncodeImage::rgba8(&img), &lightcraft_codecs::EncodeMeta::default())
.map_err(|e| e.to_string())?;
std::fs::write(&p, png).map_err(|e| e.to_string())?;
written.push(p);
}
}
_ => return Err("synth-merge: expected hdr or panorama".into()),
}
for p in written {
println!("{}", p.display());
}
Ok(())
}
fn take_value<'a>(args: &'a [String], i: &mut usize, flag: &str) -> Result<&'a str, String> {
*i += 1;
args.get(*i).map(String::as_str).ok_or_else(|| format!("{flag} needs a value"))

View File

@ -21,11 +21,12 @@ lightcraft-preview = { workspace = true }
lightcraft-codecs = { workspace = true }
lightcraft-raw = { workspace = true }
lightcraft-meta = { workspace = true }
lightcraft-merge = { workspace = true }
rayon = { workspace = true }
serde = { workspace = true }
serde_json = { workspace = true }
thiserror = { workspace = true }
log = { workspace = true }
rayon = { workspace = true }
ab_glyph = { workspace = true }
[dev-dependencies]

View File

@ -0,0 +1,72 @@
//! Photo Merge commands: `merge.hdr`, `merge.panorama`, `merge.hdrPanorama`.
//!
//! Run synchronously here (CLI, MCP, control channel); the desktop UI plans the same job with
//! [`Session::plan_merge`] and runs it on a worker thread. `preview: true` merges ≤ 1024 px frames
//! and (with `previewPath`) writes the rendered preview as PNG instead of adding a photo.
use lightcraft_catalog::PhotoId;
use serde_json::{Value, json};
use super::{CommandSpec, bad, bool_or, cmd, str_param};
use crate::{Result, Session};
fn can_merge(s: &Session) -> std::result::Result<(), String> {
if cfg!(target_arch = "wasm32") && s.media.file_bytes.is_none() {
return Err("Photo Merge needs access to the photo files".into());
}
// the photo count is checked when run (an explicit `ids` parameter may name them)
if s.catalog.photos().nth(1).is_none() { Err("the library has fewer than 2 photos".into()) } else { Ok(()) }
}
fn run(s: &mut Session, id: &str, p: &Value) -> Result<Value> {
let (kind, finish) = crate::merge::parse(id, p)?;
let ids: Vec<PhotoId> = s.targets(p);
if !bool_or(p, "preview", false) {
return s.merge_now(kind, finish, &ids);
}
let job = s.plan_merge(kind, finish, &ids, true)?;
let out = job.run(&|_, _| true).map_err(|e| bad(id, e))?;
let mut info = out.info;
if let (Some(path), Some(img)) = (str_param(p, "previewPath"), &out.preview) {
let png = lightcraft_codecs::encode_png(&lightcraft_codecs::EncodeImage::rgba8(img), &lightcraft_codecs::EncodeMeta::default())
.map_err(|e| bad(id, e.to_string()))?;
std::fs::write(path, png).map_err(|e| bad(id, format!("{path}: {e}")))?;
info["previewPath"] = json!(path);
}
if let Some(img) = &out.preview {
info["previewSize"] = json!([img.width, img.height]);
}
Ok(info)
}
pub fn specs() -> Vec<CommandSpec> {
vec![
cmd!(
"merge.hdr",
"HDR Merge",
[],
None,
"{ids?, align=true, deghost=none|low|medium|high, autoSettings=true, stack=false, preview=false, showOverlay=false, previewPath?}",
can_merge,
|s, p| run(s, "merge.hdr", p)
),
cmd!(
"merge.panorama",
"Panorama Merge",
[],
None,
"{ids?, projection=auto|spherical|cylindrical|perspective, boundaryWarp=0..100, autoCrop=false, fillEdges=false, autoSettings=true, stack=false, maxMegapixels=40, preview=false, previewPath?}",
can_merge,
|s, p| run(s, "merge.panorama", p)
),
cmd!(
"merge.hdrPanorama",
"HDR Panorama Merge",
[],
None,
"{ids?, bracket=0 (auto from EXIF), align, deghost, projection, boundaryWarp, autoCrop, fillEdges, autoSettings, stack, preview, previewPath?}",
can_merge,
|s, p| run(s, "merge.hdrPanorama", p)
),
]
}

View File

@ -8,6 +8,7 @@ mod develop;
mod edit;
mod library;
mod masks;
mod merge;
mod organize;
mod preset_files;
mod query;
@ -104,6 +105,7 @@ pub fn command_specs() -> &'static [CommandSpec] {
v.extend(develop::specs());
v.extend(masks::specs());
v.extend(organize::specs());
v.extend(merge::specs());
v.extend(query::specs());
v.extend(xmp::specs());
v.extend(preset_files::specs());

View File

@ -18,6 +18,7 @@ pub mod files;
pub mod import;
pub mod library;
pub mod media;
pub mod merge;
pub mod presets;
pub mod sidecar;
mod view;
@ -339,6 +340,8 @@ mod tests_import;
#[cfg(test)]
mod tests_library;
#[cfg(test)]
mod tests_merge;
#[cfg(test)]
mod tests_organize;
#[cfg(test)]
mod tests_xmp;

View File

@ -102,6 +102,8 @@ pub struct MediaCache {
pub file_loader: Option<FileLoader>,
pub file_probe: Option<FileProbe>,
pub preview_loader: Option<PreviewLoader>,
/// Reads a photo file's bytes (Photo Merge); `None` = the local file system.
pub file_bytes: Option<crate::merge::ByteReader>,
scenes: Vec<lightcraft_scenes::Scene>,
/// Rendered thumbnails (memory, plus disk once a library is attached).
pub rendered: Arc<PreviewCache>,
@ -117,6 +119,7 @@ impl Default for MediaCache {
file_loader: None,
file_probe: None,
preview_loader: None,
file_bytes: None,
scenes: Vec::new(),
rendered: Arc::new(PreviewCache::memory(RENDERED_MEM_BYTES)),
}

413
crates/engine/src/merge.rs Normal file
View File

@ -0,0 +1,413 @@
//! Photo Merge (HDR, Panorama, HDR Panorama) in the library.
//!
//! A merge is planned from the session ([`crate::Session::plan_merge`] → [`MergeJob`], which is
//! `Send` and self-contained), run anywhere (a worker thread in the app, inline for the CLI/MCP:
//! [`MergeJob::run`], with progress and cancellation), and finished on the session
//! ([`crate::Session::finish_merge`]): the DNG is written next to the first source photo
//! (`<name>-HDR.dng`, `-Pano.dng`, `-HDR-Pano.dng`), imported, selected, and given Auto Settings /
//! the auto-crop rectangle when asked. A *preview* run works on ≤ 1024 px frames and returns a
//! rendered image (with the deghost overlay when asked) instead of a file.
use std::sync::Arc;
use lightcraft_catalog::{PhotoId, Source};
use lightcraft_develop::DevelopSettings;
use lightcraft_merge::{Deghost, Frame, HdrOptions, PanoOptions, Projection};
use lightcraft_pipeline::{Quality, RenderRequest, SourceInfo};
use lightcraft_raster::{Rgb32f, Rgba8};
use serde::{Deserialize, Serialize};
use serde_json::{Value, json};
use crate::{EngineError, Result, Session};
/// Long edge of the frames of a preview merge.
pub const PREVIEW_EDGE: usize = 1024;
#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
#[serde(tag = "kind", rename_all = "camelCase")]
pub enum MergeKind {
Hdr {
hdr: HdrOptions,
},
Panorama {
pano: PanoOptions,
},
/// Merge each bracket (`bracket` consecutive photos; 0 = detect from EXIF), then stitch.
HdrPanorama {
hdr: HdrOptions,
pano: PanoOptions,
bracket: usize,
},
}
impl MergeKind {
fn suffix(&self) -> &'static str {
match self {
MergeKind::Hdr { .. } => "HDR",
MergeKind::Panorama { .. } => "Pano",
MergeKind::HdrPanorama { .. } => "HDR-Pano",
}
}
fn min_photos(&self) -> usize {
match self {
MergeKind::HdrPanorama { .. } => 4,
_ => 2,
}
}
}
/// Options that apply after the merge.
#[derive(Clone, Debug, Default, PartialEq, Serialize, Deserialize)]
#[serde(default, rename_all = "camelCase")]
pub struct MergeFinish {
/// Run Auto Settings on the result.
pub auto_settings: bool,
/// Stack the result on top of its sources (collapsed); the result's caption also lists them.
pub stack: bool,
/// Preview only: tint deghosted areas red.
pub show_overlay: bool,
}
/// Reads a source photo's bytes (by catalog path).
pub type ByteReader = Arc<dyn Fn(&str) -> std::result::Result<Vec<u8>, String> + Send + Sync>;
/// A planned merge, detached from the session.
#[derive(Clone)]
pub struct MergeJob {
pub kind: MergeKind,
pub finish: MergeFinish,
pub sources: Vec<(PhotoId, String)>,
/// `Some(edge)`: preview at that frame size; `None`: the full-resolution merge.
pub preview: Option<usize>,
pub read: ByteReader,
}
/// What a merge produced.
pub struct MergeOutput {
/// The DNG (also for previews, at preview size).
pub dng: Vec<u8>,
/// Rendered preview (previews only).
pub preview: Option<Rgba8>,
/// Auto-crop rectangle (normalised x, y, w, h) for panoramas.
pub crop: Option<[f64; 4]>,
/// Statistics for the caller (exposures, alignment, projection, …).
pub info: Value,
}
pub type ProgressFn<'a> = dyn Fn(f32, &str) -> bool + Sync + 'a;
fn err(e: impl std::fmt::Display) -> String {
e.to_string()
}
impl MergeJob {
fn load(&self, orient: bool, progress: &ProgressFn, span: (f32, f32)) -> std::result::Result<Vec<Frame>, String> {
use rayon::prelude::*;
let done = std::sync::atomic::AtomicUsize::new(0);
let n = self.sources.len();
let frames: std::result::Result<Vec<Frame>, String> = self
.sources
.par_iter()
.map(|(_, path)| {
let bytes = (self.read)(path)?;
let f = lightcraft_merge::load_frame(&bytes, self.preview, orient).map_err(|e| format!("{path}: {e}"));
let k = done.fetch_add(1, std::sync::atomic::Ordering::Relaxed) + 1;
progress(span.0 + (span.1 - span.0) * k as f32 / n as f32, "Reading photos");
f
})
.collect();
frames
}
/// Run the merge. `progress(fraction, stage)` returns false to cancel.
pub fn run(&self, progress: &ProgressFn) -> std::result::Result<MergeOutput, String> {
if self.sources.len() < self.kind.min_photos() {
return Err(format!("select at least {} photos", self.kind.min_photos()));
}
if !progress(0.0, "Reading photos") {
return Err("cancelled".into());
}
let sub = |a: f32, b: f32| move |f: f32, s: &str| progress(a + (b - a) * f, s);
let (image, color, orientation, metadata, baseline, samples, crop, mut info, ghost) = match &self.kind {
MergeKind::Hdr { hdr } => {
let frames = self.load(false, progress, (0.0, 0.3))?;
let p = sub(0.3, 0.9);
let r = lightcraft_merge::merge_hdr(frames, hdr, &p).map_err(err)?;
let info = json!({
"reference": self.sources[r.reference].0.0,
"ev": r.ev, "evExif": r.ev_exif,
"alignment": r.alignments.iter().map(|a| json!({"model": a.model, "inliers": a.inliers, "rmsPx": a.rms})).collect::<Vec<_>>(),
});
let ghost = (!r.ghost.data.is_empty()).then(|| r.ghost.oriented(r.orientation));
(r.radiance, r.color, r.orientation, r.metadata, r.baseline_exposure, lightcraft_merge::DngSamples::Half, None, info, ghost)
}
MergeKind::Panorama { pano } => {
let frames = self.load(true, progress, (0.0, 0.25))?;
let p = sub(0.25, 0.9);
let r = lightcraft_merge::stitch(frames, pano, &p).map_err(err)?;
let info = pano_info(&r, &self.sources);
(
r.image,
r.color,
lightcraft_geom::Orientation::Normal,
r.metadata,
r.baseline_exposure,
lightcraft_merge::DngSamples::U16,
r.crop,
info,
None,
)
}
MergeKind::HdrPanorama { hdr, pano, bracket } => {
let frames = self.load(false, progress, (0.0, 0.2))?;
let b =
if *bracket >= 2 { *bracket } else { detect_bracket(&frames).ok_or("could not detect the bracket size from EXIF; set it")? };
if !frames.len().is_multiple_of(b) || frames.len() / b < 2 {
return Err(format!("{} photos don't form brackets of {b} for at least two positions", frames.len()));
}
let groups = frames.len() / b;
let mut merged = Vec::new();
let mut it = frames.into_iter();
for g in 0..groups {
let set: Vec<Frame> = it.by_ref().take(b).collect();
let p = sub(0.2 + 0.4 * g as f32 / groups as f32, 0.2 + 0.4 * (g + 1) as f32 / groups as f32);
let r = lightcraft_merge::merge_hdr(set, hdr, &p).map_err(err)?;
let mut f = Frame {
image: r.radiance,
clip: f32::INFINITY,
exposure: None,
raw: r.raw,
color: r.color,
orientation: r.orientation,
metadata: r.metadata,
baseline_exposure: r.baseline_exposure,
};
f.orient();
merged.push(f);
}
let p = sub(0.6, 0.9);
let r = lightcraft_merge::stitch(merged, pano, &p).map_err(err)?;
let mut info = pano_info(&r, &self.sources);
info["bracket"] = json!(b);
(
r.image,
r.color,
lightcraft_geom::Orientation::Normal,
r.metadata,
r.baseline_exposure,
lightcraft_merge::DngSamples::Half,
r.crop,
info,
None,
)
}
};
if !progress(0.9, "Writing DNG") {
return Err("cancelled".into());
}
let mut metadata = metadata;
if self.finish.stack {
let names: Vec<&str> =
self.sources.iter().map(|(_, p)| std::path::Path::new(p).file_name().and_then(|n| n.to_str()).unwrap_or(p)).collect();
metadata.caption = Some(format!("Merged from {}", names.join(", ")));
}
metadata.software = Some("LightCraft Photo Merge".into());
let dng = lightcraft_merge::write_linear_dng(&image, &color, orientation, &metadata, baseline, samples).map_err(err)?;
info["width"] = json!(image.width);
info["height"] = json!(image.height);
let preview = match self.preview {
Some(edge) => Some(render_preview(&dng, edge, self.finish.auto_settings, crop, ghost.as_ref().filter(|_| self.finish.show_overlay))?),
None => None,
};
progress(1.0, "Done");
Ok(MergeOutput { dng, preview, crop, info })
}
}
fn pano_info(r: &lightcraft_merge::PanoResult, sources: &[(PhotoId, String)]) -> Value {
json!({
"projection": r.projection,
"used": r.used.iter().filter_map(|&i| sources.get(i).map(|s| s.0.0)).collect::<Vec<_>>(),
"fov": [r.fov.0, r.fov.1],
"rmsPx": r.rms,
"gains": r.gains,
"crop": r.crop,
})
}
/// Bracket size from EXIF: the smallest period of the exposure sequence (≥ 2).
pub fn detect_bracket(frames: &[Frame]) -> Option<usize> {
let ev: Vec<f64> = frames.iter().map(|f| f.exposure).collect::<Option<_>>()?;
let n = ev.len();
(2..=n / 2).find(|&b| n.is_multiple_of(b) && (0..n).all(|k| (ev[k] - ev[k % b]).abs() < 0.2) && (1..b).any(|k| (ev[k] - ev[0]).abs() > 0.5))
}
/// Render a merged DNG like the library will show it (default raw settings, optionally Auto
/// Settings and the auto crop), with the deghost overlay tinted red.
fn render_preview(
dng: &[u8],
edge: usize,
auto: bool,
crop: Option<[f64; 4]>,
ghost: Option<&lightcraft_raster::Plane>,
) -> std::result::Result<Rgba8, String> {
let (src, info) = crate::files::load_bytes(dng, edge)?;
let mut s = DevelopSettings::for_raw(info.as_shot_temp, info.as_shot_tint);
if auto {
apply_auto(&mut s, &src, &info);
}
if let Some(c) = crop {
s.crop.geometry.rect = lightcraft_geom::Rect::from_xywh(c[0], c[1], c[2], c[3]);
}
let req = RenderRequest { max_w: edge, max_h: edge, quality: Quality::Full, apply_crop: false };
let mut img = lightcraft_pipeline::render(&src, &info, &s, &req).image;
if let Some(g) = ghost {
let (w, h) = (img.width, img.height);
for y in 0..h {
for x in 0..w {
let gx = ((x as f32 + 0.5) / w as f32 * g.width as f32) as usize;
let gy = ((y as f32 + 0.5) / h as f32 * g.height as f32) as usize;
let a = g.get(gx.min(g.width - 1), gy.min(g.height - 1)).clamp(0.0, 1.0) * 0.6;
if a > 0.01 {
let p = &mut img.data[y * w + x];
p[0] = (p[0] as f32 * (1.0 - a) + 255.0 * a) as u8;
p[1] = (p[1] as f32 * (1.0 - a)) as u8;
p[2] = (p[2] as f32 * (1.0 - a)) as u8;
}
}
}
}
if let Some(c) = crop {
// draw the auto-crop frame
let (w, h) = (img.width as f64, img.height as f64);
let (x0, y0, x1, y1) = ((c[0] * w) as usize, (c[1] * h) as usize, ((c[0] + c[2]) * w) as usize, ((c[1] + c[3]) * h) as usize);
for y in y0..y1.min(img.height) {
for x in x0..x1.min(img.width) {
if y == y0 || y + 1 == y1 || x == x0 || x + 1 == x1 {
img.data[y * img.width + x] = [255, 255, 255, 255];
}
}
}
}
Ok(img)
}
fn apply_auto(s: &mut DevelopSettings, src: &Rgb32f, info: &SourceInfo) {
let a = lightcraft_pipeline::auto::auto_tone(src, info, s);
s.light.exposure = a.exposure;
s.light.contrast = a.contrast;
s.light.highlights = a.highlights;
s.light.shadows = a.shadows;
s.light.whites = a.whites;
s.light.blacks = a.blacks;
s.color.vibrance = a.vibrance;
s.color.saturation = a.saturation;
}
/// A free file name `dir/<stem>-<suffix>.dng` (`-2`, `-3`… when taken).
fn output_path(first: &str, suffix: &str) -> String {
let p = std::path::Path::new(first);
let dir = p.parent().map(|d| d.to_path_buf()).unwrap_or_default();
let stem = p.file_stem().map(|s| s.to_string_lossy().to_string()).unwrap_or_else(|| "Merge".into());
for k in 1.. {
let name = if k == 1 { format!("{stem}-{suffix}.dng") } else { format!("{stem}-{suffix}-{k}.dng") };
let cand = dir.join(name);
if !cand.exists() {
return cand.to_string_lossy().to_string();
}
}
unreachable!()
}
/// Parse `merge.*` command parameters into a kind and finishing options.
pub fn parse(id: &str, p: &Value) -> Result<(MergeKind, MergeFinish)> {
let bad = |m: String| EngineError::BadParams { cmd: id.into(), msg: m };
let deghost: Deghost = serde_json::from_value(p.get("deghost").cloned().unwrap_or(json!("none"))).map_err(|e| bad(format!("deghost: {e}")))?;
let hdr = HdrOptions { align: p.get("align").and_then(Value::as_bool).unwrap_or(true), deghost };
let projection: Projection =
serde_json::from_value(p.get("projection").cloned().unwrap_or(json!("auto"))).map_err(|e| bad(format!("projection: {e}")))?;
let pano = PanoOptions {
projection,
boundary_warp: p.get("boundaryWarp").and_then(Value::as_f64).unwrap_or(0.0).clamp(0.0, 100.0),
auto_crop: p.get("autoCrop").and_then(Value::as_bool).unwrap_or(false),
fill_edges: p.get("fillEdges").and_then(Value::as_bool).unwrap_or(false),
max_megapixels: p.get("maxMegapixels").and_then(Value::as_f64).unwrap_or(40.0),
};
let kind = match id {
"merge.hdr" => MergeKind::Hdr { hdr },
"merge.panorama" => MergeKind::Panorama { pano },
_ => MergeKind::HdrPanorama { hdr, pano, bracket: p.get("bracket").and_then(Value::as_u64).unwrap_or(0) as usize },
};
let finish = MergeFinish {
auto_settings: p.get("autoSettings").and_then(Value::as_bool).unwrap_or(true),
stack: p.get("stack").and_then(Value::as_bool).unwrap_or(false),
show_overlay: p.get("showOverlay").and_then(Value::as_bool).unwrap_or(false),
};
Ok((kind, finish))
}
impl Session {
/// Plan a merge of `ids` (catalog order of the given ids is kept).
pub fn plan_merge(&self, kind: MergeKind, finish: MergeFinish, ids: &[PhotoId], preview: bool) -> Result<MergeJob> {
let mut sources = Vec::new();
for id in ids {
let p = self.catalog.photo(*id).ok_or(lightcraft_catalog::CatalogError::NoPhoto(*id))?;
match &p.source {
Source::File { path } => sources.push((*id, path.clone())),
Source::Demo { .. } => return Err(EngineError::Other(format!("{} is a demo photo; Photo Merge needs photo files", p.file_name))),
}
}
if sources.len() < kind.min_photos() {
return Err(EngineError::Other(format!("select at least {} photos to merge", kind.min_photos())));
}
let read: ByteReader = match &self.media.file_bytes {
Some(r) => r.clone(),
None => Arc::new(|path: &str| std::fs::read(path).map_err(|e| format!("{path}: {e}"))),
};
Ok(MergeJob { kind, finish, sources, preview: preview.then_some(PREVIEW_EDGE), read })
}
/// Write a finished (non-preview) merge next to its first source, import it, select it and
/// apply Auto Settings / the auto crop. Returns `{id, path, …info}`.
pub fn finish_merge(&mut self, job: &MergeJob, out: MergeOutput) -> Result<Value> {
let first = &job.sources.first().ok_or_else(|| EngineError::Other("nothing merged".into()))?.1;
let path = output_path(first, job.kind.suffix());
std::fs::write(&path, &out.dng).map_err(|e| EngineError::Other(format!("{path}: {e}")))?;
let report = crate::import::import(self, std::slice::from_ref(&path), crate::import::ImportMode::Add)?;
let id = report
.imported
.first()
.copied()
.map(PhotoId)
.ok_or_else(|| EngineError::Other(format!("the merged file {path} could not be imported: {:?}", report.failed)))?;
if job.finish.stack {
let sources: Vec<PhotoId> = job.sources.iter().map(|(p, _)| *p).collect();
if let Some(op) = self.catalog.stack_with_ops(id, &sources) {
self.commit("Stack with Sources", op)?;
}
}
self.selection = crate::Selection::single(id);
if job.finish.auto_settings {
self.execute("develop.auto", &json!({}))?;
}
if let Some(c) = out.crop
&& let Some(d) = self.develop_of(id)
{
let mut d = (*d).clone();
d.crop.geometry.rect = lightcraft_geom::Rect::from_xywh(c[0], c[1], c[2], c[3]);
self.set_develop(id, d, "Auto Crop")?;
}
let mut info = out.info;
info["id"] = json!(id.0);
info["path"] = json!(path);
info["bytes"] = json!(out.dng.len());
Ok(info)
}
/// Plan, run and finish a merge synchronously (CLI, MCP, control channel).
pub fn merge_now(&mut self, kind: MergeKind, finish: MergeFinish, ids: &[PhotoId]) -> Result<Value> {
let job = self.plan_merge(kind, finish, ids, false)?;
let out = job.run(&|_, _| true).map_err(EngineError::Other)?;
self.finish_merge(&job, out)
}
}

View File

@ -0,0 +1,90 @@
//! Photo Merge through the command layer: synthetic brackets / panorama views on disk → merged
//! DNG next to the sources → imported, selected, Auto Settings / auto crop applied.
use serde_json::json;
use crate::Session;
fn temp_dir(tag: &str) -> std::path::PathBuf {
let d = std::env::temp_dir().join(format!("lc-merge-{tag}-{}", std::process::id()));
let _ = std::fs::remove_dir_all(&d);
std::fs::create_dir_all(&d).unwrap();
d
}
fn import(s: &mut Session, paths: Vec<String>) -> Vec<u64> {
let r = s.execute("library.import", &json!({"paths": paths})).unwrap();
r["imported"].as_array().unwrap().iter().map(|v| v.as_u64().unwrap()).collect()
}
#[test]
fn hdr_merge_command_creates_and_imports_a_dng() {
let dir = temp_dir("hdr");
let mut paths = Vec::new();
for (i, b) in lightcraft_merge::synth::bracket_dngs(480, 320, &[-2.0, 0.0, 2.0]).into_iter().enumerate() {
let p = dir.join(format!("IMG_{i}.dng"));
std::fs::write(&p, b).unwrap();
paths.push(p.to_string_lossy().to_string());
}
let mut s = Session::new().with_fs();
let ids = import(&mut s, paths);
assert_eq!(ids.len(), 3);
// preview: rendered PNG, nothing added
let png = dir.join("preview.png");
let pv = s.execute("merge.hdr", &json!({"ids": ids, "deghost": "medium", "preview": true, "showOverlay": true, "previewPath": png})).unwrap();
assert!(png.exists() && pv["previewSize"][0].as_u64().unwrap() > 100, "{pv}");
assert_eq!(s.catalog.photos().count(), 3);
let ev: Vec<f64> = pv["ev"].as_array().unwrap().iter().map(|v| v.as_f64().unwrap()).collect();
assert!((ev[0] + 2.0).abs() < 0.1 && (ev[2] - 2.0).abs() < 0.1, "{ev:?}");
// full merge
let r = s.execute("merge.hdr", &json!({"ids": ids, "autoSettings": true, "stack": true})).unwrap();
let path = r["path"].as_str().unwrap();
assert!(path.ends_with("IMG_0-HDR.dng"), "{path}");
let id = lightcraft_catalog::PhotoId(r["id"].as_u64().unwrap());
let p = s.catalog.photo(id).unwrap().clone();
assert_eq!((p.width, p.height), (480, 320));
assert_eq!(p.kind, lightcraft_catalog::MediaKind::Raw);
assert_eq!(s.selection.active, Some(id));
let d = s.develop_of(id).unwrap();
assert!(d.light.exposure != 0.0 || d.light.highlights != 0.0, "auto settings applied");
let img = s.render_now(id, 256, 256).unwrap();
assert_eq!(img.image.width, 256);
// Create Stack: the result tops a collapsed stack of its sources, so the grid shows the result
let st = s.catalog.stack_of(id).expect("merge result stacked").clone();
assert_eq!(st.top(), id);
assert_eq!(st.photos[1..].iter().map(|p| p.0).collect::<Vec<_>>(), ids);
assert!(st.collapsed);
assert_eq!(s.visible_cloned(), vec![id]);
// the second merge doesn't overwrite the first
let r2 = s.execute("merge.hdr", &json!({"ids": ids, "autoSettings": false})).unwrap();
assert!(r2["path"].as_str().unwrap().ends_with("IMG_0-HDR-2.dng"));
let id2 = lightcraft_catalog::PhotoId(r2["id"].as_u64().unwrap());
assert!(s.catalog.stack_of(id2).is_none(), "no stack without the option");
let _ = std::fs::remove_dir_all(&dir);
}
#[test]
fn panorama_command_with_auto_crop() {
let dir = temp_dir("pano");
let mut paths = Vec::new();
for (i, v) in lightcraft_merge::synth::pano_views(480, 360, 420.0, &[-25.0, 0.0, 25.0]).into_iter().enumerate() {
let img = v.to_srgb8();
let png = lightcraft_codecs::encode_png(&lightcraft_codecs::EncodeImage::rgba8(&img), &lightcraft_codecs::EncodeMeta::default()).unwrap();
let p = dir.join(format!("P_{i}.png"));
std::fs::write(&p, png).unwrap();
paths.push(p.to_string_lossy().to_string());
}
let mut s = Session::new().with_fs();
let ids = import(&mut s, paths);
let r = s.execute("merge.panorama", &json!({"ids": ids, "projection": "cylindrical", "autoCrop": true, "autoSettings": false})).unwrap();
assert_eq!(r["projection"], "cylindrical");
assert_eq!(r["used"].as_array().unwrap().len(), 3);
let id = lightcraft_catalog::PhotoId(r["id"].as_u64().unwrap());
let p = s.catalog.photo(id).unwrap().clone();
assert!(p.width > 800 && p.height > 300, "{}×{}", p.width, p.height);
let crop = s.develop_of(id).unwrap().crop.geometry.rect;
assert!(crop.width() < 1.0 && crop.width() * crop.height() > 0.5, "{crop:?}");
// too few photos
assert!(s.execute("merge.panorama", &json!({"ids": [ids[0]]})).is_err());
let _ = std::fs::remove_dir_all(&dir);
}

21
crates/merge/Cargo.toml Normal file
View File

@ -0,0 +1,21 @@
[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 }
lightcraft-scenes = { workspace = true }
serde = { workspace = true }
thiserror = { workspace = true }
rayon = { workspace = true }

142
crates/merge/src/align.rs Normal file
View File

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

View File

@ -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());
}
}

200
crates/merge/src/frame.rs Normal file
View File

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

449
crates/merge/src/hdr.rs Normal file
View File

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

58
crates/merge/src/lib.rs Normal file
View File

@ -0,0 +1,58 @@
//! 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 mod synth;
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) }
}

247
crates/merge/src/linalg.rs Normal file
View File

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

142
crates/merge/src/output.rs Normal file
View File

@ -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());
}
}

View File

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

View File

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

View File

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

View File

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

View File

@ -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());
}
}

291
crates/merge/src/ransac.rs Normal file
View File

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

81
crates/merge/src/synth.rs Normal file
View File

@ -0,0 +1,81 @@
//! Synthetic merge inputs from procedural scenes (tests, demos, UI checks — no media files):
//! exposure brackets as linear DNGs and overlapping panorama views as PNGs.
use lightcraft_geom::Orientation;
use lightcraft_meta::Metadata;
use lightcraft_raster::Rgb32f;
use crate::frame::FrameColor;
use crate::linalg::{self, M3};
use crate::output::{DngSamples, write_linear_dng};
fn scene(kind: lightcraft_scenes::Kind) -> lightcraft_scenes::Scene {
lightcraft_scenes::demo_library().into_iter().find(|s| s.kind == kind).expect("scene kind in the demo library")
}
fn srgb_color() -> FrameColor {
FrameColor::Linear { to_xyz_d50: lightcraft_color::bradford(lightcraft_color::D65, lightcraft_color::D50).mul(&lightcraft_color::SRGB.to_xyz()) }
}
/// Exposure brackets (`evs`, e.g. `[-2, 0, 2]`) of a sunset scene as 16-bit linear DNGs, each
/// shifted by a few pixels (handheld) with EXIF exposure times; `w × h` pixels.
pub fn bracket_dngs(w: usize, h: usize, evs: &[f64]) -> Vec<Vec<u8>> {
let m = 12usize;
let truth = scene(lightcraft_scenes::Kind::OceanSunset).render(w + 2 * m, h + 2 * m);
let mut l: Vec<f32> = truth.data.iter().map(|p| p[1]).collect();
l.sort_by(|a, b| a.total_cmp(b));
let s = 0.18 / l[l.len() / 2].max(1e-6);
let to_srgb = lightcraft_color::REC2020.to_space(&lightcraft_color::SRGB).to_f32();
evs.iter()
.enumerate()
.map(|(i, ev)| {
let k = 2f32.powf(*ev as f32) * s;
let (dx, dy) = ((i as f32 * 2.3) % 5.0 - 2.0, (i as f32 * 1.7) % 4.0 - 1.5);
let img = Rgb32f::from_fn(w, h, |x, y| {
let p = truth.sample_bilinear(x as f32 + 0.5 + m as f32 - dx, y as f32 + 0.5 + m as f32 - dy);
let q = [0, 1, 2].map(|r| to_srgb[r][0] * p[0] + to_srgb[r][1] * p[1] + to_srgb[r][2] * p[2]);
q.map(|v| (v * k).clamp(0.0, 1.0))
});
let meta = Metadata {
make: Some("LightCraft".into()),
model: Some("Synthetic Bracket".into()),
exposure_time: Some(2f64.powf(*ev) / 125.0),
f_number: Some(8.0),
iso: Some(100),
..Default::default()
};
write_linear_dng(&img, &srgb_color(), Orientation::Normal, &meta, 0.0, DngSamples::U16).expect("encode")
})
.collect()
}
/// `yaws.len()` overlapping views (degrees of yaw, small pitch/roll wobble) of a canyon scene used
/// as a 200° × 70° spherical environment, `w × h` pixels with focal length `f`, as linear sRGB.
pub fn pano_views(w: usize, h: usize, f: f64, yaws: &[f64]) -> Vec<Rgb32f> {
let env = scene(lightcraft_scenes::Kind::Canyon).render(2400, 840);
let mut l: Vec<f32> = env.data.iter().map(|p| p[1]).collect();
l.sort_by(|a, b| a.total_cmp(b));
let s = 0.25 / l[l.len() / 2].max(1e-6);
let (theta, phi) = (200.0f64, 70.0f64);
yaws.iter()
.enumerate()
.map(|(i, yaw)| {
let wob = (i as f64 * 1.3).sin();
let ry = linalg::rodrigues([0.0, yaw.to_radians(), 0.0]);
let rx = linalg::rodrigues([(wob * 1.5).to_radians(), 0.0, 0.0]);
let rz = linalg::rodrigues([0.0, 0.0, (wob * 0.8).to_radians()]);
let c2w: M3 = linalg::mul3(&linalg::mul3(&ry, &rx), &rz);
Rgb32f::from_fn(w, h, |x, y| {
let d = linalg::apply3(&c2w, [(x as f64 + 0.5 - w as f64 / 2.0) / f, (y as f64 + 0.5 - h as f64 / 2.0) / f, 1.0]);
let t = d[0].atan2(d[2]).to_degrees();
let p = (d[1] / linalg::norm(d)).asin().to_degrees();
if t.abs() > theta / 2.0 || p.abs() > phi / 2.0 {
return [0.0; 3];
}
let px = (t / theta + 0.5) * env.width as f64;
let py = (p / phi + 0.5) * env.height as f64;
env.sample_bilinear(px as f32, py as f32).map(|v| (v * s).min(1.0))
})
})
.collect()
}

163
crates/merge/tests/hdr.rs Normal file
View File

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

247
crates/merge/tests/pano.rs Normal file
View File

@ -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());
}

View File

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

View File

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

View File

@ -77,7 +77,7 @@ pub fn inspect(app: &LightcraftApp, ctx: &egui::Context) -> Value {
"selection": app.session.selection.ids.iter().map(|p| p.0).collect::<Vec<_>>(),
"activeMask": app.session.active_mask,
"widgetCount": app.widgets.len(),
"perf": {"frameMs": app.perf.frame_ms, "fps": app.perf.fps, "lastRenderMs": app.renderer.last_main_ms, "renderQueue": app.renderer.queued(), "rendersInFlight": app.renderer.in_flight(), "rendersDone": app.renderer.completed, "thumbTextures": app.renderer.thumb_textures(), "gpu": (lightcraft_engine::gpu::ready() && lightcraft_engine::gpu::available()).then(lightcraft_engine::gpu::adapter_name).flatten()},
"perf": {"frameMs": app.perf.frame_ms, "fps": app.perf.fps, "lastRenderMs": app.renderer.last_main_ms, "renderQueue": app.renderer.queued(), "rendersInFlight": app.renderer.in_flight(), "mergeRunning": app.merge.busy(), "lastMerge": app.merge.last_result, "rendersDone": app.renderer.completed, "thumbTextures": app.renderer.thumb_textures(), "gpu": (lightcraft_engine::gpu::ready() && lightcraft_engine::gpu::available()).then(lightcraft_engine::gpu::adapter_name).flatten()},
"loupe": app.loupe_shown.map(|(p, src)| json!({"photo": p.0, "source": src, "pending": app.renderer.is_pending(crate::render::Slot::Main)})),
"status": app.ui.status,
})

View File

@ -171,7 +171,7 @@ impl Headless {
/// Is anything still in progress (renders, queued input)?
pub fn busy(&self) -> bool {
self.app.renderer.in_flight() > 0 || !self.app.synthetic.is_empty() || !self.events.is_empty()
self.app.renderer.in_flight() > 0 || self.app.merge.busy() || !self.app.synthetic.is_empty() || !self.events.is_empty()
}
/// Run frames until nothing is pending (renders finished, input consumed) for a few frames in

View File

@ -10,6 +10,7 @@ pub mod headless;
pub mod icons;
pub mod menubar;
pub mod menus;
pub mod merge;
pub mod panels;
pub mod render;
pub mod shortcuts;
@ -91,6 +92,8 @@ pub struct LightcraftApp {
/// What the loupe drew last frame: photo and source ("render", "cached", "embedded", "small",
/// "thumb", "none").
pub loupe_shown: Option<(lightcraft_catalog::PhotoId, &'static str)>,
/// Photo Merge dialog previews and background merges.
pub merge: merge::MergeState,
}
impl LightcraftApp {
@ -120,6 +123,7 @@ impl LightcraftApp {
widgets: vec![],
gesture: None,
loupe_shown: None,
merge: merge::MergeState::default(),
}
}
@ -183,7 +187,7 @@ impl LightcraftApp {
return;
}
let now = now_ms();
let busy = self.renderer.in_flight() > 0;
let busy = self.renderer.in_flight() > 0 || self.merge.busy();
let mut shots = std::mem::take(&mut self.pending_screenshots);
let mut shadow_ticked = false;
shots.retain_mut(|s| {
@ -240,7 +244,13 @@ impl LightcraftApp {
let memory = main.memory(|m| m.clone());
view.ctx.memory_mut(|m| *m = memory);
view.run(headless::HeadlessView::raw_input(size, ppp, time, vec![]), |ui| self.ui(ui));
let img = capture.then(|| view.paint(&self.renderer.cpu_textures()));
let img = capture.then(|| {
let mut tex = self.renderer.cpu_textures();
if let (Some((t, ..)), Some(px)) = (&self.merge.preview, &self.merge.preview_pixels) {
tex.insert(t.id(), crate::softpaint::CpuTexture::linear(px.clone()));
}
view.paint(&tex)
});
self.shadow = Some(view);
img
}
@ -290,6 +300,7 @@ impl LightcraftApp {
ctx.request_repaint();
}
self.renderer.poll(ctx, &mut self.session);
merge::poll(self, ctx);
self.session.persist_if_dirty();
self.collect_screenshots(ctx);
self.issue_screenshots(ctx);

View File

@ -126,6 +126,7 @@ const LAYOUT: &[(&str, &[&str])] = &[
"version.create",
"photo.virtualCopy",
"@Stack",
"@Photo Merge",
"---",
"develop.auto",
"develop.treatment",

View File

@ -65,6 +65,9 @@ pub const UI_COMMANDS: &[UiCommand] = &[
("dialog.autoStack", "Auto-Stack by Capture Time…", None, "Photo>Stack"),
("dialog.copySettings", "Choose Edit Settings to Copy…", Some("Cmd+Shift+C"), "Edit"),
("dialog.export", "Export…", Some("Cmd+Shift+E"), "File"),
("dialog.mergeHdr", "HDR…", Some("Ctrl+H"), "Photo>Photo Merge"),
("dialog.mergePanorama", "Panorama…", Some("Ctrl+M"), "Photo>Photo Merge"),
("dialog.mergeHdrPanorama", "HDR Panorama…", None, "Photo>Photo Merge"),
("file.addPhotos", "Add Photos…", Some("Cmd+Shift+I"), "File"),
("file.importPresets", "Import Presets…", None, "File"),
("file.exportPresets", "Export Presets…", None, "File"),
@ -337,6 +340,9 @@ pub fn run_ui_command(app: &mut LightcraftApp, id: &str, p: &Value) -> Option<Re
});
Ok(Value::Null)
}
"dialog.mergeHdr" => crate::merge::open(app, "merge.hdr"),
"dialog.mergePanorama" => crate::merge::open(app, "merge.panorama"),
"dialog.mergeHdrPanorama" => crate::merge::open(app, "merge.hdrPanorama"),
"app.about" => {
app.ui.dialog = Some(Dialog::About);
Ok(Value::Null)
@ -435,6 +441,7 @@ pub fn ui_enabled(app: &LightcraftApp, id: &str) -> bool {
"file.exportPresets" => app.session.presets.iter().any(|p| !p.builtin),
"view.compare" => app.session.catalog.len() > 1,
"compare.swap" | "compare.makeSelect" => app.ui.view == ViewMode::Compare,
s if s.starts_with("dialog.merge") => app.session.targets(&serde_json::json!({})).len() >= 2 && app.merge.final_task.is_none(),
_ => true,
}
}

389
crates/ui-egui/src/merge.rs Normal file
View File

@ -0,0 +1,389 @@
//! Photo Merge in the UI: the HDR / Panorama / HDR Panorama dialog (options + live preview) and
//! the background merge jobs.
//!
//! Previews and the final merge run [`lightcraft_engine::merge::MergeJob`]s on worker threads with
//! progress and cancellation; the window stays responsive. Changing an option cancels the running
//! preview and starts a new one. The final merge keeps running after the dialog closes (progress
//! in a toast); when it finishes, the result is written, imported and selected
//! ([`lightcraft_engine::Session::finish_merge`]).
use std::sync::atomic::{AtomicBool, Ordering};
use std::sync::mpsc::{Receiver, channel};
use std::sync::{Arc, Mutex};
use lightcraft_catalog::PhotoId;
use lightcraft_engine::merge::{MergeJob, MergeOutput};
use serde::{Deserialize, Serialize};
use serde_json::{Value, json};
use crate::LightcraftApp;
/// The dialog's options (also the parameters of the `merge.*` commands).
#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
#[serde(default, rename_all = "camelCase")]
pub struct MergeDialog {
/// "merge.hdr" | "merge.panorama" | "merge.hdrPanorama".
pub command: String,
pub align: bool,
/// none | low | medium | high
pub deghost: String,
pub show_overlay: bool,
pub auto_settings: bool,
pub stack: bool,
/// auto | spherical | cylindrical | perspective
pub projection: String,
pub boundary_warp: f64,
pub auto_crop: bool,
pub fill_edges: bool,
/// HDR panorama: photos per bracket (0 = from EXIF).
pub bracket: u32,
}
impl Default for MergeDialog {
fn default() -> Self {
MergeDialog {
command: "merge.hdr".into(),
align: true,
deghost: "none".into(),
show_overlay: false,
auto_settings: true,
stack: false,
projection: "auto".into(),
boundary_warp: 0.0,
auto_crop: true,
fill_edges: false,
bracket: 0,
}
}
}
impl MergeDialog {
pub fn for_command(command: &str) -> MergeDialog {
MergeDialog { command: command.into(), ..Default::default() }
}
pub fn title(&self) -> &'static str {
match self.command.as_str() {
"merge.panorama" => "Panorama Merge Preview",
"merge.hdrPanorama" => "HDR Panorama Merge Preview",
_ => "HDR Merge Preview",
}
}
pub fn is_hdr(&self) -> bool {
self.command != "merge.panorama"
}
pub fn is_pano(&self) -> bool {
self.command != "merge.hdr"
}
pub fn params(&self) -> Value {
serde_json::to_value(self).unwrap_or_default()
}
}
/// A merge running on a worker thread.
pub struct MergeTask {
pub preview: bool,
pub job: MergeJob,
/// The dialog options it was started with.
pub options: MergeDialog,
pub progress: Arc<Mutex<(f32, String)>>,
pub cancel: Arc<AtomicBool>,
rx: Receiver<Result<MergeOutput, String>>,
}
/// Merge state of the app.
#[derive(Default)]
pub struct MergeState {
pub preview_task: Option<MergeTask>,
pub final_task: Option<MergeTask>,
/// The latest preview (texture, size, its options, stats) and the latest error.
pub preview: Option<(egui::TextureHandle, [usize; 2], MergeDialog, Value)>,
pub preview_pixels: Option<Arc<egui::ColorImage>>,
pub error: Option<String>,
/// Options whose preview failed (not retried until they change).
pub failed_options: Option<MergeDialog>,
/// The last finished merge (`{id, path, …}`).
pub last_result: Option<Value>,
/// Photos the dialog merges (the selection when it opened).
pub ids: Vec<PhotoId>,
}
impl MergeState {
/// A preview or a merge is running.
pub fn busy(&self) -> bool {
self.preview_task.is_some() || self.final_task.is_some()
}
}
fn spawn(app: &LightcraftApp, options: &MergeDialog, ids: &[PhotoId], preview: bool) -> Result<MergeTask, String> {
let (kind, finish) = lightcraft_engine::merge::parse(&options.command, &options.params()).map_err(|e| e.to_string())?;
let job = app.session.plan_merge(kind, finish, ids, preview).map_err(|e| e.to_string())?;
let progress = Arc::new(Mutex::new((0.0, "Starting".to_string())));
let cancel = Arc::new(AtomicBool::new(false));
let (tx, rx) = channel();
let (j, p, c) = (job.clone(), progress.clone(), cancel.clone());
let work = move || {
let r = j.run(&|f, stage| {
if let Ok(mut g) = p.lock() {
*g = (f, stage.to_string());
}
!c.load(Ordering::Relaxed)
});
let _ = tx.send(r);
};
#[cfg(not(target_arch = "wasm32"))]
std::thread::Builder::new().name("photo-merge".into()).spawn(work).map_err(|e| e.to_string())?;
#[cfg(target_arch = "wasm32")]
work();
Ok(MergeTask { preview, job, options: options.clone(), progress, cancel, rx })
}
/// Open the merge dialog for the selection.
pub fn open(app: &mut LightcraftApp, command: &str) -> Result<Value, String> {
let ids = app.session.targets(&json!({}));
if ids.len() < 2 {
return Err("select at least 2 photos to merge".into());
}
app.merge.ids = ids;
app.merge.error = None;
app.merge.preview = None;
app.merge.preview_pixels = None;
if let Some(t) = app.merge.preview_task.take() {
t.cancel.store(true, Ordering::Relaxed);
}
app.ui.dialog = Some(crate::state::Dialog::Merge { opts: MergeDialog::for_command(command) });
Ok(json!({"photos": app.merge.ids.len()}))
}
/// Start the full-resolution merge (the dialog's Merge button / `ui.dialog.confirm`).
pub fn start_final(app: &mut LightcraftApp, opts: &MergeDialog) -> Result<Value, String> {
if app.merge.final_task.is_some() {
return Err("a merge is already running".into());
}
if let Some(t) = app.merge.preview_task.take() {
t.cancel.store(true, Ordering::Relaxed);
}
let ids = app.merge.ids.clone();
let task = spawn(app, opts, &ids, false)?;
app.merge.final_task = Some(task);
Ok(json!({"started": true, "photos": ids.len()}))
}
/// Per frame: keep the preview in sync with the dialog's options, collect finished jobs.
pub fn poll(app: &mut LightcraftApp, ctx: &egui::Context) {
// preview for the open dialog
let dialog_opts = match &app.ui.dialog {
Some(crate::state::Dialog::Merge { opts }) => Some(opts.clone()),
_ => None,
};
match &dialog_opts {
Some(opts) => {
let current = app.merge.preview_task.as_ref().map(|t| &t.options).or(app.merge.preview.as_ref().map(|p| &p.2));
let stale = current != Some(opts) && !last_failed_matches(app, opts);
if stale {
if let Some(t) = app.merge.preview_task.take() {
t.cancel.store(true, Ordering::Relaxed);
}
let ids = app.merge.ids.clone();
match spawn(app, opts, &ids, true) {
Ok(t) => {
app.merge.error = None;
app.merge.preview_task = Some(t);
}
Err(e) => {
app.merge.error = Some(e);
app.merge.failed_options = Some(opts.clone());
}
}
}
}
None => {
if let Some(t) = app.merge.preview_task.take() {
t.cancel.store(true, Ordering::Relaxed);
}
}
}
if let Some(t) = &app.merge.preview_task {
match t.rx.try_recv() {
Ok(Ok(out)) => {
let t = app.merge.preview_task.take().expect("task");
if let Some(img) = out.preview {
let color = Arc::new(egui::ColorImage::from_rgba_unmultiplied([img.width, img.height], &img.as_bytes()));
let tex = ctx.load_texture("merge-preview", color.clone(), egui::TextureOptions::LINEAR);
app.merge.preview_pixels = Some(color);
app.merge.preview = Some((tex, [img.width, img.height], t.options, out.info));
}
}
Ok(Err(e)) => {
let t = app.merge.preview_task.take().expect("task");
if e != "cancelled" {
app.merge.error = Some(e);
app.merge.failed_options = Some(t.options);
}
}
Err(_) => ctx.request_repaint_after(std::time::Duration::from_millis(50)),
}
}
// the final merge
if let Some(t) = &app.merge.final_task {
match t.rx.try_recv() {
Ok(r) => {
let t = app.merge.final_task.take().expect("task");
let what = if t.options.command == "merge.panorama" {
"Panorama"
} else if t.options.command == "merge.hdrPanorama" {
"HDR panorama"
} else {
"HDR"
};
match r.map_err(lightcraft_engine::EngineError::Other).and_then(|out| app.session.finish_merge(&t.job, out)) {
Ok(v) => {
app.merge.last_result = Some(v);
app.toast(ctx, format!("{what} merge added"));
}
Err(e) => {
app.ui.status = e.to_string();
app.toast(ctx, format!("{what} merge failed: {e}"));
}
}
}
Err(_) => {
let (f, stage) = t.progress.lock().map(|g| g.clone()).unwrap_or_default();
let now = ctx.input(|i| i.time);
app.ui.toast = Some((format!("Merging… {stage} {:.0}%", f * 100.0), now + 0.5));
ctx.request_repaint_after(std::time::Duration::from_millis(100));
}
}
}
}
fn last_failed_matches(app: &LightcraftApp, opts: &MergeDialog) -> bool {
app.merge.failed_options.as_ref() == Some(opts)
}
/// The dialog body (options on the left, preview on the right).
pub fn body(app: &mut LightcraftApp, ui: &mut egui::Ui, opts: &mut MergeDialog) {
let t = crate::theme::Tokens::get(ui.ctx());
ui.horizontal_top(|ui| {
// preview
ui.vertical(|ui| {
let box_size = egui::vec2(520.0, 360.0);
let (rect, _) = ui.allocate_exact_size(box_size, egui::Sense::hover());
ui.painter().rect_filled(rect, 4.0, t.canvas);
crate::widgets::register(ui.ctx(), "merge.preview", rect);
if let Some((tex, size, _, _)) = &app.merge.preview {
let s = (box_size.x / size[0] as f32).min(box_size.y / size[1] as f32);
let r = egui::Rect::from_center_size(rect.center(), egui::vec2(size[0] as f32 * s, size[1] as f32 * s));
ui.painter().image(tex.id(), r, egui::Rect::from_min_max(egui::pos2(0.0, 0.0), egui::pos2(1.0, 1.0)), egui::Color32::WHITE);
}
if let Some(task) = &app.merge.preview_task {
let (f, stage) = task.progress.lock().map(|g| g.clone()).unwrap_or_default();
ui.painter().text(
rect.center_bottom() - egui::vec2(0.0, 14.0),
egui::Align2::CENTER_CENTER,
format!("{stage}… {:.0}%", f * 100.0),
t.font(12.0),
t.text,
);
}
if let Some(e) = &app.merge.error {
ui.painter().text(rect.center(), egui::Align2::CENTER_CENTER, e, t.font(13.0), t.text);
}
let info = app.merge.preview.as_ref().map(|p| p.3.clone());
if let Some(info) = info {
let txt = if opts.is_pano() {
format!(
"{} photos · {} · {:.0}° × {:.0}°",
info["used"].as_array().map_or(0, Vec::len),
info["projection"].as_str().unwrap_or(""),
info["fov"][0].as_f64().unwrap_or(0.0),
info["fov"][1].as_f64().unwrap_or(0.0)
)
} else {
let ev: Vec<String> =
info["ev"].as_array().map(|a| a.iter().filter_map(Value::as_f64).map(|v| format!("{v:+.1}")).collect()).unwrap_or_default();
format!("{} photos · exposures {} EV", ev.len(), ev.join(" / "))
};
ui.label(egui::RichText::new(txt).color(t.text_dim));
}
});
ui.add_space(12.0);
// options
ui.vertical(|ui| {
ui.set_width(220.0);
if opts.is_pano() {
ui.label(egui::RichText::new("Projection").color(t.text_label));
for (v, l) in [("auto", "Auto Select"), ("spherical", "Spherical"), ("cylindrical", "Cylindrical"), ("perspective", "Perspective")] {
if ui.radio(opts.projection == v, l).clicked() {
opts.projection = v.to_string();
}
}
ui.add_space(4.0);
let mut bw = opts.boundary_warp;
if num(ui, &BOUNDARY_WARP, &mut bw) {
opts.boundary_warp = bw.round();
}
ui.checkbox(&mut opts.fill_edges, "Fill Edges");
ui.checkbox(&mut opts.auto_crop, "Auto Crop");
ui.add_space(4.0);
}
if opts.is_hdr() {
ui.checkbox(&mut opts.align, "Auto Align");
ui.label(egui::RichText::new("Deghost Amount").color(t.text_label));
ui.horizontal(|ui| {
for (i, (v, l)) in [("none", "None"), ("low", "Low"), ("medium", "Med"), ("high", "High")].iter().enumerate() {
if crate::widgets::text_button(ui, &format!("mergeDeghost-{i}"), l, opts.deghost == *v).clicked() {
opts.deghost = v.to_string();
}
}
});
if opts.command == "merge.hdr" {
ui.add_enabled_ui(opts.deghost != "none", |ui| ui.checkbox(&mut opts.show_overlay, "Show Deghost Overlay"));
}
if opts.command == "merge.hdrPanorama" {
let mut b = opts.bracket as f64;
if num(ui, &BRACKET, &mut b) {
opts.bracket = b.round() as u32;
}
}
ui.add_space(4.0);
}
ui.checkbox(&mut opts.auto_settings, "Auto Settings");
ui.checkbox(&mut opts.stack, "Create Stack");
ui.add_space(8.0);
ui.label(egui::RichText::new(format!("{} photos", app.merge.ids.len())).color(t.text_dim));
});
});
}
const BOUNDARY_WARP: lightcraft_develop::ControlSpec = lightcraft_develop::ControlSpec {
id: "merge.boundaryWarp",
label: "Boundary Warp",
section: lightcraft_develop::Section::Light,
min: 0.0,
max: 100.0,
default: 0.0,
step: 1.0,
decimals: 0,
track: lightcraft_develop::Track::Plain,
};
const BRACKET: lightcraft_develop::ControlSpec = lightcraft_develop::ControlSpec {
id: "merge.bracket",
label: "Photos per bracket (0 = auto)",
section: lightcraft_develop::Section::Light,
min: 0.0,
max: 9.0,
default: 0.0,
step: 1.0,
decimals: 0,
track: lightcraft_develop::Track::Plain,
};
fn num(ui: &mut egui::Ui, spec: &lightcraft_develop::ControlSpec, v: &mut f64) -> bool {
match crate::widgets::slider(ui, spec, *v, true, None).value {
Some(n) => {
*v = n;
true
}
None => false,
}
}

View File

@ -27,6 +27,7 @@ pub fn show(app: &mut LightcraftApp, ctx: &egui::Context) {
Dialog::CreatePreset { .. } => "Create Preset",
Dialog::CopySettings { .. } => "Choose Edit Settings to Copy",
Dialog::Export { .. } => "Export",
Dialog::Merge { opts } => opts.title(),
Dialog::About => "About LightCraft",
Dialog::Shortcuts => "Keyboard Shortcuts",
};
@ -187,6 +188,7 @@ pub fn show(app: &mut LightcraftApp, ctx: &egui::Context) {
});
field(ui, "Folder", |ui| ui.add(egui::TextEdit::singleline(dir).desired_width(f32::INFINITY)));
}
Dialog::Merge { opts } => crate::merge::body(app, ui, opts),
Dialog::About => {
ui.label(egui::RichText::new("LightCraft").font(t.semibold(20.0)).color(t.text));
ui.label(format!("Version {} — a clean-room, pure-Rust photo library and raw developer.", env!("CARGO_PKG_VERSION")));
@ -233,7 +235,12 @@ pub fn show(app: &mut LightcraftApp, ctx: &egui::Context) {
if !informational && ui.button("Cancel").clicked() {
close = true;
}
if ui.button(if informational { "Close" } else { "OK" }).clicked() {
let ok = match dlg {
Dialog::Merge { .. } => "Merge",
_ if informational => "Close",
_ => "OK",
};
if ui.button(ok).clicked() {
if informational {
close = true;
} else {
@ -287,6 +294,7 @@ pub fn confirm_dialog(app: &mut LightcraftApp, dlg: &Dialog) -> Result<serde_jso
"metadata": opts.metadata, "removeLocation": opts.remove_location, "watermark": opts.watermark,
}),
),
Dialog::Merge { opts } => crate::merge::start_final(app, opts),
Dialog::About | Dialog::Shortcuts => Ok(serde_json::Value::Null),
}
}

View File

@ -354,6 +354,14 @@ pub fn context_menu(app: &mut LightcraftApp, ui: &mut egui::Ui, id: PhotoId) {
if ui.button("Reset Edits").clicked() {
let _ = app.run("develop.reset", json!({}));
}
ui.menu_button("Photo Merge", |ui| {
let n = app.session.targets(&json!({})).len();
for (id, label) in [("dialog.mergeHdr", "HDR…"), ("dialog.mergePanorama", "Panorama…"), ("dialog.mergeHdrPanorama", "HDR Panorama…")] {
if ui.add_enabled(n >= 2, egui::Button::new(label)).clicked() {
let _ = app.run(id, json!({}));
}
}
});
ui.separator();
if ui.button("Rotate Left").clicked() {
let _ = app.run("photo.rotateLeft", json!({}));

View File

@ -52,10 +52,11 @@ pub fn parse(s: &str) -> Option<(Modifiers, Key)> {
fn matches(i: &egui::InputState, m: Modifiers, k: Key) -> bool {
i.events.iter().any(|e| match e {
egui::Event::Key { key, pressed: true, modifiers, .. } => {
(*key == k || (k == Key::Backspace && *key == Key::Delete))
&& modifiers.command == m.command
&& modifiers.shift == m.shift
&& modifiers.alt == m.alt
// `Ctrl` is the physical Control key (on macOS distinct from Cmd; elsewhere Cmd = Ctrl);
// "Delete" (⌫ = Backspace) also matches forward-delete
let ctrl_ok = if m.ctrl { modifiers.ctrl } else { !modifiers.ctrl || modifiers.command };
let cmd_ok = m.ctrl || modifiers.command == m.command;
(*key == k || (k == Key::Backspace && *key == Key::Delete)) && ctrl_ok && cmd_ok && modifiers.shift == m.shift && modifiers.alt == m.alt
}
_ => false,
})

View File

@ -163,6 +163,10 @@ pub enum Dialog {
limit_kb: u32,
dir: String,
},
/// Photo Merge (HDR / Panorama / HDR Panorama) options; the preview lives in the app.
Merge {
opts: crate::merge::MergeDialog,
},
About,
Shortcuts,
}

64
docs/merge.md Normal file
View File

@ -0,0 +1,64 @@
# Photo Merge (HDR, Panorama, HDR Panorama)
`crates/merge` (layer L3, no UI dependencies) implements the algorithms; `crates/engine/src/merge.rs`
plans, runs and finishes merges in the library; the desktop UI shows a dialog with a live preview.
All algorithms are our own implementations from published papers (cited in the module docs).
## Using it
- **App:** select the photos, then *Photo Merge ▸ HDR…* (`Ctrl+H`), *Panorama…* (`Ctrl+M`) or
*HDR Panorama…* from the photo context menu (or the `dialog.merge*` commands). The dialog previews
the merge at ≤ 1024 px and re-runs the preview when an option changes; *Merge* runs the full merge
in the background (progress in a toast). The result is written next to the first photo as
`<name>-HDR.dng`, `-Pano.dng` or `-HDR-Pano.dng`, imported and selected.
- **Commands** (CLI, MCP, control channel):
- `merge.hdr {ids?, align=true, deghost=none|low|medium|high, autoSettings=true, stack=false, preview=false, showOverlay=false, previewPath?}`
- `merge.panorama {ids?, projection=auto|spherical|cylindrical|perspective, boundaryWarp=0..100, autoCrop=false, fillEdges=false, autoSettings=true, maxMegapixels=40, preview, previewPath?}`
- `merge.hdrPanorama {ids?, bracket=0 (from EXIF), …both sets of options}`
- **CLI:** `lightcraft-cli merge hdr|panorama|hdr-panorama [options] FILES…`;
`lightcraft-cli synth-merge hdr|panorama -o DIR` writes synthetic inputs (procedural scene).
## What it does
**HDR** (`hdr.rs`): frames ordered by EXIF exposure (or median luminance); the middle one is the
reference. Auto Align: features on exposure-normalised luminance, RANSAC translation / similarity /
homography (the simplest model that explains ≥ 97 % of the matches). Exposure ratios are measured
(median log ratio of well-exposed pixels between brightness neighbours), EXIF is the fallback.
Radiance is a weighted mean in linear light, weight `k · (1 − smoothstep(0.8·clip, 0.96·clip, max
channel))`. Deghost: a per-pixel consistency test against the reference radiance (with a noise
floor), dilated and feathered masks; the union is the deghost overlay. Output: 16-bit float DNG.
**Panorama** (`pano/`): MOPS-style features, pairwise RANSAC homographies with Brown & Lowe's
verification, focal from homographies (or EXIF 35 mm), rotations chained along a maximum spanning
tree, Levenberg–Marquardt bundle adjustment of rotations + focal lengths, straightening.
Projections: spherical, cylindrical, perspective; *Auto* picks perspective for compact fields of
view < 110°, cylindrical when the vertical field is < 90°, else spherical. Exposure: EXIF, then
gain compensation from the overlaps. Seams: each pixel goes to the image with the highest
centre-weighted feather weight; multi-band blending (Laplacian pyramids) in log space. Boundary
Warp: a separable stretch of the smoothed valid-area edges to the canvas edges. Auto Crop: the
largest rectangle of image data (set as the develop crop, so it stays adjustable). Fill Edges:
push–pull diffusion. Output: 16-bit integer DNG (16-bit float for HDR panoramas).
**DNG output** (`output.rs`): LinearRaw, 3 samples, Deflate tiles; values scaled below 1.0 with the
scale in `BaselineExposure`; raw sources keep their DNG colour tags (camera RGB), other sources get
`ColorMatrix1`/`ForwardMatrix1` of their RGB space.
## Measured (synthetic tests, `cargo test -p lightcraft-merge`)
- HDR, 3 brackets ±2 EV with 14-bit noise and handheld shifts: alignment error < 0.3 px, exposure
error < 0.02 EV, radiance median log error < 0.03 (noise-limited, bias < 0.01), recovered
highlights < 0.05; moving object: ghost error 0.2–2.0 without deghosting, < 0.06 with High.
- Panorama, 5 views over 182° with a 1.3× exposure drift: relative rotations < 0.1°, focal < 1 %,
BA RMS ≈ 0.35 px; median log error vs ground truth < 0.02; no seam steps.
## Limits (known gaps)
- No lens-distortion model in the panorama bundle adjustment (rectilinear lenses assumed); no
360° wrap-around; flat (translating-camera) scans are fitted with the rotation model.
- No camera response calibration for JPEG brackets (sRGB decoding is taken as linear light).
- Boundary Warp is a separable stretch, not a content-preserving mesh warp; Fill Edges is diffusion,
not patch synthesis.
- *Create Stack*: the result goes on top of a collapsed stack of its sources (`Catalog::stack_with_ops`); its
caption also lists them.
- Panorama output is limited to `maxMegapixels` (default 40) — blending memory is ~25 bytes/pixel.
- The browser build has no file access for merges yet.

View File

@ -35,7 +35,7 @@ and prints the summary; `cargo xtask parity --write` refreshes the summary table
| L. Presets (PRE) | 1 | 4 | 2 | 1 | 0/2 (0%) | 1/2 (50%) |
| M. Versions & history (VER) | 3 | 1 | 1 | 0 | 1/1 (100%) | 2/3 (67%) |
| N. Copy / paste / sync (SYNC) | 4 | 0 | 1 | 0 | 3/3 (100%) | 1/1 (100%) |
| O. Merge (MERGE) | 0 | 0 | 4 | 0 | — | — |
| O. Merge (MERGE) | 3 | 0 | 1 | 0 | — | — |
| P. Enhance (ENH) | 0 | 0 | 2 | 0 | — | — |
| Q. HDR (HDR) | 0 | 0 | 5 | 0 | — | — |
| R. Video (VID) | 0 | 0 | 5 | 1 | — | 0/2 (0%) |
@ -48,7 +48,7 @@ and prints the summary; `cargo xtask parity --write` refreshes the summary table
| Y. Menus | 46 | 18 | 21 | 8 | 37/47 (79%) | 7/23 (30%) |
| Z. Keyboard shortcuts (desktop) | 46 | 16 | 18 | 1 | 39/52 (75%) | 7/23 (30%) |
| Lightroom Classic extras | 4 | 19 | 56 | 9 | — | 3/21 (14%) |
| **Total** | 205 | 103 | 192 | 37 | 152/198 (77%) | 48/144 (33%) |
| **Total** | 208 | 103 | 189 | 37 | 152/198 (77%) | 48/144 (33%) |
<!-- /parity:summary -->
## Top gaps
@ -401,9 +401,9 @@ under the UI command that wraps it (a no-op error); the UI command now wins.
| Id | Feature | Tier | Status | Evidence | Missing / notes |
|---|---|---|---|---|---|
| LR-MERGE-HDR | HDR merge | P2 | ⬜ | | |
| LR-MERGE-PANO | Panorama | P2 | ⬜ | | |
| LR-MERGE-HDRPANO | HDR panorama | P2 | ⬜ | | |
| LR-MERGE-HDR | HDR merge | P2 | ✅ | `cmd:merge.hdr`, `cmd:dialog.mergeHdr` | auto align, deghost None–High + overlay, auto settings, float DNG, Create Stack; JPEG brackets treated as linear |
| LR-MERGE-PANO | Panorama | P2 | ✅ | `cmd:merge.panorama`, `cmd:dialog.mergePanorama` | spherical/cylindrical/perspective + auto, boundary warp, auto crop, fill edges (diffusion), DNG; no lens model / 360° wrap |
| LR-MERGE-HDRPANO | HDR panorama | P2 | ✅ | `cmd:merge.hdrPanorama`, `cmd:dialog.mergeHdrPanorama` | brackets grouped by EXIF |
| LR-MERGE-HEADLESS | Merge with last settings | P2 | ⬜ | | |
## P. Enhance (ENH)