diff --git a/Cargo.lock b/Cargo.lock index 6aa522a..dac2564 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1529,6 +1529,20 @@ dependencies = [ "log", ] +[[package]] +name = "dr-pano" +version = "0.12.2" +dependencies = [ + "dr-decode", + "dr-types", + "env_logger", + "log", + "ndarray", + "ort", + "ort-tract", + "thiserror 2.0.20", +] + [[package]] name = "dr-pipeline" version = "0.12.2" diff --git a/Cargo.toml b/Cargo.toml index 37fbdd3..7b3ae01 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -11,6 +11,7 @@ members = [ "core/dr-ingest", "core/dr-gpu", "core/dr-lens", + "core/dr-pano", "core/dr-pipeline", "core/dr-preset-xmp", "core/dr-segment", @@ -48,6 +49,8 @@ dr-film = { path = "core/dr-film" } dr-ingest = { path = "core/dr-ingest" } dr-gpu = { path = "core/dr-gpu" } dr-lens = { path = "core/dr-lens" } +# Optional runtime, like `dr-segment`: the geometry never needs a model. +dr-pano = { path = "core/dr-pano", default-features = false } dr-pipeline = { path = "core/dr-pipeline" } dr-preset-xmp = { path = "core/dr-preset-xmp" } # `default-features = false` belongs *here*, not on each dependant: a member diff --git a/core/dr-pano/Cargo.toml b/core/dr-pano/Cargo.toml new file mode 100644 index 0000000..a63c09c --- /dev/null +++ b/core/dr-pano/Cargo.toml @@ -0,0 +1,37 @@ +[package] +name = "dr-pano" +version.workspace = true +edition.workspace = true +rust-version.workspace = true +license.workspace = true +# Guards against a Git LFS pointer being embedded in place of the weights. +build = "build.rs" + +[dependencies] +thiserror.workspace = true +log.workspace = true + +# Inference for the learned keypoint detector, on the same footing as +# `dr-segment`: `ort` is the API, tract is the engine, and both are optional +# so that the geometry — matching, the rotation solve, the projections — is a +# dependency-free crate that tests without a model. +ort = { workspace = true, optional = true } +ort-tract = { workspace = true, optional = true } +ndarray = { workspace = true, optional = true } + +[dev-dependencies] +# The example aligns real frames from their embedded previews. +dr-decode.workspace = true +dr-types.workspace = true +env_logger.workspace = true + +[features] +default = ["xfeat", "embedded-model"] + +# The XFeat detector (FR-MRG-8). Off, the crate has no model and no runtime, +# and `Detector` has no implementation — a build that only wants the geometry. +xfeat = ["dep:ort", "dep:ort-tract", "dep:ndarray"] + +# Compile the weights into the binary, for the same reason `dr-segment` does: +# Android hands the app no path to read a model from (ARCH §6.9). +embedded-model = ["xfeat"] diff --git a/core/dr-pano/build.rs b/core/dr-pano/build.rs new file mode 100644 index 0000000..4d7de34 --- /dev/null +++ b/core/dr-pano/build.rs @@ -0,0 +1,37 @@ +//! Check the model is a model and not an LFS pointer. +//! +//! `models/keypoints/*.onnx` is stored in Git LFS (see `.gitattributes`). A +//! clone made without git-lfs, or with `GIT_LFS_SKIP_SMUDGE` set, leaves a +//! ~130-byte text pointer at that path instead of the weights, and +//! `include_bytes!` would embed it without complaint. Same guard as +//! `dr-segment`'s, for the same failure. + +use std::path::Path; + +const MODEL: &str = "../../models/keypoints/xfeat-1024.onnx"; + +fn main() { + println!("cargo:rerun-if-changed={MODEL}"); + println!("cargo:rerun-if-changed=build.rs"); + + if std::env::var_os("CARGO_FEATURE_EMBEDDED_MODEL").is_none() { + return; + } + + let path = Path::new(MODEL); + let Ok(bytes) = std::fs::read(path) else { + panic!( + "\n\n{MODEL} is missing.\n\ + It ships in Git LFS. Run `git lfs install && git lfs pull`, or build \ + with `--no-default-features` for a geometry-only build.\n" + ); + }; + + if bytes.starts_with(b"version https://git-lfs.github.com/spec/") { + panic!( + "\n\n{MODEL} is a Git LFS pointer, not the model.\n\ + Run `git lfs install && git lfs pull`, or build with \ + `--no-default-features` for a geometry-only build.\n" + ); + } +} diff --git a/core/dr-pano/src/align.rs b/core/dr-pano/src/align.rs new file mode 100644 index 0000000..e00b8b2 --- /dev/null +++ b/core/dr-pano/src/align.rs @@ -0,0 +1,421 @@ +//! 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`). +//! +//! 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, + 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>, + /// Focal length in pixels of the features' image. + pub focal: f64, + pub links: Vec, + 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, + } + } +} + +/// Align a set of frames from their features. +/// +/// 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 { + 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 = Vec::new(); + let mut matched_any = vec![false; n]; + for i in 0..n { + for j in i + 1..n { + let matches: Vec = match_features(&frames[i], &frames[j], opts.min_similarity); + if matches.len() < 4 { + continue; + } + matched_any[i] = true; + matched_any[j] = true; + 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; + 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, + }); + } + } + + // 3: the focal length. + let mut estimates: Vec = links + .iter() + .filter_map(|l| homography::focal_from_homography(&l.h)) + .filter(|f| f.is_finite() && *f > 0.0) + .collect(); + let longest = frames.iter().map(|f| f.width.max(f.height) as f64).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> = vec![None; n]; + let mut unaligned = Vec::new(); + if links.is_empty() { + for k in 0..n { + unaligned.push(( + k, + if matched_any[k] { + 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 = 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 = 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 adjusted = bundle::adjust(start, &obs, &opts.adjust)?; + 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, 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 = (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 = (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 = (0..DESCRIPTOR_LEN).map(|_| rnd() as f32).collect(); + let norm = desc.iter().map(|v| v * v).sum::().sqrt(); + let desc: Vec = desc.iter().map(|v| v / norm).collect(); + for k in 0..n { + 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 { + frames[k].keypoints.push(Keypoint { + x: (x + rnd() * 0.6) as f32, + y: (y + rnd() * 0.6) as f32, + score: 1.0, + }); + frames[k].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 one_frame_is_refused() { + let (frames, _) = synthetic_sweep(1, 0.3, 1400.0, 640, 480); + assert!(matches!( + align(&frames, &AlignOptions::default()), + Err(PanoError::Input(_)) + )); + } +} diff --git a/core/dr-pano/src/bundle.rs b/core/dr-pano/src/bundle.rs new file mode 100644 index 0000000..df3b236 --- /dev/null +++ b/core/dr-pano/src/bundle.rs @@ -0,0 +1,434 @@ +//! Bundle adjustment: every rotation and the focal length, refined together. +//! +//! The pairwise homographies (`homography.rs`) each know about two frames. +//! Chained around a loop they disagree with themselves by the accumulated +//! error, and a twelve-frame sweep chained end to end drifts by a visible +//! amount. This solves for all the rotations at once, against every inlier +//! match of every pair, so the error is spread rather than accumulated — +//! Brown & Lowe's step 4, with the camera model reduced to what a panorama +//! needs: one rotation per frame and one focal length shared by all. +//! +//! Levenberg–Marquardt with a numerical Jacobian. Analytic derivatives of a +//! rotation's projection are not hard, but they are a second place the +//! model is written down, and the model is small: forty parameters, a few +//! thousand residuals, a Jacobian that costs forty residual evaluations. +//! The whole solve is milliseconds. Correctness over cleverness, and one +//! definition of the projection to keep right. + +use crate::linalg::{DMat, Mat3, Vec3}; +use crate::PanoError; + +/// A point in one image, centred on the principal point, in pixels. +pub type Point = (f64, f64); + +/// One inlier match between two frames. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct Observation { + pub i: usize, + pub j: usize, + pub pi: Point, + pub pj: Point, +} + +/// What the adjustment starts from and returns: a rotation per frame +/// (camera to world; frame 0 is the world) and the focal length in pixels. +#[derive(Debug, Clone, PartialEq)] +pub struct Cameras { + pub rotations: Vec, + pub focal: f64, +} + +impl Cameras { + /// The unit direction, in world space, that pixel `p` of frame `i` looks + /// along. + pub fn bearing(&self, i: usize, p: Point) -> Vec3 { + self.rotations[i] * Vec3::new(p.0, p.1, self.focal).normalised() + } + + /// Where world direction `d` lands in frame `j`, or `None` if it is + /// behind the camera. + pub fn project(&self, j: usize, d: Vec3) -> Option { + let c = self.rotations[j].transpose() * d; + if c.z() <= 1e-9 { + return None; + } + Some((self.focal * c.x() / c.z(), self.focal * c.y() / c.z())) + } +} + +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct AdjustOptions { + pub max_iterations: usize, + /// Residuals beyond this many pixels are down-weighted (Huber), so a + /// mismatch RANSAC let through pulls with bounded force. + pub huber_px: f64, + /// Whether the focal length is a free parameter. Off, it is held at the + /// starting value — for a set whose rotations are all small, the focal + /// length is weakly observable and better taken from the homographies' + /// median than pulled about by noise. + pub refine_focal: bool, +} + +impl Default for AdjustOptions { + fn default() -> Self { + AdjustOptions { + max_iterations: 50, + huber_px: 3.0, + refine_focal: true, + } + } +} + +/// The adjusted cameras and the fit. +#[derive(Debug, Clone, PartialEq)] +pub struct Adjusted { + pub cameras: Cameras, + /// Root-mean-square reprojection error over all observations, in pixels + /// (unweighted, so an outlier RANSAC missed shows here rather than + /// hiding under its Huber weight). + pub rms_px: f64, + pub iterations: usize, +} + +/// Refine `start` against `observations`. +/// +/// Frame 0's rotation is held fixed: the world frame is arbitrary and +/// fixing one camera removes the freedom. Every other frame must appear in +/// at least one observation or its rotation is undetermined and the normal +/// equations are singular — the caller (`align`) guarantees it by only +/// adjusting frames a spanning tree reached. +pub fn adjust( + start: Cameras, + observations: &[Observation], + opts: &AdjustOptions, +) -> Result { + let n_frames = start.rotations.len(); + if n_frames < 2 || observations.is_empty() { + let rms = rms(&start, observations); + return Ok(Adjusted { + cameras: start, + rms_px: rms, + iterations: 0, + }); + } + // Every adjustable frame must be constrained by something, or its + // block of the normal equations is zero and the solve is meaningless — + // checked here, by name, rather than left to surface as a step that + // fails to lower the cost. + let mut seen = vec![false; n_frames]; + for o in observations { + seen[o.i] = true; + seen[o.j] = true; + } + if let Some(k) = (1..n_frames).find(|&k| !seen[k]) { + return Err(PanoError::Geometry(format!( + "frame {k} has no observations constraining it" + ))); + } + let n_rot = 3 * (n_frames - 1); + let n_params = n_rot + usize::from(opts.refine_focal); + let n_res = 2 * observations.len(); + + // Parameters are *increments* on the current cameras, re-applied each + // accepted step: rotation k ← exp(δ_k) · rotation k, focal ← f · exp(δ_f). + // Composing on the left keeps the increment in world space, where a + // small rotation means the same thing for every frame. + let apply = |base: &Cameras, x: &[f64]| -> Cameras { + let mut rotations = base.rotations.clone(); + for k in 1..n_frames { + let w = Vec3::new(x[3 * (k - 1)], x[3 * (k - 1) + 1], x[3 * (k - 1) + 2]); + rotations[k] = (Mat3::exp(w) * base.rotations[k]).orthonormalised(); + } + let focal = if opts.refine_focal { + base.focal * x[n_rot].exp() + } else { + base.focal + }; + Cameras { rotations, focal } + }; + + let residuals = |c: &Cameras, out: &mut Vec| { + out.clear(); + for o in observations { + let d = c.bearing(o.i, o.pi); + match c.project(o.j, d) { + Some((x, y)) => { + out.push(x - o.pj.0); + out.push(y - o.pj.1); + } + None => { + // Behind the camera: as wrong as a residual can be + // without being infinite. The Huber weight caps its pull. + out.push(1e4); + out.push(1e4); + } + } + } + }; + + let weights = |r: &[f64], out: &mut Vec| { + out.clear(); + for pair in r.chunks_exact(2) { + let m = (pair[0] * pair[0] + pair[1] * pair[1]).sqrt(); + let w = if m > opts.huber_px { + opts.huber_px / m + } else { + 1.0 + }; + out.push(w); + out.push(w); + } + }; + + // The robust cost itself, not the weighted sum of squares: the weights + // above are the IRLS linearisation for one step, and comparing two + // steps by sums taken under different weights would accept the wrong + // ones. Huber: quadratic within the threshold, linear beyond it. + let cost = |r: &[f64]| -> f64 { + r.chunks_exact(2) + .map(|pair| { + let m = (pair[0] * pair[0] + pair[1] * pair[1]).sqrt(); + if m <= opts.huber_px { + m * m + } else { + 2.0 * opts.huber_px * m - opts.huber_px * opts.huber_px + } + }) + .sum() + }; + + let mut cameras = start; + let mut r = Vec::with_capacity(n_res); + let mut w = Vec::with_capacity(n_res); + residuals(&cameras, &mut r); + weights(&r, &mut w); + let mut current = cost(&r); + + let mut lambda = 1e-3; + let mut jac = vec![0.0f64; n_res * n_params]; + let mut r_plus = Vec::with_capacity(n_res); + let zero = vec![0.0f64; n_params]; + let mut iterations = 0; + + for _ in 0..opts.max_iterations { + iterations += 1; + + // Numerical Jacobian about the current cameras (x = 0). + const H: f64 = 1e-6; + for p in 0..n_params { + let mut x = zero.clone(); + x[p] = H; + let c_plus = apply(&cameras, &x); + residuals(&c_plus, &mut r_plus); + for (k, (rp, r0)) in r_plus.iter().zip(&r).enumerate() { + jac[k * n_params + p] = (rp - r0) / H; + } + } + + // Normal equations, weighted: (JᵀWJ + λ·diag) δ = −JᵀWr. + let mut a = DMat::zeros(n_params); + let mut b = vec![0.0f64; n_params]; + for k in 0..n_res { + let row = &jac[k * n_params..(k + 1) * n_params]; + let wk = w[k]; + for p in 0..n_params { + b[p] -= wk * row[p] * r[k]; + for q in 0..n_params { + a[(p, q)] += wk * row[p] * row[q]; + } + } + } + + // Try steps with increasing damping until one lowers the cost. + let mut accepted = false; + for _ in 0..10 { + let mut damped = a.clone(); + for p in 0..n_params { + let d = a[(p, p)]; + damped[(p, p)] = d + lambda * d.max(1e-9); + } + let Some(delta) = damped.solve_spd(&b) else { + return Err(PanoError::Geometry( + "the adjustment's normal equations are singular: a frame has no \ + observations constraining it" + .into(), + )); + }; + let candidate = apply(&cameras, &delta); + residuals(&candidate, &mut r_plus); + let c_new = cost(&r_plus); + if c_new < current { + let improvement = (current - c_new) / current.max(1e-12); + let step: f64 = delta.iter().map(|d| d * d).sum::().sqrt(); + cameras = candidate; + std::mem::swap(&mut r, &mut r_plus); + weights(&r, &mut w); + current = c_new; + lambda = (lambda / 3.0).max(1e-9); + accepted = true; + // Converged when a *lightly damped* step no longer helps. A + // heavily damped step is small by construction and would + // pass an improvement test long before the minimum. + if step < 1e-10 || (improvement < 1e-8 && lambda < 1e-2) { + return Ok(Adjusted { + rms_px: rms(&cameras, observations), + cameras, + iterations, + }); + } + break; + } + lambda *= 5.0; + } + if !accepted { + break; + } + } + + Ok(Adjusted { + rms_px: rms(&cameras, observations), + cameras, + iterations, + }) +} + +/// Unweighted RMS reprojection error in pixels. +pub fn rms(c: &Cameras, observations: &[Observation]) -> f64 { + if observations.is_empty() { + return 0.0; + } + let sum: f64 = observations + .iter() + .map(|o| match c.project(o.j, c.bearing(o.i, o.pi)) { + Some((x, y)) => (x - o.pj.0).powi(2) + (y - o.pj.1).powi(2), + None => 1e8, + }) + .sum(); + (sum / observations.len() as f64).sqrt() +} + +#[cfg(test)] +mod tests { + use super::*; + + /// A synthetic sweep: `n` cameras panned by `step` radians each with a + /// little pitch and roll, `f` pixels, and matches between neighbours + /// from a cloud of world directions. + fn sweep(n: usize, step: f64, f: f64, noise_px: f64) -> (Cameras, Vec) { + let mut rotations = Vec::new(); + for k in 0..n { + let yaw = step * k as f64; + let pitch = 0.01 * ((k * 7) % 3) as f64; + let roll = 0.005 * ((k * 5) % 4) as f64; + let r = Mat3::rotation(Vec3::new(0.0, 1.0, 0.0), yaw) + * Mat3::rotation(Vec3::new(1.0, 0.0, 0.0), pitch) + * Mat3::rotation(Vec3::new(0.0, 0.0, 1.0), roll); + rotations.push(r); + } + let truth = Cameras { rotations, focal: f }; + + // World directions: a fan across the whole sweep. + let mut obs = Vec::new(); + let mut seed = 12345u64; + let mut rnd = || { + seed = seed + .wrapping_mul(6364136223846793005) + .wrapping_add(1442695040888963407); + ((seed >> 33) as f64 / (1u64 << 31) as f64) - 0.5 + }; + let total = step * (n as f64 - 1.0); + for _ in 0..400 * 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()) + .normalised(); + // Visible in which frames? Within ±0.35 f of centre. + let mut seen: Vec<(usize, Point)> = Vec::new(); + for k in 0..n { + if let Some(p) = truth.project(k, d) { + if p.0.abs() < 0.35 * f && p.1.abs() < 0.25 * f { + seen.push((k, (p.0 + rnd() * noise_px, p.1 + rnd() * noise_px))); + } + } + } + for a in 0..seen.len() { + for b in a + 1..seen.len() { + obs.push(Observation { + i: seen[a].0, + j: seen[b].0, + pi: seen[a].1, + pj: seen[b].1, + }); + } + } + } + (truth, obs) + } + + fn angle_between(a: Mat3, b: Mat3) -> f64 { + (a.transpose() * b).log().norm() + } + + #[test] + fn a_perturbed_start_converges_back_to_the_truth() { + let (truth, obs) = sweep(6, 0.3, 1400.0, 0.0); + assert!(obs.len() > 500); + // Perturb every rotation but the first by ~1°, and the focal by 5%. + let mut start = truth.clone(); + for k in 1..6 { + let w = Vec3::new(0.01, -0.015, 0.008) * (k as f64 / 3.0); + start.rotations[k] = Mat3::exp(w) * start.rotations[k]; + } + start.focal *= 1.05; + let before = rms(&start, &obs); + let out = adjust(start, &obs, &AdjustOptions::default()).expect("solvable"); + assert!(out.rms_px < 1e-3, "rms {} (was {before})", out.rms_px); + assert!((out.cameras.focal - 1400.0).abs() < 0.5, "focal {}", out.cameras.focal); + for k in 0..6 { + let err = angle_between(out.cameras.rotations[k], truth.rotations[k]); + assert!(err < 1e-5, "frame {k} off by {err} rad"); + } + } + + #[test] + fn noise_is_averaged_rather_than_accumulated() { + let (truth, obs) = sweep(8, 0.25, 1400.0, 1.0); + let mut start = truth.clone(); + for k in 1..8 { + start.rotations[k] = Mat3::exp(Vec3::new(0.0, 0.004 * k as f64, 0.0)) * start.rotations[k]; + } + let out = adjust(start, &obs, &AdjustOptions::default()).expect("solvable"); + // ±0.5 px of uniform noise on every coordinate has an RMS of 0.41 px + // per axis, so the fit's RMS over both axes should sit near 0.58 and + // cannot be much below it. + assert!(out.rms_px < 0.7, "rms {}", out.rms_px); + + // The focal length and the sweep are nearly degenerate for a + // single row: only the perspective inside each overlap pins the + // focal, and a pixel of noise is worth about a tenth of a percent of + // it. What that error does is scale every yaw by the same factor — + // a uniform stretch of the panorama, invisible in the result — so the + // absolute rotation error grows linearly along the sweep and is not + // the measure of the solve. The residual after removing that stretch + // is. + let f_ratio = out.cameras.focal / 1400.0; + assert!((f_ratio - 1.0).abs() < 5e-3, "focal {}", out.cameras.focal); + for k in 0..8 { + let yaw_k = 0.25 * k as f64; + let expected_stretch = (f_ratio - 1.0).abs() * yaw_k; + let err = angle_between(out.cameras.rotations[k], truth.rotations[k]); + assert!( + err < expected_stretch + 1.5e-4, + "frame {k} off by {err} rad, {expected_stretch} of it the focal's" + ); + } + } + + #[test] + fn a_frame_without_observations_is_refused() { + let (truth, mut obs) = sweep(4, 0.3, 1400.0, 0.0); + obs.retain(|o| o.i != 3 && o.j != 3); + let err = adjust(truth, &obs, &AdjustOptions::default()).unwrap_err(); + assert!(matches!(err, PanoError::Geometry(_))); + } +} diff --git a/core/dr-pano/src/features.rs b/core/dr-pano/src/features.rs new file mode 100644 index 0000000..8bd7841 --- /dev/null +++ b/core/dr-pano/src/features.rs @@ -0,0 +1,336 @@ +//! Keypoints with descriptors, and the decoder that reads them out of +//! XFeat's dense maps. +//! +//! The network (S15.2) produces three maps at an eighth of the input +//! resolution and stops; everything from there to a list of keypoints is +//! this file, in plain Rust, for the reason `dr-segment` decodes yolo26's +//! heads itself: the post-processing is cheap, shape-dependent and exactly +//! the kind of graph tract parses badly. It is a port of the reference +//! `XFeat.detectAndCompute`, step for step, so that a keypoint here is the +//! keypoint the paper's numbers were measured on. + +/// One detected point, in the pixel coordinates of the image it was +/// detected in, with the detector's confidence. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct Keypoint { + pub x: f32, + pub y: f32, + /// The reliability the detector assigned; higher is better, and the + /// scale is the detector's own — comparable within one model only. + pub score: f32, +} + +/// The keypoints of one image and their descriptors. +#[derive(Debug, Clone, PartialEq)] +pub struct Features { + pub keypoints: Vec, + /// `keypoints.len() × DESCRIPTOR_LEN`, each row L2-normalised, so that a + /// dot product between two rows is their cosine similarity. + pub descriptors: Vec, + /// The image the coordinates are in. + pub width: usize, + pub height: usize, +} + +/// The length of one descriptor. XFeat's is 64; the matcher does not care +/// what the number is, only that both sides agree. +pub const DESCRIPTOR_LEN: usize = 64; + +impl Features { + pub fn len(&self) -> usize { + self.keypoints.len() + } + + pub fn is_empty(&self) -> bool { + self.keypoints.is_empty() + } + + pub fn descriptor(&self, i: usize) -> &[f32] { + &self.descriptors[i * DESCRIPTOR_LEN..(i + 1) * DESCRIPTOR_LEN] + } +} + +/// XFeat's three output maps, as the network hands them back. +/// +/// All three are `channels × height × width` at an eighth of the input, in +/// the NCHW order the ONNX export declares (`feats [1, 64, H/8, W/8]`, +/// `keypoints [1, 65, H/8, W/8]`, `heatmap [1, 1, H/8, W/8]`). +pub struct XFeatMaps<'a> { + /// 64 channels: the dense descriptor field. + pub feats: &'a [f32], + /// 65 channels: for each 8×8 cell, a logit per position plus one for + /// "no keypoint here". + pub keypoints: &'a [f32], + /// 1 channel: reliability. + pub heatmap: &'a [f32], + /// The maps' width and height (the input's, divided by eight). + pub width: usize, + pub height: usize, +} + +/// How the decoder picks keypoints. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct DecodeOptions { + /// Keep at most this many, by score. The reference default is 4096. + pub top_k: usize, + /// A cell position's softmax probability must exceed this to be a + /// keypoint at all. The reference default is 0.05. + pub threshold: f32, + /// Ignore keypoints within this many pixels of the map's edge. A frame + /// padded into the detector's fixed input (`Gray::padded`) has a hard + /// edge where the padding starts, and the detector fires on it. + pub border: usize, +} + +impl Default for DecodeOptions { + fn default() -> Self { + DecodeOptions { + top_k: 4096, + threshold: 0.05, + border: 4, + } + } +} + +/// Decode keypoints and descriptors from the network's maps. +/// +/// The reference, step for step: +/// 1. softmax over the 65 logits of each cell, keep the 64 positions; +/// 2. pixel-shuffle those into a full-resolution keypoint heatmap — channel +/// `c` of cell `(cx, cy)` is pixel `(cx·8 + c%8, cy·8 + c/8)`; +/// 3. 5×5 non-maximum suppression over that heatmap, above `threshold`; +/// 4. score each survivor by its heatmap value times the reliability map +/// sampled bilinearly at its position; +/// 5. keep the `top_k` by score; +/// 6. sample the descriptor field bilinearly at each and L2-normalise. +/// +/// Bilinear where the reference samples the descriptor field bicubically: +/// a quarter-pixel's difference in a field that is smooth by construction, +/// and one interpolator rather than two to keep correct. +pub fn decode_xfeat(maps: &XFeatMaps<'_>, opts: &DecodeOptions) -> Features { + let (w8, h8) = (maps.width, maps.height); + let (w, h) = (w8 * 8, h8 * 8); + let cells = w8 * h8; + debug_assert_eq!(maps.keypoints.len(), 65 * cells); + debug_assert_eq!(maps.feats.len(), DESCRIPTOR_LEN * cells); + debug_assert_eq!(maps.heatmap.len(), cells); + + // 1 + 2: softmax per cell, scattered into the full-resolution heatmap. + let mut heat = vec![0.0f32; w * h]; + for cy in 0..h8 { + for cx in 0..w8 { + let cell = cy * w8 + cx; + let logit = |c: usize| maps.keypoints[c * cells + cell]; + let max = (0..65).map(logit).fold(f32::MIN, f32::max); + let mut sum = 0.0f32; + let mut exps = [0.0f32; 65]; + for (c, e) in exps.iter_mut().enumerate() { + *e = (logit(c) - max).exp(); + sum += *e; + } + for (c, e) in exps.iter().enumerate().take(64) { + let (dx, dy) = (c % 8, c / 8); + heat[(cy * 8 + dy) * w + cx * 8 + dx] = e / sum; + } + } + } + + // 3: a pixel survives if it is the maximum of its 5×5 neighbourhood and + // above threshold. Ties go to every tied pixel, as the reference's + // `x == max_pool(x)` does. + let border = opts.border.max(2); + let mut survivors: Vec<(usize, usize, f32)> = Vec::new(); + for y in border..h.saturating_sub(border) { + for x in border..w.saturating_sub(border) { + let v = heat[y * w + x]; + if v <= opts.threshold { + continue; + } + let mut is_max = true; + 'nb: for ny in y - 2..=y + 2 { + for nx in x - 2..=x + 2 { + if heat[ny * w + nx] > v { + is_max = false; + break 'nb; + } + } + } + if is_max { + survivors.push((x, y, v)); + } + } + } + + // 4: heatmap value × reliability, the latter sampled at the keypoint's + // position in map coordinates (`align_corners = False`: pixel `x` of the + // full image is `x / 8 - 0.5` in the map). + let sample = |field: &[f32], channels: usize, c: usize, x: f32, y: f32| -> f32 { + let fx = (x / 8.0 - 0.5).clamp(0.0, (w8 - 1) as f32); + let fy = (y / 8.0 - 0.5).clamp(0.0, (h8 - 1) as f32); + let x0 = fx as usize; + let y0 = fy as usize; + let x1 = (x0 + 1).min(w8 - 1); + let y1 = (y0 + 1).min(h8 - 1); + let tx = fx - x0 as f32; + let ty = fy - y0 as f32; + let at = |xx: usize, yy: usize| field[c * (w8 * h8) + yy * w8 + xx]; + let _ = channels; + let top = at(x0, y0) * (1.0 - tx) + at(x1, y0) * tx; + let bot = at(x0, y1) * (1.0 - tx) + at(x1, y1) * tx; + top * (1.0 - ty) + bot * ty + }; + let mut scored: Vec<(usize, usize, f32)> = survivors + .into_iter() + .map(|(x, y, v)| { + let r = sample(maps.heatmap, 1, 0, x as f32, y as f32); + (x, y, v * r) + }) + .collect(); + + // 5: best first, then cut. `sort_unstable_by` on a total order of the + // score; NaN cannot occur — every input is a probability or a sigmoid. + scored.sort_unstable_by(|a, b| b.2.total_cmp(&a.2)); + scored.truncate(opts.top_k); + + // 6: descriptors. + let mut keypoints = Vec::with_capacity(scored.len()); + let mut descriptors = Vec::with_capacity(scored.len() * DESCRIPTOR_LEN); + for (x, y, score) in scored { + let (xf, yf) = (x as f32, y as f32); + let start = descriptors.len(); + for c in 0..DESCRIPTOR_LEN { + descriptors.push(sample(maps.feats, DESCRIPTOR_LEN, c, xf, yf)); + } + let norm = descriptors[start..] + .iter() + .map(|v| v * v) + .sum::() + .sqrt() + .max(1e-12); + for v in &mut descriptors[start..] { + *v /= norm; + } + keypoints.push(Keypoint { + x: xf, + y: yf, + score, + }); + } + + Features { + keypoints, + descriptors, + width: w, + height: h, + } +} + +#[cfg(test)] +mod tests { + use super::*; + + /// Maps for a `w8 × h8` grid where every cell says "no keypoint" except + /// the listed ones, which put all their weight on one position. + fn maps(w8: usize, h8: usize, hot: &[(usize, usize, usize)]) -> (Vec, Vec, Vec) { + let cells = w8 * h8; + let mut kp = vec![0.0f32; 65 * cells]; + // "None" strongly preferred everywhere. + for cell in 0..cells { + kp[64 * cells + cell] = 10.0; + } + for &(cx, cy, c) in hot { + let cell = cy * w8 + cx; + kp[64 * cells + cell] = 0.0; + kp[c * cells + cell] = 10.0; + } + let heat = vec![0.5f32; cells]; + // Descriptors: channel c is constant c across the field, so any + // sampled descriptor is the same known vector. + let mut feats = vec![0.0f32; DESCRIPTOR_LEN * cells]; + for c in 0..DESCRIPTOR_LEN { + for v in &mut feats[c * cells..(c + 1) * cells] { + *v = c as f32; + } + } + (feats, kp, heat) + } + + #[test] + fn a_hot_cell_position_becomes_a_keypoint_at_the_right_pixel() { + // Cell (2, 1), channel 8*3 + 5 = 29 → pixel (2*8 + 5, 1*8 + 3). + let (f, k, h) = maps(8, 8, &[(2, 1, 29)]); + let out = decode_xfeat( + &XFeatMaps { + feats: &f, + keypoints: &k, + heatmap: &h, + width: 8, + height: 8, + }, + &DecodeOptions::default(), + ); + assert_eq!(out.len(), 1); + assert_eq!((out.keypoints[0].x, out.keypoints[0].y), (21.0, 11.0)); + assert_eq!((out.width, out.height), (64, 64)); + // Score is the softmax weight (~1) times the reliability (0.5). + assert!((out.keypoints[0].score - 0.5).abs() < 5e-3); + } + + #[test] + fn descriptors_are_unit_length() { + let (f, k, h) = maps(8, 8, &[(3, 3, 0), (5, 5, 63)]); + let out = decode_xfeat( + &XFeatMaps { + feats: &f, + keypoints: &k, + heatmap: &h, + width: 8, + height: 8, + }, + &DecodeOptions::default(), + ); + assert_eq!(out.len(), 2); + for i in 0..2 { + let n: f32 = out.descriptor(i).iter().map(|v| v * v).sum(); + assert!((n - 1.0).abs() < 1e-5); + } + } + + #[test] + fn top_k_keeps_the_best() { + let (f, k, mut h) = maps(8, 8, &[(1, 1, 0), (3, 3, 0), (5, 5, 0)]); + // Make cell (3, 3) the most reliable. + h[3 * 8 + 3] = 0.9; + let out = decode_xfeat( + &XFeatMaps { + feats: &f, + keypoints: &k, + heatmap: &h, + width: 8, + height: 8, + }, + &DecodeOptions { + top_k: 1, + ..Default::default() + }, + ); + assert_eq!(out.len(), 1); + assert_eq!((out.keypoints[0].x, out.keypoints[0].y), (24.0, 24.0)); + } + + #[test] + fn the_border_is_excluded() { + let (f, k, h) = maps(8, 8, &[(0, 0, 0)]); + let out = decode_xfeat( + &XFeatMaps { + feats: &f, + keypoints: &k, + heatmap: &h, + width: 8, + height: 8, + }, + &DecodeOptions::default(), + ); + assert!(out.is_empty()); + } +} diff --git a/core/dr-pano/src/homography.rs b/core/dr-pano/src/homography.rs new file mode 100644 index 0000000..f061a91 --- /dev/null +++ b/core/dr-pano/src/homography.rs @@ -0,0 +1,365 @@ +//! Pairwise geometry: a homography between two frames, found robustly. +//! +//! Two frames of a panorama are related by a rotation, and a rotation seen +//! through one lens is a homography of the image plane — `H = K R Kᵀ⁻¹`. The +//! homography is estimated first, from matches, because it does not need +//! the focal length; the focal length is then *read off* it (§ below), and +//! the rotation follows from both. This is the order Brown & Lowe (2007) +//! and OpenCV's stitcher use, and it is what makes the pipeline work when +//! EXIF says nothing about the lens. +//! +//! Coordinates throughout are **centred**: the principal point is the +//! origin. The focal formulae assume it, and centring before the DLT also +//! conditions the linear system — Hartley's normalisation, done once by the +//! caller rather than inside every solve. + +use crate::linalg::{DMat, Mat3, Vec3}; + +/// A point in one image, centred on the principal point. +pub type Point = (f64, f64); + +/// Apply a homography to a point. +pub fn apply(h: &Mat3, p: Point) -> Option { + let v = *h * Vec3::new(p.0, p.1, 1.0); + if v.z().abs() < 1e-12 { + return None; + } + Some((v.x() / v.z(), v.y() / v.z())) +} + +/// Least-squares homography from at least four correspondences by the +/// direct linear transform, with `h33` fixed at 1. +/// +/// Fixing `h33` turns the homogeneous 8×9 system into an ordinary 8-unknown +/// least-squares problem that the normal equations and a Cholesky +/// factorisation solve without an SVD. The one homography it cannot +/// represent — `h33 = 0`, a point at the origin mapped to infinity — does +/// not occur between overlapping frames of one scene. +/// +/// The points should be scaled to order one (divide by the focal length or +/// the image size) before calling: the normal equations square the +/// conditioning, and pixel coordinates in the thousands make them singular +/// in `f64`. +pub fn dlt(pairs: &[(Point, Point)]) -> Option { + if pairs.len() < 4 { + return None; + } + // Each pair gives two rows of A h = b with h = (h11..h32). + // x' = (h11 x + h12 y + h13) / (h31 x + h32 y + 1) + // → h11 x + h12 y + h13 - h31 x x' - h32 y x' = x' + let mut ata = DMat::zeros(8); + let mut atb = [0.0f64; 8]; + for &((x, y), (xp, yp)) in pairs { + let rows: [([f64; 8], f64); 2] = [ + ([x, y, 1.0, 0.0, 0.0, 0.0, -x * xp, -y * xp], xp), + ([0.0, 0.0, 0.0, x, y, 1.0, -x * yp, -y * yp], yp), + ]; + for (a, b) in rows { + for i in 0..8 { + atb[i] += a[i] * b; + for j in 0..8 { + ata[(i, j)] += a[i] * a[j]; + } + } + } + } + let h = ata.solve_spd(&atb)?; + Some(Mat3([ + [h[0], h[1], h[2]], + [h[3], h[4], h[5]], + [h[6], h[7], 1.0], + ])) +} + +/// A homography with the correspondences that agree with it. +#[derive(Debug, Clone, PartialEq)] +pub struct RobustHomography { + pub h: Mat3, + /// Indices into the input pairs. + pub inliers: Vec, +} + +/// RANSAC over [`dlt`] on four-point samples, then a final least-squares +/// fit over every inlier. +/// +/// `threshold` is the reprojection distance, in the same units as the +/// points, within which a pair counts as agreeing. The iteration count +/// adapts to the inlier ratio found so far in the usual way, capped at +/// `max_iterations`. `seed` makes a run reproducible (NFR-MRG-2): the +/// sampling is a small linear congruential generator, not the system's. +pub fn ransac_homography( + pairs: &[(Point, Point)], + threshold: f64, + max_iterations: usize, + seed: u64, +) -> Option { + if pairs.len() < 4 { + return None; + } + let n = pairs.len(); + let thr2 = threshold * threshold; + let mut rng = Lcg(seed); + let mut best: Option<(Vec, Mat3)> = None; + let mut iterations = max_iterations; + let mut i = 0; + while i < iterations { + i += 1; + let sample = rng.distinct4(n); + let Some(h) = dlt(&sample.map(|k| pairs[k])) else { + continue; + }; + let inliers: Vec = (0..n) + .filter(|&k| agrees(&h, pairs[k], thr2)) + .collect(); + if best.as_ref().is_none_or(|(b, _)| inliers.len() > b.len()) { + // Adapt: enough iterations to have drawn one all-inlier sample + // with probability 0.99, given the ratio seen so far. + let w = inliers.len() as f64 / n as f64; + let p_all = w.powi(4); + if p_all > 0.0 && p_all < 1.0 { + let needed = ((1.0 - 0.99f64).ln() / (1.0 - p_all).ln()).ceil() as usize; + iterations = iterations.min(needed.max(i + 1)); + } + best = Some((inliers, h)); + } + } + let (inliers, h) = best?; + if inliers.len() < 4 { + return None; + } + // Refit on every inlier, and keep the refit only if it did not lose + // support — a least-squares fit over a set with a few borderline points + // can be pulled off the consensus the sample found. + let refit: Vec<(Point, Point)> = inliers.iter().map(|&k| pairs[k]).collect(); + let h = match dlt(&refit) { + Some(r) => { + let count = (0..n).filter(|&k| agrees(&r, pairs[k], thr2)).count(); + if count >= inliers.len() { + r + } else { + h + } + } + None => h, + }; + let inliers: Vec = (0..n).filter(|&k| agrees(&h, pairs[k], thr2)).collect(); + Some(RobustHomography { h, inliers }) +} + +fn agrees(h: &Mat3, (p, q): (Point, Point), thr2: f64) -> bool { + match apply(h, p) { + Some((x, y)) => { + let (dx, dy) = (x - q.0, y - q.1); + dx * dx + dy * dy <= thr2 + } + None => false, + } +} + +/// The focal length a homography implies, if it implies one. +/// +/// For `H = K R K⁻¹` with `K = diag(f, f, 1)` and the principal point at the +/// origin, the orthonormality of `R` gives two independent estimates of `f²` +/// from the first two rows and two from the first two columns; each is +/// taken where it is positive and the better-conditioned of the pair is +/// chosen, as OpenCV's `focalsFromHomography` does. The geometric mean of +/// the row and column estimates is returned. `None` when the homography is +/// too close to a pure translation to say anything — every estimate is then +/// a ratio of small numbers. +pub fn focal_from_homography(h: &Mat3) -> Option { + let m = h.0; + let (h00, h01, h02) = (m[0][0], m[0][1], m[0][2]); + let (h10, h11, h12) = (m[1][0], m[1][1], m[1][2]); + let (h20, h21) = (m[2][0], m[2][1]); + + let pick = |mut v1: f64, mut v2: f64, d1: f64, d2: f64| -> Option { + if v1 < v2 { + std::mem::swap(&mut v1, &mut v2); + } + if v1 > 0.0 && v2 > 0.0 { + Some((if d1.abs() > d2.abs() { v1 } else { v2 }).sqrt()) + } else if v1 > 0.0 { + Some(v1.sqrt()) + } else { + None + } + }; + + // From the third row. + let d1 = h20 * h21; + let d2 = (h21 - h20) * (h21 + h20); + let f1 = if d1.abs() > 1e-12 || d2.abs() > 1e-12 { + let v1 = if d1.abs() > 1e-12 { + -(h00 * h01 + h10 * h11) / d1 + } else { + f64::NAN + }; + let v2 = if d2.abs() > 1e-12 { + (h00 * h00 + h10 * h10 - h01 * h01 - h11 * h11) / d2 + } else { + f64::NAN + }; + pick(nan_to_neg(v1), nan_to_neg(v2), d1, d2) + } else { + None + }; + + // From the third column. + let d1 = h00 * h10 + h01 * h11; + let d2 = h00 * h00 + h01 * h01 - h10 * h10 - h11 * h11; + let f0 = if d1.abs() > 1e-12 || d2.abs() > 1e-12 { + let v1 = if d1.abs() > 1e-12 { + -h02 * h12 / d1 + } else { + f64::NAN + }; + let v2 = if d2.abs() > 1e-12 { + (h12 * h12 - h02 * h02) / d2 + } else { + f64::NAN + }; + pick(nan_to_neg(v1), nan_to_neg(v2), d1, d2) + } else { + None + }; + + match (f0, f1) { + (Some(a), Some(b)) => Some((a * b).sqrt()), + (Some(a), None) | (None, Some(a)) => Some(a), + (None, None) => None, + } +} + +fn nan_to_neg(v: f64) -> f64 { + if v.is_finite() { + v + } else { + -1.0 + } +} + +/// The rotation a homography encodes for a known focal length: +/// `R = K⁻¹ H K`, re-orthonormalised, with the scale of `H` divided out. +pub fn rotation_from_homography(h: &Mat3, f: f64) -> Mat3 { + let m = h.0; + // K⁻¹ H K with K = diag(f, f, 1): scale the third row by f and the + // third column by 1/f. + let r = Mat3([ + [m[0][0], m[0][1], m[0][2] / f], + [m[1][0], m[1][1], m[1][2] / f], + [m[2][0] * f, m[2][1] * f, m[2][2]], + ]); + r.orthonormalised() +} + +/// A small deterministic generator for RANSAC's samples. +struct Lcg(u64); + +impl Lcg { + fn next(&mut self) -> u64 { + // Knuth's MMIX constants. + self.0 = self + .0 + .wrapping_mul(6364136223846793005) + .wrapping_add(1442695040888963407); + self.0 >> 33 + } + + fn below(&mut self, n: usize) -> usize { + (self.next() % n as u64) as usize + } + + fn distinct4(&mut self, n: usize) -> [usize; 4] { + let mut s = [0usize; 4]; + for i in 0..4 { + loop { + let k = self.below(n); + if !s[..i].contains(&k) { + s[i] = k; + break; + } + } + } + s + } +} + +#[cfg(test)] +mod tests { + use super::*; + + /// Points under a known rotation seen through a known focal length, + /// in centred image coordinates scaled by that focal length. + fn synthetic(f: f64, r: Mat3, n: usize, noise: f64, seed: u64) -> Vec<(Point, Point)> { + let mut rng = Lcg(seed); + let mut out = Vec::new(); + while out.len() < n { + // A point on the first image plane, within ±0.3 f of centre. + let x = (rng.below(6001) as f64 - 3000.0) / 10000.0; + let y = (rng.below(4001) as f64 - 2000.0) / 10000.0; + let b = Vec3::new(x, y, 1.0); + let v = r * b; + if v.z() <= 0.2 { + continue; + } + let nx = (rng.below(2001) as f64 - 1000.0) / 1000.0 * noise; + let ny = (rng.below(2001) as f64 - 1000.0) / 1000.0 * noise; + out.push(((x, y), (v.x() / v.z() + nx, v.y() / v.z() + ny))); + } + let _ = f; + out + } + + #[test] + fn dlt_recovers_a_known_homography_exactly() { + let r = Mat3::exp(Vec3::new(0.05, 0.3, 0.02)); + let pairs = synthetic(1.0, r, 12, 0.0, 1); + let h = dlt(&pairs).expect("solvable"); + for &(p, q) in &pairs { + let (x, y) = apply(&h, p).unwrap(); + assert!((x - q.0).abs() < 1e-9 && (y - q.1).abs() < 1e-9); + } + } + + #[test] + fn ransac_finds_the_consensus_among_outliers() { + let r = Mat3::exp(Vec3::new(-0.02, 0.25, 0.01)); + let mut pairs = synthetic(1.0, r, 60, 0.0005, 2); + // Forty outliers: wrong second point. + let mut rng = Lcg(9); + for _ in 0..40 { + let k = rng.below(60); + let (p, _) = pairs[k]; + pairs.push((p, ((rng.below(1000) as f64 - 500.0) / 1000.0, 0.1))); + } + let robust = ransac_homography(&pairs, 0.003, 500, 3).expect("found"); + assert!(robust.inliers.len() >= 55, "{} inliers", robust.inliers.len()); + assert!(robust.inliers.iter().all(|&k| k < 60)); + } + + #[test] + fn focal_is_read_off_a_rotation_homography() { + // H in *pixel* coordinates for f = 1400: K R K⁻¹. + let f = 1400.0; + let r = Mat3::exp(Vec3::new(0.03, 0.35, -0.01)); + let m = r.0; + let h = Mat3([ + [m[0][0], m[0][1], m[0][2] * f], + [m[1][0], m[1][1], m[1][2] * f], + [m[2][0] / f, m[2][1] / f, m[2][2]], + ]); + let est = focal_from_homography(&h).expect("estimable"); + assert!((est - f).abs() / f < 1e-6, "{est}"); + let back = rotation_from_homography(&h, f); + for i in 0..3 { + for j in 0..3 { + assert!((back.0[i][j] - m[i][j]).abs() < 1e-9); + } + } + } + + #[test] + fn the_identity_implies_no_focal() { + assert!(focal_from_homography(&Mat3::IDENTITY).is_none()); + } +} diff --git a/core/dr-pano/src/image.rs b/core/dr-pano/src/image.rs new file mode 100644 index 0000000..d53ab09 --- /dev/null +++ b/core/dr-pano/src/image.rs @@ -0,0 +1,259 @@ +//! The grayscale proxy a detector reads. +//! +//! Alignment runs on proxies (FR-MRG-7) — a detector at 1024 px sees +//! everything it needs, and the full-resolution frames never leave the GPU. +//! This is that proxy: one channel, `f32` in `0.0..=1.0`, upright, and no +//! larger than the detector's fixed input. + +/// A single-channel image, row-major, values in `0.0..=1.0`. +#[derive(Debug, Clone, PartialEq)] +pub struct Gray { + pub width: usize, + pub height: usize, + pub data: Vec, +} + +impl Gray { + /// From tightly packed 8-bit RGBA, by the Rec. 709 luma weights. + /// + /// The proxy is what a detector looks at, not what the photographer + /// sees, so which luma is used matters less than that it is the same one + /// for every frame — a keypoint's descriptor must not change between two + /// frames because they were converted differently. + pub fn from_rgba8(rgba: &[u8], width: usize, height: usize) -> Gray { + let n = width * height; + assert!(rgba.len() >= n * 4, "rgba buffer is short for {width}×{height}"); + let data = rgba[..n * 4] + .chunks_exact(4) + .map(|p| { + (0.2126 * f32::from(p[0]) + 0.7152 * f32::from(p[1]) + 0.0722 * f32::from(p[2])) + / 255.0 + }) + .collect(); + Gray { + width, + height, + data, + } + } + + /// Apply an EXIF orientation so the image is upright. + /// + /// Learned detectors are not rotation-invariant — a descriptor of a + /// feature seen sideways is a different descriptor — and a portrait set + /// (the 6D fixture is one) would match poorly or not at all fed as + /// stored. The camera says which way is up; the proxy is turned before + /// anything looks at it, and the composite is written upright. + /// + /// The value is the EXIF `Orientation` tag. Mirrored values (2, 4, 5, 7) + /// are not produced by any camera and are treated as their unmirrored + /// counterparts. + pub fn oriented(&self, orientation: u16) -> Gray { + match orientation { + 3 | 4 => self.rotated_180(), + 6 | 5 => self.rotated_90_cw(), + 8 | 7 => self.rotated_90_ccw(), + _ => self.clone(), + } + } + + fn rotated_90_cw(&self) -> Gray { + let (w, h) = (self.width, self.height); + let mut data = vec![0.0; w * h]; + for y in 0..h { + for x in 0..w { + // Source (x, y) lands at (h - 1 - y, x) in an h-wide image. + data[x * h + (h - 1 - y)] = self.data[y * w + x]; + } + } + Gray { + width: h, + height: w, + data, + } + } + + fn rotated_90_ccw(&self) -> Gray { + let (w, h) = (self.width, self.height); + let mut data = vec![0.0; w * h]; + for y in 0..h { + for x in 0..w { + // Source (x, y) lands at (y, w - 1 - x) in an h-wide image. + data[(w - 1 - x) * h + y] = self.data[y * w + x]; + } + } + Gray { + width: h, + height: w, + data, + } + } + + fn rotated_180(&self) -> Gray { + let mut data = self.data.clone(); + data.reverse(); + Gray { + width: self.width, + height: self.height, + data, + } + } + + /// Resample to exactly `width × height` by area averaging on the way + /// down and bilinear on the way up. + /// + /// Area averaging, not point sampling, for a reduction: a 5472 px frame + /// to 1024 is a factor of five, and picking one source pixel in + /// twenty-five aliases every edge the detector is looking for. + pub fn resampled(&self, width: usize, height: usize) -> Gray { + if width == self.width && height == self.height { + return self.clone(); + } + let mut data = vec![0.0f32; width * height]; + let sx = self.width as f64 / width as f64; + let sy = self.height as f64 / height as f64; + if sx >= 1.0 && sy >= 1.0 { + for oy in 0..height { + let y0 = (oy as f64 * sy) as usize; + let y1 = (((oy + 1) as f64 * sy) as usize).clamp(y0 + 1, self.height); + for ox in 0..width { + let x0 = (ox as f64 * sx) as usize; + let x1 = (((ox + 1) as f64 * sx) as usize).clamp(x0 + 1, self.width); + let mut sum = 0.0f32; + for y in y0..y1 { + let row = &self.data[y * self.width..(y + 1) * self.width]; + sum += row[x0..x1].iter().sum::(); + } + data[oy * width + ox] = sum / ((y1 - y0) * (x1 - x0)) as f32; + } + } + } else { + for oy in 0..height { + let fy = ((oy as f64 + 0.5) * sy - 0.5).max(0.0); + let y0 = (fy as usize).min(self.height - 1); + let y1 = (y0 + 1).min(self.height - 1); + let ty = (fy - y0 as f64) as f32; + for ox in 0..width { + let fx = ((ox as f64 + 0.5) * sx - 0.5).max(0.0); + let x0 = (fx as usize).min(self.width - 1); + let x1 = (x0 + 1).min(self.width - 1); + let tx = (fx - x0 as f64) as f32; + let p = |x: usize, y: usize| self.data[y * self.width + x]; + let top = p(x0, y0) * (1.0 - tx) + p(x1, y0) * tx; + let bot = p(x0, y1) * (1.0 - tx) + p(x1, y1) * tx; + data[oy * width + ox] = top * (1.0 - ty) + bot * ty; + } + } + } + Gray { + width, + height, + data, + } + } + + /// Scale so the image fits inside `max_width × max_height`, preserving + /// aspect, never enlarging. Returns the image and the scale applied, + /// which is what maps a proxy keypoint back to the source. + pub fn fitted(&self, max_width: usize, max_height: usize) -> (Gray, f64) { + let scale = (max_width as f64 / self.width as f64) + .min(max_height as f64 / self.height as f64) + .min(1.0); + let w = ((self.width as f64 * scale).round() as usize).max(1); + let h = ((self.height as f64 * scale).round() as usize).max(1); + (self.resampled(w, h), w as f64 / self.width as f64) + } + + /// Copy into the top-left of a `width × height` canvas, zero elsewhere. + /// + /// The detector's input is a fixed shape (S15.2), and a frame that fits + /// inside it is padded rather than stretched: stretching changes the + /// aspect and with it every descriptor. + pub fn padded(&self, width: usize, height: usize) -> Gray { + assert!(self.width <= width && self.height <= height); + let mut data = vec![0.0; width * height]; + for y in 0..self.height { + data[y * width..y * width + self.width] + .copy_from_slice(&self.data[y * self.width..(y + 1) * self.width]); + } + Gray { + width, + height, + data, + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + + fn ramp(w: usize, h: usize) -> Gray { + Gray { + width: w, + height: h, + data: (0..w * h).map(|i| i as f32).collect(), + } + } + + #[test] + fn rotating_four_quarter_turns_is_the_identity() { + let g = ramp(5, 3); + let mut r = g.clone(); + for _ in 0..4 { + r = r.rotated_90_cw(); + } + assert_eq!(r, g); + assert_eq!(g.rotated_90_cw().rotated_90_ccw(), g); + assert_eq!(g.rotated_180().rotated_180(), g); + } + + #[test] + fn a_clockwise_turn_moves_the_top_left_to_the_top_right() { + // 2×3 image, pixel values by position. + let g = ramp(2, 3); + let r = g.rotated_90_cw(); + assert_eq!((r.width, r.height), (3, 2)); + // Top-left of source (value 0) is at top-right of result. + assert_eq!(r.data[2], 0.0); + // Bottom-left of source (value 4) is at top-left of result. + assert_eq!(r.data[0], 4.0); + } + + #[test] + fn orientation_8_is_a_counter_clockwise_turn() { + let g = ramp(4, 2); + assert_eq!(g.oriented(8), g.rotated_90_ccw()); + assert_eq!(g.oriented(6), g.rotated_90_cw()); + assert_eq!(g.oriented(1), g); + } + + #[test] + fn downsampling_by_two_averages_blocks() { + let g = Gray { + width: 4, + height: 2, + data: vec![0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0], + }; + let r = g.resampled(2, 1); + assert_eq!(r.data, vec![2.5, 4.5]); + } + + #[test] + fn fitting_never_enlarges_and_reports_the_scale() { + let g = ramp(100, 50); + let (f, s) = g.fitted(1024, 768); + assert_eq!((f.width, f.height), (100, 50)); + assert_eq!(s, 1.0); + let (f, s) = g.fitted(50, 50); + assert_eq!((f.width, f.height), (50, 25)); + assert_eq!(s, 0.5); + } + + #[test] + fn padding_places_the_image_at_the_origin() { + let g = ramp(2, 2); + let p = g.padded(3, 3); + assert_eq!(p.data, vec![0.0, 1.0, 0.0, 2.0, 3.0, 0.0, 0.0, 0.0, 0.0]); + } +} diff --git a/core/dr-pano/src/lib.rs b/core/dr-pano/src/lib.rs new file mode 100644 index 0000000..76c754b --- /dev/null +++ b/core/dr-pano/src/lib.rs @@ -0,0 +1,65 @@ +//! TRACES: FR-MRG-1 | FR-MRG-10 +//! Panorama geometry — from several frames to the rotations that relate +//! them, and the projections that lay them out. +//! +//! This is the CPU half of a merge (FR-MRG-10): keypoints, matching, the +//! rotation solve and the choice of output surface. The per-pixel half — +//! rendering, warping, seams, blending — is the GPU's and lives in +//! `dr-gpu`, driven from above; nothing here touches a full-resolution +//! pixel. The split is the whole design (panorama.md §4): everything in +//! this crate is bounded by the number of frames, not the size of the +//! composite, and runs on proxies. +//! +//! # Layout +//! +//! - [`image`] — the grayscale proxy a detector reads: oriented, resampled. +//! - [`features`] — keypoints with descriptors, and the XFeat decoder. +//! - [`xfeat`] — the network under tract (feature `xfeat`). +//! - [`matching`] — mutual nearest neighbours. +//! - [`homography`] — a robust pairwise homography, the focal length read +//! off it, and the rotation it implies. +//! - [`bundle`] — every rotation and the focal length refined together. +//! - [`align`] — the whole thing, from features to cameras, honest about +//! what it could not place. +//! - [`projection`] — perspective, cylindrical, spherical. +//! - [`linalg`] — the small dense algebra all of it uses. +//! +//! # What it depends on +//! +//! Nothing, without the `xfeat` feature: the geometry is pure Rust with +//! hand-rolled linear algebra (`linalg` says why) so that it tests without +//! a model, a GPU or a device, on synthetic sets whose answer is known +//! exactly. With the feature it adds the same `ort`-over-tract runtime the +//! rest of the application already carries. + +pub mod align; +pub mod bundle; +pub mod features; +pub mod homography; +pub mod image; +pub mod linalg; +pub mod matching; +pub mod projection; +#[cfg(feature = "xfeat")] +pub mod xfeat; + +pub use align::{align, AlignOptions, Alignment, Link, Unaligned}; +pub use bundle::Cameras; +pub use features::{Features, Keypoint}; +pub use image::Gray; +pub use projection::Projection; + +#[derive(Debug, thiserror::Error)] +pub enum PanoError { + #[error("bad input: {0}")] + Input(String), + #[error("geometry: {0}")] + Geometry(String), + #[error("model: {0}")] + Model(String), + #[error("could not read the model: {0}")] + ModelRead(#[source] std::io::Error), + #[cfg(feature = "xfeat")] + #[error("inference: {0}")] + Inference(#[source] ort::Error), +} diff --git a/core/dr-pano/src/linalg.rs b/core/dr-pano/src/linalg.rs new file mode 100644 index 0000000..9db7dcf --- /dev/null +++ b/core/dr-pano/src/linalg.rs @@ -0,0 +1,366 @@ +//! The small dense linear algebra the geometry needs, and nothing more. +//! +//! Hand-rolled rather than pulled in, and the decision was made on purpose +//! (2026-09-19): the largest system this crate ever solves is a rotation +//! per frame plus one focal length — forty unknowns for a dozen frames — +//! and everything else is three-vectors. A general linear-algebra crate +//! would be the largest dependency in `dr-pano` by an order of magnitude, +//! for a Cholesky factorisation that is thirty lines. +//! +//! `f64` throughout. The geometry is solved once per merge on a few thousand +//! matches; there is no reason to give up precision for speed here, and the +//! bundle adjustment's normal equations are poorly conditioned enough near +//! convergence that `f32` would stall it. + +use std::ops::{Add, Index, IndexMut, Mul, Neg, Sub}; + +/// A vector in three dimensions. +#[derive(Debug, Clone, Copy, PartialEq, Default)] +pub struct Vec3(pub [f64; 3]); + +impl Vec3 { + pub const fn new(x: f64, y: f64, z: f64) -> Self { + Vec3([x, y, z]) + } + + pub fn dot(self, o: Vec3) -> f64 { + self.0[0] * o.0[0] + self.0[1] * o.0[1] + self.0[2] * o.0[2] + } + + pub fn cross(self, o: Vec3) -> Vec3 { + Vec3([ + self.0[1] * o.0[2] - self.0[2] * o.0[1], + self.0[2] * o.0[0] - self.0[0] * o.0[2], + self.0[0] * o.0[1] - self.0[1] * o.0[0], + ]) + } + + pub fn norm(self) -> f64 { + self.dot(self).sqrt() + } + + /// The unit vector along `self`, or `self` unchanged if it is zero. + pub fn normalised(self) -> Vec3 { + let n = self.norm(); + if n > 0.0 { + self * (1.0 / n) + } else { + self + } + } + + pub fn x(self) -> f64 { + self.0[0] + } + pub fn y(self) -> f64 { + self.0[1] + } + pub fn z(self) -> f64 { + self.0[2] + } +} + +impl Add for Vec3 { + type Output = Vec3; + fn add(self, o: Vec3) -> Vec3 { + Vec3([self.0[0] + o.0[0], self.0[1] + o.0[1], self.0[2] + o.0[2]]) + } +} + +impl Sub for Vec3 { + type Output = Vec3; + fn sub(self, o: Vec3) -> Vec3 { + Vec3([self.0[0] - o.0[0], self.0[1] - o.0[1], self.0[2] - o.0[2]]) + } +} + +impl Mul for Vec3 { + type Output = Vec3; + fn mul(self, s: f64) -> Vec3 { + Vec3([self.0[0] * s, self.0[1] * s, self.0[2] * s]) + } +} + +impl Neg for Vec3 { + type Output = Vec3; + fn neg(self) -> Vec3 { + Vec3([-self.0[0], -self.0[1], -self.0[2]]) + } +} + +/// A 3×3 matrix, row-major. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct Mat3(pub [[f64; 3]; 3]); + +impl Mat3 { + pub const IDENTITY: Mat3 = Mat3([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]); + + /// The matrix whose columns are `a`, `b`, `c`. + pub fn from_columns(a: Vec3, b: Vec3, c: Vec3) -> Mat3 { + Mat3([ + [a.0[0], b.0[0], c.0[0]], + [a.0[1], b.0[1], c.0[1]], + [a.0[2], b.0[2], c.0[2]], + ]) + } + + pub fn transpose(self) -> Mat3 { + let m = self.0; + Mat3([ + [m[0][0], m[1][0], m[2][0]], + [m[0][1], m[1][1], m[2][1]], + [m[0][2], m[1][2], m[2][2]], + ]) + } + + pub fn column(self, i: usize) -> Vec3 { + Vec3([self.0[0][i], self.0[1][i], self.0[2][i]]) + } + + pub fn trace(self) -> f64 { + self.0[0][0] + self.0[1][1] + self.0[2][2] + } + + /// The rotation about `axis` (any length) by `angle` radians — Rodrigues. + pub fn rotation(axis: Vec3, angle: f64) -> Mat3 { + let k = axis.normalised(); + let (s, c) = angle.sin_cos(); + let t = 1.0 - c; + let (x, y, z) = (k.0[0], k.0[1], k.0[2]); + Mat3([ + [t * x * x + c, t * x * y - s * z, t * x * z + s * y], + [t * x * y + s * z, t * y * y + c, t * y * z - s * x], + [t * x * z - s * y, t * y * z + s * x, t * z * z + c], + ]) + } + + /// The rotation whose axis-angle vector is `w` (direction is the axis, + /// length is the angle). The exponential map; [`Self::log`] inverts it. + pub fn exp(w: Vec3) -> Mat3 { + let angle = w.norm(); + if angle < 1e-12 { + // First-order: I + [w]×, which is what the limit is and avoids + // dividing by the angle. + let (x, y, z) = (w.0[0], w.0[1], w.0[2]); + return Mat3([[1.0, -z, y], [z, 1.0, -x], [-y, x, 1.0]]); + } + Mat3::rotation(w, angle) + } + + /// The axis-angle vector of a rotation matrix. Inverse of [`Self::exp`]. + pub fn log(self) -> Vec3 { + let m = self.0; + let cos = ((self.trace() - 1.0) * 0.5).clamp(-1.0, 1.0); + let axis = Vec3([m[2][1] - m[1][2], m[0][2] - m[2][0], m[1][0] - m[0][1]]); + if cos > 1.0 - 1e-6 { + // Small angle: `acos` near 1 loses everything below ~1e-8 to + // rounding, but the antisymmetric part is `2 sin θ · axis` and + // keeps it. First order, exact to the precision that matters. + return axis * 0.5; + } + let angle = cos.acos(); + if angle > std::f64::consts::PI - 1e-6 { + // Near π the antisymmetric part vanishes; take the axis from the + // symmetric part instead. Rare for a panorama, but the solver may + // pass through it on a bad start and must not return NaN. + let d = Vec3([ + ((m[0][0] + 1.0) * 0.5).max(0.0).sqrt(), + ((m[1][1] + 1.0) * 0.5).max(0.0).sqrt(), + ((m[2][2] + 1.0) * 0.5).max(0.0).sqrt(), + ]); + return d.normalised() * angle; + } + axis * (angle / (2.0 * angle.sin())) + } + + /// Re-orthonormalise a matrix that has drifted from a rotation through + /// accumulated products. Gram–Schmidt on the columns; cheap and adequate + /// for drift of the size floating-point products produce. + pub fn orthonormalised(self) -> Mat3 { + let a = self.column(0).normalised(); + let b = (self.column(1) - a * a.dot(self.column(1))).normalised(); + let c = a.cross(b); + Mat3::from_columns(a, b, c) + } +} + +impl Mul for Mat3 { + type Output = Vec3; + fn mul(self, v: Vec3) -> Vec3 { + let m = self.0; + Vec3([ + m[0][0] * v.0[0] + m[0][1] * v.0[1] + m[0][2] * v.0[2], + m[1][0] * v.0[0] + m[1][1] * v.0[1] + m[1][2] * v.0[2], + m[2][0] * v.0[0] + m[2][1] * v.0[1] + m[2][2] * v.0[2], + ]) + } +} + +impl Mul for Mat3 { + type Output = Mat3; + fn mul(self, o: Mat3) -> Mat3 { + let mut r = [[0.0; 3]; 3]; + for (i, row) in r.iter_mut().enumerate() { + for (j, cell) in row.iter_mut().enumerate() { + *cell = (0..3).map(|k| self.0[i][k] * o.0[k][j]).sum(); + } + } + Mat3(r) + } +} + +/// A dense square matrix, for the normal equations. +#[derive(Debug, Clone, PartialEq)] +pub struct DMat { + n: usize, + data: Vec, +} + +impl DMat { + pub fn zeros(n: usize) -> DMat { + DMat { + n, + data: vec![0.0; n * n], + } + } + + pub fn n(&self) -> usize { + self.n + } + + /// Solve `self · x = b` for a symmetric positive-definite `self` by + /// Cholesky factorisation. `None` if the matrix is not positive definite, + /// which for the normal equations means the problem is not determined by + /// the data — a frame with no matches, for instance — and the caller + /// should say so rather than proceed. + /// + /// Destroys neither input: the factor is built in a copy. The systems + /// here are at most a few dozen unknowns and the copy is nothing. + pub fn solve_spd(&self, b: &[f64]) -> Option> { + let n = self.n; + debug_assert_eq!(b.len(), n); + let mut l = vec![0.0; n * n]; + for j in 0..n { + let mut d = self[(j, j)]; + for k in 0..j { + d -= l[j * n + k] * l[j * n + k]; + } + if d <= 0.0 || !d.is_finite() { + return None; + } + let djj = d.sqrt(); + l[j * n + j] = djj; + for i in j + 1..n { + let mut s = self[(i, j)]; + for k in 0..j { + s -= l[i * n + k] * l[j * n + k]; + } + l[i * n + j] = s / djj; + } + } + // Forward: L y = b. + let mut y = vec![0.0; n]; + for i in 0..n { + let mut s = b[i]; + for k in 0..i { + s -= l[i * n + k] * y[k]; + } + y[i] = s / l[i * n + i]; + } + // Back: Lᵀ x = y. + let mut x = vec![0.0; n]; + for i in (0..n).rev() { + let mut s = y[i]; + for k in i + 1..n { + s -= l[k * n + i] * x[k]; + } + x[i] = s / l[i * n + i]; + } + Some(x) + } +} + +impl Index<(usize, usize)> for DMat { + type Output = f64; + fn index(&self, (i, j): (usize, usize)) -> &f64 { + &self.data[i * self.n + j] + } +} + +impl IndexMut<(usize, usize)> for DMat { + fn index_mut(&mut self, (i, j): (usize, usize)) -> &mut f64 { + &mut self.data[i * self.n + j] + } +} + +#[cfg(test)] +mod tests { + use super::*; + + fn close(a: f64, b: f64) -> bool { + (a - b).abs() < 1e-9 + } + + #[test] + fn exp_and_log_are_inverses() { + for w in [ + Vec3::new(0.1, -0.2, 0.3), + Vec3::new(1.0, 0.0, 0.0), + Vec3::new(0.0, 0.0, 2.5), + Vec3::new(1e-9, 0.0, 0.0), + ] { + let back = Mat3::exp(w).log(); + for i in 0..3 { + assert!(close(back.0[i], w.0[i]), "{w:?} -> {back:?}"); + } + } + } + + #[test] + fn a_rotation_is_orthonormal_and_preserves_length() { + let r = Mat3::exp(Vec3::new(0.4, 0.5, -0.6)); + let rt = r.transpose() * r; + for i in 0..3 { + for j in 0..3 { + assert!(close(rt.0[i][j], Mat3::IDENTITY.0[i][j])); + } + } + let v = Vec3::new(1.0, 2.0, 3.0); + assert!(close((r * v).norm(), v.norm())); + } + + #[test] + fn rotation_about_z_turns_x_towards_y() { + let r = Mat3::rotation(Vec3::new(0.0, 0.0, 1.0), std::f64::consts::FRAC_PI_2); + let v = r * Vec3::new(1.0, 0.0, 0.0); + assert!(close(v.x(), 0.0) && close(v.y(), 1.0) && close(v.z(), 0.0)); + } + + #[test] + fn cholesky_solves_a_small_spd_system() { + // A = Bᵀ B for a random-ish B is SPD by construction. + let b = [[2.0, 1.0, 0.0], [1.0, 3.0, 1.0], [0.0, 1.0, 4.0], [1.0, 1.0, 1.0]]; + let mut a = DMat::zeros(3); + for i in 0..3 { + for j in 0..3 { + a[(i, j)] = (0..4).map(|k| b[k][i] * b[k][j]).sum(); + } + } + let x_true = [1.0, -2.0, 0.5]; + let rhs: Vec = (0..3) + .map(|i| (0..3).map(|j| a[(i, j)] * x_true[j]).sum()) + .collect(); + let x = a.solve_spd(&rhs).expect("spd"); + for i in 0..3 { + assert!(close(x[i], x_true[i]), "{x:?}"); + } + } + + #[test] + fn cholesky_refuses_an_indefinite_matrix() { + let mut a = DMat::zeros(2); + a[(0, 0)] = 1.0; + a[(1, 1)] = -1.0; + assert!(a.solve_spd(&[1.0, 1.0]).is_none()); + } +} diff --git a/core/dr-pano/src/matching.rs b/core/dr-pano/src/matching.rs new file mode 100644 index 0000000..f716c5a --- /dev/null +++ b/core/dr-pano/src/matching.rs @@ -0,0 +1,142 @@ +//! Descriptor matching between two images. +//! +//! Mutual nearest neighbour on cosine similarity, with a floor on the +//! similarity — the reference XFeat's own matcher (`match_mkpts`, +//! `min_cossim = 0.82`). For a panorama that is enough: one lens, one +//! scene, near-pure rotation and 20–40 % overlap make the matching problem +//! easy, and what is hard — sky, repeated structure, exposure drift — is +//! handled by the detector's descriptors and by RANSAC downstream, not by a +//! cleverer matcher. A learned matcher (LightGlue) is the step after this +//! one fails on a real set, and it has not (panorama.md §6). +//! +//! Brute force. `4096 × 4096 × 64` multiply-adds is a billion, which is +//! tens of milliseconds a pair on one core, and there are at most a few +//! dozen pairs. Not worth an index. + +use crate::features::{Features, DESCRIPTOR_LEN}; + +/// A correspondence: keypoint `a` in the first image matches keypoint `b` +/// in the second, with the cosine similarity of their descriptors. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct Match { + pub a: usize, + pub b: usize, + pub similarity: f32, +} + +/// Match two sets of features. +/// +/// A pair is kept when each is the other's nearest neighbour and their +/// similarity is at least `min_similarity`. +pub fn match_features(a: &Features, b: &Features, min_similarity: f32) -> Vec { + if a.is_empty() || b.is_empty() { + return Vec::new(); + } + let best_ab = nearest(a, b); + let best_ba = nearest(b, a); + best_ab + .iter() + .enumerate() + .filter_map(|(ia, &(ib, sim))| { + (best_ba[ib].0 == ia && sim >= min_similarity).then_some(Match { + a: ia, + b: ib, + similarity: sim, + }) + }) + .collect() +} + +/// For each descriptor in `from`, the index of its nearest in `to` and the +/// similarity. +fn nearest(from: &Features, to: &Features) -> Vec<(usize, f32)> { + (0..from.len()) + .map(|i| { + let d = from.descriptor(i); + let mut best = (0usize, f32::MIN); + for j in 0..to.len() { + let s = dot(d, to.descriptor(j)); + if s > best.1 { + best = (j, s); + } + } + best + }) + .collect() +} + +#[inline] +fn dot(a: &[f32], b: &[f32]) -> f32 { + // Written as a plain loop over a fixed length so the compiler + // vectorises it; the length is a constant and the slices are exact. + let mut s = 0.0f32; + for k in 0..DESCRIPTOR_LEN { + s += a[k] * b[k]; + } + s +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::features::Keypoint; + + /// Features whose descriptors are unit vectors along the given axes. + fn along(axes: &[usize]) -> Features { + let mut descriptors = vec![0.0; axes.len() * DESCRIPTOR_LEN]; + for (i, &ax) in axes.iter().enumerate() { + descriptors[i * DESCRIPTOR_LEN + ax] = 1.0; + } + Features { + keypoints: axes + .iter() + .map(|_| Keypoint { + x: 0.0, + y: 0.0, + score: 1.0, + }) + .collect(), + descriptors, + width: 1, + height: 1, + } + } + + #[test] + fn identical_descriptors_match_mutually() { + let a = along(&[0, 1, 2]); + let b = along(&[2, 0, 1]); + let m = match_features(&a, &b, 0.8); + let mut pairs: Vec<(usize, usize)> = m.iter().map(|m| (m.a, m.b)).collect(); + pairs.sort(); + assert_eq!(pairs, vec![(0, 1), (1, 2), (2, 0)]); + assert!(m.iter().all(|m| (m.similarity - 1.0).abs() < 1e-6)); + } + + #[test] + fn a_descriptor_with_no_counterpart_is_unmatched() { + let a = along(&[0, 1, 5]); + let b = along(&[0, 1]); + let m = match_features(&a, &b, 0.8); + assert_eq!(m.len(), 2); + assert!(m.iter().all(|m| m.a != 2)); + } + + #[test] + fn mutuality_breaks_a_one_sided_match() { + // b0 is the nearest to both a0 and a1, but a0 is its nearest — a1 + // must not be matched to it. + let mut a = along(&[0, 0]); + a.descriptors[DESCRIPTOR_LEN] = 0.9; + a.descriptors[DESCRIPTOR_LEN + 1] = (1.0f32 - 0.81).sqrt(); + let b = along(&[0]); + let m = match_features(&a, &b, 0.0); + assert_eq!(m.len(), 1); + assert_eq!((m[0].a, m[0].b), (0, 0)); + } + + #[test] + fn empty_input_is_empty_output() { + assert!(match_features(&along(&[]), &along(&[1]), 0.5).is_empty()); + } +} diff --git a/core/dr-pano/src/projection.rs b/core/dr-pano/src/projection.rs new file mode 100644 index 0000000..c2bc08f --- /dev/null +++ b/core/dr-pano/src/projection.rs @@ -0,0 +1,186 @@ +//! TRACES: FR-MRG-4 +//! The surface the composite is drawn on. +//! +//! A panorama is a set of directions; a picture is a plane. The projection +//! is the map between them, and the three offered are the three every +//! stitcher offers because each is right for a different field of view: +//! perspective keeps straight lines straight and cannot reach 180°; +//! cylindrical keeps verticals vertical and stretches nothing horizontally, +//! for the wide single row; spherical for anything that also looks up. +//! +//! Every function here is the *inverse* map — output pixel to direction — +//! because that is what a gather needs (`lens.rs` in `dr-pipeline` says +//! why a warp is written that way), and it is the function the WGSL warp +//! will repeat verbatim. The forward map exists for bounds only. + +use crate::linalg::Vec3; + +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum Projection { + Perspective, + Cylindrical, + Spherical, +} + +impl Projection { + /// Which projection a field of view calls for. + /// + /// Perspective stretches the edges by `1 / cos` of the angle from the + /// centre, which is 2× at 60° and unbounded at 90°; the switch is where + /// that stretch starts to look like a mistake. Spherical is for a set + /// that spans enough vertically that a cylinder would stretch the top + /// and bottom the same way. + pub fn suggest(horizontal_fov: f64, vertical_fov: f64) -> Projection { + if horizontal_fov < 70f64.to_radians() && vertical_fov < 70f64.to_radians() { + Projection::Perspective + } else if vertical_fov < 100f64.to_radians() { + Projection::Cylindrical + } else { + Projection::Spherical + } + } + + /// The direction an output point looks along. `scale` is the output's + /// focal length in pixels: the radius of the cylinder or sphere, or the + /// plane's distance. Coordinates are centred on the projection's origin + /// (the direction `+z`). + pub fn to_direction(self, scale: f64, u: f64, v: f64) -> Vec3 { + match self { + Projection::Perspective => Vec3::new(u, v, scale).normalised(), + Projection::Cylindrical => { + let theta = u / scale; + Vec3::new(theta.sin(), v / scale, theta.cos()).normalised() + } + Projection::Spherical => { + let theta = u / scale; + let phi = v / scale; + Vec3::new(theta.sin() * phi.cos(), phi.sin(), theta.cos() * phi.cos()) + } + } + } + + /// Where a direction lands on the output, or `None` where the + /// projection cannot show it (behind a perspective plane, at a + /// cylinder's poles). + pub fn from_direction(self, scale: f64, d: Vec3) -> Option<(f64, f64)> { + let (x, y, z) = (d.x(), d.y(), d.z()); + match self { + Projection::Perspective => (z > 1e-9).then(|| (scale * x / z, scale * y / z)), + Projection::Cylindrical => { + let r = (x * x + z * z).sqrt(); + (r > 1e-9).then(|| (scale * x.atan2(z), scale * y / r)) + } + Projection::Spherical => { + let r = (x * x + z * z).sqrt(); + Some((scale * x.atan2(z), scale * y.atan2(r))) + } + } + } +} + +/// The output rectangle a set of frames covers, in centred output pixels. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct Bounds { + pub min_u: f64, + pub min_v: f64, + pub max_u: f64, + pub max_v: f64, +} + +impl Bounds { + pub fn width(&self) -> f64 { + self.max_u - self.min_u + } + pub fn height(&self) -> f64 { + self.max_v - self.min_v + } +} + +/// Bounds of the frames' footprints under `projection`, by walking each +/// frame's border. +/// +/// `frame_size` is the frames' width and height in the same pixels the +/// cameras' focal length is in. The border is sampled rather than only its +/// corners because under a cylinder the widest point of a rolled frame is +/// not a corner. +pub fn bounds( + projection: Projection, + scale: f64, + cameras: &crate::bundle::Cameras, + frame_size: (f64, f64), +) -> Option { + let (w, h) = frame_size; + let mut b: Option = None; + let steps = 64; + for k in 0..cameras.rotations.len() { + for s in 0..steps { + let t = s as f64 / steps as f64; + for p in [ + (-w / 2.0 + w * t, -h / 2.0), + (-w / 2.0 + w * t, h / 2.0), + (-w / 2.0, -h / 2.0 + h * t), + (w / 2.0, -h / 2.0 + h * t), + ] { + let d = cameras.bearing(k, p); + let Some((u, v)) = projection.from_direction(scale, d) else { + continue; + }; + b = Some(match b { + None => Bounds { + min_u: u, + min_v: v, + max_u: u, + max_v: v, + }, + Some(b) => Bounds { + min_u: b.min_u.min(u), + min_v: b.min_v.min(v), + max_u: b.max_u.max(u), + max_v: b.max_v.max(v), + }, + }); + } + } + } + b +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn to_and_from_direction_are_inverses() { + for proj in [Projection::Perspective, Projection::Cylindrical, Projection::Spherical] { + for (u, v) in [(0.0, 0.0), (300.0, -200.0), (-900.0, 450.0)] { + let d = proj.to_direction(1000.0, u, v); + let (bu, bv) = proj.from_direction(1000.0, d).expect("in front"); + assert!((bu - u).abs() < 1e-9 && (bv - v).abs() < 1e-9, "{proj:?} {u} {v}"); + } + } + } + + #[test] + fn the_origin_looks_down_z_in_every_projection() { + for proj in [Projection::Perspective, Projection::Cylindrical, Projection::Spherical] { + let d = proj.to_direction(500.0, 0.0, 0.0); + assert!((d.z() - 1.0).abs() < 1e-12); + } + } + + #[test] + fn a_cylinder_maps_ninety_degrees_to_a_quarter_turn_of_pixels() { + let d = Vec3::new(1.0, 0.0, 0.0); + let (u, v) = Projection::Cylindrical.from_direction(100.0, d).unwrap(); + assert!((u - 100.0 * std::f64::consts::FRAC_PI_2).abs() < 1e-9); + assert_eq!(v, 0.0); + assert!(Projection::Perspective.from_direction(100.0, d).is_none()); + } + + #[test] + fn suggestion_widens_with_the_field() { + assert_eq!(Projection::suggest(0.5, 0.5), Projection::Perspective); + assert_eq!(Projection::suggest(2.5, 0.8), Projection::Cylindrical); + assert_eq!(Projection::suggest(3.0, 2.5), Projection::Spherical); + } +} diff --git a/core/dr-pano/src/xfeat.rs b/core/dr-pano/src/xfeat.rs new file mode 100644 index 0000000..3cb8b49 --- /dev/null +++ b/core/dr-pano/src/xfeat.rs @@ -0,0 +1,137 @@ +//! TRACES: FR-MRG-8 +//! The XFeat detector — the network under tract, and the decoder after it. +//! +//! Apache-2.0 weights (`models/LICENCE.md`), exported at a fixed shape by +//! `tools/export-xfeat.sh` and loaded through the same `ort`-over-tract +//! backend `dr-segment` and `dr-face` use, so this adds no runtime and no C +//! to the tree. ~300 ms per frame on the reference desktop, ~400 ms on the +//! tablet (S15.2, S15.4). + +use crate::features::{decode_xfeat, DecodeOptions, Features, XFeatMaps, DESCRIPTOR_LEN}; +use crate::image::Gray; +use crate::PanoError; + +/// The input shape the shipped export was made for. A different size is a +/// different file (`tools/export-xfeat.sh`). +pub const INPUT_WIDTH: usize = 1024; +pub const INPUT_HEIGHT: usize = 768; + +#[cfg(feature = "embedded-model")] +const EMBEDDED_MODEL: &[u8] = include_bytes!("../../../models/keypoints/xfeat-1024.onnx"); + +/// A loaded detector. +pub struct XFeat { + session: ort::session::Session, + pub options: DecodeOptions, +} + +impl XFeat { + /// The weights compiled into the binary. + #[cfg(feature = "embedded-model")] + pub fn embedded() -> Result { + Self::from_bytes(EMBEDDED_MODEL) + } + + pub fn from_path(path: &std::path::Path) -> Result { + let bytes = std::fs::read(path).map_err(PanoError::ModelRead)?; + Self::from_bytes(&bytes) + } + + pub fn from_bytes(bytes: &[u8]) -> Result { + install_backend(); + let session = ort::session::Session::builder() + .map_err(PanoError::Inference)? + .commit_from_memory(bytes) + .map_err(PanoError::Inference)?; + Ok(XFeat { + session, + options: DecodeOptions::default(), + }) + } + + /// Detect keypoints in an upright grayscale image. + /// + /// The image is fitted into the network's fixed input — scaled down if + /// larger, never up, and padded to the right and bottom — and the + /// keypoints come back in the coordinates of `image` itself, so a + /// caller that already scaled a frame to a proxy maps them on with the + /// scale it used and nothing else. + pub fn detect(&mut self, image: &Gray) -> Result { + let (fitted, scale) = image.fitted(INPUT_WIDTH, INPUT_HEIGHT); + let padded = fitted.padded(INPUT_WIDTH, INPUT_HEIGHT); + + let input = ndarray::Array::from_shape_vec( + ndarray::IxDyn(&[1, 1, INPUT_HEIGHT, INPUT_WIDTH]), + padded.data, + ) + .expect("shape matches the buffer by construction"); + let tensor = ort::value::Tensor::from_array(input).map_err(PanoError::Inference)?; + let outputs = self + .session + .run(ort::inputs![tensor]) + .map_err(PanoError::Inference)?; + + let (w8, h8) = (INPUT_WIDTH / 8, INPUT_HEIGHT / 8); + let expect = |i: usize, channels: usize| -> Result, PanoError> { + let (shape, data) = outputs[i] + .try_extract_tensor::() + .map_err(PanoError::Inference)?; + let dims: Vec = shape.iter().copied().collect(); + if dims != [1, channels as i64, h8 as i64, w8 as i64] { + return Err(PanoError::Model(format!( + "output {i} is {dims:?}, expected [1, {channels}, {h8}, {w8}] — \ + not the export this decoder was written for" + ))); + } + Ok(data.to_vec()) + }; + let feats = expect(0, DESCRIPTOR_LEN)?; + let keypoints = expect(1, 65)?; + let heatmap = expect(2, 1)?; + + let mut features = decode_xfeat( + &XFeatMaps { + feats: &feats, + keypoints: &keypoints, + heatmap: &heatmap, + width: w8, + height: h8, + }, + &self.options, + ); + + // Back to the caller's image: drop anything the padding produced, + // undo the fit. + let border = self.options.border as f32; + let limit_x = fitted.width as f32 - border; + let limit_y = fitted.height as f32 - border; + let mut kept_kp = Vec::with_capacity(features.len()); + let mut kept_desc = Vec::with_capacity(features.descriptors.len()); + for (i, kp) in features.keypoints.iter().enumerate() { + if kp.x >= limit_x || kp.y >= limit_y { + continue; + } + kept_kp.push(crate::features::Keypoint { + x: (kp.x / scale as f32), + y: (kp.y / scale as f32), + score: kp.score, + }); + kept_desc.extend_from_slice(features.descriptor(i)); + } + features.keypoints = kept_kp; + features.descriptors = kept_desc; + features.width = image.width; + features.height = image.height; + Ok(features) + } +} + +fn install_backend() { + use std::sync::Once; + static ONCE: Once = Once::new(); + ONCE.call_once(|| { + // False if another crate installed it first, which is fine: there is + // one backend compiled in for it to have chosen. + let _ = ort::set_api(ort_tract::api()); + }); +} diff --git a/docs/traceability.md b/docs/traceability.md index c7f3c6a..4daa791 100644 --- a/docs/traceability.md +++ b/docs/traceability.md @@ -9,18 +9,18 @@ Denominators are parsed from [`requirements.md`](requirements.md) at run time, n | Metric | Value | |---|---| -| Source files scanned | 356 | -| TRACES tags found | 1504 | +| Source files scanned | 368 | +| TRACES tags found | 1508 | | Requirements defined | 184 | | Requirements deferred (post-v1) | 24 | -| Requirements covered | 144 | -| **Coverage** | **78.3%** (144/184) | +| Requirements covered | 148 | +| **Coverage** | **80.4%** (148/184) | ### By type | Type | Covered | Defined | |---|---|---| -| FR | 106 | 127 | +| FR | 110 | 127 | | NFR | 34 | 50 | | R | 4 | 7 | @@ -98,8 +98,12 @@ _None._ | FR-EXP-7 | [`ui/dr-ui/src/activity.rs:83`](../ui/dr-ui/src/activity.rs#L83), [`ui/dr-ui/src/export.rs:1`](../ui/dr-ui/src/export.rs#L1), [`ui/dr-ui/src/export.rs:992`](../ui/dr-ui/src/export.rs#L992), [`ui/dr-ui/src/lib.rs:262`](../ui/dr-ui/src/lib.rs#L262), [`ui/dr-ui/src/lib.rs:2752`](../ui/dr-ui/src/lib.rs#L2752), [`ui/dr-ui/src/lib.rs:601`](../ui/dr-ui/src/lib.rs#L601), [`ui/dr-ui/src/lib.rs:647`](../ui/dr-ui/src/lib.rs#L647), [`ui/dr-ui/src/lib.rs:674`](../ui/dr-ui/src/lib.rs#L674), [`ui/dr-ui/src/library_ui.rs:4275`](../ui/dr-ui/src/library_ui.rs#L4275), [`ui/dr-ui/src/library_ui.rs:719`](../ui/dr-ui/src/library_ui.rs#L719), [`ui/dr-ui/src/library_ui.rs:8026`](../ui/dr-ui/src/library_ui.rs#L8026), [`ui/dr-ui/src/library_ui.rs:8103`](../ui/dr-ui/src/library_ui.rs#L8103), [`ui/dr-ui/src/library_ui.rs:8115`](../ui/dr-ui/src/library_ui.rs#L8115), [`ui/dr-ui/src/library_ui.rs:819`](../ui/dr-ui/src/library_ui.rs#L819), [`ui/dr-ui/src/library_ui.rs:876`](../ui/dr-ui/src/library_ui.rs#L876), [`ui/dr-ui/ui/app.slint:1091`](../ui/dr-ui/ui/app.slint#L1091), [`ui/dr-ui/ui/app.slint:1630`](../ui/dr-ui/ui/app.slint#L1630), [`ui/dr-ui/ui/library.slint:1583`](../ui/dr-ui/ui/library.slint#L1583), [`ui/dr-ui/ui/library.slint:4278`](../ui/dr-ui/ui/library.slint#L4278) | | FR-EXP-8 | [`core/dr-decode/src/lib.rs:330`](../core/dr-decode/src/lib.rs#L330), [`core/dr-decode/src/lib.rs:354`](../core/dr-decode/src/lib.rs#L354), [`core/dr-decode/src/lib.rs:368`](../core/dr-decode/src/lib.rs#L368), [`core/dr-decode/src/lib.rs:71`](../core/dr-decode/src/lib.rs#L71), [`core/dr-decode/src/lib.rs:79`](../core/dr-decode/src/lib.rs#L79), [`core/dr-decode/src/lib.rs:82`](../core/dr-decode/src/lib.rs#L82), [`core/dr-decode/src/locate.rs:1164`](../core/dr-decode/src/locate.rs#L1164), [`core/dr-decode/src/locate.rs:1223`](../core/dr-decode/src/locate.rs#L1223), [`core/dr-decode/src/locate.rs:316`](../core/dr-decode/src/locate.rs#L316), [`core/dr-decode/src/locate.rs:487`](../core/dr-decode/src/locate.rs#L487), [`core/dr-decode/src/locate.rs:571`](../core/dr-decode/src/locate.rs#L571), [`core/dr-decode/src/locate.rs:584`](../core/dr-decode/src/locate.rs#L584), [`core/dr-decode/src/locate.rs:667`](../core/dr-decode/src/locate.rs#L667), [`core/dr-export/examples/export.rs:99`](../core/dr-export/examples/export.rs#L99), [`core/dr-export/src/encode.rs:117`](../core/dr-export/src/encode.rs#L117), [`core/dr-export/src/encode.rs:161`](../core/dr-export/src/encode.rs#L161), [`core/dr-export/src/encode.rs:1`](../core/dr-export/src/encode.rs#L1), [`core/dr-export/src/encode.rs:206`](../core/dr-export/src/encode.rs#L206), [`core/dr-export/src/encode.rs:235`](../core/dr-export/src/encode.rs#L235), [`core/dr-export/src/encode.rs:311`](../core/dr-export/src/encode.rs#L311), [`core/dr-export/src/encode.rs:325`](../core/dr-export/src/encode.rs#L325), [`core/dr-export/src/encode.rs:408`](../core/dr-export/src/encode.rs#L408), [`core/dr-export/src/encode.rs:456`](../core/dr-export/src/encode.rs#L456), [`core/dr-export/src/encode.rs:70`](../core/dr-export/src/encode.rs#L70), [`core/dr-export/src/encode.rs:795`](../core/dr-export/src/encode.rs#L795), [`core/dr-export/src/encode.rs:809`](../core/dr-export/src/encode.rs#L809), [`core/dr-export/src/encode.rs:850`](../core/dr-export/src/encode.rs#L850), [`core/dr-export/src/encode.rs:898`](../core/dr-export/src/encode.rs#L898), [`core/dr-export/src/exif.rs:1`](../core/dr-export/src/exif.rs#L1), [`core/dr-export/src/lib.rs:136`](../core/dr-export/src/lib.rs#L136), [`core/dr-export/src/metadata.rs:1`](../core/dr-export/src/metadata.rs#L1), [`core/dr-export/src/metadata.rs:41`](../core/dr-export/src/metadata.rs#L41), [`core/dr-export/src/metadata.rs:74`](../core/dr-export/src/metadata.rs#L74), [`core/dr-types/src/lib.rs:655`](../core/dr-types/src/lib.rs#L655), [`core/dr-types/src/settings.rs:647`](../core/dr-types/src/settings.rs#L647), [`ui/dr-ui/src/develop.rs:4539`](../ui/dr-ui/src/develop.rs#L4539), [`ui/dr-ui/src/develop.rs:4651`](../ui/dr-ui/src/develop.rs#L4651), [`ui/dr-ui/src/develop.rs:729`](../ui/dr-ui/src/develop.rs#L729), [`ui/dr-ui/src/export.rs:467`](../ui/dr-ui/src/export.rs#L467), [`ui/dr-ui/src/export.rs:634`](../ui/dr-ui/src/export.rs#L634), [`ui/dr-ui/src/export.rs:680`](../ui/dr-ui/src/export.rs#L680), [`ui/dr-ui/src/export.rs:702`](../ui/dr-ui/src/export.rs#L702), [`ui/dr-ui/src/export.rs:826`](../ui/dr-ui/src/export.rs#L826), [`ui/dr-ui/src/export.rs:844`](../ui/dr-ui/src/export.rs#L844), [`ui/dr-ui/src/lib.rs:262`](../ui/dr-ui/src/lib.rs#L262), [`ui/dr-ui/src/lib.rs:285`](../ui/dr-ui/src/lib.rs#L285), [`ui/dr-ui/src/lib.rs:633`](../ui/dr-ui/src/lib.rs#L633), [`ui/dr-ui/src/settings_ui.rs:1`](../ui/dr-ui/src/settings_ui.rs#L1) | | FR-EXP-9 | [`core/dr-decode/src/lib.rs:510`](../core/dr-decode/src/lib.rs#L510), [`core/dr-export/src/lib.rs:128`](../core/dr-export/src/lib.rs#L128), [`core/dr-export/src/lib.rs:1`](../core/dr-export/src/lib.rs#L1), [`core/dr-gpu/src/adjust.rs:1106`](../core/dr-gpu/src/adjust.rs#L1106), [`ui/dr-ui/examples/face_native.rs:1`](../ui/dr-ui/examples/face_native.rs#L1), [`ui/dr-ui/src/develop.rs:4672`](../ui/dr-ui/src/develop.rs#L4672), [`ui/dr-ui/src/develop.rs:4710`](../ui/dr-ui/src/develop.rs#L4710), [`ui/dr-ui/src/develop.rs:6878`](../ui/dr-ui/src/develop.rs#L6878), [`ui/dr-ui/src/lib.rs:601`](../ui/dr-ui/src/lib.rs#L601), [`ui/dr-ui/src/library.rs:4052`](../ui/dr-ui/src/library.rs#L4052), [`ui/dr-ui/src/library.rs:4522`](../ui/dr-ui/src/library.rs#L4522), [`ui/dr-ui/src/library.rs:4558`](../ui/dr-ui/src/library.rs#L4558), [`ui/dr-ui/tests/export_ignores_the_viewport.rs:1`](../ui/dr-ui/tests/export_ignores_the_viewport.rs#L1) | +| FR-MRG-1 | [`core/dr-pano/src/align.rs:1`](../core/dr-pano/src/align.rs#L1), [`core/dr-pano/src/lib.rs:1`](../core/dr-pano/src/lib.rs#L1) | +| FR-MRG-10 | [`core/dr-pano/src/lib.rs:1`](../core/dr-pano/src/lib.rs#L1) | | FR-MRG-3 | [`core/dr-decode/examples/linear_dng.rs:1`](../core/dr-decode/examples/linear_dng.rs#L1) | -| FR-MRG-8 | [`core/dr-segment/examples/onnx_probe.rs:1`](../core/dr-segment/examples/onnx_probe.rs#L1) | +| FR-MRG-4 | [`core/dr-pano/src/projection.rs:1`](../core/dr-pano/src/projection.rs#L1) | +| FR-MRG-5 | [`core/dr-pano/src/align.rs:1`](../core/dr-pano/src/align.rs#L1) | +| FR-MRG-8 | [`core/dr-pano/src/xfeat.rs:1`](../core/dr-pano/src/xfeat.rs#L1), [`core/dr-segment/examples/onnx_probe.rs:1`](../core/dr-segment/examples/onnx_probe.rs#L1) | | FR-NC-1 | [`core/dr-sync-nextcloud/src/auth.rs:132`](../core/dr-sync-nextcloud/src/auth.rs#L132), [`core/dr-sync-nextcloud/src/auth.rs:44`](../core/dr-sync-nextcloud/src/auth.rs#L44), [`core/dr-sync-nextcloud/src/provider.rs:1`](../core/dr-sync-nextcloud/src/provider.rs#L1), [`core/dr-sync/src/account.rs:361`](../core/dr-sync/src/account.rs#L361), [`ui/dr-ui/src/launch.rs:277`](../ui/dr-ui/src/launch.rs#L277), [`ui/dr-ui/src/launch.rs:61`](../ui/dr-ui/src/launch.rs#L61), [`ui/dr-ui/src/launch_ui.rs:417`](../ui/dr-ui/src/launch_ui.rs#L417) | | FR-NC-10 | [`core/dr-sync/src/account.rs:226`](../core/dr-sync/src/account.rs#L226), [`ui/dr-ui/src/export.rs:1`](../ui/dr-ui/src/export.rs#L1), [`ui/dr-ui/src/lib.rs:674`](../ui/dr-ui/src/lib.rs#L674), [`ui/dr-ui/src/library.rs:1170`](../ui/dr-ui/src/library.rs#L1170), [`ui/dr-ui/src/library.rs:2463`](../ui/dr-ui/src/library.rs#L2463), [`ui/dr-ui/src/library.rs:589`](../ui/dr-ui/src/library.rs#L589), [`ui/dr-ui/src/library.rs:916`](../ui/dr-ui/src/library.rs#L916), [`ui/dr-ui/src/library_ui.rs:2062`](../ui/dr-ui/src/library_ui.rs#L2062), [`ui/dr-ui/src/library_ui.rs:4275`](../ui/dr-ui/src/library_ui.rs#L4275), [`ui/dr-ui/src/library_ui.rs:660`](../ui/dr-ui/src/library_ui.rs#L660), [`ui/dr-ui/src/sidecar_cache.rs:1`](../ui/dr-ui/src/sidecar_cache.rs#L1) | | FR-NC-12 | [`core/dr-sync-folder/src/lib.rs:1`](../core/dr-sync-folder/src/lib.rs#L1), [`core/dr-sync-nextcloud/src/lib.rs:1022`](../core/dr-sync-nextcloud/src/lib.rs#L1022), [`core/dr-sync-nextcloud/src/lib.rs:40`](../core/dr-sync-nextcloud/src/lib.rs#L40), [`core/dr-sync-nextcloud/src/provider.rs:1`](../core/dr-sync-nextcloud/src/provider.rs#L1), [`core/dr-sync/src/account.rs:1`](../core/dr-sync/src/account.rs#L1), [`core/dr-sync/src/account.rs:87`](../core/dr-sync/src/account.rs#L87), [`core/dr-sync/src/lib.rs:218`](../core/dr-sync/src/lib.rs#L218), [`core/dr-sync/src/lib.rs:51`](../core/dr-sync/src/lib.rs#L51), [`core/dr-sync/src/provider.rs:106`](../core/dr-sync/src/provider.rs#L106), [`core/dr-sync/src/provider.rs:1`](../core/dr-sync/src/provider.rs#L1), [`core/dr-sync/src/provider.rs:53`](../core/dr-sync/src/provider.rs#L53), [`core/dr-sync/src/reachability.rs:1`](../core/dr-sync/src/reachability.rs#L1), [`ui/dr-ui/src/remote.rs:1`](../ui/dr-ui/src/remote.rs#L1) | @@ -210,7 +214,7 @@ Defined in `requirements.md` and marked `(post-v1)` on the defining line. Not in ## Not yet tagged -40 of 184 requirements have no implementation tag. Expected while the codebase is young; each should gain one as it is built. +36 of 184 requirements have no implementation tag. Expected while the codebase is young; each should gain one as it is built.
Show untagged requirements @@ -222,12 +226,8 @@ Defined in `requirements.md` and marked `(post-v1)` on the defining line. Not in - FR-DEV-3g - FR-DSP-2 - FR-DSP-4 -- FR-MRG-1 -- FR-MRG-10 - FR-MRG-11 - FR-MRG-2 -- FR-MRG-4 -- FR-MRG-5 - FR-MRG-6 - FR-MRG-7 - FR-MRG-9