Files
DarkRoom/core/dr-pano/src/align.rs
T
dtourolle 621a6b8313 Leave a panorama frame out without leaving the page
A frame that did not fit ended the job with its name, and the only way
on was Back, a smaller selection and every frame read, demosaiced and
searched for keypoints again. Each row on the page now has a box. An
unticked frame is left out and the rest are solved again from what the
first pass measured.

dr_pano::align is split for it: match_pairs does the matching and the
pairwise RANSAC once over every frame (about 5 s for twelve), and solve
takes a subset of the frames and uses only the links among them (about
0.1 s). Solving a subset by re-aligning it also moved every RANSAC seed,
which are keyed on frame position, and on the fixture set that was
enough to lose a marginal link and strand a neighbour of the frame left
out.

A frame that cannot be placed no longer stops the job either. Its row
names it and Merge stays off until it is unticked; the headless example
leaves such frames out the same way, and takes --leave-out N to try it.
2026-09-27 16:29:39 -04:00

596 lines
21 KiB
Rust

//! TRACES: FR-MRG-1 | FR-MRG-5
//! From features to cameras: the alignment of a whole set.
//!
//! 1. Match every pair of frames (`matching`).
//! 2. For each pair with enough matches, a robust homography
//! (`homography::ransac_homography`); a pair is a *link* when its inliers
//! pass Brown & Lowe's test, `n_inliers > 8 + 0.3 · n_matches`, which
//! is what separates a real overlap from a coincidence of descriptors.
//! 3. The focal length: the median of what the links' homographies imply,
//! or the caller's hint if none of them implies anything.
//! 4. A spanning tree over the links, strongest first, from the
//! best-connected frame; rotations chained along it.
//! 5. Bundle adjustment over every link's inliers (`bundle`).
//!
//! Steps 1 and 2 are [`match_pairs`] and most of the time; 3 to 5 are
//! [`solve`], which takes a subset of the frames. Leaving a frame out is
//! then a solve over the pairs already measured — the same links, not a
//! fresh RANSAC whose seeds would move with the frames' positions.
//!
//! What it refuses to do is guess. A frame the tree does not reach is
//! reported by index with the reason (FR-MRG-5) and left out of the
//! cameras; the caller decides whether a set with a hole is worth
//! stitching, and the requirement says it is not.
use crate::bundle::{self, AdjustOptions, Cameras, Observation};
use crate::features::Features;
use crate::homography::{self, RobustHomography};
use crate::linalg::Mat3;
use crate::matching::{match_features, Match};
use crate::PanoError;
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct AlignOptions {
/// Descriptor similarity floor for a match (`matching`).
pub min_similarity: f32,
/// RANSAC agreement distance, in pixels of the features' image.
pub ransac_px: f64,
pub ransac_iterations: usize,
/// A pair needs at least this many inliers to be a link, on top of
/// Brown & Lowe's ratio test.
pub min_inliers: usize,
/// Focal length in pixels of the features' image, if the caller knows
/// it (EXIF and a sensor width). Used only when the homographies do not
/// determine one.
pub focal_hint: Option<f64>,
pub adjust: AdjustOptions,
/// For RANSAC's sampling: the same seed gives the same alignment
/// (NFR-MRG-2).
pub seed: u64,
}
impl Default for AlignOptions {
fn default() -> Self {
AlignOptions {
min_similarity: 0.82,
ransac_px: 3.0,
ransac_iterations: 1000,
min_inliers: 12,
focal_hint: None,
adjust: AdjustOptions::default(),
seed: 0x5eed,
}
}
}
/// An overlap the alignment trusts.
#[derive(Debug, Clone, PartialEq)]
pub struct Link {
pub i: usize,
pub j: usize,
pub matches: usize,
pub inliers: usize,
/// Maps centred points of `i` to centred points of `j`.
pub h: Mat3,
}
/// Why a frame is not in the alignment.
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum Unaligned {
/// Not enough matches with any other frame to try a geometry.
NoMatches,
/// Matches existed but none survived RANSAC as a real overlap.
NoOverlap,
/// Overlaps existed but only with frames that are themselves unaligned.
Disconnected,
}
impl std::fmt::Display for Unaligned {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
f.write_str(match self {
Unaligned::NoMatches => "too few matching features with any other frame",
Unaligned::NoOverlap => "no consistent overlap with any other frame",
Unaligned::Disconnected => "overlaps only with frames that could not be aligned",
})
}
}
/// The result: cameras for the aligned frames, and the rest named.
#[derive(Debug, Clone, PartialEq)]
pub struct Alignment {
/// One rotation per input frame, camera to world, for aligned frames;
/// `None` for the unaligned. The reference frame is the best-connected
/// one and has the identity.
pub rotations: Vec<Option<Mat3>>,
/// Focal length in pixels of the features' image.
pub focal: f64,
pub links: Vec<Link>,
pub unaligned: Vec<(usize, Unaligned)>,
/// Bundle adjustment's RMS reprojection error, in pixels.
pub rms_px: f64,
}
impl Alignment {
pub fn is_complete(&self) -> bool {
self.unaligned.is_empty()
}
/// The cameras of the aligned frames, indexed as the input — a frame
/// that is not aligned is given the identity, so this is only useful
/// when [`Self::is_complete`].
pub fn cameras(&self) -> Cameras {
Cameras {
rotations: self
.rotations
.iter()
.map(|r| r.unwrap_or(Mat3::IDENTITY))
.collect(),
focal: self.focal,
}
}
}
/// Every pair of a set measured: steps 1 and 2, the expensive part, kept
/// so that a solve over a subset reuses it.
#[derive(Debug, Clone, PartialEq)]
pub struct Pairs {
/// Each frame's long edge, for the focal length's clamp.
long_edges: Vec<f64>,
/// Pairs with enough matches to try a geometry, whether or not it held.
matched: Vec<(usize, usize)>,
links: Vec<Link>,
/// Every link's inliers, in pixels, centred.
observations: Vec<Observation>,
}
impl Pairs {
/// How many frames were measured.
pub fn len(&self) -> usize {
self.long_edges.len()
}
pub fn is_empty(&self) -> bool {
self.long_edges.is_empty()
}
}
/// Align a set of frames from their features: [`match_pairs`], then
/// [`solve`] over all of them.
///
/// Every `Features` must be in its own frame's pixel coordinates with the
/// image size filled in; points are centred on the image centre here. The
/// frames must all come from the same lens at the same focal length, which
/// is the panorama assumption and not checked — the caller has the EXIF.
pub fn align(frames: &[Features], opts: &AlignOptions) -> Result<Alignment, PanoError> {
let pairs = match_pairs(frames, opts)?;
solve(&pairs, &vec![true; frames.len()], opts)
}
/// Steps 1 and 2: every pair matched, and a robust homography for each
/// pair with enough matches.
pub fn match_pairs(frames: &[Features], opts: &AlignOptions) -> Result<Pairs, PanoError> {
let n = frames.len();
if n < 2 {
return Err(PanoError::Input(
"a panorama needs at least two frames".into(),
));
}
let centre = |k: usize, i: usize| -> (f64, f64) {
let kp = frames[k].keypoints[i];
(
f64::from(kp.x) - frames[k].width as f64 / 2.0,
f64::from(kp.y) - frames[k].height as f64 / 2.0,
)
};
// Scale for the DLT's conditioning: points of order one.
let scale = 1.0
/ frames
.iter()
.map(|f| f.width.max(f.height) as f64)
.fold(1.0, f64::max);
// 1 + 2: every pair.
let mut links = Vec::new();
let mut observations: Vec<Observation> = Vec::new();
let mut matched = Vec::new();
let t_match = std::time::Instant::now();
for i in 0..n {
for j in i + 1..n {
let matches: Vec<Match> = match_features(&frames[i], &frames[j], opts.min_similarity);
log::debug!("pair {i}-{j}: {} matches", matches.len());
if matches.len() < 4 {
continue;
}
matched.push((i, j));
let pairs: Vec<((f64, f64), (f64, f64))> = matches
.iter()
.map(|m| {
let (a, b) = (centre(i, m.a), centre(j, m.b));
((a.0 * scale, a.1 * scale), (b.0 * scale, b.1 * scale))
})
.collect();
let Some(RobustHomography { h, inliers }) = homography::ransac_homography(
&pairs,
opts.ransac_px * scale,
opts.ransac_iterations,
opts.seed ^ ((i as u64) << 32 | j as u64),
) else {
continue;
};
let needed = (8.0 + 0.3 * matches.len() as f64).ceil() as usize;
log::debug!("pair {i}-{j}: {} inliers, {needed} needed", inliers.len());
if inliers.len() <= needed || inliers.len() < opts.min_inliers {
continue;
}
// Back to pixels: H_px = S⁻¹ H S.
let m = h.0;
let h_px = Mat3([
[m[0][0], m[0][1], m[0][2] / scale],
[m[1][0], m[1][1], m[1][2] / scale],
[m[2][0] * scale, m[2][1] * scale, m[2][2]],
]);
for &k in &inliers {
let (a, b) = pairs[k];
observations.push(Observation {
i,
j,
pi: (a.0 / scale, a.1 / scale),
pj: (b.0 / scale, b.1 / scale),
});
}
links.push(Link {
i,
j,
matches: matches.len(),
inliers: inliers.len(),
h: h_px,
});
}
}
log::debug!("matching and pairwise geometry in {:?}", t_match.elapsed());
Ok(Pairs {
long_edges: frames
.iter()
.map(|f| f.width.max(f.height) as f64)
.collect(),
matched,
links,
observations,
})
}
/// Steps 3 to 5 over the frames `keep` marks, from pairs already measured.
///
/// The result is indexed by the kept frames in order: its frame `k` is the
/// `k`-th frame `keep` marks. Only pairs whose frames are both kept take
/// part, so a frame whose only overlap was with one left out is reported
/// as unaligned, as it would be had it never been measured with it.
pub fn solve(pairs: &Pairs, keep: &[bool], opts: &AlignOptions) -> Result<Alignment, PanoError> {
if keep.len() != pairs.len() {
return Err(PanoError::Input(format!(
"{} flags for {} frames",
keep.len(),
pairs.len()
)));
}
// Input index to the solve's.
let mut slot = vec![None; keep.len()];
let mut n = 0usize;
for (k, &kept) in keep.iter().enumerate() {
if kept {
slot[k] = Some(n);
n += 1;
}
}
if n < 2 {
return Err(PanoError::Input(
"a panorama needs at least two frames".into(),
));
}
let both = |i: usize, j: usize| Some((slot[i]?, slot[j]?));
let mut matched_any = vec![false; n];
for &(i, j) in &pairs.matched {
if let Some((i, j)) = both(i, j) {
matched_any[i] = true;
matched_any[j] = true;
}
}
let links: Vec<Link> = pairs
.links
.iter()
.filter_map(|l| {
let (i, j) = both(l.i, l.j)?;
Some(Link { i, j, ..l.clone() })
})
.collect();
let observations: Vec<Observation> = pairs
.observations
.iter()
.filter_map(|o| {
let (i, j) = both(o.i, o.j)?;
Some(Observation { i, j, ..*o })
})
.collect();
// 3: the focal length.
let mut estimates: Vec<f64> = links
.iter()
.filter_map(|l| homography::focal_from_homography(&l.h))
.filter(|f| f.is_finite() && *f > 0.0)
.collect();
let longest = pairs
.long_edges
.iter()
.zip(keep)
.filter(|(_, &kept)| kept)
.map(|(&e, _)| e)
.fold(0.0, f64::max);
let focal = if !estimates.is_empty() {
estimates.sort_by(f64::total_cmp);
let median = estimates[estimates.len() / 2];
// A homography of a nearly pure pan can imply almost anything;
// clamp to the range a real lens on this sensor can reach.
median.clamp(0.3 * longest, 6.0 * longest)
} else if let Some(hint) = opts.focal_hint {
hint
} else {
// No overlap said anything and nobody told us: a normal lens.
longest
};
// 4: spanning tree, strongest link first, from the best-connected frame.
let mut rotations: Vec<Option<Mat3>> = vec![None; n];
let mut unaligned = Vec::new();
if links.is_empty() {
for (k, &matched) in matched_any.iter().enumerate() {
unaligned.push((
k,
if matched {
Unaligned::NoOverlap
} else {
Unaligned::NoMatches
},
));
}
return Ok(Alignment {
rotations,
focal,
links,
unaligned,
rms_px: 0.0,
});
}
let mut degree = vec![0usize; n];
for l in &links {
degree[l.i] += l.inliers;
degree[l.j] += l.inliers;
}
let root = (0..n).max_by_key(|&k| degree[k]).unwrap_or(0);
rotations[root] = Some(Mat3::IDENTITY);
loop {
// The strongest link from an aligned frame to an unaligned one.
let best = links
.iter()
.filter(|l| rotations[l.i].is_some() != rotations[l.j].is_some())
.max_by_key(|l| l.inliers);
let Some(l) = best else { break };
let r_ij = homography::rotation_from_homography(&l.h, focal);
// H_ij takes points of i to j, so bearings b_j = R_ij b_i, and with
// world = R_i · cam_i: R_j = R_i · R_ijᵀ.
if let Some(ri) = rotations[l.i] {
rotations[l.j] = Some((ri * r_ij.transpose()).orthonormalised());
} else if let Some(rj) = rotations[l.j] {
rotations[l.i] = Some((rj * r_ij).orthonormalised());
}
}
for k in 0..n {
if rotations[k].is_none() {
let reason = if !matched_any[k] {
Unaligned::NoMatches
} else if links.iter().any(|l| l.i == k || l.j == k) {
Unaligned::Disconnected
} else {
Unaligned::NoOverlap
};
unaligned.push((k, reason));
}
}
// 5: adjust the aligned frames together. The reference frame must be
// index 0 of the adjustment (it holds frame 0 fixed), so the aligned
// frames are renumbered with the root first.
let aligned: Vec<usize> = std::iter::once(root)
.chain((0..n).filter(|&k| k != root && rotations[k].is_some()))
.collect();
let index_of = |k: usize| aligned.iter().position(|&a| a == k);
let start = Cameras {
rotations: aligned.iter().map(|&k| rotations[k].unwrap()).collect(),
focal,
};
let obs: Vec<Observation> = observations
.iter()
.filter_map(|o| {
Some(Observation {
i: index_of(o.i)?,
j: index_of(o.j)?,
pi: o.pi,
pj: o.pj,
})
})
.collect();
let t_adjust = std::time::Instant::now();
let adjusted = bundle::adjust(start, &obs, &opts.adjust)?;
log::debug!(
"bundle adjustment: {} observations, {} iterations in {:?}",
obs.len(),
adjusted.iterations,
t_adjust.elapsed()
);
for (slot, &k) in aligned.iter().enumerate() {
rotations[k] = Some(adjusted.cameras.rotations[slot]);
}
Ok(Alignment {
rotations,
focal: adjusted.cameras.focal,
links,
unaligned,
rms_px: adjusted.rms_px,
})
}
#[cfg(test)]
mod tests {
use super::*;
use crate::features::{Keypoint, DESCRIPTOR_LEN};
use crate::linalg::Vec3;
/// Frames of a synthetic sweep: world directions with random unit
/// descriptors, each frame seeing the ones in its field of view.
fn synthetic_sweep(
n: usize,
step: f64,
f: f64,
w: usize,
h: usize,
) -> (Vec<Features>, Cameras) {
let mut seed = 777u64;
let mut rnd = || {
seed = seed
.wrapping_mul(6364136223846793005)
.wrapping_add(1442695040888963407);
((seed >> 33) as f64 / (1u64 << 31) as f64) - 0.5
};
let rotations: Vec<Mat3> = (0..n)
.map(|k| {
Mat3::rotation(Vec3::new(0.0, 1.0, 0.0), step * k as f64)
* Mat3::rotation(Vec3::new(1.0, 0.0, 0.0), 0.02 * ((k % 3) as f64 - 1.0))
})
.collect();
let truth = Cameras {
rotations,
focal: f,
};
let total = step * (n as f64 - 1.0);
let mut frames: Vec<Features> = (0..n)
.map(|_| Features {
keypoints: Vec::new(),
descriptors: Vec::new(),
width: w,
height: h,
})
.collect();
for _ in 0..600 * n {
let yaw = rnd() * (total + 0.8) + total / 2.0;
let pitch = rnd() * 0.5;
let d = Vec3::new(
yaw.sin() * pitch.cos(),
pitch.sin(),
yaw.cos() * pitch.cos(),
);
let desc: Vec<f32> = (0..DESCRIPTOR_LEN).map(|_| rnd() as f32).collect();
let norm = desc.iter().map(|v| v * v).sum::<f32>().sqrt();
let desc: Vec<f32> = desc.iter().map(|v| v / norm).collect();
for (k, frame) in frames.iter_mut().enumerate() {
if let Some(p) = truth.project(k, d) {
let (x, y) = (p.0 + w as f64 / 2.0, p.1 + h as f64 / 2.0);
if x >= 0.0 && x < w as f64 && y >= 0.0 && y < h as f64 {
frame.keypoints.push(Keypoint {
x: (x + rnd() * 0.6) as f32,
y: (y + rnd() * 0.6) as f32,
score: 1.0,
});
frame.descriptors.extend_from_slice(&desc);
}
}
}
}
(frames, truth)
}
fn angle_between(a: Mat3, b: Mat3) -> f64 {
(a.transpose() * b).log().norm()
}
#[test]
fn a_synthetic_sweep_is_aligned_to_its_truth() {
let (frames, truth) = synthetic_sweep(6, 0.3, 1400.0, 1024, 768);
let out = align(&frames, &AlignOptions::default()).expect("aligned");
assert!(out.is_complete(), "unaligned: {:?}", out.unaligned);
assert_eq!(out.links.len(), 5 + 4, "links: {}", out.links.len());
assert!((out.focal - 1400.0).abs() < 15.0, "focal {}", out.focal);
assert!(out.rms_px < 1.0, "rms {}", out.rms_px);
// Relative rotations match the truth's, whichever frame is the root.
let root = out
.rotations
.iter()
.position(|r| *r == Some(Mat3::IDENTITY))
.unwrap();
for k in 0..6 {
let rel_truth = truth.rotations[root].transpose() * truth.rotations[k];
let rel_out = out.rotations[k].unwrap();
let err = angle_between(rel_truth, rel_out);
assert!(err < 2e-3, "frame {k} off by {err} rad");
}
}
#[test]
fn a_frame_from_nowhere_is_named_not_guessed() {
let (mut frames, _) = synthetic_sweep(4, 0.3, 1400.0, 1024, 768);
// Frame 3 gets descriptors nobody else has.
for v in &mut frames[3].descriptors {
*v = -*v;
}
let out = align(&frames, &AlignOptions::default()).expect("aligned");
assert_eq!(out.unaligned.len(), 1);
assert_eq!(out.unaligned[0].0, 3);
assert!(out.rotations[3].is_none());
assert!(out.rotations[..3].iter().all(Option::is_some));
}
#[test]
fn a_frame_left_out_is_solved_without_measuring_again() {
let (frames, truth) = synthetic_sweep(6, 0.3, 1400.0, 1024, 768);
let opts = AlignOptions::default();
let pairs = match_pairs(&frames, &opts).expect("measured");
// The first frame left out: five cameras, indexed as the kept
// frames, and the links among them only.
let keep = [false, true, true, true, true, true];
let out = solve(&pairs, &keep, &opts).expect("solved");
assert!(out.is_complete(), "unaligned: {:?}", out.unaligned);
assert_eq!(out.rotations.len(), 5);
assert_eq!(out.links.len(), 4 + 3, "links: {}", out.links.len());
let root = out
.rotations
.iter()
.position(|r| *r == Some(Mat3::IDENTITY))
.unwrap();
for k in 0..5 {
let rel_truth = truth.rotations[root + 1].transpose() * truth.rotations[k + 1];
let err = angle_between(rel_truth, out.rotations[k].unwrap());
assert!(err < 2e-3, "frame {k} off by {err} rad");
}
// A frame in the middle left out splits the sweep only if nothing
// spans the gap; at 0.3 rad steps its neighbours still overlap.
let keep = [true, true, false, true, true, true];
let out = solve(&pairs, &keep, &opts).expect("solved");
assert!(out.is_complete(), "unaligned: {:?}", out.unaligned);
// And the whole set solved from the pairs is `align`'s answer.
assert_eq!(
solve(&pairs, &[true; 6], &opts).expect("solved"),
align(&frames, &opts).expect("aligned")
);
}
#[test]
fn one_frame_is_refused() {
let (frames, _) = synthetic_sweep(1, 0.3, 1400.0, 640, 480);
assert!(matches!(
align(&frames, &AlignOptions::default()),
Err(PanoError::Input(_))
));
}
}