diff --git a/crates/track/src/camtrack/bundle.rs b/crates/track/src/camtrack/bundle.rs index 077f04b71..7f8dc5bcc 100644 --- a/crates/track/src/camtrack/bundle.rs +++ b/crates/track/src/camtrack/bundle.rs @@ -12,7 +12,7 @@ use rayon::prelude::*; use super::Distortion; use super::geometry::Pose; -use super::linalg::*; +use super::linalg::{M3, V3, add, cholesky_solve, dot, inv3, mm, mtv, mv, rodrigues, scale, skew, sub}; /// A camera in the adjustment. #[derive(Clone, Copy, Debug, PartialEq)] @@ -174,6 +174,29 @@ fn cam_jac(l: &Layout, ci: usize, c: &BaCam, xc: V3, a: &[[f64; 3]; 2], p: [f64; } } +/// Smooth per-camera log focal lengths in the normal equations. +fn add_focal_prior(cams: &[BaCam], o: &BaOpts, l: &Layout, s_mat: &mut [f64], bc: &mut [f64]) { + let n = l.n; + // Focal prior. + if let FocalMode::PerCamera(k) = o.focal { + for i in 0..cams.len().saturating_sub(1) { + let r = k * (cams[i].f.ln() - cams[i + 1].f.ln()); + let (a, b) = (l.fidx[i], l.fidx[i + 1]); + let ia = (a != usize::MAX).then_some(a); + let ib = (b != usize::MAX).then_some(b); + for (idx, g) in [(ia, k), (ib, -k)] { + let Some(x) = idx else { continue }; + bc[x] -= g * r; + for (idy, h) in [(ia, k), (ib, -k)] { + if let Some(y) = idy { + s_mat[x * n + y] += g * h; + } + } + } + } + } +} + /// Run Levenberg–Marquardt. `cams` and `pts` are updated in place. pub fn adjust(cams: &mut [BaCam], pts: &mut [V3], obs: &[BaObs], o: &BaOpts) -> BaReport { let l = layout(cams.len(), o); @@ -298,24 +321,7 @@ pub fn adjust(cams: &mut [BaCam], pts: &mut [V3], obs: &[BaObs], o: &BaOpts) -> let sys = (0..npts).map(|j| build_point(j, &mut s, &mut b)).collect(); (s, b, sys) }; - // Focal prior. - if let FocalMode::PerCamera(k) = o.focal { - for i in 0..cams.len().saturating_sub(1) { - let r = k * (cams[i].f.ln() - cams[i + 1].f.ln()); - let (a, b) = (l.fidx[i], l.fidx[i + 1]); - let ia = (a != usize::MAX).then_some(a); - let ib = (b != usize::MAX).then_some(b); - for (idx, g) in [(ia, k), (ib, -k)] { - let Some(x) = idx else { continue }; - bc[x] -= g * r; - for (idy, h) in [(ia, k), (ib, -k)] { - if let Some(y) = idy { - s_mat[x * n + y] += g * h; - } - } - } - } - } + add_focal_prior(cams, o, &l, &mut s_mat, &mut bc); let base_s = s_mat.clone(); let base_bc = bc.clone(); let mean_diag = ((0..n).map(|i| base_s[i * n + i]).sum::() / n.max(1) as f64).max(1e-12); diff --git a/crates/track/src/camtrack/geometry.rs b/crates/track/src/camtrack/geometry.rs index 6956ac2b5..5d74ad9b3 100644 --- a/crates/track/src/camtrack/geometry.rs +++ b/crates/track/src/camtrack/geometry.rs @@ -5,7 +5,10 @@ //! All image points here are **normalised** (calibrated) coordinates `x = (u − c) / f`. //! A camera is `(R, C)` with `X_cam = R (X − C)`. -use super::linalg::*; +use super::linalg::{ + I3, M3, V3, add, col, cross, det3, dot, from_cols, madd, mm, mscale, mt, mtv, mv, nearest_rotation, normalize, null_vector, outer, scale, sub, svd3, + sym_eigen, +}; /// A deterministic xorshift generator for RANSAC samples. pub struct Rng(pub u64); @@ -361,6 +364,7 @@ pub fn procrustes(a: &[V3], b: &[V3]) -> Option { #[cfg(test)] mod tests { + use super::super::linalg::{log_rot, norm, rodrigues}; use super::*; fn scene() -> (Vec, Pose) { diff --git a/crates/track/src/camtrack/solver.rs b/crates/track/src/camtrack/solver.rs index a4e2c21e6..d00f7106f 100644 --- a/crates/track/src/camtrack/solver.rs +++ b/crates/track/src/camtrack/solver.rs @@ -7,7 +7,7 @@ use rayon::prelude::*; use super::bundle::{self, BaCam, BaObs, BaOpts, FocalMode}; use super::geometry::{self, Pose}; -use super::linalg::*; +use super::linalg::{V3, add, log_rot, mm, mt, mtv, mv, norm, normalize, rodrigues, scale, slerp, sub, sym_eigen}; use super::{CameraSolve, CameraTracks, Distortion, ShotType, SolveMethod, SolveSettings, SolvedFrame, SolvedPoint}; use crate::fit::{Model, ransac}; diff --git a/crates/track/src/face.rs b/crates/track/src/face.rs index 45c7c4baa..b65dc3d81 100644 --- a/crates/track/src/face.rs +++ b/crates/track/src/face.rs @@ -547,6 +547,32 @@ fn fit_from_face(face: &Face, topo: &Topology, offset: [f64; 2]) -> Option Option<(usize, usize)> { + let (w, h) = size; + let (ci, cj) = center; + for r in 0..radius { + for dj in -r..=r { + for di in -r..=r { + if di.abs() != r && dj.abs() != r { + continue; + } + let (i, j) = (ci as i64 + di, cj as i64 + dj); + if i >= 0 && j >= 0 && (i as usize) < w && (j as usize) < h && skin[j as usize * w + i as usize] { + return Some((i as usize, j as usize)); + } + } + } + } + None +} + +struct FaceGeometry { + center: [f64; 2], + height: f64, + roll: f64, +} + /// A pixel window of a frame in layer coordinates. struct Window { x0: i64, @@ -752,6 +778,59 @@ impl FaceTracker { Some(fit) } + /// Connected dark features inside the fitted face, in normalized face coordinates. + fn feature_blobs(&self, frame: &Frame, win: &Window, skin: &[bool], luma: &[f32], outline: &[[f64; 2]], geometry: FaceGeometry) -> Vec { + let FaceGeometry { center, height: b, roll } = geometry; + let (w, h) = (win.w, win.h); + // Non-skin components inside the face outline (shrunk a little: anti-aliased edges). + let inner: Vec<[f64; 2]> = outline.iter().map(|p| [center[0] + 0.9 * (p[0] - center[0]), center[1] + 0.9 * (p[1] - center[1])]).collect(); + let mut feat = vec![false; w * h]; + for j in 0..h { + for i in 0..w { + if !skin[j * w + i] && crate::mask::inside_polygon(&inner, win.layer(frame, i, j)) { + feat[j * w + i] = true; + } + } + } + let min_area = (0.0006 * b * b).max(2.0) as usize; + let mut seen = vec![false; w * h]; + let mut blobs: Vec = vec![]; + let skin_luma = self.skin.med[0] as f64; + for j0 in 0..h { + for i0 in 0..w { + if !feat[j0 * w + i0] || seen[j0 * w + i0] { + continue; + } + let mut st = vec![(i0, j0)]; + seen[j0 * w + i0] = true; + let mut pts = vec![]; + while let Some((i, j)) = st.pop() { + let q = to_face(center, b, roll, win.layer(frame, i, j)); + pts.push((q, skin_luma - luma[j * w + i] as f64)); + for dj in -1i64..=1 { + for di in -1i64..=1 { + let (a, c) = (i as i64 + di, j as i64 + dj); + if a < 0 || c < 0 || a as usize >= w || c as usize >= h { + continue; + } + let (a, c) = (a as usize, c as usize); + if feat[c * w + a] && !seen[c * w + a] { + seen[c * w + a] = true; + st.push((a, c)); + } + } + } + } + if pts.len() >= min_area { + let n = pts.len() as f64; + let centroid = [pts.iter().map(|p| p.0[0]).sum::() / n, pts.iter().map(|p| p.0[1]).sum::() / n]; + blobs.push(Blob { pts, centroid }); + } + } + } + blobs + } + fn fit(&self, frame: &Frame) -> Option { let pv = &self.prev; let b0 = pv.height; @@ -769,22 +848,7 @@ impl FaceTracker { } // Seed: the skin pixel nearest the predicted centre. let (ci, cj) = win.cell(frame, pv.center)?; - let mut seed = None; - 'search: for r in 0..(0.6 * b0).ceil() as i64 { - for dj in -r..=r { - for di in -r..=r { - if di.abs() != r && dj.abs() != r { - continue; - } - let (i, j) = (ci as i64 + di, cj as i64 + dj); - if i >= 0 && j >= 0 && (i as usize) < w && (j as usize) < h && skin[j as usize * w + i as usize] { - seed = Some((i as usize, j as usize)); - break 'search; - } - } - } - } - let (si, sj) = seed?; + let (si, sj) = nearest_skin(&skin, (w, h), (ci, cj), (0.6 * b0).ceil() as i64)?; // Flood fill (4-connected). let mut comp = vec![false; w * h]; let mut stack = vec![(si, sj)]; @@ -855,52 +919,7 @@ impl FaceTracker { roll += 180.0; } let b = major; - // Non-skin components inside the face outline (shrunk a little: anti-aliased edges). - let inner: Vec<[f64; 2]> = outline.iter().map(|p| [center[0] + 0.9 * (p[0] - center[0]), center[1] + 0.9 * (p[1] - center[1])]).collect(); - let mut feat = vec![false; w * h]; - for j in 0..h { - for i in 0..w { - if !skin[j * w + i] && crate::mask::inside_polygon(&inner, win.layer(frame, i, j)) { - feat[j * w + i] = true; - } - } - } - let min_area = (0.0006 * b * b).max(2.0) as usize; - let mut seen = vec![false; w * h]; - let mut blobs: Vec = vec![]; - let skin_luma = self.skin.med[0] as f64; - for j0 in 0..h { - for i0 in 0..w { - if !feat[j0 * w + i0] || seen[j0 * w + i0] { - continue; - } - let mut st = vec![(i0, j0)]; - seen[j0 * w + i0] = true; - let mut pts = vec![]; - while let Some((i, j)) = st.pop() { - let q = to_face(center, b, roll, win.layer(frame, i, j)); - pts.push((q, skin_luma - luma[j * w + i] as f64)); - for dj in -1i64..=1 { - for di in -1i64..=1 { - let (a, c) = (i as i64 + di, j as i64 + dj); - if a < 0 || c < 0 || a as usize >= w || c as usize >= h { - continue; - } - let (a, c) = (a as usize, c as usize); - if feat[c * w + a] && !seen[c * w + a] { - seen[c * w + a] = true; - st.push((a, c)); - } - } - } - } - if pts.len() >= min_area { - let n = pts.len() as f64; - let centroid = [pts.iter().map(|p| p.0[0]).sum::() / n, pts.iter().map(|p| p.0[1]).sum::() / n]; - blobs.push(Blob { pts, centroid }); - } - } - } + let blobs = self.feature_blobs(frame, &win, &skin, &luma, &outline, FaceGeometry { center, height: b, roll }); // Match components to features near their expected positions. let ex = &self.expect; let mean_of = |ids: &[usize]| { diff --git a/crates/track/src/roto.rs b/crates/track/src/roto.rs index 305065d2c..0948d3fe0 100644 --- a/crates/track/src/roto.rs +++ b/crates/track/src/roto.rs @@ -407,7 +407,9 @@ fn beta(img: &[[f32; 3]], free: &[bool], w: usize, h: usize) -> f64 { /// One min-cut over the `free` pixels; others keep `labels`. `prior`: P(foreground) per pixel /// and how much it weighs against the colour models. -fn cut(img: &[[f32; 3]], w: usize, h: usize, labels: &mut [u8], free: &[bool], hard: &[u8], fg: &Gmm, bg: &Gmm, prior: Option<(&[f32], f64)>) { +fn cut(img: &[[f32; 3]], size: (usize, usize), labels: &mut [u8], free: &[bool], hard: &[u8], models: (&Gmm, &Gmm), prior: Option<(&[f32], f64)>) { + let (w, h) = size; + let (fg, bg) = models; let mut idx = vec![u32::MAX; w * h]; let mut nodes = vec![]; for i in 0..w * h { @@ -622,6 +624,36 @@ pub struct Prior { pub weight: f64, } +/// Fit foreground/background colour models from strokes and the propagated matte. +fn initial_models(rgb: &[[f32; 3]], hard: &[u8], fgs: &[bool], strokes: &[&Stroke], prior: Option<&Prior>, opts: &SegOpts, size: (usize, usize)) -> (Gmm, Gmm) { + let (w, h) = size; + // Initial colour models: the strokes (and the propagated matte). + match prior { + Some(p) => { + let fgc = samples(rgb, |i| (hard[i] == HARD_FG) || (hard[i] == 0 && p.soft[i] >= 0.5)); + let bgc = samples(rgb, |i| (hard[i] == HARD_BG) || (hard[i] == 0 && p.soft[i] < 0.5)); + let f = Gmm::fit(&fgc, COMPONENTS, MAX_SAMPLES); + let b = Gmm::fit(&bgc, COMPONENTS, MAX_SAMPLES); + (if f.is_empty() { p.fg.clone() } else { f }, if b.is_empty() { p.bg.clone() } else { b }) + } + None => { + let fgc = samples(rgb, |i| hard[i] == HARD_FG); + let mut bgc = samples(rgb, |i| hard[i] == HARD_BG); + if bgc.is_empty() { + // No background strokes: what is far from the foreground strokes. + let d = matting::distance(fgs, w, h); + let r = strokes.iter().filter(|s| s.kind == StrokeKind::Fg).map(|s| s.radius).fold(1.0, f64::max) * opts.scale; + let far = (4.0 * r).max(0.25 * w.max(h) as f64) as f32; + bgc = samples(rgb, |i| d[i] > far); + if bgc.is_empty() { + bgc = samples(rgb, |i| hard[i] != HARD_FG); + } + } + (Gmm::fit(&fgc, COMPONENTS, MAX_SAMPLES), Gmm::fit(&bgc, COMPONENTS, MAX_SAMPLES)) + } + } +} + /// Segment a frame from its strokes (and a propagated prior): GrabCut iterations on a reduced /// copy, then a full-resolution cut in a band around the boundary. pub fn segment(img: &Image, strokes: &[&Stroke], prior: Option<&Prior>, opts: &SegOpts, refine_radius: f64) -> FrameSeg { @@ -649,31 +681,7 @@ pub fn segment(img: &Image, strokes: &[&Stroke], prior: Option<&Prior>, opts: &S let soft_prior: Option> = prior.map(|p| if p.weight > PRIOR { p.soft.clone() } else { matting::gauss(&p.soft, w, h, (opts.search_radius * opts.scale / 4.0).max(1.0)) }); let pw = prior.map_or(PRIOR, |p| p.weight); - // Initial colour models: the strokes (and the propagated matte). - let (fg0, bg0) = match prior { - Some(p) => { - let fgc = samples(&rgb, |i| (hard[i] == HARD_FG) || (hard[i] == 0 && p.soft[i] >= 0.5)); - let bgc = samples(&rgb, |i| (hard[i] == HARD_BG) || (hard[i] == 0 && p.soft[i] < 0.5)); - let f = Gmm::fit(&fgc, COMPONENTS, MAX_SAMPLES); - let b = Gmm::fit(&bgc, COMPONENTS, MAX_SAMPLES); - (if f.is_empty() { p.fg.clone() } else { f }, if b.is_empty() { p.bg.clone() } else { b }) - } - None => { - let fgc = samples(&rgb, |i| hard[i] == HARD_FG); - let mut bgc = samples(&rgb, |i| hard[i] == HARD_BG); - if bgc.is_empty() { - // No background strokes: what is far from the foreground strokes. - let d = matting::distance(&fgs, w, h); - let r = strokes.iter().filter(|s| s.kind == StrokeKind::Fg).map(|s| s.radius).fold(1.0, f64::max) * opts.scale; - let far = (4.0 * r).max(0.25 * w.max(h) as f64) as f32; - bgc = samples(&rgb, |i| d[i] > far); - if bgc.is_empty() { - bgc = samples(&rgb, |i| hard[i] != HARD_FG); - } - } - (Gmm::fit(&fgc, COMPONENTS, MAX_SAMPLES), Gmm::fit(&bgc, COMPONENTS, MAX_SAMPLES)) - } - }; + let (fg0, bg0) = initial_models(&rgb, &hard, &fgs, strokes, prior, opts, (w, h)); let iters = if opts.best { 4 } else { 3 }; let work = if opts.best { 960 } else { 480 }; let f = w.max(h).div_ceil(work).max(1); @@ -721,7 +729,7 @@ pub fn segment(img: &Image, strokes: &[&Stroke], prior: Option<&Prior>, opts: &S }; let free = vec![true; cw * ch]; for it in 0..iters { - cut(&crgb, cw, ch, &mut cl, &free, &chard, &fgm, &bgm, cprior.as_deref().map(|c| (c, pw))); + cut(&crgb, (cw, ch), &mut cl, &free, &chard, (&fgm, &bgm), cprior.as_deref().map(|c| (c, pw))); if it + 1 < iters { let (a, b) = fit_models(&crgb, &cl); if !a.is_empty() && !b.is_empty() { @@ -743,11 +751,11 @@ pub fn segment(img: &Image, strokes: &[&Stroke], prior: Option<&Prior>, opts: &S let sd = matting::signed_distance(&labels, w, h); let band = f as f32 + 1.5; let free: Vec = sd.iter().map(|d| d.abs() <= band).collect(); - cut(&rgb, w, h, &mut labels, &free, &hard, &fgm, &bgm, soft_prior.as_deref().map(|c| (c, pw))); + cut(&rgb, (w, h), &mut labels, &free, &hard, (&fgm, &bgm), soft_prior.as_deref().map(|c| (c, pw))); } else { let free = vec![true; w * h]; for it in 0..iters { - cut(&rgb, w, h, &mut labels, &free, &hard, &fgm, &bgm, soft_prior.as_deref().map(|c| (c, pw))); + cut(&rgb, (w, h), &mut labels, &free, &hard, (&fgm, &bgm), soft_prior.as_deref().map(|c| (c, pw))); if it + 1 < iters { let (a, b) = fit_models(&rgb, &labels); if !a.is_empty() && !b.is_empty() { @@ -843,7 +851,7 @@ pub fn propagate(prev_img: &Image, prev: &FrameSeg, next_img: &Image, strokes: & let hard = vec![0u8; w * h]; let soft = matting::gauss(&prior.soft, w, h, (r as f64 / 4.0).max(1.0)); let (fg, bg) = if prev.fg.is_empty() || prev.bg.is_empty() { fit_models(&rgb, &labels) } else { (prev.fg.clone(), prev.bg.clone()) }; - cut(&rgb, w, h, &mut labels, &free, &hard, &fg, &bg, Some((&soft, PRIOR))); + cut(&rgb, (w, h), &mut labels, &free, &hard, (&fg, &bg), Some((&soft, PRIOR))); let seeds = core_of(&prior.soft, w, h, r as f64); clean(&mut labels, &seeds, &vec![false; w * h], w, h); finish(next_img, &rgb, labels, strokes, Some(&prior.refine), opts, refine_radius) diff --git a/crates/track/src/roto/matting.rs b/crates/track/src/roto/matting.rs index 439d9c43e..cb22b680a 100644 --- a/crates/track/src/roto/matting.rs +++ b/crates/track/src/roto/matting.rs @@ -200,7 +200,8 @@ pub fn guided_filter(img: &[[f32; 3]], p: &[f32], w: usize, h: usize, r: usize, /// (K. He, J. Sun and X. Tang, "Fast Matting Using Large Kernel Matting Laplacian Matrices", /// CVPR 2010): `(L·x)ᵢ = Σ_{k∋i} (xᵢ − a_kᵀIᵢ − b_k)` with `(a_k, b_k)` the ridge-regression /// fit of `x` to the colours in window k. -pub fn closed_form(img: &[[f32; 3]], init: &[f32], known: &[bool], w: usize, h: usize, r: usize, eps: f32, iters: usize) -> Vec { +pub fn closed_form(img: &[[f32; 3]], init: &[f32], known: &[bool], size: (usize, usize), r: usize, eps: f32, iters: usize) -> Vec { + let (w, h) = size; let n = w * h; if n == 0 || known.iter().all(|k| *k) { return init.to_vec(); diff --git a/crates/track/src/roto/refine.rs b/crates/track/src/roto/refine.rs index 2a498e286..71f9261b2 100644 --- a/crates/track/src/roto/refine.rs +++ b/crates/track/src/roto/refine.rs @@ -123,7 +123,7 @@ pub fn matte_alpha( q[j] = if band[j] { q[j].clamp(0.0, 1.0) } else { cp[j] }; } if p.higher_quality { - q = matting::closed_form(&cimg, &q, &outside, cw, ch, 1, 1e-5 * (1.0 + p.smooth.max(0.0)) as f32, 50); + q = matting::closed_form(&cimg, &q, &outside, (cw, ch), 1, 1e-5 * (1.0 + p.smooth.max(0.0)) as f32, 50); } // Refine Edge Matte adjustments. if p.smooth > 1.0 { @@ -177,7 +177,8 @@ pub fn matte_alpha( /// Decontaminate edge colours of `rgb` in place for the matte `alpha` (see /// [`matting::decontaminate`]); returns the decontamination map. `band`: the Refine Edge band /// (pixels where decontamination always applies). -pub fn decontaminate(rgb: &mut [[f32; 3]], alpha: &[f32], band: &[u8], w: usize, h: usize, p: &MatteParams, band_radius: f64, scale: f64) -> Vec { +pub fn decontaminate(rgb: &mut [[f32; 3]], alpha: &[f32], band: &[u8], size: (usize, usize), p: &MatteParams, band_radius: f64, scale: f64) -> Vec { + let (w, h) = size; if !p.decontaminate { return vec![0.0; w * h]; } diff --git a/crates/track/src/roto/tests.rs b/crates/track/src/roto/tests.rs index 3c4645cb2..92c7bec04 100644 --- a/crates/track/src/roto/tests.rs +++ b/crates/track/src/roto/tests.rs @@ -246,7 +246,7 @@ fn refine_edge_recovers_a_soft_edge_and_decontaminates() { assert!(frac as f64 > 0.6 * n as f64, "{frac}/{n} fractional"); // Decontamination: edge pixels' colour is the foreground's, not a blue blend. let mut col = rgb.clone(); - let map = decontaminate(&mut col, &a, &seg.refine, W as usize, H as usize, &p, 9.0, 1.0); + let map = decontaminate(&mut col, &a, &seg.refine, (W as usize, H as usize), &p, 9.0, 1.0); let (mut blue_before, mut blue_after, mut m) = (0.0f64, 0.0f64, 0usize); for i in 0..a.len() { if a[i] > 0.2 && a[i] < 0.8 {