Merge pull request #2341 from vangel76/camtrack-validate-input

Camera tracking: validate solver input and bound workspaces
This commit is contained in:
BFlat
2026-10-29 10:08:05 +07:00
committed by GitHub
4 changed files with 300 additions and 43 deletions
+51 -1
View File
@@ -67,8 +67,9 @@ pub struct Track2D {
}
impl Track2D {
/// Exclusive end frame, saturated when the track exceeds the frame index range.
pub fn end(&self) -> u32 {
self.start + self.pts.len() as u32
u32::try_from(self.pts.len()).ok().and_then(|n| self.start.checked_add(n)).unwrap_or(u32::MAX)
}
/// Position on frame `k` when the track covers it.
pub fn at(&self, k: u32) -> Option<[f64; 2]> {
@@ -483,3 +484,52 @@ pub fn point_color(id: u32) -> [f32; 3] {
// Bright, slightly desaturated.
[0.25 + 0.75 * r, 0.25 + 0.75 * g, 0.25 + 0.75 * b]
}
#[cfg(test)]
mod validation_tests {
use super::*;
fn tracks() -> CameraTracks {
CameraTracks { size: [1920.0, 1080.0], factor: 1.0, frame_duration: 1.0 / 30.0, frames: 3, ..Default::default() }
}
#[test]
fn track_end_saturates_instead_of_wrapping() {
assert_eq!(Track2D { start: u32::MAX, pts: vec![[0.0; 2]], ..Default::default() }.end(), u32::MAX);
assert_eq!(Track2D { start: 3, pts: vec![[0.0; 2]; 2], ..Default::default() }.end(), 5);
}
#[test]
fn solve_rejects_large_frame_count_before_allocating() {
let t = CameraTracks { frames: u32::MAX, ..tracks() };
let error = solve(&t, &SolveSettings::default(), None).unwrap_err();
assert!(error.0.contains("frame or track count"));
}
#[test]
fn solve_rejects_overflowing_and_out_of_clip_tracks() {
for start in [u32::MAX, 2] {
let t = CameraTracks { tracks: vec![Track2D { start, pts: vec![[0.0; 2]; 3], ..Default::default() }], ..tracks() };
let error = solve(&t, &SolveSettings::default(), None).unwrap_err();
assert!(error.0.contains("frame range"));
}
}
#[test]
fn solve_rejects_non_finite_track_positions() {
let t = CameraTracks { tracks: vec![Track2D { pts: vec![[f32::NAN, 0.0]; 3], ..Default::default() }], ..tracks() };
let error = solve(&t, &SolveSettings::default(), None).unwrap_err();
assert!(error.0.contains("finite"));
}
#[test]
fn solve_rejects_invalid_dimensions_timing_and_scale() {
for t in [
CameraTracks { size: [f64::NAN, 1080.0], ..tracks() },
CameraTracks { frame_duration: 0.0, ..tracks() },
CameraTracks { factor: f64::INFINITY, ..tracks() },
] {
assert!(solve(&t, &SolveSettings::default(), None).is_err());
}
}
}
+166 -18
View File
@@ -10,9 +10,9 @@
use rayon::prelude::*;
use super::Distortion;
use super::geometry::Pose;
use super::linalg::{M3, V3, add, cholesky_solve, dot, inv3, mm, mtv, mv, rodrigues, scale, skew, sub};
use super::{Distortion, SolveError};
/// A camera in the adjustment.
#[derive(Clone, Copy, Debug, PartialEq)]
@@ -52,6 +52,7 @@ pub struct BaOpts {
pub fixed: Vec<bool>,
/// Huber threshold (pixels).
pub huber: f64,
/// Maximum LM iterations (at most 1000).
pub max_iter: usize,
/// Adjust the shared radial distortion `(k1, k2)` too.
pub distortion: bool,
@@ -76,21 +77,89 @@ fn huber(e2: f64, k: f64) -> f64 {
if e2 <= k * k { e2 } else { 2.0 * k * e2.sqrt() - k * k }
}
/// Residual (pixels) of one observation, `None` behind the camera.
// Bound the temporary systems before any input-sized allocation. Includes dense matrices
// for the main system, trial copies and both Rayon phases, plus per-point/observation data.
const MAX_BA_BYTES: usize = 256 * 1024 * 1024;
const MAX_BA_CAMERAS: usize = 4096;
const MAX_BA_POINTS: usize = 200_000;
const MAX_BA_OBSERVATIONS: usize = 1_000_000;
fn valid_camera(c: &BaCam) -> bool {
c.f.is_finite()
&& c.f > 0.0
&& c.pose.r.iter().flatten().chain(&c.pose.c).all(|v| v.is_finite())
&& [c.dist.k1, c.dist.k2, c.dist.radius].iter().all(|v| v.is_finite())
}
fn validate(cams: &[BaCam], pts: &[V3], obs: &[BaObs], o: &BaOpts) -> Result<(), SolveError> {
if cams.len() > MAX_BA_CAMERAS || pts.len() > MAX_BA_POINTS || obs.len() > MAX_BA_OBSERVATIONS {
return Err(SolveError("bundle adjustment input exceeds the supported size".into()));
}
if !o.huber.is_finite() || o.huber <= 0.0 || o.max_iter > 1000 || matches!(o.focal, FocalMode::PerCamera(k) if !k.is_finite() || k < 0.0) {
return Err(SolveError("invalid bundle adjustment options".into()));
}
if cams.iter().any(|c| !valid_camera(c)) || pts.iter().flatten().any(|v| !v.is_finite()) {
return Err(SolveError("bundle adjustment cameras and points must be finite, with positive focal lengths".into()));
}
for ob in obs {
if cams.get(ob.cam as usize).is_none() || pts.get(ob.pt as usize).is_none() {
return Err(SolveError("bundle adjustment observation index is out of range".into()));
}
if ob.uv.iter().any(|v| !v.is_finite()) {
return Err(SolveError("bundle adjustment observations must be finite".into()));
}
}
Ok(())
}
fn matrix_len(n: usize, npts: usize, nobs: usize) -> Result<usize, SolveError> {
let threads = rayon::current_num_threads();
let build_workers = if n <= 700 && npts > 64 { threads.clamp(1, 16) } else { 0 };
let schur_workers = if n <= 1200 && npts >= 64 && threads > 1 { threads.clamp(1, 8) } else { 0 };
// The two worker phases do not coexist; allow the larger phase plus main/trial buffers.
let copies = 4 + build_workers.max(schur_workers);
let len = n.checked_mul(n);
let bytes = len
.and_then(|v| v.checked_add(n))
.and_then(|v| v.checked_mul(size_of::<f64>()))
.and_then(|v| v.checked_mul(copies))
.and_then(|v| npts.checked_mul(320).and_then(|p| v.checked_add(p)))
.and_then(|v| nobs.checked_mul(2 * 9 * size_of::<(usize, V3)>() + size_of::<BaObs>() + 2 * size_of::<u32>()).and_then(|p| v.checked_add(p)));
match (len, bytes) {
(Some(len), Some(bytes)) if bytes <= MAX_BA_BYTES => Ok(len),
_ => Err(SolveError("bundle adjustment exceeds its 256 MiB workspace budget; use fewer cameras or observations".into())),
}
}
/// Residual (pixels), or `None` for an invalid or behind-camera projection.
pub fn residual(c: &BaCam, x: V3, uv: [f64; 2]) -> Option<[f64; 2]> {
if !valid_camera(c) || x.iter().chain(&uv).any(|v| !v.is_finite()) {
return None;
}
project_residual(c, x, uv)
}
fn project_residual(c: &BaCam, x: V3, uv: [f64; 2]) -> Option<[f64; 2]> {
let p = c.pose.to_cam(x);
if p[2] <= 1e-9 {
if p[2] <= 1e-9 || p.iter().any(|v| !v.is_finite()) {
return None;
}
let q = c.dist.distort([c.f * p[0] / p[2], c.f * p[1] / p[2]]);
Some([q[0] - uv[0], q[1] - uv[1]])
let r = [q[0] - uv[0], q[1] - uv[1]];
r.iter().all(|v| v.is_finite()).then_some(r)
}
/// Total robust cost.
pub fn cost(cams: &[BaCam], pts: &[V3], obs: &[BaObs], o: &BaOpts) -> f64 {
/// Total robust cost. Rejects invalid inputs before starting parallel evaluation.
pub fn cost(cams: &[BaCam], pts: &[V3], obs: &[BaObs], o: &BaOpts) -> Result<f64, SolveError> {
validate(cams, pts, obs, o)?;
let c = validated_cost(cams, pts, obs, o);
if c.is_finite() { Ok(c) } else { Err(SolveError("bundle adjustment cost is not finite".into())) }
}
fn validated_cost(cams: &[BaCam], pts: &[V3], obs: &[BaObs], o: &BaOpts) -> f64 {
let c: f64 = obs
.par_iter()
.map(|ob| match residual(&cams[ob.cam as usize], pts[ob.pt as usize], ob.uv) {
.map(|ob| match cams.get(ob.cam as usize).zip(pts.get(ob.pt as usize)).and_then(|(c, x)| project_residual(c, *x, ob.uv)) {
Some(r) => huber(r[0] * r[0] + r[1] * r[1], o.huber),
None => BEHIND,
})
@@ -198,15 +267,21 @@ fn add_focal_prior(cams: &[BaCam], o: &BaOpts, l: &Layout, s_mat: &mut [f64], bc
}
/// Run Levenberg–Marquardt. `cams` and `pts` are updated in place.
pub fn adjust(cams: &mut [BaCam], pts: &mut [V3], obs: &[BaObs], o: &BaOpts) -> BaReport {
/// Invalid inputs or workspaces exceeding 256 MiB return an error before mutation.
pub fn adjust(cams: &mut [BaCam], pts: &mut [V3], obs: &[BaObs], o: &BaOpts) -> Result<BaReport, SolveError> {
validate(cams, pts, obs, o)?;
let l = layout(cams.len(), o);
let matrix_len = matrix_len(l.n, pts.len(), obs.len())?;
let npts = pts.len();
// Observations by point.
let mut by_pt: Vec<Vec<u32>> = vec![vec![]; npts];
for (i, ob) in obs.iter().enumerate() {
by_pt[ob.pt as usize].push(i as u32);
}
let mut cur = cost(cams, pts, obs, o);
let mut cur = validated_cost(cams, pts, obs, o);
if !cur.is_finite() {
return Err(SolveError("bundle adjustment cost is not finite".into()));
}
let initial = cur;
let t0 = web_time::Instant::now();
let mut tries = 0usize;
@@ -214,7 +289,7 @@ pub fn adjust(cams: &mut [BaCam], pts: &mut [V3], obs: &[BaObs], o: &BaOpts) ->
let mut lambda = 1e-4;
let mut iters = 0;
if l.n == 0 && o.fix_points {
return BaReport { initial_cost: initial, cost: cur, iterations: 0 };
return Ok(BaReport { initial_cost: initial, cost: cur, iterations: 0 });
}
let n = l.n;
let parallel = n <= 700;
@@ -297,12 +372,12 @@ pub fn adjust(cams: &mut [BaCam], pts: &mut [V3], obs: &[BaObs], o: &BaOpts) ->
let parts: Vec<(Vec<f64>, Vec<f64>, Vec<PtSys>)> = (0..chunks)
.into_par_iter()
.map(|t| {
let (mut s, mut b) = (vec![0.0; n * n], vec![0.0; n]);
let (mut s, mut b) = (vec![0.0; matrix_len], vec![0.0; n]);
let v = (t * per..((t + 1) * per).min(npts)).map(|j| build_point(j, &mut s, &mut b)).collect();
(s, b, v)
})
.collect();
let mut s = vec![0.0; n * n];
let mut s = vec![0.0; matrix_len];
let mut b = vec![0.0; n];
let mut sys = Vec::with_capacity(npts);
for (ps, pb, pv) in parts {
@@ -316,7 +391,7 @@ pub fn adjust(cams: &mut [BaCam], pts: &mut [V3], obs: &[BaObs], o: &BaOpts) ->
}
(s, b, sys)
} else {
let mut s = vec![0.0; n * n];
let mut s = vec![0.0; matrix_len];
let mut b = vec![0.0; n];
let sys = (0..npts).map(|j| build_point(j, &mut s, &mut b)).collect();
(s, b, sys)
@@ -376,7 +451,7 @@ pub fn adjust(cams: &mut [BaCam], pts: &mut [V3], obs: &[BaObs], o: &BaOpts) ->
let parts: Vec<(Vec<f64>, Vec<f64>)> = (0..chunks)
.into_par_iter()
.map(|t| {
let (mut ds, mut dr) = (vec![0.0; n * n], vec![0.0; n]);
let (mut ds, mut dr) = (vec![0.0; matrix_len], vec![0.0; n]);
reduce(t * per..((t + 1) * per).min(npts), &mut ds, &mut dr);
(ds, dr)
})
@@ -430,7 +505,11 @@ pub fn adjust(cams: &mut [BaCam], pts: &mut [V3], obs: &[BaObs], o: &BaOpts) ->
*x = add(*x, mv(vi, b));
});
}
let c2 = cost(&nc, &np, obs, o);
if nc.iter().any(|c| !valid_camera(c)) || np.iter().flatten().any(|v| !v.is_finite()) {
lambda *= 4.0;
continue;
}
let c2 = validated_cost(&nc, &np, obs, o);
if c2.is_finite() && c2 < cur {
cams.copy_from_slice(&nc);
pts.copy_from_slice(&np);
@@ -449,7 +528,7 @@ pub fn adjust(cams: &mut [BaCam], pts: &mut [V3], obs: &[BaObs], o: &BaOpts) ->
t0.elapsed().as_secs_f64()
);
}
return BaReport { initial_cost: initial, cost: cur, iterations: iters };
return Ok(BaReport { initial_cost: initial, cost: cur, iterations: iters });
}
break;
}
@@ -462,7 +541,7 @@ pub fn adjust(cams: &mut [BaCam], pts: &mut [V3], obs: &[BaObs], o: &BaOpts) ->
if dbg {
eprintln!("ba n={} pts={} obs={} iters={iters} tries={tries} {:.2}s cost {initial:.1} -> {cur:.1}", l.n, npts, obs.len(), t0.elapsed().as_secs_f64());
}
BaReport { initial_cost: initial, cost: cur, iterations: iters }
Ok(BaReport { initial_cost: initial, cost: cur, iterations: iters })
}
#[cfg(test)]
@@ -470,6 +549,75 @@ mod tests {
use super::*;
use crate::camtrack::geometry::Rng;
fn test_camera() -> BaCam {
BaCam { pose: Pose::IDENTITY, f: 900.0, dist: Distortion::NONE }
}
#[test]
fn cost_rejects_invalid_observation_indices() {
for (cam, pt) in [(1, 0), (0, 1)] {
let obs = [BaObs { cam, pt, uv: [0.0; 2] }];
assert!(cost(&[test_camera()], &[[0.0, 0.0, 5.0]], &obs, &BaOpts::default()).is_err());
}
}
#[test]
fn adjust_rejects_invalid_observation_without_mutating_inputs() {
for (cam, pt) in [(1, 0), (0, 1)] {
let mut cams = [test_camera()];
let mut pts = [[0.0, 0.0, 5.0]];
let before = (cams, pts);
let obs = [BaObs { cam, pt, uv: [0.0; 2] }];
assert!(adjust(&mut cams, &mut pts, &obs, &BaOpts::default()).is_err());
assert_eq!((cams, pts), before);
}
}
#[test]
fn residual_rejects_non_finite_inputs_and_projection_overflow() {
let cam = test_camera();
assert!(residual(&cam, [0.0, 0.0, f64::NAN], [0.0; 2]).is_none());
assert!(residual(&cam, [0.0, 0.0, 5.0], [f64::INFINITY, 0.0]).is_none());
for bad in [
BaCam { f: f64::NAN, ..cam },
BaCam { dist: Distortion { k1: f64::INFINITY, ..Distortion::NONE }, ..cam },
BaCam { pose: Pose { c: [f64::NAN; 3], ..Pose::IDENTITY }, ..cam },
BaCam { pose: Pose { r: [[f64::INFINITY; 3]; 3], ..Pose::IDENTITY }, ..cam },
] {
assert!(residual(&bad, [0.0, 0.0, 5.0], [0.0; 2]).is_none());
}
assert!(residual(&cam, [f64::MAX, 0.0, 1.0], [0.0; 2]).is_none());
}
#[test]
fn adjustment_rejects_non_finite_points_observations_and_options() {
let cam = test_camera();
let pts = [[0.0, 0.0, 5.0]];
let obs = [BaObs { cam: 0, pt: 0, uv: [0.0; 2] }];
assert!(cost(&[cam], &[[f64::NAN, 0.0, 5.0]], &obs, &BaOpts::default()).is_err());
assert!(cost(&[cam], &pts, &[BaObs { uv: [f64::NAN; 2], ..obs[0] }], &BaOpts::default()).is_err());
for opts in [
BaOpts { huber: f64::NAN, ..Default::default() },
BaOpts { huber: 0.0, ..Default::default() },
BaOpts { focal: FocalMode::PerCamera(f64::INFINITY), ..Default::default() },
BaOpts { max_iter: usize::MAX, ..Default::default() },
] {
let mut trial_pts = pts;
assert!(adjust(&mut [cam], &mut trial_pts, &obs, &opts).is_err());
}
}
#[test]
fn matrix_workspace_is_checked_and_bounded_before_allocation() {
assert!(matrix_len(usize::MAX, 0, 0).is_err());
assert!(matrix_len(1, usize::MAX, 0).is_err());
assert!(matrix_len(1, 0, usize::MAX).is_err());
let mut cams = vec![test_camera(); 1000];
assert!(adjust(&mut cams, &mut [], &[], &BaOpts::default()).is_err());
assert!(matrix_len(6, 1000, 1000).is_ok());
assert!(matrix_len(955, 900, 160 * 900).is_ok());
}
#[test]
fn bundle_adjustment_converges_from_noise() {
let mut rng = Rng(99);
@@ -499,7 +647,7 @@ mod tests {
let mut fixed = vec![false; 6];
fixed[0] = true;
fixed[1] = true;
let rep = adjust(&mut cams, &mut pts, &obs, &BaOpts { fixed, max_iter: 50, ..Default::default() });
let rep = adjust(&mut cams, &mut pts, &obs, &BaOpts { fixed, max_iter: 50, ..Default::default() }).unwrap();
let rms = (rep.cost / obs.len() as f64).sqrt();
assert!(rms < 1e-4, "{rep:?}");
}
+11
View File
@@ -247,6 +247,9 @@ pub fn reproj2(pose: &Pose, x: V3, o: [f64; 2]) -> Option<f64> {
/// Pick the essential decomposition with the most points in front of both cameras; returns
/// the second camera (the first is the identity) and that count.
pub fn select_pose(cands: &[(M3, V3)], x1: &[[f64; 2]], x2: &[[f64; 2]]) -> Option<(Pose, usize)> {
if x1.len() != x2.len() || x1.is_empty() || x1.iter().chain(x2).flatten().any(|v| !v.is_finite()) {
return None;
}
let mut best: Option<(Pose, usize)> = None;
for (r, t) in cands {
let p2 = Pose::from_rt(*r, *t);
@@ -384,6 +387,14 @@ mod tests {
[c[0] / c[2], c[1] / c[2]]
}
#[test]
fn pose_selection_rejects_mismatched_or_non_finite_correspondences() {
let candidates = [(I3, [0.1, 0.0, 0.0])];
assert!(select_pose(&candidates, &[[0.0; 2]], &[]).is_none());
assert!(select_pose(&candidates, &[], &[[0.0; 2]]).is_none());
assert!(select_pose(&candidates, &[[f64::NAN, 0.0]], &[[0.0; 2]]).is_none());
}
#[test]
fn essential_recovers_relative_pose() {
let (pts, p2) = scene();
+72 -24
View File
@@ -42,13 +42,54 @@ struct Data {
noise_bits: std::sync::atomic::AtomicU64,
}
// Input budgets apply before building the duplicated per-frame/per-track lookup tables.
const MAX_SOLVE_FRAMES: usize = 100_000;
const MAX_SOLVE_TRACKS: usize = 200_000;
const MAX_SOLVE_OBSERVATIONS: usize = 2_000_000;
fn validate_tracks(t: &CameraTracks) -> Result<(), SolveError> {
if t.frames as usize > MAX_SOLVE_FRAMES || t.tracks.len() > MAX_SOLVE_TRACKS {
return err("camera tracks exceed the supported frame or track count");
}
if t.size.iter().any(|v| !v.is_finite() || *v <= 0.0)
|| !t.factor.is_finite()
|| t.factor < 0.0
|| !t.start.is_finite()
|| !t.frame_duration.is_finite()
|| t.frame_duration <= 0.0
{
return err("camera tracks have invalid dimensions, timing or analysis scale");
}
let mut count = 0usize;
for tr in &t.tracks {
count = count
.checked_add(tr.pts.len())
.filter(|n| *n <= MAX_SOLVE_OBSERVATIONS)
.ok_or_else(|| SolveError("camera tracks exceed the supported observation count".into()))?;
let end = u32::try_from(tr.pts.len()).ok().and_then(|n| tr.start.checked_add(n));
if end.is_none_or(|end| end > t.frames) {
return err("camera track extends beyond the clip's frame range");
}
if tr.pts.iter().flatten().any(|v| !v.is_finite()) {
return err("camera track positions must be finite");
}
}
Ok(())
}
impl Data {
fn new(t: &CameraTracks, deleted: &[u32]) -> Data {
fn new(t: &CameraTracks, deleted: &[u32]) -> Result<Data, SolveError> {
validate_tracks(t)?;
if deleted.len() > MAX_SOLVE_TRACKS {
return err("too many deleted camera track ids");
}
let n = t.frames as usize;
let c = [t.size[0] * 0.5, t.size[1] * 0.5];
let del: HashSet<u32> = deleted.iter().copied().collect();
let mut tracks = vec![];
let mut by_frame = vec![vec![]; n];
let mut by_frame = Vec::new();
by_frame.try_reserve_exact(n).map_err(|_| SolveError("not enough memory for camera track frames".into()))?;
by_frame.resize_with(n, Vec::new);
for tr in t.tracks.iter().filter(|tr| !del.contains(&tr.id) && tr.pts.len() >= 3) {
let ti = tracks.len() as u32;
let pts: Vec<[f64; 2]> = tr.pts.iter().map(|p| [p[0] as f64 - c[0], p[1] as f64 - c[1]]).collect();
@@ -60,7 +101,7 @@ impl Data {
tracks.push((tr.id, tr.start as usize, pts));
}
let noise = (0.5 * t.factor).max(0.5);
Data { n, w: t.size[0], tracks, by_frame, noise_bits: std::sync::atomic::AtomicU64::new(noise.to_bits()) }
Ok(Data { n, w: t.size[0], tracks, by_frame, noise_bits: std::sync::atomic::AtomicU64::new(noise.to_bits()) })
}
fn at(&self, ti: u32, k: usize) -> Option<[f64; 2]> {
let (_, s, p) = &self.tracks[ti as usize];
@@ -174,10 +215,10 @@ impl Recon {
}
/// Bundle-adjust the cameras of `frames` (the first stays fixed) and their points.
fn adjust(&mut self, d: &Data, frames: &[usize], cfg: &Cfg, focal: FocalMode, iters: usize) -> f64 {
fn adjust(&mut self, d: &Data, frames: &[usize], cfg: &Cfg, focal: FocalMode, iters: usize) -> Result<f64, SolveError> {
let (mut cams, mut pts, obs, ids) = self.observations(d, frames);
if obs.is_empty() {
return f64::INFINITY;
return Ok(f64::INFINITY);
}
let mut fixed = vec![false; cams.len()];
fixed[0] = true;
@@ -186,14 +227,14 @@ impl Recon {
&mut pts,
&obs,
&BaOpts { focal, fix_points: false, fix_centers: cfg.kind == Kind::Tripod, fixed, huber: d.huber(), max_iter: iters, distortion: cfg.distortion },
);
)?;
for (k, c) in frames.iter().zip(&cams) {
self.cams.insert(*k, *c);
}
for (ti, x) in ids.iter().zip(&pts) {
self.pts.insert(*ti, *x);
}
rep.cost
Ok(rep.cost)
}
/// Mark observations above the inlier threshold as outliers; drop points left with fewer
@@ -302,8 +343,12 @@ impl Recon {
break;
}
let mut cams = vec![start];
bundle::adjust(&mut cams, &mut pts.clone(), &obs, &opts(d.inlier_thr() * 4.0));
bundle::adjust(&mut cams, &mut pts.clone(), &obs, &opts(d.huber()));
if bundle::adjust(&mut cams, &mut pts.clone(), &obs, &opts(d.inlier_thr() * 4.0)).is_err() {
return false;
}
if bundle::adjust(&mut cams, &mut pts.clone(), &obs, &opts(d.huber())).is_err() {
return false;
}
let inl = obs.iter().filter(|o| bundle::residual(&cams[0], pts[o.pt as usize], o.uv).is_some_and(|r| r[0] * r[0] + r[1] * r[1] <= thr2)).count();
if best.as_ref().is_none_or(|b| inl > b.0) {
best = Some((inl, cams[0]));
@@ -317,7 +362,9 @@ impl Recon {
let keep: Vec<BaObs> =
obs.iter().filter(|o| bundle::residual(&cam, pts[o.pt as usize], o.uv).is_some_and(|r| r[0] * r[0] + r[1] * r[1] <= thr2)).copied().collect();
let mut cams = vec![cam];
bundle::adjust(&mut cams, &mut pts.clone(), &keep, &opts(d.huber()));
if bundle::adjust(&mut cams, &mut pts.clone(), &keep, &opts(d.huber())).is_err() {
return false;
}
for o in &obs {
let ok = bundle::residual(&cams[0], pts[o.pt as usize], o.uv).is_some_and(|r| r[0] * r[0] + r[1] * r[1] <= thr2);
if !ok {
@@ -471,10 +518,11 @@ fn grow(d: &Data, kfs: &[usize], j: usize, init: &(Pose, HashMap<u32, V3>), f: f
r.pts = init.1.clone();
let inner = |r: &mut Recon, frames: &[usize]| {
let focal = if cfg.per_frame_f && frames.len() > 2 { FocalMode::PerCamera(F_PRIOR) } else { FocalMode::Fixed };
r.adjust(d, frames, cfg, focal, 10);
r.adjust(d, frames, cfg, focal, 10)?;
r.reject(d);
Ok::<(), SolveError>(())
};
inner(&mut r, &[a, b]);
inner(&mut r, &[a, b]).ok()?;
let order: Vec<usize> = kfs[1..j].iter().chain(&kfs[j + 1..]).copied().collect();
let mut done = 0;
for k in order {
@@ -491,11 +539,11 @@ fn grow(d: &Data, kfs: &[usize], j: usize, init: &(Pose, HashMap<u32, V3>), f: f
done += 1;
let frames: Vec<usize> = r.cams.keys().copied().collect();
if done % 3 == 0 || frames.len() <= 6 {
inner(&mut r, &frames);
inner(&mut r, &frames).ok()?;
}
}
let frames: Vec<usize> = r.cams.keys().copied().collect();
inner(&mut r, &frames);
inner(&mut r, &frames).ok()?;
(r.cams.len() >= 2 && r.pts.len() >= 8).then_some(r)
}
@@ -531,10 +579,10 @@ fn grow_tripod(d: &Data, kfs: &[usize], f: f64, cfg: &Cfg, cancel: Option<&Atomi
r.triangulate(d, Kind::Tripod, 0.0, 20.0);
let frames: Vec<usize> = r.cams.keys().copied().collect();
let focal = if cfg.per_frame_f && frames.len() > 2 { FocalMode::PerCamera(F_PRIOR) } else { FocalMode::Fixed };
r.adjust(d, &frames, cfg, focal, 15);
r.adjust(d, &frames, cfg, focal, 15).ok()?;
r.reject(d);
r.triangulate(d, Kind::Tripod, 0.0, 1.0);
r.adjust(d, &frames, cfg, focal, 15);
r.adjust(d, &frames, cfg, focal, 15).ok()?;
(r.cams.len() >= 2 && r.pts.len() >= 8).then_some(r)
}
@@ -684,7 +732,7 @@ pub fn solve(tracks: &CameraTracks, s: &SolveSettings, cancel: Option<&AtomicBoo
for pass in 0..2 {
let und = (pass > 0).then(|| undistorted(tracks, &dist));
let c = solve_core(und.as_ref().unwrap_or(tracks), s, cancel)?;
let d = Data::new(tracks, &s.deleted);
let d = Data::new(tracks, &s.deleted)?;
d.set_noise(c.d.noise());
let mut r = c.r;
for cam in r.cams.values_mut() {
@@ -699,10 +747,10 @@ pub fn solve(tracks: &CameraTracks, s: &SolveSettings, cancel: Option<&AtomicBoo
// The pinhole solve's outliers include the distorted frame edges: start from all
// observations (the Huber loss keeps real outliers in check), then reject again.
r.outliers.clear();
r.adjust(&d, &frames, &cfg, c.focal, 30);
r.adjust(&d, &frames, &cfg, c.focal, 30)?;
r.outliers.clear();
r.reject(&d);
r.adjust(&d, &frames, &cfg, c.focal, 30);
r.adjust(&d, &frames, &cfg, c.focal, 30)?;
dist = r.cams[&frames[0]].dist;
for cam in r.cams.values_mut() {
cam.dist = dist;
@@ -742,7 +790,7 @@ pub fn solve(tracks: &CameraTracks, s: &SolveSettings, cancel: Option<&AtomicBoo
}
fn solve_core(tracks: &CameraTracks, s: &SolveSettings, cancel: Option<&AtomicBool>) -> Result<Core, SolveError> {
let d = Data::new(tracks, &s.deleted);
let d = Data::new(tracks, &s.deleted)?;
if d.n < 2 {
return err("the clip must be at least two frames long");
}
@@ -848,7 +896,7 @@ fn solve_core(tracks: &CameraTracks, s: &SolveSettings, cancel: Option<&AtomicBo
lap("tripod");
// Keyframes: adjust with the focal length free, reject, re-triangulate, adjust again.
let kf_frames: Vec<usize> = r.cams.keys().copied().collect();
r.adjust(&d, &kf_frames, &cfg, focal, 40);
r.adjust(&d, &kf_frames, &cfg, focal, 40)?;
r.reject(&d);
// Re-estimate the tracking noise from the inlier residuals (robustly, from their median:
// for 2D Gaussian noise of deviation σ the median distance is 1.1774 σ), then tighten the
@@ -876,9 +924,9 @@ fn solve_core(tracks: &CameraTracks, s: &SolveSettings, cancel: Option<&AtomicBo
r.reject(&d);
}
r.triangulate(&d, kind, 1.0, 1.0);
r.adjust(&d, &kf_frames, &cfg, focal, 40);
r.adjust(&d, &kf_frames, &cfg, focal, 40)?;
r.reject(&d);
r.adjust(&d, &kf_frames, &cfg, focal, 20);
r.adjust(&d, &kf_frames, &cfg, focal, 20)?;
if cancelled() {
return err("cancelled");
}
@@ -910,7 +958,7 @@ fn solve_core(tracks: &CameraTracks, s: &SolveSettings, cancel: Option<&AtomicBo
let all: Vec<usize> = r.cams.keys().copied().collect();
r.reject(&d);
if all.len() <= 160 {
r.adjust(&d, &all, &cfg, focal, 25);
r.adjust(&d, &all, &cfg, focal, 25)?;
r.reject(&d);
} else {
// Points only (cameras fixed) would be cheap but the keyframe solution is already