M13.3: Subspace Warp content-preserving mesh warp, face tracker library

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

crates/effects/src/warp_stab.rs

@@ -17,8 +17,11 @@
 //!   Edge Cropping ignores bad edge pixels of the input frames;
 //! - **Show Track Points** draws the analysed background features.
 //!
-//! *Subspace Warp* is approximated by a per-frame perspective warp (no content-preserving
-//! mesh warp yet), and *Rolling Shutter Ripple* is stored but has no effect.
+//! *Subspace Warp* warps each frame with a content-preserving mesh fitted to the analysis'
+//! subspace-smoothed feature trajectories ([`effectcraft_track::subspace`]); analyses made
+//! before trajectories were stored fall back to a perspective warp. *Rolling Shutter Ripple*
+//! picks the mesh: Automatic Reduction (12 × 8, stiffer) or Enhanced Reduction (twice the rows,
+//! softer), which follows the row-wise wobble of rolling-shutter shake more closely.
 //!
 //! Without an analysis (or with one made for another layer size) the input passes through.
 
@@ -30,7 +33,7 @@ use effectcraft_keyframe::Value;
 use effectcraft_project::ParamUi;
 use effectcraft_raster::{Image, Px};
 use effectcraft_track::Homography;
-use effectcraft_track::stabilize::{Framing, Method, Plan, StabResult, StabSettings, WarpAnalysis, plan};
+use effectcraft_track::stabilize::{Framing, Method, Plan, Ripple, StabResult, StabSettings, WarpAnalysis, plan};
 use rayon::prelude::*;
 
 use crate::{Buf, EffectCtx, EffectSpec, Params, num, p, popup, slider};
@@ -145,6 +148,7 @@ pub fn settings(params: &Params, fps: f64) -> StabSettings {
         // AE's 0…100 slider (50 = balanced) onto the planner's −100…100.
         crop_less_smooth_more: (f(params, "advanced/cropLessSmoothMore") - 50.0) * 2.0,
         fps,
+        ripple: if e(params, "advanced/rollingShutterRipple") == 1 { Ripple::Enhanced } else { Ripple::Automatic },
     }
 }
 
@@ -175,7 +179,7 @@ fn settings_key(s: &StabSettings) -> u64 {
     for v in [s.smoothness, s.max_scale, s.action_safe, s.additional_scale, s.crop_less_smooth_more, s.fps] {
         h.write_u64(v.to_bits());
     }
-    h.write(format!("{:?}{:?}{:?}{:?}", s.result, s.method, s.preserve_scale, s.framing).as_bytes());
+    h.write(format!("{:?}{:?}{:?}{:?}{:?}", s.result, s.method, s.preserve_scale, s.framing, s.ripple).as_bytes());
     h.finish()
 }
 
@@ -275,6 +279,13 @@ fn render(ctx: &EffectCtx, buf: Buf) -> Buf {
     }
     let sc = if buf.scale > 0.0 { buf.scale } else { 1.0 };
     let Some(inv) = pl.warps[k].inverse() else { return buf };
+    // Output point → source point of this frame (a mesh lookup for Subspace Warp).
+    let mesh = pl.meshes.get(k);
+    let finv = pl.framing.inverse().unwrap_or(Homography::IDENTITY);
+    let src_of = |p: [f64; 2]| match mesh {
+        Some(m) => m.sample(finv.apply(p)),
+        None => inv.apply(p),
+    };
     let (ow, oh) = (((size[0] * sc).ceil() as u32).max(1), ((size[1] * sc).ceil() as u32).max(1));
     let synth = matches!(settings(ctx.params, 30.0).framing, Framing::SynthesizeEdges);
     let crop = [
@@ -306,7 +317,7 @@ fn render(ctx: &EffectCtx, buf: Buf) -> Buf {
             {
                 continue;
             }
-            let q = inv.apply(p);
+            let q = src_of(p);
             if !synth {
                 *px = sample(src, q);
                 continue;
@@ -337,12 +348,12 @@ fn render(ctx: &EffectCtx, buf: Buf) -> Buf {
     });
     let need = need_m.into_inner().unwrap_or_else(|e| e.into_inner());
     if synth && need.iter().any(|w| *w > 0.0) {
-        synthesize(ctx, &a, k, &inv, &crop, sc, &mut out, &need);
+        synthesize(ctx, &a, k, &src_of, &crop, sc, &mut out, &need);
     }
     if b(ctx.params, "advanced/showTrackPoints") {
         let r = (2.5 * f(ctx.params, "advanced/trackPointSize") / 100.0 * sc).max(0.5);
         for q in &a.frames[k].points {
-            let o = pl.warps[k].apply([q[0] as f64, q[1] as f64]);
+            let Some(o) = pl.output_of(k, [q[0] as f64, q[1] as f64]) else { continue };
             dot(&mut out, [o[0] * sc, o[1] * sc], r, [0.25, 0.75, 1.0, 1.0]);
         }
     }
@@ -352,7 +363,16 @@ fn render(ctx: &EffectCtx, buf: Buf) -> Buf {
 /// Fill the pixels flagged in `need` (weight = how much synthesized content they take) from
 /// neighbouring frames, nearest first.
 #[allow(clippy::too_many_arguments)]
-fn synthesize(ctx: &EffectCtx, a: &WarpAnalysis, k: usize, inv: &Homography, crop: &[f64; 4], sc: f64, out: &mut Image, need: &[f32]) {
+fn synthesize(
+    ctx: &EffectCtx,
+    a: &WarpAnalysis,
+    k: usize,
+    src_of: &(dyn Fn([f64; 2]) -> [f64; 2] + Sync),
+    crop: &[f64; 4],
+    sc: f64,
+    out: &mut Image,
+    need: &[f32],
+) {
     let Some(host) = ctx.env.host else { return };
     let n = a.frames.len();
     let range = ((f(ctx.params, "advanced/synthesizeInputRange").max(0.0) / a.frame_duration.max(1e-6)).round() as usize).min(60);
@@ -376,7 +396,7 @@ fn synthesize(ctx: &EffectCtx, a: &WarpAnalysis, k: usize, inv: &Homography, cro
                 .filter(|i| !filled[**i])
                 .filter_map(|i| {
                     let (x, y) = (i % ow, i / ow);
-                    let q = inv.apply([(x as f64 + 0.5) / sc, (y as f64 + 0.5) / sc]);
+                    let q = src_of([(x as f64 + 0.5) / sc, (y as f64 + 0.5) / sc]);
                     let qj = m.apply(q);
                     inside(crop, qj).then(|| (*i, nb.img.sample_bilinear_clamped(qj[0] * nb.scale + nb.offset[0], qj[1] * nb.scale + nb.offset[1])))
                 })
@@ -412,7 +432,7 @@ mod tests {
     use effectcraft_track::stabilize::FrameMotion;
 
     fn analysis_json(n: usize, w: f64, h: f64, shift: impl Fn(usize) -> [f64; 2]) -> String {
-        let mut a = WarpAnalysis { version: 1, start: 0.0, frame_duration: 1.0 / 25.0, size: [w, h], detailed: false, frames: vec![] };
+        let mut a = WarpAnalysis { version: 1, start: 0.0, frame_duration: 1.0 / 25.0, size: [w, h], detailed: false, frames: vec![], tracks: vec![] };
         for k in 0..n {
             let mut fm = FrameMotion::default();
             if k > 0 {

crates/track/src/camtrack/features.rs

@@ -28,11 +28,19 @@ pub struct TrackAnalyzer {
     frame: u32,
     size: [f64; 2],
     factor: f64,
+    /// (analysis long side, wanted features) overriding the options' defaults.
+    budget: Option<(u32, usize)>,
 }
 
 impl TrackAnalyzer {
     pub fn new(size: [f64; 2], opts: AnalyzeOpts) -> TrackAnalyzer {
-        TrackAnalyzer { opts, prev: None, active: vec![], tracks: vec![], frame: 0, size, factor: 1.0 }
+        TrackAnalyzer { opts, prev: None, active: vec![], tracks: vec![], frame: 0, size, factor: 1.0, budget: None }
+    }
+
+    /// An analyzer working at most at `max_side` pixels with about `want` live features (the
+    /// Warp Stabilizer's trajectories for Subspace Warp).
+    pub fn with_budget(size: [f64; 2], max_side: u32, want: usize) -> TrackAnalyzer {
+        TrackAnalyzer { budget: Some((max_side, want.max(8))), ..TrackAnalyzer::new(size, AnalyzeOpts::default()) }
     }
 
     pub fn frames(&self) -> u32 {
@@ -45,6 +53,9 @@ impl TrackAnalyzer {
     }
 
     fn params(&self) -> (u32, usize) {
+        if let Some(b) = self.budget {
+            return b;
+        }
         if self.opts.detailed { (1280, 900) } else { (720, 500) }
     }
 

crates/track/src/lib.rs

@@ -19,6 +19,7 @@
 //! after subtracting the frame's `offset` (see [`Frame`]).
 
 pub mod camtrack;
+pub mod face;
 pub mod fit;
 pub mod klt;
 pub mod mask;
@@ -26,6 +27,7 @@ pub mod plane;
 pub mod roto;
 pub mod solve;
 pub mod stabilize;
+pub mod subspace;
 
 use effectcraft_raster::Image;
 use rayon::prelude::*;

crates/track/src/solve.rs

@@ -37,6 +37,12 @@ pub fn solve_linear<const N: usize>(mut a: [[f64; N]; N], mut b: [f64; N]) -> Op
 #[derive(Clone, Copy, Debug, PartialEq)]
 pub struct Homography(pub [[f64; 3]; 3]);
 
+impl Default for Homography {
+    fn default() -> Self {
+        Homography::IDENTITY
+    }
+}
+
 impl Homography {
     pub const IDENTITY: Homography = Homography([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]);
 

crates/track/src/stabilize.rs

@@ -18,14 +18,23 @@
 //!    Video Stabilization with Robust L1 Optimal Camera Paths" (CVPR 2011), applied here
 //!    per frame rather than through a linear program.
 //!
+//!
+//! *Subspace Warp* additionally follows long feature trajectories through the clip
+//! ([`crate::camtrack::TrackAnalyzer`]) and stabilises each frame with a content-preserving mesh
+//! warp fitted to the subspace-smoothed trajectories ([`crate::subspace`]): scenes with depth
+//! (parallax) that one homography per frame cannot hold still are stabilised locally. The
+//! framing (crop, auto-scale, Maximum Scale relaxation) works the same on meshes.
+//!
 //! Coordinates are layer pixels; transforms map *source* frame pixels to *stabilized* pixels.
 
 use serde::{Deserialize, Serialize};
 
 use crate::Frame;
-use crate::fit::{Model, ransac};
+use crate::camtrack::{Track2D, TrackAnalyzer};
+use crate::fit::{Model, least_squares, ransac};
 use crate::klt::{CornerOpts, GrayPyramid, LkOpts, analysis_factor, good_features, track};
 use crate::solve::Homography;
+use crate::subspace::{Mesh, MeshOpts, smooth_trajectories, solve_mesh};
 
 /// Inter-frame motion of one frame (from the previous frame to this one) for each model.
 #[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
@@ -82,6 +91,10 @@ pub struct WarpAnalysis {
     pub detailed: bool,
     /// Per frame; frame 0's motion is the identity.
     pub frames: Vec<FrameMotion>,
+    /// Long feature trajectories (layer pixels, 0.1 px) for Subspace Warp; empty in analyses
+    /// made before it (which then stabilise like Perspective).
+    #[serde(skip_serializing_if = "Vec::is_empty")]
+    pub tracks: Vec<Track2D>,
 }
 
 impl WarpAnalysis {
@@ -120,6 +133,7 @@ pub struct Analyzer {
     prev: Option<GrayPyramid>,
     frames: Vec<FrameMotion>,
     size: [f64; 2],
+    trajectories: TrackAnalyzer,
 }
 
 fn round_pts(p: &[[f64; 2]], n: usize) -> Vec<[f32; 2]> {
@@ -129,7 +143,8 @@ fn round_pts(p: &[[f64; 2]], n: usize) -> Vec<[f32; 2]> {
 
 impl Analyzer {
     pub fn new(size: [f64; 2], opts: AnalyzeOpts) -> Analyzer {
-        Analyzer { opts, prev: None, frames: vec![], size }
+        let (side, want) = if opts.detailed { (640, 260) } else { (480, 160) };
+        Analyzer { opts, prev: None, frames: vec![], size, trajectories: TrackAnalyzer::with_budget(size, side, want) }
     }
 
     pub fn len(&self) -> usize {
@@ -144,6 +159,7 @@ impl Analyzer {
         let (max_side, max_feat) = if self.opts.detailed { (960, 700) } else { (480, 350) };
         let factor = analysis_factor(frame.img.width, frame.img.height, max_side);
         let pyr = GrayPyramid::from_image(frame.img, frame.offset, factor, 3);
+        self.trajectories.push(frame);
         let mut fm = FrameMotion::default();
         if let Some(prev) = &self.prev
             && prev.width() == pyr.width()
@@ -184,7 +200,17 @@ impl Analyzer {
     }
 
     pub fn finish(self, start: f64, frame_duration: f64) -> WarpAnalysis {
-        WarpAnalysis { version: 1, start, frame_duration, size: self.size, detailed: self.opts.detailed, frames: self.frames }
+        let mut tracks = self.trajectories.finish(start, frame_duration).tracks;
+        // Short trajectories hardly constrain the subspace: keep those of at least 6 frames, at
+        // 0.1 px (smaller JSON).
+        tracks.retain(|t| t.pts.len() >= 6);
+        for (i, t) in tracks.iter_mut().enumerate() {
+            t.id = i as u32;
+            for p in t.pts.iter_mut() {
+                *p = [(p[0] * 10.0).round() / 10.0, (p[1] * 10.0).round() / 10.0];
+            }
+        }
+        WarpAnalysis { version: 2, start, frame_duration, size: self.size, detailed: self.opts.detailed, frames: self.frames, tracks }
     }
 }
 
@@ -205,11 +231,31 @@ pub enum Method {
     /// Position, Scale, Rotation.
     Similarity,
     Perspective,
-    /// Approximated by Perspective (see the module docs of the effect).
+    /// Content-preserving mesh warps fitted to subspace-smoothed trajectories.
     #[default]
     SubspaceWarp,
 }
 
+/// Advanced ▸ Rolling Shutter Ripple: how freely Subspace Warp's mesh may bend row by row (a
+/// rolling-shutter camera exposes rows at different times, so shake wobbles the frame).
+#[derive(Clone, Copy, Debug, Default, PartialEq, Eq, Hash)]
+pub enum Ripple {
+    /// A 12 × 8 mesh.
+    #[default]
+    Automatic,
+    /// Twice the rows and a softer shape term: follows row-wise wobble more closely.
+    Enhanced,
+}
+
+impl Ripple {
+    pub fn mesh(self) -> MeshOpts {
+        match self {
+            Ripple::Automatic => MeshOpts { cols: 12, rows: 8, alpha: 1.0, prior: 0.02 },
+            Ripple::Enhanced => MeshOpts { cols: 12, rows: 16, alpha: 0.5, prior: 0.02 },
+        }
+    }
+}
+
 /// Borders ▸ Framing.
 #[derive(Clone, Copy, Debug, Default, PartialEq, Eq, Hash)]
 pub enum Framing {
@@ -239,6 +285,7 @@ pub struct StabSettings {
     pub crop_less_smooth_more: f64,
     /// Frames per second (for Smoothness).
     pub fps: f64,
+    pub ripple: Ripple,
 }
 
 impl Default for StabSettings {
@@ -254,6 +301,7 @@ impl Default for StabSettings {
             additional_scale: 100.0,
             crop_less_smooth_more: 0.0,
             fps: 30.0,
+            ripple: Ripple::Automatic,
         }
     }
 }
@@ -269,8 +317,15 @@ impl StabSettings {
 /// The per-frame stabilization of a clip.
 #[derive(Clone, Debug, Default)]
 pub struct Plan {
-    /// Final source → output transform per frame (layer pixels), framing scale included.
+    /// Final source → output transform per frame (layer pixels), framing scale included. For
+    /// mesh warps, the closest homography (for overlays and edge synthesis).
     pub warps: Vec<Homography>,
+    /// Subspace Warp: per frame, the mesh (stabilised pixels before framing → source pixels);
+    /// empty for the other methods.
+    pub meshes: Vec<Mesh>,
+    /// The framing transform (auto-scale and Additional Scale about the centre) after the
+    /// stabilisation.
+    pub framing: Homography,
     /// The visible region in output pixels `[x0, y0, x1, y1]` (Stabilize, Crop / Auto-scale);
     /// `None` = everything the warp produces.
     pub crop: Option<[f64; 4]>,
@@ -280,6 +335,24 @@ pub struct Plan {
     pub valid_fraction: f64,
 }
 
+impl Plan {
+    /// The source position (frame `k`'s layer pixels) shown at output point `p`.
+    pub fn source_of(&self, k: usize, p: [f64; 2]) -> Option<[f64; 2]> {
+        if let Some(m) = self.meshes.get(k) {
+            let inv = self.framing.inverse()?;
+            return Some(m.sample(inv.apply(p)));
+        }
+        Some(self.warps.get(k)?.inverse()?.apply(p))
+    }
+    /// Where source point `q` of frame `k` lands in the output.
+    pub fn output_of(&self, k: usize, q: [f64; 2]) -> Option<[f64; 2]> {
+        if let Some(m) = self.meshes.get(k) {
+            return m.forward(q).map(|p| self.framing.apply(p));
+        }
+        Some(self.warps.get(k)?.apply(q))
+    }
+}
+
 fn convex_contains(q: &[[f64; 2]; 4], p: [f64; 2]) -> bool {
     let mut sign = 0.0;
     for i in 0..4 {
@@ -418,14 +491,100 @@ pub fn corrections(a: &WarpAnalysis, s: &StabSettings) -> Vec<Homography> {
     out
 }
 
+/// One frame's stabilisation before framing.
+#[derive(Clone, Debug)]
+enum Warp {
+    H(Homography),
+    M(Mesh),
+}
+
+impl Warp {
+    fn valid(&self, size: [f64; 2]) -> f64 {
+        match self {
+            Warp::H(h) => valid_fraction(h, size),
+            Warp::M(m) => m.valid_fraction(),
+        }
+    }
+    /// Blend towards the identity (`l` = 1 keeps the warp).
+    fn relax(&self, l: f64) -> Warp {
+        match self {
+            Warp::H(h) => Warp::H(Homography::IDENTITY.lerp(h, l)),
+            Warp::M(m) => Warp::M(m.lerp_identity(l)),
+        }
+    }
+    fn homography(&self) -> Homography {
+        match self {
+            Warp::H(h) => *h,
+            Warp::M(m) => m.homography(),
+        }
+    }
+}
+
+/// Remove the uniform scale of the correspondences' best similarity (Preserve Scale): the
+/// wanted positions are scaled back about their centroid.
+fn unscale(pairs: &mut [([f64; 2], [f64; 2])]) {
+    let (src, dst): (Vec<[f64; 2]>, Vec<[f64; 2]>) = pairs.iter().copied().unzip();
+    let Some(h) = least_squares(Model::Similarity, &src, &dst) else { return };
+    let m = h.0;
+    let k = (m[0][0] * m[1][1] - m[0][1] * m[1][0]).abs().sqrt();
+    if k < 1e-6 {
+        return;
+    }
+    let n = dst.len() as f64;
+    let c = [dst.iter().map(|q| q[0]).sum::<f64>() / n, dst.iter().map(|q| q[1]).sum::<f64>() / n];
+    for (_, d) in pairs.iter_mut() {
+        *d = [c[0] + (d[0] - c[0]) / k, c[1] + (d[1] - c[1]) / k];
+    }
+}
+
+/// Subspace Warp's meshes (stabilised → source, before framing) for every frame; `None` for the
+/// other methods and for analyses without trajectories. `base`: the perspective corrections
+/// ([`corrections`]), the meshes' prior.
+pub fn meshes(a: &WarpAnalysis, s: &StabSettings, base: &[Homography]) -> Option<Vec<Mesh>> {
+    if s.method != Method::SubspaceWarp || a.tracks.is_empty() || base.len() != a.frames.len() {
+        return None;
+    }
+    let n = a.frames.len();
+    let o = s.ripple.mesh();
+    let size = a.size;
+    let pairs: Vec<Vec<([f64; 2], [f64; 2])>> = match s.result {
+        StabResult::NoMotion => (0..n).map(|t| a.tracks.iter().filter_map(|tr| tr.at(t as u32)).map(|p| (p, base[t].apply(p))).collect()).collect(),
+        StabResult::SmoothMotion => {
+            let sigma = s.sigma();
+            if sigma < 0.3 {
+                return Some((0..n).map(|_| Mesh::identity(o.cols, o.rows, size)).collect());
+            }
+            smooth_trajectories(&a.tracks, n, sigma, 9)
+        }
+    };
+    use rayon::prelude::*;
+    Some(
+        pairs
+            .into_par_iter()
+            .enumerate()
+            .map(|(t, mut pr)| {
+                if s.preserve_scale && pr.len() >= 2 {
+                    unscale(&mut pr);
+                }
+                if pr.len() < 6 { Mesh::from_homography(o.cols, o.rows, size, &base[t]) } else { solve_mesh(size, &o, &pr, &base[t]) }
+            })
+            .collect(),
+    )
+}
+
 /// Stabilization plan: corrections plus framing.
 pub fn plan(a: &WarpAnalysis, s: &StabSettings) -> Plan {
-    let mut b = corrections(a, s);
+    use rayon::prelude::*;
+    let base = corrections(a, s);
+    let mut b: Vec<Warp> = match meshes(a, s, &base) {
+        Some(m) => m.into_iter().map(Warp::M).collect(),
+        None => base.into_iter().map(Warp::H).collect(),
+    };
     let size = a.size;
     let c = [size[0] / 2.0, size[1] / 2.0];
     let add = (s.additional_scale / 100.0).max(0.01);
     let z_add = Homography::scale_about(add, c);
-    let fracs = |b: &[Homography]| b.iter().map(|h| valid_fraction(h, size)).collect::<Vec<f64>>();
+    let fracs = |b: &[Warp]| b.par_iter().map(|w| w.valid(size)).collect::<Vec<f64>>();
     let mut fr = fracs(&b);
     let mut vmin = fr.iter().copied().fold(1.0, f64::min);
     let margin = (s.action_safe / 100.0).clamp(0.0, 0.45);
@@ -435,22 +594,24 @@ pub fn plan(a: &WarpAnalysis, s: &StabSettings) -> Plan {
         // Relax the correction on frames that need more than Maximum Scale.
         let target = need_cover / max_scale;
         let n = b.len();
-        let mut lam = vec![1.0f64; n];
-        for t in 0..n {
-            if fr[t] >= target {
-                continue;
-            }
-            let (mut lo, mut hi) = (0.0, 1.0);
-            for _ in 0..24 {
-                let mid = 0.5 * (lo + hi);
-                if valid_fraction(&Homography::IDENTITY.lerp(&b[t], mid), size) >= target {
-                    lo = mid;
-                } else {
-                    hi = mid;
+        let lam: Vec<f64> = (0..n)
+            .into_par_iter()
+            .map(|t| {
+                if fr[t] >= target {
+                    return 1.0;
                 }
-            }
-            lam[t] = lo;
-        }
+                let (mut lo, mut hi) = (0.0, 1.0);
+                for _ in 0..24 {
+                    let mid = 0.5 * (lo + hi);
+                    if b[t].relax(mid).valid(size) >= target {
+                        lo = mid;
+                    } else {
+                        hi = mid;
+                    }
+                }
+                lo
+            })
+            .collect();
         // Spread the relaxation smoothly (min filter, then a Gaussian of the same radius) so
         // the camera doesn't jump.
         let r = ((s.sigma() * 0.5).ceil() as usize).max(1);
@@ -468,7 +629,7 @@ pub fn plan(a: &WarpAnalysis, s: &StabSettings) -> Plan {
                 (acc / ws).min(lam[t])
             })
             .collect();
-        b = b.iter().zip(&smooth).map(|(h, l)| if *l >= 1.0 { *h } else { Homography::IDENTITY.lerp(h, *l) }).collect();
+        b = b.iter().zip(&smooth).map(|(w, l)| if *l >= 1.0 { w.clone() } else { w.relax(*l) }).collect();
         fr = fracs(&b);
         vmin = fr.iter().copied().fold(1.0, f64::min);
     }
@@ -483,13 +644,23 @@ pub fn plan(a: &WarpAnalysis, s: &StabSettings) -> Plan {
         }
     };
     let z = Homography::scale_about(auto, c);
-    let warps = b.iter().map(|h| z_add.then_after(&z.then_after(h))).collect();
-    Plan { warps, crop, auto_scale: auto, valid_fraction: vmin }
+    let framing = z_add.then_after(&z);
+    let warps = b.par_iter().map(|w| framing.then_after(&w.homography())).collect();
+    let meshes = b
+        .into_iter()
+        .filter_map(|w| match w {
+            Warp::M(m) => Some(m),
+            Warp::H(_) => None,
+        })
+        .collect();
+    Plan { warps, meshes, framing, crop, auto_scale: auto, valid_fraction: vmin }
 }
 
 #[cfg(test)]
 mod tests {
     use super::*;
+    use effectcraft_raster::Image;
+    use rayon::prelude::*;
 
     /// A synthetic analysis: a slow pan plus per-frame jitter.
     fn shaky(n: usize) -> (WarpAnalysis, Vec<[f64; 2]>) {
@@ -497,7 +668,7 @@ mod tests {
             let j = [((k * 7919) % 13) as f64 - 6.0, ((k * 104729) % 11) as f64 - 5.0];
             [2.0 * k as f64 + j[0], 0.5 * k as f64 + j[1]]
         };
-        let mut a = WarpAnalysis { version: 1, start: 0.0, frame_duration: 1.0 / 30.0, size: [640.0, 360.0], detailed: false, frames: vec![] };
+        let mut a = WarpAnalysis { version: 1, start: 0.0, frame_duration: 1.0 / 30.0, size: [640.0, 360.0], detailed: false, frames: vec![], tracks: vec![] };
         let path: Vec<[f64; 2]> = (0..n).map(pos).collect();
         for k in 0..n {
             let mut f = FrameMotion::default();
@@ -573,4 +744,139 @@ mod tests {
         let j = WarpAnalysis::from_json(&a.to_json()).unwrap();
         assert_eq!(j, a);
     }
+
+    // ------------------------------------------------------------ synthetic shaky footage
+
+    fn hash(i: i64, j: i64, k: i64) -> f64 {
+        let mut v = (i.wrapping_mul(73_856_093) ^ j.wrapping_mul(19_349_663) ^ k.wrapping_mul(83_492_791)) as u64;
+        v ^= v >> 13;
+        v = v.wrapping_mul(0x5bd1_e995);
+        v ^= v >> 15;
+        (v % 10_000) as f64 / 10_000.0
+    }
+
+    fn value_noise(x: f64, y: f64, cell: f64, seed: i64) -> f64 {
+        let (fx, fy) = (x / cell, y / cell);
+        let (i, j) = (fx.floor() as i64, fy.floor() as i64);
+        let (u, v) = (fx - i as f64, fy - j as f64);
+        let (u, v) = (u * u * (3.0 - 2.0 * u), v * v * (3.0 - 2.0 * v));
+        let a = hash(i, j, seed) * (1.0 - u) + hash(i + 1, j, seed) * u;
+        let b = hash(i, j + 1, seed) * (1.0 - u) + hash(i + 1, j + 1, seed) * u;
+        a * (1.0 - v) + b * v
+    }
+
+    fn texture(x: f64, y: f64) -> f32 {
+        (0.55 * value_noise(x, y, 5.0, 1) + 0.45 * value_noise(x, y, 13.0, 2)) as f32
+    }
+
+    /// A scene with depth: a near bump in the middle moves `1 + depth` times as much as the
+    /// background when the camera shakes (parallax a homography cannot model).
+    struct Scene {
+        size: [f64; 2],
+        cams: Vec<[f64; 2]>,
+    }
+
+    impl Scene {
+        fn new(n: usize) -> Scene {
+            let cams = (0..n)
+                .map(|t| {
+                    let j = [5.0 * (hash(t as i64, 7, 3) - 0.5), 5.0 * (hash(t as i64, 11, 5) - 0.5)];
+                    [0.4 * t as f64 + j[0], 0.15 * t as f64 + j[1]]
+                })
+                .collect();
+            Scene { size: [320.0, 180.0], cams }
+        }
+        fn depth(&self, s: [f64; 2]) -> f64 {
+            1.0 + 1.2 * (-((s[0] - 160.0).powi(2) + (s[1] - 90.0).powi(2)) / (2.0 * 38.0f64.powi(2))).exp()
+        }
+        /// Where scene point `s` appears on frame `t`.
+        fn image_of(&self, t: usize, s: [f64; 2]) -> [f64; 2] {
+            let (c, d) = (self.cams[t], self.depth(s));
+            [s[0] + c[0] * d, s[1] + c[1] * d]
+        }
+        fn frame(&self, t: usize) -> Image {
+            let (w, h) = (self.size[0] as u32, self.size[1] as u32);
+            let mut im = Image::new(w, h);
+            let c = self.cams[t];
+            im.rows_mut().for_each(|(y, row)| {
+                for (x, px) in row.iter_mut().enumerate() {
+                    let p = [x as f64 + 0.5, y as f64 + 0.5];
+                    let mut s = p;
+                    for _ in 0..10 {
+                        let d = self.depth(s);
+                        s = [p[0] - c[0] * d, p[1] - c[1] * d];
+                    }
+                    let v = texture(s[0], s[1]);
+                    *px = [v, v, v, 1.0];
+                }
+            });
+            im
+        }
+    }
+
+    fn residual_jitter(scene: &Scene, out: impl Fn(usize, [f64; 2]) -> Option<[f64; 2]>) -> f64 {
+        let n = scene.cams.len();
+        let mut tot = 0.0;
+        let mut cnt = 0;
+        for gy in 0..5 {
+            for gx in 0..7 {
+                let s = [70.0 + gx as f64 * 30.0, 45.0 + gy as f64 * 22.0];
+                let ys: Vec<Option<[f64; 2]>> = (0..n).map(|t| out(t, scene.image_of(t, s))).collect();
+                for t in 6..n - 6 {
+                    if let (Some(a), Some(b), Some(c)) = (ys[t - 1], ys[t], ys[t + 1]) {
+                        tot += (c[0] - 2.0 * b[0] + a[0]).hypot(c[1] - 2.0 * b[1] + a[1]);
+                        cnt += 1;
+                    }
+                }
+            }
+        }
+        tot / cnt.max(1) as f64
+    }
+
+    #[test]
+    fn subspace_warp_beats_perspective_on_parallax() {
+        let n = 48;
+        let scene = Scene::new(n);
+        let mut an = Analyzer::new(scene.size, AnalyzeOpts::default());
+        for t in 0..n {
+            let img = scene.frame(t);
+            an.push(&Frame { img: &img, offset: [0.0; 2] });
+        }
+        let a = an.finish(0.0, 1.0 / 24.0);
+        assert!(a.tracks.len() >= 40, "trajectories: {}", a.tracks.len());
+        // The analysis round-trips with its trajectories.
+        assert_eq!(WarpAnalysis::from_json(&a.to_json()).unwrap(), a);
+        let base = StabSettings { framing: Framing::StabilizeOnly, fps: 24.0, ..Default::default() };
+        let before = residual_jitter(&scene, |_, x| Some(x));
+        let ppl = plan(&a, &StabSettings { method: Method::Perspective, ..base });
+        assert!(ppl.meshes.is_empty());
+        let persp = residual_jitter(&scene, |t, x| ppl.output_of(t, x));
+        let spl = plan(&a, &StabSettings { method: Method::SubspaceWarp, ..base });
+        assert_eq!(spl.meshes.len(), n);
+        let sub = residual_jitter(&scene, |t, x| spl.output_of(t, x));
+        let epl = plan(&a, &StabSettings { method: Method::SubspaceWarp, ripple: Ripple::Enhanced, ..base });
+        let enh = residual_jitter(&scene, |t, x| epl.output_of(t, x));
+        eprintln!("jitter: before {before:.3}, perspective {persp:.3}, subspace {sub:.3}, enhanced {enh:.3}");
+        assert!(persp < 0.8 * before, "perspective {persp} vs {before}");
+        assert!(sub < 0.6 * persp, "subspace {sub} vs perspective {persp}");
+        assert!(enh < 0.6 * persp, "enhanced {enh} vs perspective {persp}");
+        // source_of and output_of agree.
+        let q = [150.0, 80.0];
+        let o = spl.output_of(10, q).unwrap();
+        let back = spl.source_of(10, o).unwrap();
+        assert!((back[0] - q[0]).abs() < 1e-3 && (back[1] - q[1]).abs() < 1e-3);
+        // Auto-scale framing covers the frame with mesh warps too.
+        let fpl = plan(&a, &StabSettings { method: Method::SubspaceWarp, fps: 24.0, ..Default::default() });
+        assert!(fpl.auto_scale >= 1.0);
+        for k in [0, n / 2, n - 1] {
+            for p in [[1.0, 1.0], [319.0, 1.0], [319.0, 179.0], [1.0, 179.0]] {
+                let q = fpl.source_of(k, p).unwrap();
+                let inside = q[0] >= -0.5 && q[1] >= -0.5 && q[0] <= 320.5 && q[1] <= 180.5;
+                assert!(inside || fpl.crop.is_some(), "frame {k}: {p:?} → {q:?}");
+            }
+        }
+        // Preserve Scale keeps working with meshes.
+        let psp = plan(&a, &StabSettings { method: Method::SubspaceWarp, preserve_scale: true, ..base });
+        assert_eq!(psp.meshes.len(), n);
+    }
 }

crates/track/src/subspace.rs

@@ -0,0 +1,447 @@
+//! Warp Stabilizer ▸ Method ▸ **Subspace Warp**: content-preserving warps fitted to smoothed
+//! feature trajectories, clean-room from the published method:
+//!
+//! - Liu, Gleicher, Wang, Jin & Agarwala, "Subspace Video Stabilization", ACM TOG 30(1), 2011:
+//!   the feature trajectories of a short window form a matrix `M` (two rows per trajectory, one
+//!   column per frame) that is close to low rank (rank 9 for a moving camera over a static
+//!   scene); its factorisation `M ≈ W C` gives *eigen-trajectories* `C`, which are smoothed (here
+//!   with a Gaussian) instead of the trajectories themselves. Every trajectory's smoothed
+//!   position is `W_i Ĉ_t`; trajectories that do not span the whole window get their
+//!   coefficients `W_i` by least squares over the frames they cover (the paper's handling of
+//!   incomplete trajectories in a moving factorisation). Smoothing in the subspace keeps the
+//!   smoothed trajectories consistent with one another, so the warp does not tear the scene.
+//! - Liu, Gleicher, Jin & Agarwala, "Content-Preserving Warps for 3D Video Stabilization", ACM
+//!   TOG 28(3), 2009: each frame is warped by a grid mesh fitted to the feature correspondences
+//!   (original → smoothed position) in least squares, with a **similarity term** asking every
+//!   grid triangle to keep its shape (`V1 = V2 + u (V3 − V2) + v R90 (V3 − V2)` with `u, v` from
+//!   the undeformed grid) so the warp bends where the scene's depth demands it but stays locally
+//!   rigid elsewhere.
+//!
+//! The mesh is stored *inversely*: a regular grid over the stabilised frame whose vertices hold
+//! the source position they show ([`Mesh`]), so rendering is a bilinear lookup per output pixel.
+//! The data term is the same linear constraint read the other way (the smoothed position's grid
+//! cell, bilinearly, must land on the original position). A weak prior ties every vertex to the
+//! frame's perspective correction, which fills feature-less regions and keeps the system well
+//! posed.
+
+use rayon::prelude::*;
+
+use crate::camtrack::Track2D;
+use crate::camtrack::linalg::{cholesky_solve, sym_eigen};
+use crate::fit::{Model, least_squares};
+use crate::solve::Homography;
+
+/// A warp mesh: a regular `cols × rows` grid over the stabilised frame (`size`), each vertex
+/// holding the source-frame position it shows.
+#[derive(Clone, Debug, PartialEq)]
+pub struct Mesh {
+    pub cols: usize,
+    pub rows: usize,
+    pub size: [f64; 2],
+    /// `(rows + 1) × (cols + 1)` source positions, row-major.
+    pub verts: Vec<[f64; 2]>,
+}
+
+impl Mesh {
+    /// The identity mesh.
+    pub fn identity(cols: usize, rows: usize, size: [f64; 2]) -> Mesh {
+        let mut m = Mesh { cols: cols.max(1), rows: rows.max(1), size, verts: vec![] };
+        m.verts = (0..=m.rows).flat_map(|j| (0..=m.cols).map(move |i| (i, j))).map(|(i, j)| m.grid(i, j)).collect();
+        m
+    }
+
+    /// The mesh of a homography (`h`: source → stabilised).
+    pub fn from_homography(cols: usize, rows: usize, size: [f64; 2], h: &Homography) -> Mesh {
+        let inv = h.inverse().unwrap_or(Homography::IDENTITY);
+        let mut m = Mesh::identity(cols, rows, size);
+        for v in m.verts.iter_mut() {
+            *v = inv.apply(*v);
+        }
+        m
+    }
+
+    fn cell(&self) -> [f64; 2] {
+        [self.size[0] / self.cols as f64, self.size[1] / self.rows as f64]
+    }
+
+    /// The undeformed position of vertex `(i, j)` (column, row).
+    pub fn grid(&self, i: usize, j: usize) -> [f64; 2] {
+        let c = self.cell();
+        [i as f64 * c[0], j as f64 * c[1]]
+    }
+
+    fn idx(&self, i: usize, j: usize) -> usize {
+        j * (self.cols + 1) + i
+    }
+
+    /// Cell and bilinear coordinates of a stabilised point (the outer cells extend beyond the
+    /// frame, so points outside it extrapolate).
+    fn locate(&self, p: [f64; 2]) -> (usize, usize, f64, f64) {
+        let c = self.cell();
+        let fx = p[0] / c[0];
+        let fy = p[1] / c[1];
+        let i = (fx.floor().max(0.0) as usize).min(self.cols - 1);
+        let j = (fy.floor().max(0.0) as usize).min(self.rows - 1);
+        (i, j, fx - i as f64, fy - j as f64)
+    }
+
+    /// The source position shown at stabilised point `p`.
+    pub fn sample(&self, p: [f64; 2]) -> [f64; 2] {
+        let (i, j, u, v) = self.locate(p);
+        let a = self.verts[self.idx(i, j)];
+        let b = self.verts[self.idx(i + 1, j)];
+        let c = self.verts[self.idx(i, j + 1)];
+        let d = self.verts[self.idx(i + 1, j + 1)];
+        let w = [(1.0 - u) * (1.0 - v), u * (1.0 - v), (1.0 - u) * v, u * v];
+        [0, 1].map(|k| w[0] * a[k] + w[1] * b[k] + w[2] * c[k] + w[3] * d[k])
+    }
+
+    /// Where source point `q` lands in the stabilised frame (Newton iterations on
+    /// [`Mesh::sample`]); `None` when it does not converge (a folded mesh).
+    pub fn forward(&self, q: [f64; 2]) -> Option<[f64; 2]> {
+        // Start from the mesh's own offset at q.
+        let s = self.sample(q);
+        let mut p = [2.0 * q[0] - s[0], 2.0 * q[1] - s[1]];
+        for _ in 0..30 {
+            let f = self.sample(p);
+            let e = [f[0] - q[0], f[1] - q[1]];
+            if e[0].hypot(e[1]) < 1e-6 {
+                return Some(p);
+            }
+            let h = 0.5;
+            let fx = self.sample([p[0] + h, p[1]]);
+            let fy = self.sample([p[0], p[1] + h]);
+            let j = [[(fx[0] - f[0]) / h, (fy[0] - f[0]) / h], [(fx[1] - f[1]) / h, (fy[1] - f[1]) / h]];
+            let det = j[0][0] * j[1][1] - j[0][1] * j[1][0];
+            if det.abs() < 1e-12 {
+                return None;
+            }
+            p[0] -= (j[1][1] * e[0] - j[0][1] * e[1]) / det;
+            p[1] -= (-j[1][0] * e[0] + j[0][0] * e[1]) / det;
+        }
+        let f = self.sample(p);
+        ((f[0] - q[0]).hypot(f[1] - q[1]) < 1e-3).then_some(p)
+    }
+
+    /// `(1 − t) · identity + t · self` per vertex.
+    pub fn lerp_identity(&self, t: f64) -> Mesh {
+        let mut m = self.clone();
+        for j in 0..=self.rows {
+            for i in 0..=self.cols {
+                let g = self.grid(i, j);
+                let k = self.idx(i, j);
+                m.verts[k] = [g[0] + (self.verts[k][0] - g[0]) * t, g[1] + (self.verts[k][1] - g[1]) * t];
+            }
+        }
+        m
+    }
+
+    /// Largest `s ∈ [0, 1]` such that the centred `s·W × s·H` rectangle of the stabilised frame
+    /// shows only source pixels (its outline, sampled densely, maps inside the source frame).
+    pub fn valid_fraction(&self) -> f64 {
+        let [w, h] = self.size;
+        let c = [w / 2.0, h / 2.0];
+        let tol = 1e-6;
+        let inside = |q: [f64; 2]| q[0] >= -tol && q[1] >= -tol && q[0] <= w + tol && q[1] <= h + tol;
+        let fits = |s: f64| {
+            let (hw, hh) = (s * w / 2.0, s * h / 2.0);
+            let n = 48;
+            (0..=n).all(|k| {
+                let t = k as f64 / n as f64;
+                let x = c[0] - hw + 2.0 * hw * t;
+                let y = c[1] - hh + 2.0 * hh * t;
+                inside(self.sample([x, c[1] - hh]))
+                    && inside(self.sample([x, c[1] + hh]))
+                    && inside(self.sample([c[0] - hw, y]))
+                    && inside(self.sample([c[0] + hw, y]))
+            })
+        };
+        if fits(1.0) {
+            return 1.0;
+        }
+        if !fits(0.0) {
+            return 0.0;
+        }
+        let (mut lo, mut hi) = (0.0, 1.0);
+        for _ in 0..24 {
+            let mid = 0.5 * (lo + hi);
+            if fits(mid) {
+                lo = mid;
+            } else {
+                hi = mid;
+            }
+        }
+        lo
+    }
+
+    /// The homography (source → stabilised) closest to the mesh in least squares.
+    pub fn homography(&self) -> Homography {
+        let mut src = vec![];
+        let mut dst = vec![];
+        for j in 0..=self.rows {
+            for i in 0..=self.cols {
+                src.push(self.verts[self.idx(i, j)]);
+                dst.push(self.grid(i, j));
+            }
+        }
+        least_squares(Model::Homography, &src, &dst).unwrap_or(Homography::IDENTITY)
+    }
+}
+
+/// Mesh resolution and rigidity.
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub struct MeshOpts {
+    pub cols: usize,
+    pub rows: usize,
+    /// Weight of the similarity (shape-preserving) term relative to one feature.
+    pub alpha: f64,
+    /// Weight of the prior towards the frame's perspective correction, per vertex.
+    pub prior: f64,
+}
+
+impl Default for MeshOpts {
+    fn default() -> Self {
+        MeshOpts { cols: 12, rows: 8, alpha: 1.0, prior: 0.02 }
+    }
+}
+
+/// Normal equations `AᵀA x = Aᵀb`, accumulated row by row.
+struct Normal {
+    n: usize,
+    a: Vec<f64>,
+    b: Vec<f64>,
+}
+
+impl Normal {
+    fn row(&mut self, coefs: &[(usize, f64)], rhs: f64, w: f64) {
+        for (ia, ca) in coefs {
+            self.b[*ia] += w * ca * rhs;
+            for (ib, cb) in coefs {
+                self.a[ia * self.n + ib] += w * ca * cb;
+            }
+        }
+    }
+}
+
+/// Fit the content-preserving warp of one frame: `pairs` are (source position, wanted
+/// stabilised position) of its features; `base` is the frame's perspective correction (source →
+/// stabilised), the prior.
+pub fn solve_mesh(size: [f64; 2], o: &MeshOpts, pairs: &[([f64; 2], [f64; 2])], base: &Homography) -> Mesh {
+    let mut m = Mesh::from_homography(o.cols, o.rows, size, base);
+    let nv = (m.cols + 1) * (m.rows + 1);
+    let n = 2 * nv;
+    let mut ne = Normal { n, a: vec![0.0; n * n], b: vec![0.0; n] };
+    // Prior.
+    for k in 0..nv {
+        let v = m.verts[k];
+        ne.row(&[(2 * k, 1.0)], v[0], o.prior);
+        ne.row(&[(2 * k + 1, 1.0)], v[1], o.prior);
+    }
+    // Data: the bilinear combination of the wanted position's cell lands on the source position.
+    let mut used = 0;
+    for (src, dst) in pairs {
+        if !(dst[0] >= 0.0 && dst[1] >= 0.0 && dst[0] <= size[0] && dst[1] <= size[1] && src[0].is_finite() && src[1].is_finite()) {
+            continue;
+        }
+        let (i, j, u, v) = m.locate(*dst);
+        let ks = [m.idx(i, j), m.idx(i + 1, j), m.idx(i, j + 1), m.idx(i + 1, j + 1)];
+        let w = [(1.0 - u) * (1.0 - v), u * (1.0 - v), (1.0 - u) * v, u * v];
+        for d in 0..2 {
+            let coefs: Vec<(usize, f64)> = ks.iter().zip(w).map(|(k, w)| (2 * k + d, w)).collect();
+            ne.row(&coefs, src[d], 1.0);
+        }
+        used += 1;
+    }
+    if used == 0 {
+        return m;
+    }
+    // Similarity: four triangles per cell, each corner against its two cell neighbours.
+    let a2 = o.alpha * o.alpha;
+    for j in 0..m.rows {
+        for i in 0..m.cols {
+            let (a, b, c, d) = ((i, j), (i + 1, j), (i + 1, j + 1), (i, j + 1));
+            for (v1, v2, v3) in [(a, b, d), (b, c, a), (c, d, b), (d, a, c)] {
+                let (g1, g2, g3) = (m.grid(v1.0, v1.1), m.grid(v2.0, v2.1), m.grid(v3.0, v3.1));
+                let dv = [g3[0] - g2[0], g3[1] - g2[1]];
+                let e = [g1[0] - g2[0], g1[1] - g2[1]];
+                let l2 = dv[0] * dv[0] + dv[1] * dv[1];
+                // R90(d) = (−d.y, d.x).
+                let u = (e[0] * dv[0] + e[1] * dv[1]) / l2;
+                let v = (e[0] * -dv[1] + e[1] * dv[0]) / l2;
+                let (k1, k2, k3) = (m.idx(v1.0, v1.1), m.idx(v2.0, v2.1), m.idx(v3.0, v3.1));
+                // x: V1x − V2x − u (V3x − V2x) + v (V3y − V2y) = 0
+                ne.row(&[(2 * k1, 1.0), (2 * k2, u - 1.0), (2 * k3, -u), (2 * k3 + 1, v), (2 * k2 + 1, -v)], 0.0, a2);
+                // y: V1y − V2y − u (V3y − V2y) − v (V3x − V2x) = 0
+                ne.row(&[(2 * k1 + 1, 1.0), (2 * k2 + 1, u - 1.0), (2 * k3 + 1, -u), (2 * k3, -v), (2 * k2, v)], 0.0, a2);
+            }
+        }
+    }
+    let mut x = ne.b;
+    if cholesky_solve(&mut ne.a, n, &mut x) {
+        for k in 0..nv {
+            m.verts[k] = [x[2 * k], x[2 * k + 1]];
+        }
+    }
+    m
+}
+
+/// Gaussian weights of the frames `lo..=hi` around `t`.
+fn gauss(lo: usize, hi: usize, t: usize, sigma: f64) -> Vec<f64> {
+    let g: Vec<f64> = (lo..=hi).map(|j| (-((j as f64 - t as f64).powi(2)) / (2.0 * sigma * sigma)).exp()).collect();
+    let s: f64 = g.iter().sum();
+    g.into_iter().map(|v| v / s.max(1e-300)).collect()
+}
+
+/// Subspace-smoothed feature positions: for every frame `t < frames`, the (original, smoothed)
+/// positions of the trajectories seen on it. `sigma`: Gaussian smoothing in frames; `rank`: the
+/// subspace dimension (9 in the paper). Frames without enough trajectories get no pairs.
+pub fn smooth_trajectories(tracks: &[Track2D], frames: usize, sigma: f64, rank: usize) -> Vec<Vec<([f64; 2], [f64; 2])>> {
+    let max_r = ((3.0 * sigma).ceil() as usize).clamp(1, 30);
+    (0..frames)
+        .into_par_iter()
+        .map(|t| {
+            let t32 = t as u32;
+            let mut r = max_r;
+            // The window around t: shrink it until enough trajectories span it.
+            let (lo, hi, complete) = loop {
+                let lo = t.saturating_sub(r);
+                let hi = (t + r).min(frames.saturating_sub(1));
+                let complete: Vec<&Track2D> = tracks.iter().filter(|tr| tr.start as usize <= lo && tr.end() as usize > hi).collect();
+                if complete.len() >= 8 || r <= 2 {
+                    break (lo, hi, complete);
+                }
+                r = (r * 2 / 3).max(2);
+            };
+            if complete.len() < 4 {
+                return vec![];
+            }
+            let m = hi - lo + 1;
+            let rows = 2 * complete.len();
+            let k = rank.min(m).min(rows).max(1);
+            // M: rows × m (row-major).
+            let mut mat = vec![0.0; rows * m];
+            for (i, tr) in complete.iter().enumerate() {
+                for c in 0..m {
+                    let p = tr.at((lo + c) as u32).unwrap_or([0.0; 2]);
+                    mat[(2 * i) * m + c] = p[0];
+                    mat[(2 * i + 1) * m + c] = p[1];
+                }
+            }
+            // Gram MᵀM and its top eigenvectors: the eigen-trajectories (orthonormal rows of C).
+            let mut gram = vec![0.0; m * m];
+            for a in 0..m {
+                for b in a..m {
+                    let s: f64 = (0..rows).map(|r| mat[r * m + a] * mat[r * m + b]).sum();
+                    gram[a * m + b] = s;
+                    gram[b * m + a] = s;
+                }
+            }
+            let (_, vecs) = sym_eigen(&gram, m);
+            // C[q][c] = vecs[c][m − 1 − q] (largest first).
+            let cmat: Vec<Vec<f64>> = (0..k).map(|q| (0..m).map(|c| vecs[c * m + (m - 1 - q)]).collect()).collect();
+            let g = gauss(lo, hi, t, sigma.max(1e-6));
+            let chat: Vec<f64> = cmat.iter().map(|row| row.iter().zip(&g).map(|(a, b)| a * b).sum()).collect();
+            let tc = t - lo;
+            let mut out = Vec::new();
+            // Complete trajectories: W = M Cᵀ.
+            for (i, _) in complete.iter().enumerate() {
+                let src = [mat[(2 * i) * m + tc], mat[(2 * i + 1) * m + tc]];
+                let mut dst = [0.0; 2];
+                for (d, slot) in dst.iter_mut().enumerate() {
+                    let row = &mat[(2 * i + d) * m..(2 * i + d + 1) * m];
+                    for q in 0..k {
+                        let w: f64 = row.iter().zip(&cmat[q]).map(|(a, b)| a * b).sum();
+                        *slot += w * chat[q];
+                    }
+                }
+                out.push((src, dst));
+            }
+            // Incomplete trajectories seen on t: coefficients by least squares over the frames
+            // they cover in the window.
+            let need = (k + 3).max(m / 2);
+            for tr in tracks.iter().filter(|tr| tr.at(t32).is_some() && !(tr.start as usize <= lo && tr.end() as usize > hi)) {
+                let cols: Vec<usize> = (0..m).filter(|c| tr.at((lo + c) as u32).is_some()).collect();
+                if cols.len() < need {
+                    continue;
+                }
+                let mut ata = vec![0.0; k * k];
+                let mut atb = [vec![0.0; k], vec![0.0; k]];
+                for &c in &cols {
+                    let p = tr.at((lo + c) as u32).expect("covered");
+                    for a in 0..k {
+                        for d in 0..2 {
+                            atb[d][a] += cmat[a][c] * p[d];
+                        }
+                        for b in 0..k {
+                            ata[a * k + b] += cmat[a][c] * cmat[b][c];
+                        }
+                    }
+                }
+                for a in 0..k {
+                    ata[a * k + a] += 1e-9;
+                }
+                let mut dst = [0.0; 2];
+                let mut ok = true;
+                for d in 0..2 {
+                    let mut x = atb[d].clone();
+                    let mut aa = ata.clone();
+                    if !cholesky_solve(&mut aa, k, &mut x) {
+                        ok = false;
+                        break;
+                    }
+                    dst[d] = x.iter().zip(&chat).map(|(a, b)| a * b).sum();
+                }
+                if ok {
+                    out.push((tr.at(t32).expect("seen"), dst));
+                }
+            }
+            out
+        })
+        .collect()
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn mesh_identity_sample_and_forward() {
+        let m = Mesh::identity(4, 3, [100.0, 60.0]);
+        assert_eq!(m.sample([37.0, 21.0]), [37.0, 21.0]);
+        assert_eq!(m.valid_fraction(), 1.0);
+        let h = Homography([[1.0, 0.02, 5.0], [-0.01, 1.0, -3.0], [0.0, 0.0, 1.0]]);
+        let mh = Mesh::from_homography(4, 3, [100.0, 60.0], &h);
+        // Affine maps are exact on a bilinear mesh.
+        let q = [40.0, 30.0];
+        let p = mh.forward(q).unwrap();
+        let e = h.apply(q);
+        assert!((p[0] - e[0]).abs() < 1e-6 && (p[1] - e[1]).abs() < 1e-6, "{p:?} {e:?}");
+        assert!(mh.homography().max_diff(&h) < 1e-6);
+        assert!(mh.valid_fraction() < 1.0);
+    }
+
+    #[test]
+    fn mesh_fits_a_translation_and_bends_for_local_motion() {
+        let size = [200.0, 120.0];
+        let mut pairs = vec![];
+        for y in 0..12 {
+            for x in 0..20 {
+                let p = [5.0 + x as f64 * 10.0, 5.0 + y as f64 * 10.0];
+                // Everything moves by (3, −2); a blob near the centre by (8, −2).
+                let bump = (-((p[0] - 100.0).powi(2) + (p[1] - 60.0).powi(2)) / (2.0 * 20.0f64.powi(2))).exp();
+                pairs.push((p, [p[0] + 3.0 + 5.0 * bump, p[1] - 2.0]));
+            }
+        }
+        let base = Homography::translation([3.0, -2.0]);
+        let m = solve_mesh(size, &MeshOpts { cols: 10, rows: 6, alpha: 0.5, prior: 0.01 }, &pairs, &base);
+        let mut err_mesh = 0.0;
+        let mut err_h = 0.0;
+        for (s, d) in &pairs {
+            let q = m.sample(*d);
+            err_mesh += (q[0] - s[0]).hypot(q[1] - s[1]);
+            let qh = base.inverse().unwrap().apply(*d);
+            err_h += (qh[0] - s[0]).hypot(qh[1] - s[1]);
+        }
+        assert!(err_mesh < 0.35 * err_h, "{err_mesh} vs {err_h}");
+    }
+}