diff --git a/crates/effects/src/camera_tracker.rs b/crates/effects/src/camera_tracker.rs
index dd97b585..ff4e51fe 100644
--- a/crates/effects/src/camera_tracker.rs
+++ b/crates/effects/src/camera_tracker.rs
@@ -8,8 +8,12 @@
 //! Shot Type, Angle of View, Solve Method and deleted points it was solved with, so changing
 //! those re-solves without re-tracking). **Method Used** and **Average Error** report the solve.
 //!
-//! The effect passes its input through; with **Render Track Points** on it draws the solved
-//! points (as the viewer shows them) into the frame. The viewer overlay, target and the Create
+//! **Solve Lens Distortion** (Advanced; also on with Detailed Analysis) solves radial lens
+//! distortion `k1, k2` jointly with the camera; **Undistort Footage** then renders the input
+//! undistorted (the pinhole image the solved camera sees) so 3D layers line up with it.
+//!
+//! The effect passes its input through (undistorted when asked); with **Render Track Points**
+//! on it draws the solved points (as the viewer shows them) into the frame. The viewer overlay, target and the Create
 //! commands live in the engine and the UI.
 
 use std::collections::HashMap;
@@ -19,6 +23,7 @@ use effectcraft_keyframe::Value;
 use effectcraft_project::ParamUi;
 use effectcraft_raster::{Image, Px};
 use effectcraft_track::camtrack::{CameraSolve, CameraTracks, ShotType, SolveMethod, SolveSettings, linalg, point_color};
+use rayon::prelude::*;
 
 use crate::warp_stab::param;
 use crate::{Buf, EffectCtx, EffectSpec, Params, num, p, popup, slider};
@@ -58,6 +63,8 @@ pub fn specs() -> Vec<EffectSpec> {
             p(METHOD_USED, "Method Used", Value::Str(String::new()), ParamUi::Hidden),
             p(AVERAGE_ERROR, "Average Error", num(0.0), ParamUi::Hidden),
             p("advanced/detailedAnalysis", "Detailed Analysis", Value::Bool(false), ParamUi::Checkbox),
+            p("advanced/lensDistortion", "Solve Lens Distortion", Value::Bool(false), ParamUi::Checkbox),
+            p("advanced/undistort", "Undistort Footage", Value::Bool(false), ParamUi::Checkbox),
             p("advanced/autoDeletePoints", "Auto-delete Points Across Time", Value::Bool(true), ParamUi::Checkbox),
             p("advanced/hideWarningBanner", "Hide Warning Banner", Value::Bool(false), ParamUi::Checkbox),
             p(TRACKS, "Tracks", Value::Str(String::new()), ParamUi::Hidden),
@@ -104,13 +111,15 @@ pub fn settings(params: &Params) -> SolveSettings {
         hfov: f(params, "horizontalAngleOfView"),
         method: SolveMethod::ALL.get(e(params, "advanced/solveMethod") as usize).copied().unwrap_or_default(),
         deleted: del,
+        lens_distortion: b(params, "advanced/lensDistortion") || b(params, "advanced/detailedAnalysis"),
     }
 }
 
 /// What a solve depends on besides the tracks.
 pub fn solve_key(st: &SolveSettings) -> String {
     let hfov = if st.shot == ShotType::SpecifyAngle { format!("{:.4}", st.hfov) } else { String::new() };
-    format!("{:?}|{hfov}|{:?}|{:?}", st.shot, st.method, st.deleted)
+    let lens = if st.lens_distortion { "|lens" } else { "" };
+    format!("{:?}|{hfov}|{:?}|{:?}{lens}", st.shot, st.method, st.deleted)
 }
 
 fn fnv(bytes: &[u8]) -> u64 {
@@ -165,18 +174,49 @@ pub fn target_radius(size_pct: f64, z: f64, median_z: f64) -> f64 {
     (6.0 * size_pct / 100.0 * (median_z / z.max(1e-9))).clamp(1.5, 40.0)
 }
 
-/// The solved points visible on frame `k`: (id, layer pixel position, camera depth).
+/// The solved points visible on frame `k`: (id, layer pixel position in the footage, camera
+/// depth).
 pub fn projected(solve: &CameraSolve, k: usize) -> Vec<(u32, [f64; 2], f64)> {
+    projected_with(solve, k, false)
+}
+
+/// [`projected`] into the footage, or (`undistorted`) into the undistorted image.
+pub fn projected_with(solve: &CameraSolve, k: usize, undistorted: bool) -> Vec<(u32, [f64; 2], f64)> {
     let Some(cam) = solve.frames.get(k) else { return vec![] };
     solve
         .visible(k)
         .filter_map(|p| {
             let c = cam.to_cam(p.pos);
-            cam.project(solve.size, p.pos).map(|q| (p.id, q, c[2]))
+            solve.project(k, p.pos, undistorted).map(|q| (p.id, q, c[2]))
         })
         .collect()
 }
 
+/// Whether the instance renders its input undistorted (Undistort Footage with a solved
+/// distortion).
+pub fn undistorts(params: &Params, solve: &CameraSolve) -> bool {
+    b(params, "advanced/undistort") && solve.distortion.is_some_and(|d| !d.is_none())
+}
+
+/// The input resampled through the solved distortion: every output pixel shows the footage
+/// at its distorted position (transparent where that falls outside the layer).
+fn undistort(buf: &Buf, sol: &CameraSolve) -> Buf {
+    let sc = if buf.scale > 0.0 { buf.scale } else { 1.0 };
+    let mut out = Image::new(buf.img.width, buf.img.height);
+    let size = sol.size;
+    out.rows_mut().for_each(|(y, row)| {
+        for (x, px) in row.iter_mut().enumerate() {
+            let p = [(x as f64 + 0.5 - buf.offset[0]) / sc, (y as f64 + 0.5 - buf.offset[1]) / sc];
+            let q = sol.to_image(p);
+            if q[0] < 0.0 || q[1] < 0.0 || q[0] > size[0] || q[1] > size[1] {
+                continue;
+            }
+            *px = buf.img.sample_bilinear_clamped(q[0] * sc + buf.offset[0], q[1] * sc + buf.offset[1]);
+        }
+    });
+    Buf { img: out, offset: buf.offset, scale: buf.scale }
+}
+
 /// Median of the depths.
 pub fn median_depth(v: &[(u32, [f64; 2], f64)]) -> f64 {
     let mut z: Vec<f64> = v.iter().map(|x| x.2).collect();
@@ -188,15 +228,23 @@ pub fn median_depth(v: &[(u32, [f64; 2], f64)]) -> f64 {
 }
 
 fn render(ctx: &EffectCtx, mut buf: Buf) -> Buf {
-    if !b(ctx.params, "renderTrackPoints") {
+    let points = b(ctx.params, "renderTrackPoints");
+    if !points && !b(ctx.params, "advanced/undistort") {
         return buf;
     }
     let Some(sol) = solve(ctx.params) else { return buf };
     if (sol.size[0] - ctx.layer_size[0]).abs() > 0.5 || (sol.size[1] - ctx.layer_size[1]).abs() > 0.5 {
         return buf;
     }
+    let und = undistorts(ctx.params, &sol);
+    if und {
+        buf = undistort(&buf, &sol);
+    }
+    if !points {
+        return buf;
+    }
     let Some(k) = sol.frame_at(ctx.time) else { return buf };
-    let pts = projected(&sol, k);
+    let pts = projected_with(&sol, k, und);
     let med = median_depth(&pts);
     let size = f(ctx.params, "trackPointSize");
     let sc = if buf.scale > 0.0 { buf.scale } else { 1.0 };
@@ -287,5 +335,37 @@ mod tests {
         let st = settings(&Params::default());
         assert_eq!(st.shot, ShotType::FixedAngle);
         assert!(solve_key(&st).contains("FixedAngle"));
+        assert!(!st.lens_distortion && !solve_key(&st).contains("lens"));
+    }
+
+    #[test]
+    fn undistort_footage_remaps_through_the_solved_lens() {
+        use effectcraft_track::camtrack::Distortion;
+        let sol = CameraSolve {
+            version: 1,
+            size: [64.0, 48.0],
+            start: 0.0,
+            frame_duration: 1.0 / 30.0,
+            frames: vec![SolvedFrame { rot: [0.0; 3], center: [0.0; 3], focal: 50.0, solved: true }],
+            distortion: Some(Distortion { k1: -0.2, k2: 0.0, radius: 40.0 }),
+            ..Default::default()
+        };
+        let mut params = params_with(&sol, false);
+        params.values.insert("advanced/undistort".into(), Value::Bool(true));
+        assert!(undistorts(&params, &sol));
+        // A horizontal ramp: undistortion keeps the centre and pulls edge content outwards.
+        let mut img = Image::new(64, 48);
+        for (i, px) in img.data.iter_mut().enumerate() {
+            let x = (i % 64) as f32 / 64.0;
+            *px = [x, x, x, 1.0];
+        }
+        let ctx = EffectCtx { params: &params, time: 0.0, layer_size: [64.0, 48.0], seed: 0, adjustment: false, env: EffectEnv::default() };
+        let out = render(&ctx, Buf { img: img.clone(), offset: [0.0; 2], scale: 1.0 });
+        let at = |im: &Image, x: u32, y: u32| im.data[im.idx(x, y)][0];
+        assert!((at(&out.img, 32, 24) - at(&img, 32, 24)).abs() < 0.02);
+        // Barrel distortion compresses the edges: the ideal pixel at x = 56 shows footage nearer the centre.
+        assert!(at(&out.img, 56, 24) < at(&img, 56, 24) - 0.01, "{} vs {}", at(&out.img, 56, 24), at(&img, 56, 24));
+        let sp = sol.to_image([56.0, 24.0]);
+        assert!(sp[0] < 56.0 && (sol.to_ideal(sp)[0] - 56.0).abs() < 1e-6);
     }
 }
diff --git a/crates/effects/src/face_track.rs b/crates/effects/src/face_track.rs
new file mode 100644
index 00000000..7b925165
--- /dev/null
+++ b/crates/effects/src/face_track.rs
@@ -0,0 +1,45 @@
+//! Face tracking results (Tracker panel ▸ Method ▸ Face Tracking (Detailed Features)):
+//!
+//! - **Face Track Points** holds one point control per facial landmark (eyebrows, eyes and
+//!   pupils, nose, nostrils, mouth, chin and jaw), keyed on every tracked frame by the engine's
+//!   mask tracker ([`effectcraft_track::face`]); expressions and other layers can follow them.
+//! - **Face Measurements** holds the measurements derived from those points (head position,
+//!   scale and orientation, eye openness, eyebrow raise, mouth openness, width and offset), keyed
+//!   by Extract & Copy Face Measurements (`track.extractFaceMeasurements`).
+//!
+//! Both pass their input through unchanged: they are data, like expression controls.
+
+use effectcraft_keyframe::Value;
+use effectcraft_project::ParamUi;
+use effectcraft_track::face::{LANDMARKS, MEASUREMENTS};
+
+use crate::{Buf, EffectCtx, EffectSpec, ParamSpec, num, p, slider};
+
+/// Effect ids.
+pub const POINTS_ID: &str = "ec.utility.facetrackpoints";
+pub const MEASUREMENTS_ID: &str = "ec.utility.facemeasurements";
+
+fn passthrough(_: &EffectCtx, b: Buf) -> Buf {
+    b
+}
+
+fn measurement_ui(id: &str) -> ParamUi {
+    match id {
+        "headPositionX" | "headPositionY" => slider(-100_000.0, 100_000.0, 0.0, 1920.0, 1),
+        "headOrientationX" | "headOrientationY" | "headOrientationZ" => slider(-180.0, 180.0, -90.0, 90.0, 1),
+        _ => slider(-100_000.0, 100_000.0, 0.0, 200.0, 1),
+    }
+}
+
+pub fn specs() -> Vec<EffectSpec> {
+    let points: Vec<ParamSpec> = LANDMARKS
+        .iter()
+        // Point defaults are fractions of the layer size: the canonical face in the middle.
+        .map(|(id, name, q)| p(id, name, Value::Vec2([0.5 + 0.25 * q[0], 0.5 + 0.25 * q[1]]), ParamUi::Point))
+        .collect();
+    let measurements: Vec<ParamSpec> = MEASUREMENTS.iter().map(|(id, name)| p(id, name, num(0.0), measurement_ui(id))).collect();
+    vec![
+        EffectSpec { id: POINTS_ID, name: "Face Track Points", category: "Utility", params: points, render: passthrough, gpu: false, float: true },
+        EffectSpec { id: MEASUREMENTS_ID, name: "Face Measurements", category: "Utility", params: measurements, render: passthrough, gpu: false, float: true },
+    ]
+}
diff --git a/crates/effects/src/lib.rs b/crates/effects/src/lib.rs
index c4a47508..febdd9da 100644
--- a/crates/effects/src/lib.rs
+++ b/crates/effects/src/lib.rs
@@ -22,6 +22,7 @@ mod distort;
 mod distort2;
 mod distort3;
 pub mod distort4;
+pub mod face_track;
 mod generate;
 mod generate2;
 mod generate3;
@@ -412,6 +413,7 @@ pub fn registry() -> &'static [EffectSpec] {
         v.extend(puppet::specs());
         v.extend(warp_stab::specs());
         v.extend(camera_tracker::specs());
+        v.extend(face_track::specs());
         v.extend(roto::specs());
         v.sort_by(|a, b| a.category.cmp(b.category).then(a.name.cmp(b.name)));
         for s in v.iter_mut() {
diff --git a/crates/effects/src/warp_stab.rs b/crates/effects/src/warp_stab.rs
index c0bbced8..ff535c41 100644
--- a/crates/effects/src/warp_stab.rs
+++ b/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 {
diff --git a/crates/engine/src/commands/camera_cmds.rs b/crates/engine/src/commands/camera_cmds.rs
index 5d454aed..68cd14fb 100644
--- a/crates/engine/src/commands/camera_cmds.rs
+++ b/crates/engine/src/commands/camera_cmds.rs
@@ -23,7 +23,7 @@ pub fn specs() -> Vec<CommandSpec> {
             "Track Camera",
             ["Animation"],
             None,
-            "{layer?, shotType?: fixed|variable|specify, aov?: degrees, solveMethod?: auto|typical|flat|tripod, detailed?, wait?}",
+            "{layer?, shotType?: fixed|variable|specify, aov?: degrees, solveMethod?: auto|typical|flat|tripod, detailed?, lensDistortion?: solve k1/k2, undistort?: render undistorted, wait?}",
             can_track,
             track_camera
         ),
@@ -140,6 +140,8 @@ fn track_camera(s: &mut Session, p: &Value) -> Result<Value> {
         return Err(bad("track.camera", "aov: 1..170 degrees"));
     }
     let detailed = b_p(p, "detailed");
+    let lens = b_p(p, "lensDistortion");
+    let undistort = b_p(p, "undistort");
     let existing = s.project.comp(cid).and_then(|c| c.layer(lid)).and_then(Layer::effects).and_then(|fx| fx.groups().find(|g| is_tracker(g)).map(|g| g.uid));
     let uid = match existing {
         Some(u) => u,
@@ -148,7 +150,7 @@ fn track_camera(s: &mut Session, p: &Value) -> Result<Value> {
             r["effects"].get(0).and_then(Value::as_u64).ok_or_else(|| EngineError::Other("the effect could not be applied".into()))?
         }
     };
-    if shot.is_some() || method.is_some() || aov.is_some() || detailed.is_some() {
+    if shot.is_some() || method.is_some() || aov.is_some() || detailed.is_some() || lens.is_some() || undistort.is_some() {
         s.edit("3D Camera Tracker Settings", None, |proj, _| {
             let g = proj
                 .comp_mut(cid)
@@ -176,6 +178,12 @@ fn track_camera(s: &mut Session, p: &Value) -> Result<Value> {
             if let Some(d) = detailed {
                 set("advanced/detailedAnalysis", KV::Bool(d));
             }
+            if let Some(d) = lens {
+                set("advanced/lensDistortion", KV::Bool(d));
+            }
+            if let Some(d) = undistort {
+                set("advanced/undistort", KV::Bool(d));
+            }
             Ok(())
         })?;
     }
@@ -222,6 +230,7 @@ fn status(s: &mut Session, p: &Value) -> Result<Value> {
     out["shotType"] = json!(settings.shot.label());
     out["solveMethod"] = json!(settings.method.label());
     out["deleted"] = json!(settings.deleted.len());
+    out["lensDistortionSolved"] = json!(settings.lens_distortion);
     let solve = ct::solve(&params);
     out["solved"] = json!(solve.is_some());
     if let Some(sv) = &solve {
@@ -233,6 +242,11 @@ fn status(s: &mut Session, p: &Value) -> Result<Value> {
         out["focalLength"] = json!(f);
         out["horizontalAngleOfView"] = json!(sv.hfov(f));
         out["groundPlane"] = json!(sv.ground.is_some());
+        out["lensDistortion"] = match &sv.distortion {
+            Some(d) => json!({"k1": d.k1, "k2": d.k2, "radius": d.radius}),
+            None => Value::Null,
+        };
+        out["undistort"] = json!(ct::undistorts(&params, sv));
         if let Some(pl) = placed(&s.project, cid, lid, uid) {
             let k = pl.frame_at(s.time_of(cid));
             let (pos, ori, zoom) = pl.camera(k);
diff --git a/crates/engine/src/commands/mask_interp.rs b/crates/engine/src/commands/mask_interp.rs
index 5df40c6d..a4548944 100644
--- a/crates/engine/src/commands/mask_interp.rs
+++ b/crates/engine/src/commands/mask_interp.rs
@@ -1,5 +1,6 @@
-//! Mask tracking (Layer ▸ Mask ▸ Track Mask, the Tracker panel in mask mode) and Smart Mask
-//! Interpolation (Window ▸ Mask Interpolation).
+//! Mask tracking (Layer ▸ Mask ▸ Track Mask, the Tracker panel in mask mode, including the Face
+//! Tracking methods and Extract & Copy Face Measurements) and Smart Mask Interpolation
+//! (Window ▸ Mask Interpolation).
 
 use effectcraft_keyframe::{Keyframe, ShapePath, Value as KV};
 use effectcraft_path::interp::{self, AddVertices, InterpOpts, Matching};
@@ -20,7 +21,7 @@ pub fn specs() -> Vec<CommandSpec> {
             "Track Mask",
             ["Animation"],
             None,
-            "{layer?, mask?: uid|index|name, method?: position|positionScale|positionScaleRotation|positionScaleRotationSkew|perspective, direction?: forward|backward|frameForward|frameBackward, start? (s), end? (s), wait?: block until done}",
+            "{layer?, mask?: uid|index|name, method?: position|positionScale|positionScaleRotation|positionScaleRotationSkew|perspective|faceOutline|faceDetailed, direction?: forward|backward|frameForward|frameBackward, start? (s), end? (s), wait?: block until done}",
             has_mask,
             track_mask
         ),
@@ -29,10 +30,19 @@ pub fn specs() -> Vec<CommandSpec> {
             "Mask Tracking Method",
             [],
             None,
-            "{method: position|positionScale|positionScaleRotation|positionScaleRotationSkew|perspective}",
+            "{method: position|positionScale|positionScaleRotation|positionScaleRotationSkew|perspective|faceOutline|faceDetailed}",
             super::always,
             mask_method
         ),
+        cmd!(
+            "track.extractFaceMeasurements",
+            "Extract & Copy Face Measurements",
+            [],
+            None,
+            "{layer?} → keys a Face Measurements effect from the layer's Face Track Points (one key per tracked frame) and copies those keys",
+            has_face_points,
+            extract_face
+        ),
         cmd!(
             "mask.interpolate",
             "Apply Mask Interpolation",
@@ -116,9 +126,12 @@ fn has_mask(s: &Session) -> std::result::Result<(), String> {
 
 fn method_p(p: &Value, cmd: &str) -> Result<Option<MaskMethod>> {
     match str_p(p, "method") {
-        Some(m) => MaskMethod::from_name(m)
-            .map(Some)
-            .ok_or_else(|| bad(cmd, format!("unknown method `{m}` (position|positionScale|positionScaleRotation|positionScaleRotationSkew|perspective)"))),
+        Some(m) => MaskMethod::from_name(m).map(Some).ok_or_else(|| {
+            bad(
+                cmd,
+                format!("unknown method `{m}` (position|positionScale|positionScaleRotation|positionScaleRotationSkew|perspective|faceOutline|faceDetailed)"),
+            )
+        }),
         None => Ok(None),
     }
 }
@@ -190,6 +203,112 @@ fn track_mask(s: &mut Session, p: &Value) -> Result<Value> {
     }))
 }
 
+// ---------------------------------------------------------------- face measurements
+
+fn face_points_layer(s: &Session, p: &Value) -> Option<(effectcraft_project::ItemId, LayerId)> {
+    use effectcraft_effects::face_track::POINTS_ID;
+    let has = |l: &Layer| l.effects().is_some_and(|fx| fx.groups().any(|g| g.match_id == POINTS_ID));
+    if p.get("layer").is_some() {
+        return layer_p(s, p, "track.extractFaceMeasurements").ok();
+    }
+    let cid = s.active_comp_id()?;
+    let comp = s.project.comp(cid)?;
+    let sel = s.state.selected_layers.iter().filter_map(|l| comp.layer(*l)).find(|l| has(l));
+    sel.or_else(|| comp.layers.iter().find(|l| has(l))).map(|l| (cid, l.id))
+}
+
+fn has_face_points(s: &Session) -> std::result::Result<(), String> {
+    has_comp(s)?;
+    face_points_layer(s, &json!({})).map(|_| ()).ok_or_else(|| "no Face Track Points: track a mask with Face Tracking (Detailed Features) first".into())
+}
+
+/// Extract & Copy Face Measurements: measure the Face Track Points at every keyed time, key them
+/// into the layer's Face Measurements effect and put those keys on the keyframe clipboard.
+fn extract_face(s: &mut Session, p: &Value) -> Result<Value> {
+    use effectcraft_effects::face_track::{MEASUREMENTS_ID, POINTS_ID};
+    use effectcraft_track::face::{LANDMARKS, MEASUREMENTS, N, face_frame, measure};
+    const CMD: &str = "track.extractFaceMeasurements";
+    let (cid, lid) = face_points_layer(s, p).ok_or_else(|| bad(CMD, "no Face Track Points: track a mask with Face Tracking (Detailed Features) first"))?;
+    let comp = s.project.comp(cid).ok_or(EngineError::NoComp)?;
+    let layer = comp.layer(lid).ok_or(EngineError::NoComp)?;
+    let g = layer.effects().and_then(|fx| fx.groups().find(|g| g.match_id == POINTS_ID)).ok_or_else(|| bad(CMD, "the layer has no Face Track Points"))?;
+    let mut times: Vec<Tick> = LANDMARKS.iter().filter_map(|l| g.get(l.0)).flat_map(|pr| pr.keys.iter().map(|k| k.time)).collect();
+    times.sort();
+    times.dedup();
+    if times.is_empty() {
+        times.push(layer.layer_time(s.time_of(cid)));
+    }
+    let rows: Vec<(Tick, [f64; 14])> = {
+        let pts_at = |t: Tick| {
+            let mut l = [[0.0; 2]; N];
+            for (i, (id, _, _)) in LANDMARKS.iter().enumerate() {
+                if let Some(KV::Vec2(v)) = g.get(id).map(|pr| pr.value_at(t)) {
+                    l[i] = v;
+                }
+            }
+            l
+        };
+        let reference = face_frame(&pts_at(times[0])).1;
+        times.iter().map(|t| (*t, measure(&pts_at(*t), reference))).collect()
+    };
+    let comp_times: Vec<Tick> = times.iter().map(|t| layer.comp_time(*t)).collect();
+    let fd = comp.frame_duration();
+    let uid = s.edit("Extract Face Measurements", None, |proj, _| {
+        let uid = crate::mask_track::face_effect(proj, cid, lid, MEASUREMENTS_ID).ok_or_else(|| bad(CMD, "Face Measurements is not available"))?;
+        let g = proj.comp_mut(cid).and_then(|c| c.layer_mut(lid)).and_then(|l| l.props.find_group_mut(uid)).ok_or(EngineError::NoComp)?;
+        for (k, (id, _)) in MEASUREMENTS.iter().enumerate() {
+            let Some(pr) = g.get_mut(id) else { continue };
+            pr.keys.clear();
+            for (t, m) in &rows {
+                effectcraft_keyframe::set_key(&mut pr.keys, Keyframe::new(*t, KV::Scalar(m[k])));
+            }
+        }
+        Ok(uid)
+    })?;
+    // Copy the keys (Edit ▸ Paste pastes them at the CTI, like copied keyframes).
+    let comp = s.project.comp(cid).ok_or(EngineError::NoComp)?;
+    let layer = comp.layer(lid).ok_or(EngineError::NoComp)?;
+    let mg = layer.props.find_group(uid).ok_or(EngineError::NoComp)?;
+    let earliest = comp_times.iter().copied().min().unwrap_or(Tick::ZERO);
+    let mut clip = vec![];
+    for (id, _) in MEASUREMENTS {
+        let Some(pr) = mg.get(id) else { continue };
+        let Some(path) = layer.props.match_path_of(pr.uid) else { continue };
+        let keys = pr
+            .keys
+            .iter()
+            .map(|k| {
+                let mut k = k.clone();
+                k.time = layer.comp_time(k.time) - earliest;
+                k
+            })
+            .collect();
+        clip.push(crate::KeyClip { layer: lid, path, keys });
+    }
+    s.state.key_clipboard = clip;
+    s.state.clip_is_keys = true;
+    // The same data as text (tab-separated, one row per frame) for other applications.
+    let mut text = String::from("EffectCraft Face Measurements\n\nFrame\tTime");
+    for (_, name) in MEASUREMENTS {
+        text.push('\t');
+        text.push_str(name);
+    }
+    text.push('\n');
+    let mut table = vec![];
+    for ((_, m), ct) in rows.iter().zip(&comp_times) {
+        let frame = if fd > Tick::ZERO { (ct.seconds() / fd.seconds()).round() as i64 } else { 0 };
+        text.push_str(&format!("{frame}\t{:.4}", ct.seconds()));
+        for v in m {
+            text.push_str(&format!("\t{v:.3}"));
+        }
+        text.push('\n');
+        let obj: serde_json::Map<String, Value> = MEASUREMENTS.iter().zip(m).map(|((id, _), v)| (id.to_string(), json!(v))).collect();
+        table.push(json!({"time": ct.seconds(), "frame": frame, "values": obj}));
+    }
+    s.bump();
+    Ok(json!({"layer": lid.0, "effect": uid, "frames": rows.len(), "measurements": table, "clipboard": text}))
+}
+
 // ---------------------------------------------------------------- mask interpolation
 
 /// Mask Interpolation panel options (kept in the editor state, serde for agents).
diff --git a/crates/engine/src/lib.rs b/crates/engine/src/lib.rs
index 17d9b235..0436f8fa 100644
--- a/crates/engine/src/lib.rs
+++ b/crates/engine/src/lib.rs
@@ -774,6 +774,8 @@ mod tests_effects;
 #[cfg(test)]
 mod tests_essential;
 #[cfg(test)]
+mod tests_face;
+#[cfg(test)]
 mod tests_fidelity;
 #[cfg(test)]
 mod tests_lottie;
diff --git a/crates/engine/src/mask_track.rs b/crates/engine/src/mask_track.rs
index 00315198..a5741052 100644
--- a/crates/engine/src/mask_track.rs
+++ b/crates/engine/src/mask_track.rs
@@ -6,6 +6,12 @@
 //! algorithm is `effectcraft_track::mask` (KLT features inside the mask + a RANSAC fit of the
 //! chosen motion model), and every frame's motion is applied to the mask's vertices and
 //! tangents, so the shape's Bezier structure is kept.
+//!
+//! The two **Face Tracking** methods use `effectcraft_track::face` instead: the mask (drawn
+//! around a face) seeds the tracker, and every frame's face outline becomes the Mask Path
+//! (*Outline Only*); *Detailed Features* also keys the facial landmarks into a **Face Track
+//! Points** effect on the layer (created on first use), from which Extract & Copy Face
+//! Measurements (`track.extractFaceMeasurements`) derives a keyed **Face Measurements** effect.
 
 use std::sync::atomic::{AtomicBool, Ordering};
 use std::sync::{Arc, Mutex};
@@ -31,11 +37,22 @@ pub enum MaskMethod {
     PositionScaleRotation,
     PositionScaleRotationSkew,
     Perspective,
+    /// Face Tracking (Outline Only).
+    FaceOutline,
+    /// Face Tracking (Detailed Features).
+    FaceDetailed,
 }
 
 impl MaskMethod {
-    pub const ALL: [MaskMethod; 5] =
-        [MaskMethod::Position, MaskMethod::PositionScale, MaskMethod::PositionScaleRotation, MaskMethod::PositionScaleRotationSkew, MaskMethod::Perspective];
+    pub const ALL: [MaskMethod; 7] = [
+        MaskMethod::Position,
+        MaskMethod::PositionScale,
+        MaskMethod::PositionScaleRotation,
+        MaskMethod::PositionScaleRotationSkew,
+        MaskMethod::Perspective,
+        MaskMethod::FaceOutline,
+        MaskMethod::FaceDetailed,
+    ];
     pub fn label(self) -> &'static str {
         match self {
             MaskMethod::Position => "Position",
@@ -43,8 +60,14 @@ impl MaskMethod {
             MaskMethod::PositionScaleRotation => "Position, Scale & Rotation",
             MaskMethod::PositionScaleRotationSkew => "Position, Scale, Rotation & Skew",
             MaskMethod::Perspective => "Perspective",
+            MaskMethod::FaceOutline => "Face Tracking (Outline Only)",
+            MaskMethod::FaceDetailed => "Face Tracking (Detailed Features)",
         }
     }
+    /// One of the Face Tracking methods.
+    pub fn is_face(self) -> bool {
+        matches!(self, MaskMethod::FaceOutline | MaskMethod::FaceDetailed)
+    }
     pub fn id(self) -> &'static str {
         match self {
             MaskMethod::Position => "position",
@@ -52,6 +75,8 @@ impl MaskMethod {
             MaskMethod::PositionScaleRotation => "positionScaleRotation",
             MaskMethod::PositionScaleRotationSkew => "positionScaleRotationSkew",
             MaskMethod::Perspective => "perspective",
+            MaskMethod::FaceOutline => "faceOutline",
+            MaskMethod::FaceDetailed => "faceDetailed",
         }
     }
     pub fn from_name(s: &str) -> Option<MaskMethod> {
@@ -67,6 +92,7 @@ impl MaskMethod {
             MaskMethod::PositionScaleRotation => Model::Similarity,
             MaskMethod::PositionScaleRotationSkew => Model::Affine,
             MaskMethod::Perspective => Model::Homography,
+            MaskMethod::FaceOutline | MaskMethod::FaceDetailed => Model::Similarity,
         }
     }
 }
@@ -109,8 +135,29 @@ pub fn transform_path(p: &ShapePath, h: &trk::Homography) -> ShapePath {
     }
 }
 
-/// One tracked frame: (layer time, comp time, mask path).
-type MaskFrame = (Tick, Tick, ShapePath);
+/// A closed, smooth mask path through a face outline polygon: `n` vertices with Catmull-Rom
+/// tangents. `like` lends its feather points.
+pub fn outline_path(poly: &[[f64; 2]], n: usize, like: &ShapePath) -> ShapePath {
+    let m = poly.len();
+    let n = n.clamp(3, m.max(3));
+    let v: Vec<[f64; 2]> = (0..n).map(|i| poly[(i * m / n) % m.max(1)]).collect();
+    let tan: Vec<[f64; 2]> = (0..n)
+        .map(|i| {
+            let (a, b) = (v[(i + n - 1) % n], v[(i + 1) % n]);
+            [(b[0] - a[0]) / 6.0, (b[1] - a[1]) / 6.0]
+        })
+        .collect();
+    ShapePath {
+        vertices: v,
+        in_tangents: tan.iter().map(|t| [-t[0], -t[1]]).collect(),
+        out_tangents: tan,
+        closed: true,
+        feather: if like.len() == n { like.feather.clone() } else { vec![] },
+    }
+}
+
+/// One tracked frame: (layer time, comp time, mask path, face landmarks for Detailed Features).
+type MaskFrame = (Tick, Tick, ShapePath, Option<[[f64; 2]; trk::face::N]>);
 
 #[derive(Default)]
 pub struct MaskTrackShared {
@@ -169,9 +216,12 @@ fn run_work(w: MaskWork, shared: &MaskTrackShared) {
     let frame_at = |t: Tick| source_frame(r, &w.project, w.comp, comp, layer, t, w.expr.as_deref());
     lock(&shared.state).total = w.times.len().saturating_sub(1) as u64;
     let Some(first) = frame_at(w.times[0]) else { return fail("the layer has no pixels to track") };
+    if w.method.is_face() {
+        return run_face(&w, shared, first, &frame_at, layer, t0);
+    }
     let mut tracker = trk::mask::MaskTracker::new(w.method.model(), flatten(&w.path, 8), &trk::Frame { img: &first.0, offset: first.1 });
     let mut path = w.path.clone();
-    lock(&shared.frames).push((layer.layer_time(w.times[0]), w.times[0], path.clone()));
+    lock(&shared.frames).push((layer.layer_time(w.times[0]), w.times[0], path.clone(), None));
     let mut next = if w.times.len() > 1 { frame_at(w.times[1]) } else { None };
     for k in 1..w.times.len() {
         if shared.cancel.load(Ordering::Relaxed) {
@@ -190,7 +240,7 @@ fn run_work(w: MaskWork, shared: &MaskTrackShared) {
             break;
         };
         path = transform_path(&path, &step.motion);
-        lock(&shared.frames).push((layer.layer_time(w.times[k]), w.times[k], path.clone()));
+        lock(&shared.frames).push((layer.layer_time(w.times[k]), w.times[k], path.clone(), None));
         let mut s = lock(&shared.state);
         s.done = k as u64;
         s.elapsed = t0.elapsed().as_secs_f64();
@@ -203,12 +253,99 @@ fn run_work(w: MaskWork, shared: &MaskTrackShared) {
     s.finished = true;
 }
 
+/// Face tracking: the face outline keys the mask; Detailed Features adds the landmarks.
+fn run_face(
+    w: &MaskWork,
+    shared: &MaskTrackShared,
+    first: (Arc<effectcraft_raster::Image>, [f64; 2]),
+    frame_at: &(dyn Fn(Tick) -> Option<(Arc<effectcraft_raster::Image>, [f64; 2])> + Sync),
+    layer: &effectcraft_project::Layer,
+    t0: web_time::Instant,
+) {
+    let detailed = w.method == MaskMethod::FaceDetailed;
+    let b = trk::mask::bounds(&flatten(&w.path, 8));
+    let Some((mut tracker, fit)) = trk::face::FaceTracker::new(&trk::Frame { img: &first.0, offset: first.1 }, b) else {
+        let mut s = lock(&shared.state);
+        s.error = Some("no face found inside the mask (draw the mask around a face)".into());
+        s.finished = true;
+        return;
+    };
+    let key = |fit: &trk::face::FaceFit| (outline_path(&fit.outline, 16, &w.path), detailed.then_some(fit.landmarks));
+    let (path, pts) = key(&fit);
+    lock(&shared.frames).push((layer.layer_time(w.times[0]), w.times[0], path, pts));
+    let mut next = if w.times.len() > 1 { frame_at(w.times[1]) } else { None };
+    for k in 1..w.times.len() {
+        if shared.cancel.load(Ordering::Relaxed) {
+            lock(&shared.state).cancelled = true;
+            break;
+        }
+        let Some(cur) = next.take() else {
+            lock(&shared.state).error = Some(format!("no frame at {:.3} s", w.times[k].seconds()));
+            break;
+        };
+        let (res, nf) =
+            rayon::join(|| tracker.step(&trk::Frame { img: &cur.0, offset: cur.1 }), || if k + 1 < w.times.len() { frame_at(w.times[k + 1]) } else { None });
+        next = nf;
+        let Some(fit) = res else {
+            lock(&shared.state).error = Some(format!("lost the face at {:.3} s", w.times[k].seconds()));
+            break;
+        };
+        let (path, pts) = key(&fit);
+        lock(&shared.frames).push((layer.layer_time(w.times[k]), w.times[k], path, pts));
+        let mut s = lock(&shared.state);
+        s.done = k as u64;
+        s.elapsed = t0.elapsed().as_secs_f64();
+        s.fps = if s.elapsed > 0.0 { k as f64 / s.elapsed } else { 0.0 };
+        s.time = w.times[k].seconds();
+        crate::offload::report(s.done, s.total);
+    }
+    let mut s = lock(&shared.state);
+    s.elapsed = t0.elapsed().as_secs_f64();
+    s.finished = true;
+}
+
+/// The layer's Face Track Points effect (created when missing). Returns its uid.
+pub(crate) fn face_points_group(p: &mut Project, comp: ItemId, layer: LayerId) -> Option<Uid> {
+    face_effect(p, comp, layer, effectcraft_effects::face_track::POINTS_ID)
+}
+
+/// The layer's effect `id`, applied when missing (data effects of face tracking).
+pub(crate) fn face_effect(p: &mut Project, comp: ItemId, layer: LayerId, id: &str) -> Option<Uid> {
+    let spec = effectcraft_effects::find(id)?;
+    let l = p.comp(comp)?.layer(layer)?;
+    if let Some(g) = l.effects().and_then(|fx| fx.groups().find(|g| g.match_id == spec.id)) {
+        return Some(g.uid);
+    }
+    let (w, h) = effectcraft_render::source_size(p, l);
+    let size = if w == 0 { p.comp(comp).map(|c| [c.width as f64, c.height as f64]).unwrap_or([1920.0, 1080.0]) } else { [w as f64, h as f64] };
+    let mut next = p.next_id;
+    let g = effectcraft_effects::instantiate(spec, &mut effectcraft_project::build::Ids(&mut next), spec.name, size);
+    let uid = g.uid;
+    let fx = p.comp_mut(comp)?.layer_mut(layer)?.props.sub_mut("effects")?;
+    fx.children.push(g.into());
+    p.next_id = next;
+    Some(uid)
+}
+
 fn write_frames(p: &mut Project, job: &MaskTrackJob, frames: &[MaskFrame]) {
+    let points = if frames.iter().any(|f| f.3.is_some()) { face_points_group(p, job.comp, job.layer) } else { None };
     let Some(layer) = p.comp_mut(job.comp).and_then(|c| c.layer_mut(job.layer)) else { return };
-    let Some(g) = layer.props.find_group_mut(job.mask) else { return };
-    let Some(pr) = g.get_mut("path") else { return };
-    for (lt, _, path) in frames {
-        effectcraft_keyframe::set_key(&mut pr.keys, Keyframe::new(*lt, Value::Path(path.clone())));
+    if let Some(g) = layer.props.find_group_mut(job.mask)
+        && let Some(pr) = g.get_mut("path")
+    {
+        for (lt, _, path, _) in frames {
+            effectcraft_keyframe::set_key(&mut pr.keys, Keyframe::new(*lt, Value::Path(path.clone())));
+        }
+    }
+    if let Some(g) = points.and_then(|u| layer.props.find_group_mut(u)) {
+        for (lt, _, _, pts) in frames {
+            let Some(pts) = pts else { continue };
+            for (i, (id, _, _)) in trk::face::LANDMARKS.iter().enumerate() {
+                if let Some(pr) = g.get_mut(id) {
+                    effectcraft_keyframe::set_key(&mut pr.keys, Keyframe::new(*lt, Value::Vec2(pts[i])));
+                }
+            }
+        }
     }
 }
 
@@ -242,7 +379,7 @@ impl Session {
             return Err("a track analysis is already running".into());
         }
         self.poll_mask_track();
-        if !wait && self.offloads() {
+        if !wait && self.offloads() && !work.method.is_face() {
             let MaskWork { comp, layer, path, method, times, .. } = work;
             return self.offload_analysis(crate::offload::WorkerJob::MaskTrack { comp, layer, mask, direction, path, method, times });
         }
diff --git a/crates/engine/src/tests_face.rs b/crates/engine/src/tests_face.rs
new file mode 100644
index 00000000..dad97935
--- /dev/null
+++ b/crates/engine/src/tests_face.rs
@@ -0,0 +1,106 @@
+//! Face tracking end to end through the commands, on synthetic face footage generated in-test
+//! (a moving, rotating face that opens its mouth and blinks).
+
+use effectcraft_keyframe::{ShapePath, Value as KV};
+use effectcraft_track::face::{FaceParams, LANDMARKS, MEASUREMENTS, render_synthetic, synth_landmarks};
+use serde_json::json;
+
+use crate::mask_track::MaskMethod;
+use crate::tests_track::{H, W, frame_time, setup};
+
+fn face(f: u32) -> FaceParams {
+    let t = f as f64;
+    FaceParams {
+        center: [160.0 + 25.0 * (t * 0.2).sin(), 120.0 + 8.0 * (t * 0.15).cos()],
+        height: 70.0 * (1.0 + 0.05 * (t * 0.1).sin()),
+        roll: 12.0 * (t * 0.18).sin(),
+        yaw: 10.0 * (t * 0.12).sin(),
+        mouth_open: 0.5 + 0.45 * (t * 0.3).sin(),
+        eye_open: [0.85, 0.85],
+        ..Default::default()
+    }
+}
+
+#[test]
+fn face_tracking_keys_outline_points_and_measurements() {
+    let (mut s, clip, _) = setup(|f| render_synthetic(W, H, &face(f), 5));
+    assert_eq!(MaskMethod::from_name("faceDetailed"), Some(MaskMethod::FaceDetailed));
+    assert_eq!(MaskMethod::from_name("Face Tracking (Outline Only)"), Some(MaskMethod::FaceOutline));
+    // A rough mask around the face.
+    let c = face(0).center;
+    let r = s.execute("mask.new", json!({"layer": clip.0, "vertices": ShapePath::rect(c, 120.0, 160.0).vertices, "closed": true})).unwrap();
+    let mask = r["mask"].as_u64().unwrap();
+    s.execute("time.set", json!({"time": 0})).unwrap();
+    let r = s.execute_checked("track.mask", json!({"layer": clip.0, "mask": mask, "method": "faceDetailed", "direction": "forward", "wait": true})).unwrap();
+    assert_eq!(r["frames"], 24, "{r}");
+    assert_eq!(s.state.mask_track_method, MaskMethod::FaceDetailed);
+    let comp = s.active_comp().unwrap();
+    let layer = comp.layer(clip).unwrap();
+    // The mask follows the face outline.
+    let path = layer.props.find_group(mask).unwrap().get("path").unwrap();
+    assert_eq!(path.keys.len(), 25);
+    for f in [0, 12, 24] {
+        let KV::Path(p) = path.value_at(frame_time(f)) else { panic!() };
+        let fp = face(f);
+        let (sn, cs) = fp.roll.to_radians().sin_cos();
+        for v in &p.vertices {
+            let d = [(v[0] - fp.center[0]) / fp.height, (v[1] - fp.center[1]) / fp.height];
+            let (u, w) = (cs * d[0] + sn * d[1], -sn * d[0] + cs * d[1]);
+            let r = ((u / 0.75).powi(2) + w * w).sqrt();
+            assert!((r - 1.0).abs() < 0.06, "frame {f}: vertex {v:?} at {r}");
+        }
+    }
+    // Face Track Points: one keyed point per landmark, on the face.
+    let fx = layer.effects().unwrap();
+    let pts = fx.groups().find(|g| g.match_id == effectcraft_effects::face_track::POINTS_ID).expect("Face Track Points");
+    let mut worst: f64 = 0.0;
+    for f in [0, 6, 12, 18, 24] {
+        let truth = synth_landmarks(&face(f));
+        for (i, (id, _, _)) in LANDMARKS.iter().enumerate() {
+            let pr = pts.get(id).unwrap();
+            assert_eq!(pr.keys.len(), 25, "{id}");
+            let KV::Vec2(v) = pr.value_at(frame_time(f)) else { panic!() };
+            worst = worst.max((v[0] - truth[i][0]).hypot(v[1] - truth[i][1]) / face(f).height);
+        }
+    }
+    assert!(worst < 0.1, "worst landmark error {worst} face heights");
+    // Extract & Copy Face Measurements.
+    assert!(s.is_enabled("track.extractFaceMeasurements"));
+    let r = s.execute_checked("track.extractFaceMeasurements", json!({"layer": clip.0})).unwrap();
+    assert_eq!(r["frames"], 25);
+    assert!(r["clipboard"].as_str().unwrap().contains("Mouth Openness"));
+    assert_eq!(s.state.key_clipboard.len(), MEASUREMENTS.len());
+    assert!(s.state.clip_is_keys);
+    let layer = s.active_comp().unwrap().layer(clip).unwrap();
+    let m = layer.effects().unwrap().groups().find(|g| g.match_id == effectcraft_effects::face_track::MEASUREMENTS_ID).expect("Face Measurements");
+    let mo = m.get("mouthOpenness").unwrap();
+    assert_eq!(mo.keys.len(), 25);
+    // Mouth openness follows the footage: widest open where the generator opens it most.
+    let open = |f: u32| mo.value_at(frame_time(f)).as_f64();
+    let (fmax, fmin) = (
+        (0..25).max_by(|a, b| face(*a).mouth_open.total_cmp(&face(*b).mouth_open)).unwrap(),
+        (0..25).min_by(|a, b| face(*a).mouth_open.total_cmp(&face(*b).mouth_open)).unwrap(),
+    );
+    assert!(open(fmax) > open(fmin) + 20.0, "{} vs {}", open(fmax), open(fmin));
+    let roll = m.get("headOrientationZ").unwrap();
+    for f in [0, 10, 20] {
+        assert!((roll.value_at(frame_time(f)).as_f64() - face(f).roll).abs() < 3.0);
+    }
+    // One undo step for the extraction.
+    assert!(s.undo());
+    assert!(!s.active_comp().unwrap().layer(clip).unwrap().effects().unwrap().groups().any(|g| g.match_id == effectcraft_effects::face_track::MEASUREMENTS_ID));
+    // Outline Only keys just the mask.
+    let r = s.execute("mask.new", json!({"layer": clip.0, "vertices": ShapePath::rect(c, 120.0, 160.0).vertices, "closed": true})).unwrap();
+    let m2 = r["mask"].as_u64().unwrap();
+    s.execute("time.set", json!({"time": 0})).unwrap();
+    s.execute_checked("track.mask", json!({"layer": clip.0, "mask": m2, "method": "faceOutline", "direction": "frameForward", "wait": true})).unwrap();
+    let layer = s.active_comp().unwrap().layer(clip).unwrap();
+    assert_eq!(layer.props.find_group(m2).unwrap().get("path").unwrap().keys.len(), 2);
+    // A mask without a face reports it.
+    let r = s.execute("mask.new", json!({"layer": clip.0, "vertices": ShapePath::rect([20.0, 20.0], 30.0, 30.0).vertices, "closed": true})).unwrap();
+    let m3 = r["mask"].as_u64().unwrap();
+    s.execute("time.set", json!({"time": 0})).unwrap();
+    let _ = s.execute("track.mask", json!({"layer": clip.0, "mask": m3, "method": "faceOutline", "direction": "forward", "wait": true}));
+    let layer = s.active_comp().unwrap().layer(clip).unwrap();
+    assert!(layer.props.find_group(m3).unwrap().get("path").unwrap().keys.is_empty());
+}
diff --git a/crates/track/src/camtrack/bundle.rs b/crates/track/src/camtrack/bundle.rs
index 72cb2f17..077f04b7 100644
--- a/crates/track/src/camtrack/bundle.rs
+++ b/crates/track/src/camtrack/bundle.rs
@@ -10,6 +10,7 @@
 
 use rayon::prelude::*;
 
+use super::Distortion;
 use super::geometry::Pose;
 use super::linalg::*;
 
@@ -19,6 +20,8 @@ pub struct BaCam {
     pub pose: Pose,
     /// Focal length (pixels).
     pub f: f64,
+    /// Lens distortion (shared by all cameras; adjusted when [`BaOpts::distortion`]).
+    pub dist: Distortion,
 }
 
 /// One observation: camera `cam` sees point `pt` at `uv` (pixels relative to the principal
@@ -50,11 +53,13 @@ pub struct BaOpts {
     /// Huber threshold (pixels).
     pub huber: f64,
     pub max_iter: usize,
+    /// Adjust the shared radial distortion `(k1, k2)` too.
+    pub distortion: bool,
 }
 
 impl Default for BaOpts {
     fn default() -> Self {
-        BaOpts { focal: FocalMode::Fixed, fix_points: false, fix_centers: false, fixed: vec![], huber: 2.0, max_iter: 30 }
+        BaOpts { focal: FocalMode::Fixed, fix_points: false, fix_centers: false, fixed: vec![], huber: 2.0, max_iter: 30, distortion: false }
     }
 }
 
@@ -77,7 +82,8 @@ pub fn residual(c: &BaCam, x: V3, uv: [f64; 2]) -> Option<[f64; 2]> {
     if p[2] <= 1e-9 {
         return None;
     }
-    Some([c.f * p[0] / p[2] - uv[0], c.f * p[1] / p[2] - uv[1]])
+    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]])
 }
 
 /// Total robust cost.
@@ -108,6 +114,8 @@ struct Layout {
     fidx: Vec<usize>,
     /// Index of the shared focal parameter.
     shared: Option<usize>,
+    /// Index of `k1` (then `k2`) when the distortion is adjusted.
+    dist: Option<usize>,
     n: usize,
     rot_only: bool,
 }
@@ -132,7 +140,11 @@ fn layout(ncam: usize, o: &BaOpts) -> Layout {
         n += 1;
         n - 1
     });
-    Layout { start, fidx, shared, n, rot_only: o.fix_centers }
+    let dist = o.distortion.then(|| {
+        n += 2;
+        n - 2
+    });
+    Layout { start, fidx, shared, dist, n, rot_only: o.fix_centers }
 }
 
 /// Camera-side Jacobian entries of one observation: (parameter index, d r / d param).
@@ -204,12 +216,21 @@ pub fn adjust(cams: &mut [BaCam], pts: &mut [V3], obs: &[BaObs], o: &BaOpts) ->
                     continue;
                 }
                 let iz = 1.0 / xc[2];
-                let p = [c.f * xc[0] * iz, c.f * xc[1] * iz];
-                let r = [p[0] - ob.uv[0], p[1] - ob.uv[1]];
+                let p0 = [c.f * xc[0] * iz, c.f * xc[1] * iz];
+                let pd = c.dist.distort(p0);
+                let r = [pd[0] - ob.uv[0], pd[1] - ob.uv[1]];
                 let e = (r[0] * r[0] + r[1] * r[1]).sqrt();
                 let w = if e <= o.huber { 1.0 } else { o.huber / e };
-                let a = [[c.f * iz, 0.0, -c.f * xc[0] * iz * iz], [0.0, c.f * iz, -c.f * xc[1] * iz * iz]];
+                let a0 = [[c.f * iz, 0.0, -c.f * xc[0] * iz * iz], [0.0, c.f * iz, -c.f * xc[1] * iz * iz]];
+                // Through the distortion: d r = D d p.
+                let (dj, dk) = c.dist.jacobian(p0);
+                let a = [0, 1].map(|i| [0, 1, 2].map(|k| dj[i][0] * a0[0][k] + dj[i][1] * a0[1][k]));
+                let p = [dj[0][0] * p0[0] + dj[0][1] * p0[1], dj[1][0] * p0[0] + dj[1][1] * p0[1]];
                 cam_jac(&l, ci, c, xc, &a, p, &mut jc);
+                if let Some(g) = l.dist {
+                    jc.push((g, [dk[0][0], dk[1][0]]));
+                    jc.push((g + 1, [dk[0][1], dk[1][1]]));
+                }
                 // Point Jacobian: A R.
                 let jx = if o.fix_points {
                     [[0.0; 3]; 2]
@@ -236,7 +257,8 @@ pub fn adjust(cams: &mut [BaCam], pts: &mut [V3], obs: &[BaObs], o: &BaOpts) ->
                             [w * (ja[0] * jx[0][0] + ja[1] * jx[1][0]), w * (ja[0] * jx[0][1] + ja[1] * jx[1][1]), w * (ja[0] * jx[0][2] + ja[1] * jx[1][2])];
                         // Each observation of a point is on another camera: only the shared focal
                         // length repeats.
-                        let found = if Some(*ia) == l.shared { ps.w.iter_mut().find(|e| e.0 == *ia) } else { None };
+                        let repeats = Some(*ia) == l.shared || l.dist.is_some_and(|g| *ia == g || *ia == g + 1);
+                        let found = if repeats { ps.w.iter_mut().find(|e| e.0 == *ia) } else { None };
                         match found {
                             Some(e) => e.1 = add(e.1, wv),
                             None => ps.w.push((*ia, wv)),
@@ -385,6 +407,10 @@ pub fn adjust(cams: &mut [BaCam], pts: &mut [V3], obs: &[BaObs], o: &BaOpts) ->
                 if let Some(g) = l.shared {
                     c.f *= dc[g].clamp(-0.5, 0.5).exp();
                 }
+                if let Some(g) = l.dist {
+                    c.dist.k1 += dc[g].clamp(-0.2, 0.2);
+                    c.dist.k2 += dc[g + 1].clamp(-0.2, 0.2);
+                }
             }
             let mut np: Vec<V3> = pts.to_vec();
             if !o.fix_points {
@@ -445,7 +471,11 @@ mod tests {
         let truth: Vec<V3> = (0..80).map(|_| [u() * 6.0, u() * 4.0, 6.0 + u() * 4.0]).collect();
         let f = 900.0;
         let cams_true: Vec<BaCam> = (0..6)
-            .map(|i| BaCam { pose: Pose { r: rodrigues([0.0, -0.05 * i as f64, 0.01 * i as f64]), c: [0.3 * i as f64, 0.05 * i as f64, 0.0] }, f })
+            .map(|i| BaCam {
+                pose: Pose { r: rodrigues([0.0, -0.05 * i as f64, 0.01 * i as f64]), c: [0.3 * i as f64, 0.05 * i as f64, 0.0] },
+                f,
+                dist: Distortion::NONE,
+            })
             .collect();
         let mut obs = vec![];
         for (ci, c) in cams_true.iter().enumerate() {
diff --git a/crates/track/src/camtrack/features.rs b/crates/track/src/camtrack/features.rs
index 30b019e5..1775a304 100644
--- a/crates/track/src/camtrack/features.rs
+++ b/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) }
     }
 
diff --git a/crates/track/src/camtrack/mod.rs b/crates/track/src/camtrack/mod.rs
index 76c0637f..e8cd28da 100644
--- a/crates/track/src/camtrack/mod.rs
+++ b/crates/track/src/camtrack/mod.rs
@@ -28,8 +28,15 @@
 //!    tripod solution when it explains the tracks as well as the general one (a quick rotation fit
 //!    at the general solution's focal length skips the tripod search when it clearly fails).
 //!
-//! The camera model has square pixels, the principal point at the layer centre and no lens
-//! distortion. Camera space is x right, y down, z forward (After Effects' convention); a camera
+//! 5. **Lens distortion** (Solve Lens Distortion, or Detailed Analysis): radial distortion
+//!    `p_d = p (1 + k1 ρ² + k2 ρ⁴)` (ρ = distance from the principal point over half the frame
+//!    diagonal; Brown's model without the tangential terms) is solved jointly with the cameras
+//!    and points in a final bundle adjustment on the raw tracks; the whole solve is then repeated
+//!    on tracks undistorted with that estimate (the linear initialisations assume a pinhole
+//!    camera) and adjusted jointly again. The solve stores `k1, k2`; the effect can undistort
+//!    its output so composited 3D layers line up.
+//!
+//! The camera model has square pixels and the principal point at the layer centre. Camera space is x right, y down, z forward (After Effects' convention); a camera
 //! is `(R, C, f)` with `X_cam = R (X − C)` and pixel `u = f x / z + w / 2`.
 //!
 //! The solve is stored in a canonical frame: the first solved frame's camera is at the origin
@@ -105,6 +112,69 @@ impl CameraTracks {
     }
 }
 
+/// Radial lens distortion around the principal point (see the module docs).
+#[derive(Clone, Copy, Debug, Default, PartialEq, Serialize, Deserialize)]
+#[serde(rename_all = "camelCase", default)]
+pub struct Distortion {
+    pub k1: f64,
+    pub k2: f64,
+    /// Normalisation radius ρ = 1 (pixels): half the frame diagonal; 0 = no distortion.
+    pub radius: f64,
+}
+
+impl Distortion {
+    pub const NONE: Distortion = Distortion { k1: 0.0, k2: 0.0, radius: 0.0 };
+
+    /// No distortion for frames of `size` (the radius set, coefficients zero).
+    pub fn for_size(size: [f64; 2]) -> Distortion {
+        Distortion { k1: 0.0, k2: 0.0, radius: 0.5 * size[0].hypot(size[1]) }
+    }
+    pub fn is_none(&self) -> bool {
+        self.radius <= 0.0 || (self.k1 == 0.0 && self.k2 == 0.0)
+    }
+    fn rho2(&self, p: [f64; 2]) -> f64 {
+        (p[0] * p[0] + p[1] * p[1]) / (self.radius * self.radius)
+    }
+    /// Ideal → distorted (pixels relative to the principal point).
+    pub fn distort(&self, p: [f64; 2]) -> [f64; 2] {
+        if self.radius <= 0.0 {
+            return p;
+        }
+        let r2 = self.rho2(p);
+        let d = 1.0 + self.k1 * r2 + self.k2 * r2 * r2;
+        [p[0] * d, p[1] * d]
+    }
+    /// Distorted → ideal (fixed-point iteration).
+    pub fn undistort(&self, q: [f64; 2]) -> [f64; 2] {
+        if self.is_none() {
+            return q;
+        }
+        let mut p = q;
+        for _ in 0..30 {
+            let r2 = self.rho2(p);
+            let d = 1.0 + self.k1 * r2 + self.k2 * r2 * r2;
+            if d.abs() < 1e-6 {
+                break;
+            }
+            p = [q[0] / d, q[1] / d];
+        }
+        p
+    }
+    /// `d distort / d p` (2 × 2, row-major) and `d distort / d (k1, k2)` (rows x, y).
+    pub fn jacobian(&self, p: [f64; 2]) -> ([[f64; 2]; 2], [[f64; 2]; 2]) {
+        if self.radius <= 0.0 {
+            return ([[1.0, 0.0], [0.0, 1.0]], [[0.0; 2]; 2]);
+        }
+        let r2 = self.rho2(p);
+        let d = 1.0 + self.k1 * r2 + self.k2 * r2 * r2;
+        // d(r2)/dp = 2 p / R².
+        let dd = (self.k1 + 2.0 * self.k2 * r2) * 2.0 / (self.radius * self.radius);
+        let j = [[d + p[0] * dd * p[0], p[0] * dd * p[1]], [p[1] * dd * p[0], d + p[1] * dd * p[1]]];
+        let k = [[p[0] * r2, p[0] * r2 * r2], [p[1] * r2, p[1] * r2 * r2]];
+        (j, k)
+    }
+}
+
 /// 3D Camera Tracker ▸ Shot Type.
 #[derive(Clone, Copy, Debug, Default, PartialEq, Eq, Hash, Serialize, Deserialize)]
 #[serde(rename_all = "camelCase")]
@@ -178,6 +248,8 @@ pub struct SolveSettings {
     pub method: SolveMethod,
     /// Track ids the user deleted.
     pub deleted: Vec<u32>,
+    /// Solve radial lens distortion jointly with the cameras.
+    pub lens_distortion: bool,
 }
 
 /// The solved camera of one frame.
@@ -261,6 +333,9 @@ pub struct CameraSolve {
     pub points: Vec<SolvedPoint>,
     #[serde(skip_serializing_if = "Option::is_none")]
     pub ground: Option<Ground>,
+    /// The solved lens distortion (Solve Lens Distortion).
+    #[serde(skip_serializing_if = "Option::is_none")]
+    pub distortion: Option<Distortion>,
 }
 
 impl CameraSolve {
@@ -355,6 +430,36 @@ impl CameraSolve {
         }
     }
 
+    /// An ideal (pinhole) image point → where it appears in the footage (layer pixels), through
+    /// the solved lens distortion.
+    pub fn to_image(&self, p: [f64; 2]) -> [f64; 2] {
+        match &self.distortion {
+            Some(d) if !d.is_none() => {
+                let c = [self.size[0] * 0.5, self.size[1] * 0.5];
+                let q = d.distort([p[0] - c[0], p[1] - c[1]]);
+                [q[0] + c[0], q[1] + c[1]]
+            }
+            _ => p,
+        }
+    }
+    /// A footage point → its ideal (undistorted) position.
+    pub fn to_ideal(&self, p: [f64; 2]) -> [f64; 2] {
+        match &self.distortion {
+            Some(d) if !d.is_none() => {
+                let c = [self.size[0] * 0.5, self.size[1] * 0.5];
+                let q = d.undistort([p[0] - c[0], p[1] - c[1]]);
+                [q[0] + c[0], q[1] + c[1]]
+            }
+            _ => p,
+        }
+    }
+    /// Project world point `x` on frame `k`: through the lens distortion into the footage, or
+    /// (`undistorted`) into the undistorted image.
+    pub fn project(&self, k: usize, x: V3, undistorted: bool) -> Option<[f64; 2]> {
+        let q = self.frames.get(k)?.project(self.size, x)?;
+        Some(if undistorted { q } else { self.to_image(q) })
+    }
+
     /// The first solved frame.
     pub fn first_solved(&self) -> Option<usize> {
         self.frames.iter().position(|c| c.solved)
diff --git a/crates/track/src/camtrack/solver.rs b/crates/track/src/camtrack/solver.rs
index 85edd2bd..5ed7b81b 100644
--- a/crates/track/src/camtrack/solver.rs
+++ b/crates/track/src/camtrack/solver.rs
@@ -8,7 +8,7 @@ use rayon::prelude::*;
 use super::bundle::{self, BaCam, BaObs, BaOpts, FocalMode};
 use super::geometry::{self, Pose};
 use super::linalg::*;
-use super::{CameraSolve, CameraTracks, ShotType, SolveMethod, SolveSettings, SolvedFrame, SolvedPoint};
+use super::{CameraSolve, CameraTracks, Distortion, ShotType, SolveMethod, SolveSettings, SolvedFrame, SolvedPoint};
 use crate::fit::{Model, ransac};
 
 fn debug() -> bool {
@@ -142,6 +142,8 @@ struct Cfg {
     focal: FocalMode,
     /// Registration estimates a focal length per frame.
     per_frame_f: bool,
+    /// The global adjustments also adjust the lens distortion.
+    distortion: bool,
 }
 
 const F_PRIOR: f64 = 30.0;
@@ -183,7 +185,7 @@ impl Recon {
             &mut cams,
             &mut pts,
             &obs,
-            &BaOpts { focal, fix_points: false, fix_centers: cfg.kind == Kind::Tripod, fixed, huber: d.huber(), max_iter: iters },
+            &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);
@@ -289,6 +291,7 @@ impl Recon {
             fixed: vec![false],
             huber,
             max_iter: 40,
+            distortion: false,
         };
         let mut best: Option<(usize, BaCam)> = None;
         // Robust start (large Huber threshold), then tighten; from the prediction and, failing
@@ -371,7 +374,11 @@ fn predict(r: &Recon, k: usize) -> Option<BaCam> {
     match (before, after) {
         (Some((ka, a)), Some((kb, b))) => {
             let t = (k - ka) as f64 / (kb - ka) as f64;
-            Some(BaCam { pose: Pose { r: slerp(&a.pose.r, &b.pose.r, t), c: add(scale(a.pose.c, 1.0 - t), scale(b.pose.c, t)) }, f: a.f * (b.f / a.f).powf(t) })
+            Some(BaCam {
+                pose: Pose { r: slerp(&a.pose.r, &b.pose.r, t), c: add(scale(a.pose.c, 1.0 - t), scale(b.pose.c, t)) },
+                f: a.f * (b.f / a.f).powf(t),
+                dist: a.dist,
+            })
         }
         (Some((ka, a)), None) => {
             // Constant velocity from the two previous cameras.
@@ -381,7 +388,7 @@ fn predict(r: &Recon, k: usize) -> Option<BaCam> {
                     let dr = log_rot(&mm(&a.pose.r, &mt(&p.pose.r)));
                     let rr = mm(&rodrigues(scale(dr, t.min(2.0))), &a.pose.r);
                     let c = add(a.pose.c, scale(sub(a.pose.c, p.pose.c), t.min(2.0)));
-                    Some(BaCam { pose: Pose { r: rr, c }, f: a.f })
+                    Some(BaCam { pose: Pose { r: rr, c }, f: a.f, dist: a.dist })
                 }
                 None => Some(a),
             }
@@ -392,7 +399,7 @@ fn predict(r: &Recon, k: usize) -> Option<BaCam> {
                 let dr = log_rot(&mm(&b.pose.r, &mt(&nx.pose.r)));
                 let rr = mm(&rodrigues(scale(dr, t.min(2.0))), &b.pose.r);
                 let c = add(b.pose.c, scale(sub(b.pose.c, nx.pose.c), t.min(2.0)));
-                Some(BaCam { pose: Pose { r: rr, c }, f: b.f })
+                Some(BaCam { pose: Pose { r: rr, c }, f: b.f, dist: b.dist })
             }
             None => Some(b),
         },
@@ -459,8 +466,8 @@ fn init_pair(d: &Data, a: usize, b: usize, f: f64, use_e: bool, use_h: bool) ->
 fn grow(d: &Data, kfs: &[usize], j: usize, init: &(Pose, HashMap<u32, V3>), f: f64, cfg: &Cfg, cancel: Option<&AtomicBool>) -> Option<Recon> {
     let (a, b) = (kfs[0], kfs[j]);
     let mut r = Recon::default();
-    r.cams.insert(a, BaCam { pose: Pose::IDENTITY, f });
-    r.cams.insert(b, BaCam { pose: init.0, f });
+    r.cams.insert(a, BaCam { pose: Pose::IDENTITY, f, dist: Distortion::NONE });
+    r.cams.insert(b, BaCam { pose: init.0, f, dist: Distortion::NONE });
     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 };
@@ -495,7 +502,7 @@ fn grow(d: &Data, kfs: &[usize], j: usize, init: &(Pose, HashMap<u32, V3>), f: f
 /// Tripod reconstruction over keyframes `kfs` at focal length `f`.
 fn grow_tripod(d: &Data, kfs: &[usize], f: f64, cfg: &Cfg, cancel: Option<&AtomicBool>) -> Option<Recon> {
     let mut r = Recon::default();
-    r.cams.insert(kfs[0], BaCam { pose: Pose::IDENTITY, f });
+    r.cams.insert(kfs[0], BaCam { pose: Pose::IDENTITY, f, dist: Distortion::NONE });
     let ray = |p: [f64; 2], f: f64| normalize([p[0] / f, p[1] / f, 1.0]);
     for w in kfs.windows(2) {
         if cancel.is_some_and(|c| c.load(Ordering::Relaxed)) {
@@ -519,7 +526,7 @@ fn grow_tripod(d: &Data, kfs: &[usize], f: f64, cfg: &Cfg, cancel: Option<&Atomi
             rel = geometry::procrustes(&ia, &ib)?;
         }
         let pa = r.cams[&a].pose.r;
-        r.cams.insert(b, BaCam { pose: Pose { r: mm(&rel, &pa), c: [0.0; 3] }, f });
+        r.cams.insert(b, BaCam { pose: Pose { r: mm(&rel, &pa), c: [0.0; 3] }, f, dist: Distortion::NONE });
     }
     r.triangulate(d, Kind::Tripod, 0.0, 20.0);
     let frames: Vec<usize> = r.cams.keys().copied().collect();
@@ -640,8 +647,101 @@ fn search_focal(w: f64, eval: &(dyn Fn(f64) -> f64 + Sync)) -> f64 {
     if best.1 <= errs[i] { best.0.exp() } else { fs[i] }
 }
 
+/// A solve before its canonical output.
+struct Core {
+    d: Data,
+    r: Recon,
+    method: SolveMethod,
+    kind: Kind,
+    cfg: Cfg,
+    focal: FocalMode,
+    kfs: Vec<usize>,
+}
+
+/// `tracks` with every position undistorted by `dist`.
+fn undistorted(tracks: &CameraTracks, dist: &Distortion) -> CameraTracks {
+    let c = [tracks.size[0] * 0.5, tracks.size[1] * 0.5];
+    let mut t = tracks.clone();
+    for tr in t.tracks.iter_mut() {
+        for p in tr.pts.iter_mut() {
+            let q = dist.undistort([p[0] as f64 - c[0], p[1] as f64 - c[1]]);
+            *p = [(q[0] + c[0]) as f32, (q[1] + c[1]) as f32];
+        }
+    }
+    t
+}
+
 /// Solve the camera of a tracked clip.
 pub fn solve(tracks: &CameraTracks, s: &SolveSettings, cancel: Option<&AtomicBool>) -> Result<CameraSolve, SolveError> {
+    if !s.lens_distortion {
+        let c = solve_core(tracks, s, cancel)?;
+        return Ok(finish(&c.d, tracks, s, c.r, c.method, c.kind));
+    }
+    // Lens distortion: solve, adjust (k1, k2) jointly on the raw tracks, then solve again on
+    // tracks undistorted with the estimate and adjust jointly once more.
+    let mut dist = Distortion::for_size(tracks.size);
+    let mut out = None;
+    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);
+        d.set_noise(c.d.noise());
+        let mut r = c.r;
+        for cam in r.cams.values_mut() {
+            cam.dist = dist;
+        }
+        let cfg = Cfg { distortion: true, ..c.cfg };
+        let all: Vec<usize> = r.cams.keys().copied().collect();
+        let frames: Vec<usize> = if all.len() <= 160 { all.clone() } else { c.kfs.iter().copied().filter(|k| r.cams.contains_key(k)).collect() };
+        if frames.len() < 2 {
+            return err("the camera could not be solved (not enough parallax or features)");
+        }
+        // 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.outliers.clear();
+        r.reject(&d);
+        r.adjust(&d, &frames, &cfg, c.focal, 30);
+        dist = r.cams[&frames[0]].dist;
+        for cam in r.cams.values_mut() {
+            cam.dist = dist;
+        }
+        if frames.len() < all.len() {
+            // Refine the other frames against the adjusted points and distortion.
+            let base = r.clone();
+            let rcfg = Cfg { distortion: false, ..cfg };
+            let refined: Vec<(usize, Option<BaCam>)> = all
+                .par_iter()
+                .filter(|k| !frames.contains(k))
+                .map(|k| {
+                    let mut local = Recon { cams: base.cams.clone(), pts: base.pts.clone(), outliers: HashSet::new() };
+                    let start = local.cams[k];
+                    let ok = local.register(&d, *k, start, &rcfg);
+                    (*k, if ok { local.cams.get(k).copied() } else { None })
+                })
+                .collect();
+            for (k, cam) in refined {
+                if let Some(cam) = cam {
+                    r.cams.insert(k, cam);
+                }
+            }
+            r.reject(&d);
+        }
+        if cancel.is_some_and(|c| c.load(Ordering::Relaxed)) {
+            return err("cancelled");
+        }
+        if debug() {
+            eprintln!("lens distortion pass {pass}: k1 {:.5} k2 {:.5}", dist.k1, dist.k2);
+        }
+        let mut sol = finish(&d, tracks, s, r, c.method, c.kind);
+        sol.distortion = Some(dist);
+        out = Some(sol);
+    }
+    out.ok_or_else(|| SolveError("the camera could not be solved".into()))
+}
+
+fn solve_core(tracks: &CameraTracks, s: &SolveSettings, cancel: Option<&AtomicBool>) -> Result<Core, SolveError> {
     let d = Data::new(tracks, &s.deleted);
     if d.n < 2 {
         return err("the clip must be at least two frames long");
@@ -653,7 +753,7 @@ pub fn solve(tracks: &CameraTracks, s: &SolveSettings, cancel: Option<&AtomicBoo
     let cancelled = || cancel.is_some_and(|c| c.load(Ordering::Relaxed));
     let fixed_f = (s.shot == ShotType::SpecifyAngle).then(|| focal_of(d.w, if s.hfov > 0.0 { s.hfov } else { 40.0 }));
     let variable = s.shot == ShotType::VariableZoom;
-    let gcfg = Cfg { kind: Kind::General, focal: FocalMode::Fixed, per_frame_f: false };
+    let gcfg = Cfg { kind: Kind::General, focal: FocalMode::Fixed, per_frame_f: false, distortion: false };
     let tcfg = Cfg { kind: Kind::Tripod, ..gcfg };
     // The focal-length search runs on the first keyframes only.
     let mini: Vec<usize> = kfs.iter().copied().take(10).collect();
@@ -738,7 +838,7 @@ pub fn solve(tracks: &CameraTracks, s: &SolveSettings, cancel: Option<&AtomicBoo
     }
     let (att, _) = results.swap_remove(pick);
     let kind = if att.method == SolveMethod::TripodPan { Kind::Tripod } else { Kind::General };
-    let cfg = Cfg { kind, focal: FocalMode::Fixed, per_frame_f: variable };
+    let cfg = Cfg { kind, focal: FocalMode::Fixed, per_frame_f: variable, distortion: false };
     let mut r = att.recon;
     let focal = match (fixed_f, variable) {
         (Some(_), _) => FocalMode::Fixed,
@@ -845,7 +945,7 @@ pub fn solve(tracks: &CameraTracks, s: &SolveSettings, cancel: Option<&AtomicBoo
         (m, _) => m,
     };
     lap("final adjustment");
-    Ok(finish(&d, tracks, s, r, method, kind))
+    Ok(Core { d, r, method, kind, cfg, focal, kfs: kf_frames })
 }
 
 /// Thickness of a point cloud: √(smallest / middle eigenvalue of its covariance).
@@ -949,5 +1049,6 @@ fn finish(d: &Data, tracks: &CameraTracks, s: &SolveSettings, r: Recon, method:
         frames,
         points,
         ground: None,
+        distortion: None,
     }
 }
diff --git a/crates/track/src/camtrack/tests.rs b/crates/track/src/camtrack/tests.rs
index 4b5a8346..ff92cfe4 100644
--- a/crates/track/src/camtrack/tests.rs
+++ b/crates/track/src/camtrack/tests.rs
@@ -105,7 +105,7 @@ fn path_error(s: &CameraSolve, t: &Truth) -> (f64, f64) {
 }
 
 fn settings(shot: ShotType, method: SolveMethod) -> SolveSettings {
-    SolveSettings { shot, hfov: 0.0, method, deleted: vec![] }
+    SolveSettings { shot, hfov: 0.0, method, deleted: vec![], lens_distortion: false }
 }
 
 #[test]
@@ -149,7 +149,8 @@ fn specified_angle_of_view_is_kept() {
     let pts = random_points(300, 5, [0.0, 0.0, 10.0], [14.0, 8.0, 8.0]);
     let tr = tracks_of(&pts, &t, 6, true);
     let hfov = 2.0 * (W * 0.5 / f).atan().to_degrees();
-    let s = solve(&tr, &SolveSettings { shot: ShotType::SpecifyAngle, hfov, method: SolveMethod::Typical, deleted: vec![] }, None).unwrap();
+    let s =
+        solve(&tr, &SolveSettings { shot: ShotType::SpecifyAngle, hfov, method: SolveMethod::Typical, deleted: vec![], lens_distortion: false }, None).unwrap();
     assert!((s.focal() - f).abs() < 1e-6);
     assert!(s.average_error < 0.5);
     let (pe, _) = path_error(&s, &t);
@@ -226,7 +227,9 @@ fn mostly_flat_scene_and_ground_plane() {
     }
     // Ground plane from the solved (coplanar) points: every point lies on it.
     let hfov = 2.0 * (W * 0.5 / f).atan().to_degrees();
-    let mut s = solve(&tr, &SolveSettings { shot: ShotType::SpecifyAngle, hfov, method: SolveMethod::MostlyFlat, deleted: vec![] }, None).expect("solve");
+    let mut s =
+        solve(&tr, &SolveSettings { shot: ShotType::SpecifyAngle, hfov, method: SolveMethod::MostlyFlat, deleted: vec![], lens_distortion: false }, None)
+            .expect("solve");
     let ids: Vec<u32> = s.visible(0).map(|p| p.id).take(12).collect();
     let g = s.ground_from(&ids, 0).expect("ground");
     let spread = s.points.iter().map(|p| dot(sub(p.pos, g.origin), g.normal).abs()).fold(0.0, f64::max);
@@ -251,7 +254,9 @@ fn deleted_tracks_are_ignored_and_cancel_stops() {
     let tr = tracks_of(&pts, &t, 14, true);
     let deleted: Vec<u32> = (0..40).collect();
     let hfov = 2.0 * (W * 0.5 / 1000.0f64).atan().to_degrees();
-    let s = solve(&tr, &SolveSettings { shot: ShotType::SpecifyAngle, hfov, method: SolveMethod::Typical, deleted: deleted.clone() }, None).unwrap();
+    let s =
+        solve(&tr, &SolveSettings { shot: ShotType::SpecifyAngle, hfov, method: SolveMethod::Typical, deleted: deleted.clone(), lens_distortion: false }, None)
+            .unwrap();
     assert!(s.points.iter().all(|p| !deleted.contains(&p.id)));
     let stop = std::sync::atomic::AtomicBool::new(true);
     let r = solve(&tr, &settings(ShotType::FixedAngle, SolveMethod::Typical), Some(&stop));
@@ -436,3 +441,55 @@ fn solve_saved_tracks() {
         t0.elapsed()
     );
 }
+
+/// Pass every track position through a lens distortion (around the frame centre).
+fn distort_tracks(t: &CameraTracks, d: &Distortion) -> CameraTracks {
+    let mut t = t.clone();
+    for tr in t.tracks.iter_mut() {
+        for p in tr.pts.iter_mut() {
+            let q = d.distort([p[0] as f64 - W / 2.0, p[1] as f64 - H / 2.0]);
+            *p = [(q[0] + W / 2.0) as f32, (q[1] + H / 2.0) as f32];
+        }
+    }
+    t
+}
+
+#[test]
+fn lens_distortion_is_recovered_in_the_bundle_adjustment() {
+    let f = 1000.0;
+    let t = moving_camera(36, |_| f);
+    let pts = random_points(450, 7, [0.0, 0.0, 10.0], [16.0, 9.0, 8.0]);
+    let truth = Distortion { k1: -0.12, k2: 0.03, radius: 0.5 * W.hypot(H) };
+    let tr = distort_tracks(&tracks_of(&pts, &t, 8, true), &truth);
+    let plain = solve(&tr, &settings(ShotType::FixedAngle, SolveMethod::Typical), None).expect("pinhole solve");
+    let lens = solve(&tr, &SolveSettings { lens_distortion: true, ..settings(ShotType::FixedAngle, SolveMethod::Typical) }, None).expect("lens solve");
+    let d = lens.distortion.expect("distortion solved");
+    eprintln!(
+        "pinhole error {:.3} px (focal {:.1}); with distortion {:.3} px (focal {:.1}), k1 {:.4} k2 {:.4}",
+        plain.average_error,
+        plain.focal(),
+        lens.average_error,
+        lens.focal(),
+        d.k1,
+        d.k2
+    );
+    assert!(plain.distortion.is_none());
+    assert!((d.k1 - truth.k1).abs() < 0.005, "k1 {}", d.k1);
+    assert!((d.k2 - truth.k2).abs() < 0.01, "k2 {}", d.k2);
+    // The distortion model explains the tracks: about the noise level, far below the pinhole fit.
+    assert!(lens.average_error < 0.3, "{}", lens.average_error);
+    assert!(lens.average_error < 0.5 * plain.average_error, "{} vs {}", lens.average_error, plain.average_error);
+    assert!((lens.focal() - f).abs() < 0.02 * f, "focal {}", lens.focal());
+    let (path, rot) = path_error(&lens, &t);
+    assert!(path < 0.01 && rot < 0.3, "path {path} rot {rot}");
+    // Projection through the lens lands on the (distorted) tracks; undistorting inverts it.
+    let p = [100.0, 80.0];
+    let q = lens.to_image(p);
+    let back = lens.to_ideal(q);
+    assert!((back[0] - p[0]).abs() < 1e-6 && (back[1] - p[1]).abs() < 1e-6);
+    let sp = &lens.points[0];
+    let k = sp.first as usize;
+    let obs = tr.track(sp.id).and_then(|t| t.at(k as u32)).expect("observed");
+    let proj = lens.project(k, sp.pos, false).expect("in front");
+    assert!((proj[0] - obs[0]).hypot(proj[1] - obs[1]) < 2.0, "{proj:?} vs {obs:?}");
+}
diff --git a/crates/track/src/face.rs b/crates/track/src/face.rs
new file mode 100644
index 00000000..62073d0a
--- /dev/null
+++ b/crates/track/src/face.rs
@@ -0,0 +1,1147 @@
+//! Face tracking (Tracker panel ▸ Method ▸ Face Tracking (Outline Only) / (Detailed Features)),
+//! a classical pipeline with no learned weights:
+//!
+//! 1. **Skin model**: on the first frame, the colour of the middle of the mask (luma and
+//!    normalised r, g chromaticity: median and median absolute deviation) defines skin; every
+//!    later frame is classified against it inside a search window around the previous face.
+//! 2. **Face region and outline**: the skin component containing the face centre is flood-filled;
+//!    its outer boundary is found along 48 rays from its centroid (so the non-skin holes of eyes,
+//!    brows, nostrils and mouth don't matter), median-smoothed, and an ellipse is fitted to the
+//!    outline polygon from its area moments (centre, half-height, roll). The outline is what
+//!    *Outline Only* keys into the mask.
+//! 3. **Facial features**: the non-skin connected components inside the face are matched to the
+//!    eyes, brows, nostrils and mouth by their position in the face frame (predicted from the
+//!    previous frame); landmark points are measured on each component (extremes along the face
+//!    axes for corners and lids, a darkness-weighted centroid for the pupil). Chin and jaw points
+//!    come from the outline.
+//! 4. **Shape model**: an active-shape point distribution model (Cootes, Taylor, Cooper & Graham,
+//!    "Active Shape Models — Their Training and Application", CVIU 61(1), 1995) regularises the
+//!    landmarks: a weighted similarity alignment and shape-parameter fit (parameters clamped to
+//!    ±3 standard deviations) fills landmarks whose feature was not found (a closed eye, a
+//!    merged component) and replaces detections that disagree with the model. The model is
+//!    *trained on synthetic shapes* produced by our own parametric face generator
+//!    ([`synth_shape`]: eye and mouth openness, mouth width, brow raise, head yaw and pitch) by
+//!    principal component analysis — no external data or weights.
+//!
+//! Coordinates are layer pixels; the face frame is centred on the face ellipse, with the
+//! half-height as unit, `u` right and `v` down (image-left is "Left").
+
+use std::sync::OnceLock;
+
+use effectcraft_raster::Image;
+
+use crate::Frame;
+use crate::camtrack::linalg::{cholesky_solve, sym_eigen};
+
+/// Landmarks: (id, display name, canonical position in the face frame).
+pub const LANDMARKS: [(&str, &str, [f64; 2]); 26] = [
+    ("leftEyebrowInner", "Left Eyebrow Inner", [-0.16, -0.44]),
+    ("leftEyebrowMiddle", "Left Eyebrow Middle", [-0.32, -0.50]),
+    ("leftEyebrowOuter", "Left Eyebrow Outer", [-0.48, -0.44]),
+    ("rightEyebrowInner", "Right Eyebrow Inner", [0.16, -0.44]),
+    ("rightEyebrowMiddle", "Right Eyebrow Middle", [0.32, -0.50]),
+    ("rightEyebrowOuter", "Right Eyebrow Outer", [0.48, -0.44]),
+    ("leftEyeInner", "Left Eye Inner", [-0.17, -0.22]),
+    ("leftEyeOuter", "Left Eye Outer", [-0.45, -0.22]),
+    ("leftEyeTop", "Left Eye Top", [-0.31, -0.292]),
+    ("leftEyeBottom", "Left Eye Bottom", [-0.31, -0.148]),
+    ("leftPupil", "Left Pupil", [-0.31, -0.22]),
+    ("rightEyeInner", "Right Eye Inner", [0.17, -0.22]),
+    ("rightEyeOuter", "Right Eye Outer", [0.45, -0.22]),
+    ("rightEyeTop", "Right Eye Top", [0.31, -0.292]),
+    ("rightEyeBottom", "Right Eye Bottom", [0.31, -0.148]),
+    ("rightPupil", "Right Pupil", [0.31, -0.22]),
+    ("noseTip", "Nose Tip", [0.0, 0.14]),
+    ("leftNostril", "Left Nostril", [-0.08, 0.20]),
+    ("rightNostril", "Right Nostril", [0.08, 0.20]),
+    ("mouthLeft", "Mouth Left", [-0.24, 0.50]),
+    ("mouthRight", "Mouth Right", [0.24, 0.50]),
+    ("mouthTop", "Mouth Top", [0.0, 0.435]),
+    ("mouthBottom", "Mouth Bottom", [0.0, 0.565]),
+    ("chin", "Chin", [0.0, 1.0]),
+    ("leftJaw", "Left Jaw", [-0.614, 0.574]),
+    ("rightJaw", "Right Jaw", [0.614, 0.574]),
+];
+
+/// Number of landmarks.
+pub const N: usize = LANDMARKS.len();
+
+/// Landmark index by id.
+pub fn landmark(id: &str) -> Option<usize> {
+    LANDMARKS.iter().position(|l| l.0 == id)
+}
+
+const BROW: [[usize; 3]; 2] = [[0, 1, 2], [3, 4, 5]];
+/// Eye landmarks: inner, outer, top, bottom, pupil.
+const EYE: [[usize; 5]; 2] = [[6, 7, 8, 9, 10], [11, 12, 13, 14, 15]];
+const NOSE_TIP: usize = 16;
+const NOSTRIL: [usize; 2] = [17, 18];
+/// Mouth: left, right, top, bottom.
+const MOUTH: [usize; 4] = [19, 20, 21, 22];
+const CHIN: usize = 23;
+const JAW: [usize; 2] = [24, 25];
+/// Half-width of the face relative to its half-height.
+const ASPECT: f64 = 0.75;
+
+// ---------------------------------------------------------------- synthetic faces
+
+/// Parameters of a synthetic face (the shape model's training data and the tests' footage).
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub struct FaceParams {
+    /// Face centre (layer pixels) and half-height (pixels).
+    pub center: [f64; 2],
+    pub height: f64,
+    /// Head rotation in degrees: roll (in the image plane), yaw (turn), pitch (nod).
+    pub roll: f64,
+    pub yaw: f64,
+    pub pitch: f64,
+    /// Eye openness (0 closed … 1 open), left and right.
+    pub eye_open: [f64; 2],
+    /// Mouth openness (0 … 1) and width (1 = neutral).
+    pub mouth_open: f64,
+    pub mouth_width: f64,
+    /// Eyebrow raise (−1 … 1), left and right.
+    pub brow_raise: [f64; 2],
+}
+
+impl Default for FaceParams {
+    fn default() -> Self {
+        FaceParams {
+            center: [0.0; 2],
+            height: 60.0,
+            roll: 0.0,
+            yaw: 0.0,
+            pitch: 0.0,
+            eye_open: [1.0; 2],
+            mouth_open: 0.2,
+            mouth_width: 1.0,
+            brow_raise: [0.0; 2],
+        }
+    }
+}
+
+fn eye_half_height(open: f64) -> f64 {
+    0.012 + 0.06 * open.clamp(0.0, 1.0)
+}
+fn mouth_half_height(open: f64) -> f64 {
+    0.015 + 0.11 * open.clamp(0.0, 1.0)
+}
+
+/// The landmarks of a synthetic face in the face frame (before roll, translation and scale).
+pub fn synth_shape(p: &FaceParams) -> [[f64; 2]; N] {
+    let mut s = [[0.0; 2]; N];
+    let (sy, sp) = (p.yaw.to_radians().sin(), p.pitch.to_radians().sin());
+    // Inner features slide with the head's turn; the nose (closer to the camera) slides more.
+    let shift = |depth: f64| [0.35 * sy * depth, 0.3 * sp * depth];
+    for (i, l) in LANDMARKS.iter().enumerate() {
+        s[i] = l.2;
+    }
+    for side in 0..2 {
+        let sg = if side == 0 { -1.0 } else { 1.0 };
+        let raise = -0.06 * p.brow_raise[side].clamp(-1.0, 1.0);
+        for &i in &BROW[side] {
+            s[i][1] += raise;
+        }
+        let h = eye_half_height(p.eye_open[side]);
+        let e = EYE[side];
+        s[e[2]] = [sg * 0.31, -0.22 - h];
+        s[e[3]] = [sg * 0.31, -0.22 + h];
+    }
+    let mw = 0.24 * p.mouth_width.clamp(0.5, 1.6);
+    let mh = mouth_half_height(p.mouth_open);
+    s[MOUTH[0]] = [-mw, 0.5];
+    s[MOUTH[1]] = [mw, 0.5];
+    s[MOUTH[2]] = [0.0, 0.5 - mh];
+    s[MOUTH[3]] = [0.0, 0.5 + mh];
+    for (i, q) in s.iter_mut().enumerate() {
+        if i == CHIN || JAW.contains(&i) {
+            continue;
+        }
+        let d = if i == NOSE_TIP {
+            1.3
+        } else if NOSTRIL.contains(&i) {
+            1.15
+        } else {
+            1.0
+        };
+        let o = shift(d);
+        q[0] += o[0];
+        q[1] += o[1];
+    }
+    s
+}
+
+/// Face frame → layer pixels.
+fn to_layer(c: [f64; 2], b: f64, roll_deg: f64, q: [f64; 2]) -> [f64; 2] {
+    let (sn, cs) = roll_deg.to_radians().sin_cos();
+    [c[0] + b * (cs * q[0] - sn * q[1]), c[1] + b * (sn * q[0] + cs * q[1])]
+}
+
+/// Layer pixels → face frame.
+fn to_face(c: [f64; 2], b: f64, roll_deg: f64, p: [f64; 2]) -> [f64; 2] {
+    let (sn, cs) = roll_deg.to_radians().sin_cos();
+    let d = [(p[0] - c[0]) / b, (p[1] - c[1]) / b];
+    [cs * d[0] + sn * d[1], -sn * d[0] + cs * d[1]]
+}
+
+/// A synthetic face's landmarks in layer pixels.
+pub fn synth_landmarks(p: &FaceParams) -> [[f64; 2]; N] {
+    synth_shape(p).map(|q| to_layer(p.center, p.height, p.roll, q))
+}
+
+fn hash(i: i64, j: i64) -> f64 {
+    let mut v = (i.wrapping_mul(73_856_093) ^ j.wrapping_mul(19_349_663)) as u64;
+    v ^= v >> 13;
+    v = v.wrapping_mul(0x5bd1_e995);
+    v ^= v >> 15;
+    (v % 10_000) as f64 / 10_000.0
+}
+
+fn seg_dist(p: [f64; 2], a: [f64; 2], b: [f64; 2]) -> f64 {
+    let d = [b[0] - a[0], b[1] - a[1]];
+    let l2 = d[0] * d[0] + d[1] * d[1];
+    let t = if l2 > 0.0 { (((p[0] - a[0]) * d[0] + (p[1] - a[1]) * d[1]) / l2).clamp(0.0, 1.0) } else { 0.0 };
+    (p[0] - a[0] - t * d[0]).hypot(p[1] - a[1] - t * d[1])
+}
+
+/// Colour of a synthetic face at face-frame point `q` (`None` outside the face).
+fn face_color(p: &FaceParams, s: &[[f64; 2]; N], q: [f64; 2]) -> Option<[f32; 3]> {
+    if (q[0] / ASPECT).powi(2) + q[1] * q[1] > 1.0 {
+        return None;
+    }
+    let shade = (0.95 - 0.08 * q[1]) as f32;
+    let mut col = [0.86 * shade, 0.64 * shade, 0.52 * shade];
+    for side in 0..2 {
+        let b = BROW[side];
+        let d = seg_dist(q, s[b[0]], s[b[1]]).min(seg_dist(q, s[b[1]], s[b[2]]));
+        if d < 0.035 {
+            col = [0.28, 0.18, 0.12];
+        }
+        let e = EYE[side];
+        let ec = [(s[e[0]][0] + s[e[1]][0]) / 2.0, (s[e[2]][1] + s[e[3]][1]) / 2.0];
+        let (hw, hh) = ((s[e[1]][0] - s[e[0]][0]).abs() / 2.0, (s[e[3]][1] - s[e[2]][1]).abs() / 2.0);
+        if ((q[0] - ec[0]) / hw).powi(2) + ((q[1] - ec[1]) / hh).powi(2) <= 1.0 {
+            let pupil = s[e[4]];
+            col = if (q[0] - pupil[0]).hypot(q[1] - pupil[1]) < 0.055 { [0.12, 0.08, 0.06] } else { [0.95, 0.95, 0.93] };
+        }
+    }
+    let _ = p;
+    for &n in &NOSTRIL {
+        if ((q[0] - s[n][0]) / 0.035).powi(2) + ((q[1] - s[n][1]) / 0.022).powi(2) <= 1.0 {
+            col = [0.3, 0.13, 0.1];
+        }
+    }
+    let m = MOUTH;
+    let mc = [(s[m[0]][0] + s[m[1]][0]) / 2.0, (s[m[2]][1] + s[m[3]][1]) / 2.0];
+    let (mw, mh) = ((s[m[1]][0] - s[m[0]][0]) / 2.0, (s[m[3]][1] - s[m[2]][1]) / 2.0);
+    if ((q[0] - mc[0]) / mw).powi(2) + ((q[1] - mc[1]) / mh).powi(2) <= 1.0 {
+        col = [0.42, 0.08, 0.1];
+    }
+    Some(col)
+}
+
+/// Render a synthetic face over a textured blue-green background (2 × 2 supersampled): the
+/// test footage of the face tracker.
+pub fn render_synthetic(w: u32, h: u32, p: &FaceParams, background_seed: i64) -> Image {
+    use rayon::prelude::*;
+    let s = synth_shape(p);
+    let mut img = Image::new(w, h);
+    img.rows_mut().for_each(|(y, row)| {
+        for (x, px) in row.iter_mut().enumerate() {
+            let mut acc = [0.0f32; 3];
+            for sy in 0..2 {
+                for sx in 0..2 {
+                    let l = [x as f64 + 0.25 + 0.5 * sx as f64, y as f64 + 0.25 + 0.5 * sy as f64];
+                    let q = to_face(p.center, p.height, p.roll, l);
+                    let c = face_color(p, &s, q).unwrap_or_else(|| {
+                        let n = hash((l[0] / 9.0) as i64 + background_seed, (l[1] / 9.0) as i64) as f32;
+                        [0.15 + 0.1 * n, 0.35 + 0.2 * n, 0.45 + 0.25 * n]
+                    });
+                    for k in 0..3 {
+                        acc[k] += c[k] * 0.25;
+                    }
+                }
+            }
+            *px = [acc[0], acc[1], acc[2], 1.0];
+        }
+    });
+    img
+}
+
+// ---------------------------------------------------------------- shape model
+
+/// Point distribution model: mean shape and principal modes (face frame).
+pub struct Pdm {
+    pub mean: Vec<f64>,
+    /// Modes as rows of length `2N`.
+    pub modes: Vec<Vec<f64>>,
+    /// Variance of each mode.
+    pub var: Vec<f64>,
+}
+
+/// The shape model, trained once on synthetic shapes.
+pub fn pdm() -> &'static Pdm {
+    static P: OnceLock<Pdm> = OnceLock::new();
+    P.get_or_init(|| {
+        let samples = 600;
+        let dim = 2 * N;
+        let mut rng = 0x2545_f491_4f6c_dd1du64;
+        let mut u = || {
+            rng ^= rng << 13;
+            rng ^= rng >> 7;
+            rng ^= rng << 17;
+            (rng % 1_000_000) as f64 / 1_000_000.0
+        };
+        let mut data = Vec::with_capacity(samples);
+        for _ in 0..samples {
+            let p = FaceParams {
+                yaw: (u() - 0.5) * 60.0,
+                pitch: (u() - 0.5) * 40.0,
+                eye_open: [u(), u()],
+                mouth_open: u(),
+                mouth_width: 0.8 + 0.45 * u(),
+                brow_raise: [u() * 2.0 - 1.0, u() * 2.0 - 1.0],
+                ..Default::default()
+            };
+            data.push(synth_shape(&p).iter().flat_map(|q| [q[0], q[1]]).collect::<Vec<f64>>());
+        }
+        let mut mean = vec![0.0; dim];
+        for d in &data {
+            for (m, v) in mean.iter_mut().zip(d) {
+                *m += v / samples as f64;
+            }
+        }
+        let mut cov = vec![0.0; dim * dim];
+        for d in &data {
+            for a in 0..dim {
+                let da = d[a] - mean[a];
+                for b in 0..dim {
+                    cov[a * dim + b] += da * (d[b] - mean[b]) / samples as f64;
+                }
+            }
+        }
+        let (vals, vecs) = sym_eigen(&cov, dim);
+        let total: f64 = vals.iter().map(|v| v.max(0.0)).sum();
+        let mut modes = vec![];
+        let mut var = vec![];
+        let mut acc = 0.0;
+        for k in (0..dim).rev() {
+            if acc >= 0.995 * total || vals[k] <= 1e-12 {
+                break;
+            }
+            acc += vals[k];
+            modes.push((0..dim).map(|r| vecs[r * dim + k]).collect());
+            var.push(vals[k]);
+        }
+        Pdm { mean, modes, var }
+    })
+}
+
+/// A similarity `p ↦ s R q + t` as `[a, b, tx, ty]` (`x = a u − b v + tx`, `y = b u + a v + ty`).
+fn fit_similarity(src: &[[f64; 2]], dst: &[[f64; 2]], w: &[f64]) -> Option<[f64; 4]> {
+    let ws: f64 = w.iter().sum();
+    if ws < 2.0 - 1e-9 {
+        return None;
+    }
+    let mut cs = [0.0; 2];
+    let mut cd = [0.0; 2];
+    for i in 0..src.len() {
+        for k in 0..2 {
+            cs[k] += w[i] * src[i][k] / ws;
+            cd[k] += w[i] * dst[i][k] / ws;
+        }
+    }
+    let (mut ss, mut dt, mut cr) = (0.0, 0.0, 0.0);
+    for i in 0..src.len() {
+        let a = [src[i][0] - cs[0], src[i][1] - cs[1]];
+        let b = [dst[i][0] - cd[0], dst[i][1] - cd[1]];
+        ss += w[i] * (a[0] * a[0] + a[1] * a[1]);
+        dt += w[i] * (a[0] * b[0] + a[1] * b[1]);
+        cr += w[i] * (a[0] * b[1] - a[1] * b[0]);
+    }
+    if ss < 1e-12 {
+        return None;
+    }
+    let (a, b) = (dt / ss, cr / ss);
+    Some([a, b, cd[0] - (a * cs[0] - b * cs[1]), cd[1] - (b * cs[0] + a * cs[1])])
+}
+
+fn sim_apply(t: &[f64; 4], q: [f64; 2]) -> [f64; 2] {
+    [t[0] * q[0] - t[1] * q[1] + t[2], t[1] * q[0] + t[0] * q[1] + t[3]]
+}
+
+fn sim_inverse(t: &[f64; 4], p: [f64; 2]) -> [f64; 2] {
+    let d = [p[0] - t[2], p[1] - t[3]];
+    let s2 = t[0] * t[0] + t[1] * t[1];
+    [(t[0] * d[0] + t[1] * d[1]) / s2, (-t[1] * d[0] + t[0] * d[1]) / s2]
+}
+
+/// Fit the shape model to detected landmarks (`w` = 0 for missing ones), starting from the
+/// face-frame similarity `init`. Returns (model landmarks in layer pixels, similarity).
+pub fn fit_shape(detected: &[[f64; 2]; N], w: &[f64; N], init: [f64; 4]) -> ([[f64; 2]; N], [f64; 4]) {
+    let m = pdm();
+    let k = m.modes.len();
+    let mut b = vec![0.0; k];
+    let mut t = init;
+    let shape = |b: &[f64]| -> [[f64; 2]; N] {
+        let mut s = [[0.0; 2]; N];
+        for i in 0..N {
+            for d in 0..2 {
+                let j = 2 * i + d;
+                s[i][d] = m.mean[j] + m.modes.iter().zip(b).map(|(row, bk)| row[j] * bk).sum::<f64>();
+            }
+        }
+        s
+    };
+    for _ in 0..10 {
+        let x = shape(&b);
+        if let Some(nt) = fit_similarity(&x, detected, w) {
+            t = nt;
+        }
+        // Detections in the model frame.
+        let y: Vec<[f64; 2]> = detected.iter().map(|p| sim_inverse(&t, *p)).collect();
+        let mut ata = vec![0.0; k * k];
+        let mut atb = vec![0.0; k];
+        for i in 0..N {
+            if w[i] <= 0.0 {
+                continue;
+            }
+            for d in 0..2 {
+                let j = 2 * i + d;
+                let r = y[i][d] - m.mean[j];
+                for a in 0..k {
+                    atb[a] += w[i] * m.modes[a][j] * r;
+                    for c in 0..k {
+                        ata[a * k + c] += w[i] * m.modes[a][j] * m.modes[c][j];
+                    }
+                }
+            }
+        }
+        for a in 0..k {
+            // A weak prior towards the mean (Tikhonov, scaled by the mode's variance).
+            ata[a * k + a] += 1e-3 / m.var[a].max(1e-9) * 1e-3 + 1e-9;
+        }
+        if cholesky_solve(&mut ata, k, &mut atb) {
+            for a in 0..k {
+                let lim = 3.0 * m.var[a].sqrt();
+                b[a] = atb[a].clamp(-lim, lim);
+            }
+        }
+    }
+    let x = shape(&b);
+    (x.map(|q| sim_apply(&t, q)), t)
+}
+
+// ---------------------------------------------------------------- tracking
+
+/// The skin colour model (luma, r and g chromaticity: median and spread).
+#[derive(Clone, Copy, Debug, PartialEq)]
+struct Skin {
+    med: [f32; 3],
+    tol: [f32; 3],
+}
+
+fn features(px: [f32; 4]) -> [f32; 3] {
+    let s = (px[0] + px[1] + px[2]).max(1e-4);
+    [0.2126 * px[0] + 0.7152 * px[1] + 0.0722 * px[2], px[0] / s, px[1] / s]
+}
+
+impl Skin {
+    fn is_skin(&self, px: [f32; 4]) -> bool {
+        let f = features(px);
+        (0..3).all(|k| (f[k] - self.med[k]).abs() <= self.tol[k])
+    }
+}
+
+/// One frame's face fit.
+#[derive(Clone, Debug, PartialEq)]
+pub struct FaceFit {
+    /// Landmarks ([`LANDMARKS`] order), layer pixels.
+    pub landmarks: [[f64; 2]; N],
+    /// Which landmarks were measured on this frame (the others come from the shape model).
+    pub detected: [bool; N],
+    /// The face outline: a closed polygon (layer pixels).
+    pub outline: Vec<[f64; 2]>,
+    /// Face ellipse centre, half-height (pixels) and roll (degrees).
+    pub center: [f64; 2],
+    pub height: f64,
+    pub roll: f64,
+}
+
+/// Follows one face through a clip.
+#[derive(Clone, Debug)]
+pub struct FaceTracker {
+    skin: Skin,
+    prev: FaceFit,
+    /// Previous landmarks in the face frame (the expected feature positions).
+    expect: [[f64; 2]; N],
+}
+
+/// A pixel window of a frame in layer coordinates.
+struct Window {
+    x0: i64,
+    y0: i64,
+    w: usize,
+    h: usize,
+}
+
+impl Window {
+    fn around(frame: &Frame, c: [f64; 2], rx: f64, ry: f64) -> Option<Window> {
+        let (ox, oy) = (frame.offset[0], frame.offset[1]);
+        let x0 = ((c[0] - rx + ox).floor() as i64).max(0);
+        let y0 = ((c[1] - ry + oy).floor() as i64).max(0);
+        let x1 = ((c[0] + rx + ox).ceil() as i64).min(frame.img.width as i64);
+        let y1 = ((c[1] + ry + oy).ceil() as i64).min(frame.img.height as i64);
+        (x1 > x0 + 4 && y1 > y0 + 4).then(|| Window { x0, y0, w: (x1 - x0) as usize, h: (y1 - y0) as usize })
+    }
+    fn px(&self, frame: &Frame, i: usize, j: usize) -> [f32; 4] {
+        frame.img.get(self.x0 + i as i64, self.y0 + j as i64)
+    }
+    /// Layer position of the centre of window pixel `(i, j)`.
+    fn layer(&self, frame: &Frame, i: usize, j: usize) -> [f64; 2] {
+        [(self.x0 + i as i64) as f64 + 0.5 - frame.offset[0], (self.y0 + j as i64) as f64 + 0.5 - frame.offset[1]]
+    }
+    fn cell(&self, frame: &Frame, p: [f64; 2]) -> Option<(usize, usize)> {
+        let x = (p[0] + frame.offset[0]).floor() as i64 - self.x0;
+        let y = (p[1] + frame.offset[1]).floor() as i64 - self.y0;
+        (x >= 0 && y >= 0 && (x as usize) < self.w && (y as usize) < self.h).then_some((x as usize, y as usize))
+    }
+}
+
+fn median(v: &mut [f32]) -> f32 {
+    v.sort_by(f32::total_cmp);
+    v[v.len() / 2]
+}
+
+/// Area-moment ellipse of a closed polygon: (centroid, half-height along the major axis,
+/// half-width, major axis angle from vertical in degrees).
+fn polygon_ellipse(poly: &[[f64; 2]]) -> Option<([f64; 2], f64, f64, f64)> {
+    let n = poly.len();
+    let (mut a, mut cx, mut cy) = (0.0, 0.0, 0.0);
+    let (mut ixx, mut iyy, mut ixy) = (0.0, 0.0, 0.0);
+    for i in 0..n {
+        let (p, q) = (poly[i], poly[(i + 1) % n]);
+        let cr = p[0] * q[1] - q[0] * p[1];
+        a += cr;
+        cx += (p[0] + q[0]) * cr;
+        cy += (p[1] + q[1]) * cr;
+        ixx += (p[1] * p[1] + p[1] * q[1] + q[1] * q[1]) * cr;
+        iyy += (p[0] * p[0] + p[0] * q[0] + q[0] * q[0]) * cr;
+        ixy += (p[0] * q[1] + 2.0 * p[0] * p[1] + 2.0 * q[0] * q[1] + q[0] * p[1]) * cr;
+    }
+    a *= 0.5;
+    if a.abs() < 1e-9 {
+        return None;
+    }
+    cx /= 6.0 * a;
+    cy /= 6.0 * a;
+    // Central second moments (per unit area): var x, var y, cov xy.
+    let vxx = iyy / (12.0 * a) - cx * cx;
+    let vyy = ixx / (12.0 * a) - cy * cy;
+    let vxy = ixy / (24.0 * a) - cx * cy;
+    let tr = vxx + vyy;
+    let det = vxx * vyy - vxy * vxy;
+    let disc = (tr * tr / 4.0 - det).max(0.0).sqrt();
+    let (l1, l2) = (tr / 2.0 + disc, (tr / 2.0 - disc).max(0.0));
+    // Major axis direction.
+    let ang = 0.5 * (2.0 * vxy).atan2(vxx - vyy); // from +x
+    let mut from_vertical = ang.to_degrees() - 90.0;
+    while from_vertical > 90.0 {
+        from_vertical -= 180.0;
+    }
+    while from_vertical < -90.0 {
+        from_vertical += 180.0;
+    }
+    Some(([cx, cy], 2.0 * l1.sqrt(), 2.0 * l2.sqrt(), from_vertical))
+}
+
+/// A non-skin component inside the face.
+struct Blob {
+    /// Face-frame positions of its pixels and their darkness (skin luma − luma).
+    pts: Vec<([f64; 2], f64)>,
+    centroid: [f64; 2],
+}
+
+impl Blob {
+    fn extreme(&self, dir: [f64; 2]) -> [f64; 2] {
+        self.pts.iter().map(|p| p.0).max_by(|a, b| (a[0] * dir[0] + a[1] * dir[1]).total_cmp(&(b[0] * dir[0] + b[1] * dir[1]))).unwrap_or(self.centroid)
+    }
+    /// Extreme along `v` (`sign` −1 = top) among the pixels within `band` of the column `u`.
+    fn lid(&self, u: f64, band: f64, sign: f64) -> [f64; 2] {
+        let c = self.pts.iter().filter(|p| (p.0[0] - u).abs() <= band).map(|p| p.0);
+        let best = c.max_by(|a, b| (sign * a[1]).total_cmp(&(sign * b[1])));
+        best.map(|p| [u, p[1]]).unwrap_or(self.centroid)
+    }
+    fn dark_centroid(&self) -> [f64; 2] {
+        let (mut s, mut w) = ([0.0; 2], 0.0);
+        for (q, d) in &self.pts {
+            let k = d.max(0.0).powi(2);
+            s[0] += k * q[0];
+            s[1] += k * q[1];
+            w += k;
+        }
+        if w > 1e-12 { [s[0] / w, s[1] / w] } else { self.centroid }
+    }
+}
+
+impl FaceTracker {
+    /// Start on `frame` with the face inside `region` (`[x0, y0, x1, y1]`, layer pixels: the
+    /// mask's bounds). `None` when no face-like skin region is found.
+    pub fn new(frame: &Frame, region: [f64; 4]) -> Option<(FaceTracker, FaceFit)> {
+        let c = [(region[0] + region[2]) / 2.0, (region[1] + region[3]) / 2.0];
+        let (hw, hh) = ((region[2] - region[0]) / 2.0, (region[3] - region[1]) / 2.0);
+        if !(hw > 4.0 && hh > 4.0) {
+            return None;
+        }
+        // Skin colour: the middle of the region (cheeks and nose).
+        let win = Window::around(frame, c, hw * 0.35, hh * 0.3)?;
+        let mut f: [Vec<f32>; 3] = [vec![], vec![], vec![]];
+        for j in 0..win.h {
+            for i in 0..win.w {
+                let v = features(win.px(frame, i, j));
+                for k in 0..3 {
+                    f[k].push(v[k]);
+                }
+            }
+        }
+        let mut med = [0.0f32; 3];
+        let mut tol = [0.0f32; 3];
+        let floors = [0.12f32, 0.025, 0.025];
+        for k in 0..3 {
+            med[k] = median(&mut f[k]);
+            let mut dev: Vec<f32> = f[k].iter().map(|v| (v - med[k]).abs()).collect();
+            tol[k] = (4.0 * median(&mut dev)).max(floors[k]);
+        }
+        let skin = Skin { med, tol };
+        let canon: [[f64; 2]; N] = LANDMARKS.map(|l| l.2);
+        let guess = FaceFit {
+            landmarks: canon.map(|q| to_layer(c, hh, 0.0, q)),
+            detected: [false; N],
+            outline: vec![],
+            center: c,
+            height: hh.max(hw / ASPECT),
+            roll: 0.0,
+        };
+        let mut t = FaceTracker { skin, prev: guess, expect: canon };
+        let fit = t.fit(frame)?;
+        // A face shows at least two of: left eye, right eye, mouth.
+        let found = [EYE[0][4], EYE[1][4], MOUTH[0]].iter().filter(|i| fit.detected[**i]).count();
+        if found < 2 {
+            return None;
+        }
+        t.accept(&fit);
+        Some((t, fit))
+    }
+
+    /// The fit of the last tracked frame.
+    pub fn last(&self) -> &FaceFit {
+        &self.prev
+    }
+
+    fn accept(&mut self, fit: &FaceFit) {
+        self.expect = fit.landmarks.map(|p| to_face(fit.center, fit.height, fit.roll, p));
+        self.prev = fit.clone();
+    }
+
+    /// Track into the next frame (either direction). `None` when the face is lost.
+    pub fn step(&mut self, frame: &Frame) -> Option<FaceFit> {
+        let fit = self.fit(frame)?;
+        self.accept(&fit);
+        Some(fit)
+    }
+
+    fn fit(&self, frame: &Frame) -> Option<FaceFit> {
+        let pv = &self.prev;
+        let b0 = pv.height;
+        let reach = 1.7 * b0;
+        let win = Window::around(frame, pv.center, reach, reach)?;
+        let (w, h) = (win.w, win.h);
+        let mut skin = vec![false; w * h];
+        let mut luma = vec![0.0f32; w * h];
+        for j in 0..h {
+            for i in 0..w {
+                let px = win.px(frame, i, j);
+                skin[j * w + i] = self.skin.is_skin(px);
+                luma[j * w + i] = features(px)[0];
+            }
+        }
+        // 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?;
+        // Flood fill (4-connected).
+        let mut comp = vec![false; w * h];
+        let mut stack = vec![(si, sj)];
+        comp[sj * w + si] = true;
+        let mut count = 0usize;
+        let (mut mx, mut my) = (0.0, 0.0);
+        while let Some((i, j)) = stack.pop() {
+            count += 1;
+            let p = win.layer(frame, i, j);
+            mx += p[0];
+            my += p[1];
+            let nb = [(i.wrapping_sub(1), j), (i + 1, j), (i, j.wrapping_sub(1)), (i, j + 1)];
+            for (a, b) in nb {
+                if a < w && b < h && !comp[b * w + a] && skin[b * w + a] {
+                    comp[b * w + a] = true;
+                    stack.push((a, b));
+                }
+            }
+        }
+        if (count as f64) < 0.15 * std::f64::consts::PI * ASPECT * b0 * b0 {
+            return None;
+        }
+        let c0 = [mx / count as f64, my / count as f64];
+        // Outline along rays: the farthest component pixel.
+        let rays = 48;
+        let mut radii = vec![0.0f64; rays];
+        for (k, r) in radii.iter_mut().enumerate() {
+            let a = k as f64 / rays as f64 * std::f64::consts::TAU;
+            let d = [a.cos(), a.sin()];
+            let mut t = 0.0;
+            while t < reach * 1.4 {
+                let p = [c0[0] + d[0] * t, c0[1] + d[1] * t];
+                match win.cell(frame, p) {
+                    Some((i, j)) => {
+                        if comp[j * w + i] {
+                            *r = t + 0.5;
+                        }
+                    }
+                    None => break,
+                }
+                t += 0.5;
+            }
+        }
+        // Median of 5, circularly.
+        let sm: Vec<f64> = (0..rays)
+            .map(|k| {
+                let mut v: Vec<f64> = (0..5).map(|d| radii[(k + rays + d - 2) % rays]).collect();
+                v.sort_by(f64::total_cmp);
+                v[2]
+            })
+            .collect();
+        let outline: Vec<[f64; 2]> = (0..rays)
+            .map(|k| {
+                let a = k as f64 / rays as f64 * std::f64::consts::TAU;
+                [c0[0] + a.cos() * sm[k], c0[1] + a.sin() * sm[k]]
+            })
+            .collect();
+        let (center, major, _minor, ang) = polygon_ellipse(&outline)?;
+        if !(major > 0.3 * b0 && major < 2.5 * b0) {
+            return None;
+        }
+        // Roll: the ellipse axis nearest the previous roll (it is ambiguous by 180°).
+        let mut roll = ang;
+        while roll - pv.roll > 90.0 {
+            roll -= 180.0;
+        }
+        while roll - pv.roll < -90.0 {
+            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<Blob> = 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::<f64>() / n, pts.iter().map(|p| p.0[1]).sum::<f64>() / n];
+                    blobs.push(Blob { pts, centroid });
+                }
+            }
+        }
+        // Match components to features near their expected positions.
+        let ex = &self.expect;
+        let mean_of = |ids: &[usize]| {
+            let n = ids.len() as f64;
+            [ids.iter().map(|i| ex[*i][0]).sum::<f64>() / n, ids.iter().map(|i| ex[*i][1]).sum::<f64>() / n]
+        };
+        let mut used = vec![false; blobs.len()];
+        let take = |target: [f64; 2], radius: f64, prefer_large: bool, used: &mut Vec<bool>| -> Option<usize> {
+            let mut best: Option<(usize, f64)> = None;
+            for (k, bl) in blobs.iter().enumerate() {
+                if used[k] {
+                    continue;
+                }
+                let d = (bl.centroid[0] - target[0]).hypot(bl.centroid[1] - target[1]);
+                if d > radius {
+                    continue;
+                }
+                let score = if prefer_large { -(bl.pts.len() as f64) } else { d };
+                if best.is_none_or(|b| score < b.1) {
+                    best = Some((k, score));
+                }
+            }
+            if let Some((k, _)) = best {
+                used[k] = true;
+            }
+            best.map(|b| b.0)
+        };
+        let mut det = [[0.0; 2]; N];
+        let mut wgt = [0.0; N];
+        let set = |i: usize, q: [f64; 2], det: &mut [[f64; 2]; N], wgt: &mut [f64; N]| {
+            det[i] = to_layer(center, b, roll, q);
+            wgt[i] = 1.0;
+        };
+        let eye_ids = [EYE[0].to_vec(), EYE[1].to_vec()];
+        let mouth = take(mean_of(&MOUTH), 0.25, true, &mut used);
+        let eyes = [take(mean_of(&eye_ids[0]), 0.15, false, &mut used), take(mean_of(&eye_ids[1]), 0.15, false, &mut used)];
+        let brows = [take(mean_of(&BROW[0]), 0.15, false, &mut used), take(mean_of(&BROW[1]), 0.15, false, &mut used)];
+        let nostrils = [take(ex[NOSTRIL[0]], 0.07, false, &mut used), take(ex[NOSTRIL[1]], 0.07, false, &mut used)];
+        for side in 0..2 {
+            let sg = if side == 0 { -1.0 } else { 1.0 };
+            if let Some(k) = eyes[side] {
+                let bl = &blobs[k];
+                let e = EYE[side];
+                let inner = bl.extreme([-sg, 0.0]);
+                let outer = bl.extreme([sg, 0.0]);
+                let mid = 0.5 * (inner[0] + outer[0]);
+                let band = 0.15 * (outer[0] - inner[0]).abs();
+                let px = 0.5 / b;
+                set(e[0], [inner[0] - sg * px, inner[1]], &mut det, &mut wgt);
+                set(e[1], [outer[0] + sg * px, outer[1]], &mut det, &mut wgt);
+                let top = bl.lid(mid, band, -1.0);
+                let bot = bl.lid(mid, band, 1.0);
+                set(e[2], [top[0], top[1] - px], &mut det, &mut wgt);
+                set(e[3], [bot[0], bot[1] + px], &mut det, &mut wgt);
+                set(e[4], bl.dark_centroid(), &mut det, &mut wgt);
+            }
+            if let Some(k) = brows[side] {
+                let bl = &blobs[k];
+                let br = BROW[side];
+                let inner = bl.extreme([-sg, 0.0]);
+                let outer = bl.extreme([sg, 0.0]);
+                set(br[0], inner, &mut det, &mut wgt);
+                set(br[2], outer, &mut det, &mut wgt);
+                // The middle: the brow's centre line at its mid column.
+                let mid = 0.5 * (inner[0] + outer[0]);
+                let band = 0.1 * (outer[0] - inner[0]).abs().max(0.05);
+                let col: Vec<f64> = bl.pts.iter().filter(|p| (p.0[0] - mid).abs() <= band).map(|p| p.0[1]).collect();
+                if !col.is_empty() {
+                    set(br[1], [mid, col.iter().sum::<f64>() / col.len() as f64], &mut det, &mut wgt);
+                }
+            }
+            if let Some(k) = nostrils[side] {
+                set(NOSTRIL[side], blobs[k].centroid, &mut det, &mut wgt);
+            }
+        }
+        if let (Some(a), Some(c)) = (nostrils[0], nostrils[1]) {
+            let (p, q) = (blobs[a].centroid, blobs[c].centroid);
+            set(NOSE_TIP, [0.5 * (p[0] + q[0]), 0.5 * (p[1] + q[1]) - 0.06], &mut det, &mut wgt);
+        }
+        if let Some(k) = mouth {
+            let bl = &blobs[k];
+            let l = bl.extreme([-1.0, 0.0]);
+            let r = bl.extreme([1.0, 0.0]);
+            let px = 0.5 / b;
+            set(MOUTH[0], [l[0] - px, l[1]], &mut det, &mut wgt);
+            set(MOUTH[1], [r[0] + px, r[1]], &mut det, &mut wgt);
+            let mid = 0.5 * (l[0] + r[0]);
+            let band = 0.12 * (r[0] - l[0]).abs();
+            let top = bl.lid(mid, band, -1.0);
+            let bot = bl.lid(mid, band, 1.0);
+            set(MOUTH[2], [top[0], top[1] - px], &mut det, &mut wgt);
+            set(MOUTH[3], [bot[0], bot[1] + px], &mut det, &mut wgt);
+        }
+        // Chin and jaw from the outline, along their canonical directions.
+        for i in [CHIN, JAW[0], JAW[1]] {
+            let q = LANDMARKS[i].2;
+            let dir = to_layer([0.0; 2], 1.0, roll, q);
+            let ang = dir[1].atan2(dir[0]);
+            let rad = ray_radius(&outline, center, ang);
+            let len = dir[0].hypot(dir[1]);
+            let p = [center[0] + dir[0] / len * rad, center[1] + dir[1] / len * rad];
+            det[i] = p;
+            wgt[i] = 1.0;
+        }
+        // Shape model: fill the missing landmarks and replace inconsistent detections.
+        let (sn, cs) = roll.to_radians().sin_cos();
+        let init = [b * cs, b * sn, center[0], center[1]];
+        let (model, _) = fit_shape(&det, &wgt, init);
+        let mut w2 = wgt;
+        for i in 0..N {
+            if wgt[i] > 0.0 && (det[i][0] - model[i][0]).hypot(det[i][1] - model[i][1]) > 0.12 * b {
+                w2[i] = 0.0;
+            }
+        }
+        let (model, sim) = if w2 != wgt { fit_shape(&det, &w2, init) } else { (model, init) };
+        let _ = sim;
+        let mut landmarks = model;
+        let mut detected = [false; N];
+        for i in 0..N {
+            if w2[i] > 0.0 {
+                landmarks[i] = det[i];
+                detected[i] = true;
+            }
+        }
+        // Roll from the jaw line (it does not move with yaw); fall back to the ellipse.
+        let jl = landmarks[JAW[0]];
+        let jr = landmarks[JAW[1]];
+        let roll_jaw = (jr[1] - jl[1]).atan2(jr[0] - jl[0]).to_degrees();
+        let roll = if (roll_jaw - roll).abs() < 20.0 { roll_jaw } else { roll };
+        Some(FaceFit { landmarks, detected, outline, center, height: b, roll })
+    }
+}
+
+/// Distance from `c` to the polygon `poly` (star-shaped around `c`) along angle `ang`.
+fn ray_radius(poly: &[[f64; 2]], c: [f64; 2], ang: f64) -> f64 {
+    let d = [ang.cos(), ang.sin()];
+    let mut best = 0.0f64;
+    let n = poly.len();
+    for i in 0..n {
+        let (p, q) = (poly[i], poly[(i + 1) % n]);
+        let e = [q[0] - p[0], q[1] - p[1]];
+        let den = d[0] * e[1] - d[1] * e[0];
+        if den.abs() < 1e-12 {
+            continue;
+        }
+        let w = [p[0] - c[0], p[1] - c[1]];
+        let t = (w[0] * e[1] - w[1] * e[0]) / den;
+        let s = (w[0] * d[1] - w[1] * d[0]) / den;
+        if t > 0.0 && (-1e-9..=1.0 + 1e-9).contains(&s) {
+            best = best.max(t);
+        }
+    }
+    best
+}
+
+// ---------------------------------------------------------------- measurements
+
+/// Face measurements (Extract & Copy Face Measurements): (id, display name).
+pub const MEASUREMENTS: [(&str, &str); 14] = [
+    ("headPositionX", "Head Position X"),
+    ("headPositionY", "Head Position Y"),
+    ("headScale", "Head Scale"),
+    ("headOrientationX", "Head Orientation X"),
+    ("headOrientationY", "Head Orientation Y"),
+    ("headOrientationZ", "Head Orientation Z"),
+    ("leftEyeOpenness", "Left Eye Openness"),
+    ("rightEyeOpenness", "Right Eye Openness"),
+    ("leftEyebrowRaise", "Left Eyebrow Distance from Eye"),
+    ("rightEyebrowRaise", "Right Eyebrow Distance from Eye"),
+    ("mouthOpenness", "Mouth Openness"),
+    ("mouthWidth", "Mouth Width"),
+    ("mouthOffsetX", "Mouth Offset X"),
+    ("mouthOffsetY", "Mouth Offset Y"),
+];
+
+/// Face centre, half-height and roll from the chin and jaw landmarks.
+pub fn face_frame(l: &[[f64; 2]; N]) -> ([f64; 2], f64, f64) {
+    let (ch, jl, jr) = (l[CHIN], l[JAW[0]], l[JAW[1]]);
+    let jm = [0.5 * (jl[0] + jr[0]), 0.5 * (jl[1] + jr[1])];
+    let k = 1.0 - LANDMARKS[JAW[0]].2[1];
+    let c = [ch[0] + (jm[0] - ch[0]) / k, ch[1] + (jm[1] - ch[1]) / k];
+    // Half-height from the chin and from the jaw width (averaged: less noise).
+    let jw = (jr[0] - jl[0]).hypot(jr[1] - jl[1]) / (2.0 * LANDMARKS[JAW[1]].2[0]);
+    let b = (0.5 * ((ch[0] - c[0]).hypot(ch[1] - c[1]) + jw)).max(1e-9);
+    let roll = (jr[1] - jl[1]).atan2(jr[0] - jl[0]).to_degrees();
+    (c, b, roll)
+}
+
+/// The measurements of a frame's landmarks ([`MEASUREMENTS`] order); `reference_height`: the
+/// face half-height that is Head Scale 100 %.
+pub fn measure(l: &[[f64; 2]; N], reference_height: f64) -> [f64; 14] {
+    let (c, b, roll) = face_frame(l);
+    let q = l.map(|p| to_face(c, b, roll, p));
+    let canon = LANDMARKS.map(|x| x.2);
+    // Inner features slide with yaw and pitch (see `synth_shape`).
+    let inner = [EYE[0][4], EYE[1][4], MOUTH[0], MOUTH[1], EYE[0][0], EYE[1][0], EYE[0][1], EYE[1][1]];
+    let n = inner.len() as f64;
+    let du = inner.iter().map(|i| q[*i][0] - canon[*i][0]).sum::<f64>() / n;
+    let dv = inner.iter().map(|i| q[*i][1] - canon[*i][1]).sum::<f64>() / n;
+    let yaw = (du / 0.35).clamp(-1.0, 1.0).asin().to_degrees();
+    let pitch = (dv / 0.3).clamp(-1.0, 1.0).asin().to_degrees();
+    let open = |e: [usize; 5]| {
+        let wdt = (q[e[1]][0] - q[e[0]][0]).abs().max(1e-9);
+        (q[e[3]][1] - q[e[2]][1]).max(0.0) / wdt * 100.0
+    };
+    let brow = |side: usize| {
+        let e = EYE[side];
+        let eye_v = 0.5 * (q[e[2]][1] + q[e[3]][1]);
+        (eye_v - q[BROW[side][1]][1]) * 100.0
+    };
+    let mw = (q[MOUTH[1]][0] - q[MOUTH[0]][0]).abs();
+    let mc = [0.5 * (q[MOUTH[0]][0] + q[MOUTH[1]][0]), 0.5 * (q[MOUTH[2]][1] + q[MOUTH[3]][1])];
+    [
+        c[0],
+        c[1],
+        b / reference_height.max(1e-9) * 100.0,
+        pitch,
+        yaw,
+        roll,
+        open(EYE[0]),
+        open(EYE[1]),
+        brow(0),
+        brow(1),
+        (q[MOUTH[3]][1] - q[MOUTH[2]][1]).max(0.0) / mw.max(1e-9) * 100.0,
+        mw / (2.0 * ASPECT) * 100.0,
+        (mc[0] - q[NOSE_TIP][0]) * 100.0,
+        (mc[1] - q[NOSE_TIP][1]) * 100.0,
+    ]
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    fn frame_params(t: usize) -> FaceParams {
+        let tf = t as f64;
+        FaceParams {
+            center: [160.0 + 30.0 * (tf * 0.15).sin(), 120.0 + 12.0 * (tf * 0.1).cos()],
+            height: 62.0 * (1.0 + 0.08 * (tf * 0.07).sin()),
+            roll: 14.0 * (tf * 0.12).sin(),
+            yaw: 12.0 * (tf * 0.09).sin(),
+            pitch: 6.0 * (tf * 0.11).cos(),
+            eye_open: [0.8 + 0.2 * (tf * 0.3).sin(), 0.8 + 0.2 * (tf * 0.3).sin()],
+            mouth_open: 0.5 + 0.45 * (tf * 0.2).sin(),
+            mouth_width: 1.0 + 0.1 * (tf * 0.13).cos(),
+            brow_raise: [0.4 * (tf * 0.17).sin(), 0.4 * (tf * 0.17).sin()],
+        }
+    }
+
+    #[test]
+    fn shape_model_spans_the_synthetic_shapes() {
+        let m = pdm();
+        assert!(m.modes.len() >= 5 && m.modes.len() <= 14, "{} modes", m.modes.len());
+        // A synthetic shape is reproduced from half its landmarks.
+        let p = FaceParams { yaw: 15.0, mouth_open: 0.8, eye_open: [0.3, 0.9], center: [100.0, 80.0], height: 50.0, roll: 10.0, ..Default::default() };
+        let truth = synth_landmarks(&p);
+        let mut w = [0.0; N];
+        for (i, wi) in w.iter_mut().enumerate() {
+            *wi = if i % 2 == 0 || i >= 19 { 1.0 } else { 0.0 };
+        }
+        let (sn, cs) = 10f64.to_radians().sin_cos();
+        let (fit, _) = fit_shape(&truth, &w, [50.0 * cs, 50.0 * sn, 100.0, 80.0]);
+        let err = (0..N).map(|i| (fit[i][0] - truth[i][0]).hypot(fit[i][1] - truth[i][1])).fold(0.0, f64::max);
+        assert!(err < 0.06 * 50.0, "max error {err}");
+    }
+
+    #[test]
+    fn tracks_a_moving_rotating_synthetic_face() {
+        let (w, h) = (320, 240);
+        let p0 = frame_params(0);
+        let img0 = render_synthetic(w, h, &p0, 3);
+        // The user's mask: a loose box around the face.
+        let region = [p0.center[0] - 60.0, p0.center[1] - 75.0, p0.center[0] + 60.0, p0.center[1] + 75.0];
+        let (mut tr, fit0) = FaceTracker::new(&Frame { img: &img0, offset: [0.0; 2] }, region).expect("face found");
+        let mut fits = vec![fit0];
+        let n = 40;
+        for t in 1..n {
+            let img = render_synthetic(w, h, &frame_params(t), 3);
+            fits.push(tr.step(&Frame { img: &img, offset: [0.0; 2] }).unwrap_or_else(|| panic!("lost at {t}")));
+        }
+        let mut per = [0.0f64; N];
+        let mut worst = 0.0f64;
+        let mut outline_err = 0.0;
+        let mut meas = vec![];
+        for (t, f) in fits.iter().enumerate() {
+            let p = frame_params(t);
+            let truth = synth_landmarks(&p);
+            for i in 0..N {
+                let e = (f.landmarks[i][0] - truth[i][0]).hypot(f.landmarks[i][1] - truth[i][1]) / p.height;
+                per[i] += e / n as f64;
+                worst = worst.max(e);
+            }
+            // Outline: on the face ellipse.
+            let mut oe = 0.0;
+            for q in &f.outline {
+                let u = to_face(p.center, p.height, p.roll, *q);
+                oe += ((u[0] / ASPECT).powi(2) + u[1] * u[1]).sqrt() - 1.0;
+            }
+            outline_err += (oe / f.outline.len() as f64).abs() / n as f64;
+            meas.push((p, measure(&f.landmarks, fits[0].height)));
+            assert!((f.roll - p.roll).abs() < 3.0, "frame {t}: roll {} vs {}", f.roll, p.roll);
+        }
+        let mean = per.iter().sum::<f64>() / N as f64;
+        eprintln!("face landmarks: mean error {mean:.4}, worst {worst:.4} face heights; outline {outline_err:.4}");
+        for (i, e) in per.iter().enumerate() {
+            assert!(*e < 0.04, "{}: mean error {e:.3} (face heights)", LANDMARKS[i].0);
+        }
+        assert!(mean < 0.02, "mean landmark error {mean:.4}");
+        assert!(worst < 0.08, "worst landmark error {worst:.3}");
+        assert!(outline_err < 0.04, "outline error {outline_err:.3}");
+        // Measurements follow the generator: mouth openness and yaw correlate strongly.
+        let corr = |a: &[f64], b: &[f64]| {
+            let n = a.len() as f64;
+            let (ma, mb) = (a.iter().sum::<f64>() / n, b.iter().sum::<f64>() / n);
+            let (mut sab, mut saa, mut sbb) = (0.0, 0.0, 0.0);
+            for (x, y) in a.iter().zip(b) {
+                sab += (x - ma) * (y - mb);
+                saa += (x - ma).powi(2);
+                sbb += (y - mb).powi(2);
+            }
+            sab / (saa * sbb).sqrt()
+        };
+        let mo: Vec<f64> = meas.iter().map(|m| m.1[10]).collect();
+        let mt: Vec<f64> = meas.iter().map(|m| m.0.mouth_open).collect();
+        assert!(corr(&mo, &mt) > 0.95, "mouth openness r = {}", corr(&mo, &mt));
+        let yw: Vec<f64> = meas.iter().map(|m| m.1[4]).collect();
+        let yt: Vec<f64> = meas.iter().map(|m| m.0.yaw).collect();
+        assert!(corr(&yw, &yt) > 0.9, "yaw r = {}", corr(&yw, &yt));
+        let rl: Vec<f64> = meas.iter().map(|m| m.1[5]).collect();
+        let rt: Vec<f64> = meas.iter().map(|m| m.0.roll).collect();
+        assert!(corr(&rl, &rt) > 0.98, "roll r = {}", corr(&rl, &rt));
+        let sc: Vec<f64> = meas.iter().map(|m| m.1[2]).collect();
+        let st: Vec<f64> = meas.iter().map(|m| m.0.height).collect();
+        assert!(corr(&sc, &st) > 0.95, "scale r = {}", corr(&sc, &st));
+    }
+
+    #[test]
+    fn measurements_of_the_synthetic_model_are_exact() {
+        let p = FaceParams { center: [50.0, 60.0], height: 40.0, roll: 20.0, yaw: 10.0, pitch: -5.0, ..Default::default() };
+        let m = measure(&synth_landmarks(&p), 40.0);
+        assert!((m[0] - 50.0).abs() < 1e-6 && (m[1] - 60.0).abs() < 1e-6, "{m:?}");
+        assert!((m[2] - 100.0).abs() < 1e-6);
+        assert!((m[5] - 20.0).abs() < 1e-6);
+        assert!((m[4] - 10.0).abs() < 0.5 && (m[3] + 5.0).abs() < 0.5, "{m:?}");
+    }
+}
diff --git a/crates/track/src/lib.rs b/crates/track/src/lib.rs
index 1c6e6a15..f4a7b88f 100644
--- a/crates/track/src/lib.rs
+++ b/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::*;
diff --git a/crates/track/src/solve.rs b/crates/track/src/solve.rs
index 6c267721..693c52f7 100644
--- a/crates/track/src/solve.rs
+++ b/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]]);
 
diff --git a/crates/track/src/stabilize.rs b/crates/track/src/stabilize.rs
index 1abe79f9..9769fcbd 100644
--- a/crates/track/src/stabilize.rs
+++ b/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);
+    }
 }
diff --git a/crates/track/src/subspace.rs b/crates/track/src/subspace.rs
new file mode 100644
index 00000000..339f9708
--- /dev/null
+++ b/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}");
+    }
+}
diff --git a/crates/ui-egui/src/panels/camera_tracker_ui.rs b/crates/ui-egui/src/panels/camera_tracker_ui.rs
index 41fb7204..c9c40f05 100644
--- a/crates/ui-egui/src/panels/camera_tracker_ui.rs
+++ b/crates/ui-egui/src/panels/camera_tracker_ui.rs
@@ -38,6 +38,8 @@ struct Active {
     show_3d: bool,
     point_size: f64,
     target_size: f64,
+    /// Undistort Footage is on: points project without the lens distortion.
+    undistort: bool,
     /// Layer pixels → comp pixels.
     l2c: Mat3,
     frame: usize,
@@ -66,6 +68,7 @@ fn active(s: &Session, comp: &Comp, ectx: &EvalCtx, l2c: &dyn Fn(&EvalCtx, &Laye
             let t = tracks.as_ref()?;
             (t.frame_duration > 0.0).then(|| ((lt - t.start) / t.frame_duration).round().clamp(0.0, t.frames.saturating_sub(1) as f64) as usize)
         });
+        let undistort = solve.as_ref().is_some_and(|sv| ct::undistorts(&params, sv));
         return Some(Active {
             layer: *lid,
             uid: *uid,
@@ -74,6 +77,7 @@ fn active(s: &Session, comp: &Comp, ectx: &EvalCtx, l2c: &dyn Fn(&EvalCtx, &Laye
             show_3d: val("showTrackPoints").map(|v| v.as_enum()).unwrap_or(1) == 1,
             point_size: val("trackPointSize").map(|v| v.as_f64()).unwrap_or(100.0),
             target_size: val("targetSize").map(|v| v.as_f64()).unwrap_or(100.0),
+            undistort,
             l2c: l2c(ectx, l),
             frame: frame.unwrap_or(0),
         });
@@ -98,18 +102,20 @@ fn l2c_scale(m: &Mat3) -> f32 {
 
 /// Draw the target disc (concentric rings in perspective) of `t` on frame `k`.
 fn draw_target(p: &egui::Painter, a: &Active, map: &ViewerMap, sv: &CameraSolve, t: &Target) {
-    let Some(cam) = sv.frames.get(a.frame) else { return };
+    if sv.frames.get(a.frame).is_none() {
+        return;
+    }
     let r = t.size * 0.5 * a.target_size / 100.0;
     for (i, frac) in [1.0, 0.66, 0.33].into_iter().enumerate() {
         let ring = ct::target_circle(t.center, t.normal, r * frac, 48);
-        let pts: Vec<Pos2> = ring.iter().filter_map(|x| cam.project(sv.size, *x)).map(|q| to_screen(a, map, q)).collect();
+        let pts: Vec<Pos2> = ring.iter().filter_map(|x| sv.project(a.frame, *x, a.undistort)).map(|q| to_screen(a, map, q)).collect();
         if pts.len() < 3 {
             return;
         }
         let fill = if i % 2 == 0 { Color32::from_rgba_unmultiplied(0xd0, 0x30, 0x30, 70) } else { Color32::from_rgba_unmultiplied(0xff, 0xff, 0xff, 50) };
         p.add(egui::Shape::convex_polygon(pts.clone(), fill, Stroke::new(1.0, Color32::from_rgb(0xe0, 0x40, 0x40))));
     }
-    if let Some(c) = cam.project(sv.size, t.center) {
+    if let Some(c) = sv.project(a.frame, t.center, a.undistort) {
         let c = to_screen(a, map, c);
         p.circle_filled(c, 3.0, Color32::from_rgb(0xff, 0xd0, 0x40));
     }
@@ -168,7 +174,7 @@ pub fn viewer_hook(
     let mut hits: Vec<Hit> = vec![];
     match (&a.solve, a.show_3d) {
         (Some(sv), true) => {
-            let pts = ct::projected(sv, a.frame);
+            let pts = ct::projected_with(sv, a.frame, a.undistort);
             let med = ct::median_depth(&pts);
             for (id, q, z) in pts {
                 let sp = to_screen(&a, map, q);
diff --git a/crates/ui-egui/src/panels/tracker.rs b/crates/ui-egui/src/panels/tracker.rs
index 364e114e..02c59636 100644
--- a/crates/ui-egui/src/panels/tracker.rs
+++ b/crates/ui-egui/src/panels/tracker.rs
@@ -404,6 +404,27 @@ fn mask_mode(app: &mut EffectcraftApp, ui: &mut egui::Ui, ctx: &egui::Context, p
         bx += 34.0;
     }
     y += 34.0;
+    // Face Tracking (Detailed Features): Extract & Copy Face Measurements.
+    let has_points = app
+        .session
+        .active_comp()
+        .and_then(|c| c.layer(layer))
+        .and_then(|l| l.effects())
+        .is_some_and(|fx| fx.groups().any(|g| g.match_id == effectcraft_engine::effects::face_track::POINTS_ID));
+    if method == MaskMethod::FaceDetailed || has_points {
+        let r = Rect::from_min_size(pos2(x0, y), vec2(w.min(260.0), 24.0));
+        if button(app, ui, r, "Extract & Copy Face Measurements", has_points && !running, "tracker.mask.extractFace") {
+            match crate::menus::invoke(app, ctx, "track.extractFaceMeasurements", json!({"layer": layer.0})) {
+                Ok(v) => {
+                    if let Some(text) = v.get("clipboard").and_then(|t| t.as_str()) {
+                        ctx.copy_text(text.to_string());
+                    }
+                }
+                Err(e) => app.ui.status = e,
+            }
+        }
+        y += 32.0;
+    }
     if let Some(pr) = app.session.mask_track_progress() {
         let bar = Rect::from_min_size(pos2(x0, y), vec2(w, 6.0));
         p.rect_filled(bar, 3.0, t.field_bg);
@@ -419,7 +440,12 @@ fn mask_mode(app: &mut EffectcraftApp, ui: &mut egui::Ui, ctx: &egui::Context, p
         );
         ctx.request_repaint_after(std::time::Duration::from_millis(50));
     } else {
-        p.text(pos2(x0, y + 8.0), Align2::LEFT_CENTER, "Tracks the pixels inside the mask and keys its Mask Path.", Tokens::ui(11.0), t.text_faint);
+        let hint = if method.is_face() {
+            "Finds the face inside the mask and keys its outline (and, with Detailed Features, Face Track Points)."
+        } else {
+            "Tracks the pixels inside the mask and keys its Mask Path."
+        };
+        p.text(pos2(x0, y + 8.0), Align2::LEFT_CENTER, hint, Tokens::ui(11.0), t.text_faint);
     }
 }
 
diff --git a/crates/ui-egui/tests/ui_camera_tracker.rs b/crates/ui-egui/tests/ui_camera_tracker.rs
index 45374fba..963ce3b6 100644
--- a/crates/ui-egui/tests/ui_camera_tracker.rs
+++ b/crates/ui-egui/tests/ui_camera_tracker.rs
@@ -41,6 +41,7 @@ fn app() -> (EffectcraftApp, u64, u64) {
             })
             .collect(),
         ground: None,
+        distortion: None,
     };
     s.execute("prop.set", json!({"layer": plate, "path": "effects/#1/solve", "value": solve.to_json()})).unwrap();
     s.camera_pending.clear();
diff --git a/docs/agents.md b/docs/agents.md
index 8689b9da..c25aa399 100644
--- a/docs/agents.md
+++ b/docs/agents.md
@@ -105,11 +105,18 @@ effectcraft-cli run clip.ecproj track.motion '{"layer":"#2"}' \
 
 `track.mask` follows the pixels inside a mask and keys its Mask Path on every frame (vertices and
 tangents move with the fitted motion). `method` is `position`, `positionScale`,
-`positionScaleRotation` (the default, also kept as the Tracker panel's Method), `positionScaleRotationSkew`
-or `perspective`; `direction` is `forward|backward|frameForward|frameBackward` from the current
+`positionScaleRotation` (the default, also kept as the Tracker panel's Method), `positionScaleRotationSkew`,
+`perspective`, `faceOutline` (Face Tracking (Outline Only): the mask, drawn around a face, becomes
+the face outline on every frame) or `faceDetailed` (Face Tracking (Detailed Features): also keys a
+Face Track Points effect with one point per landmark — eyebrows, eyes, pupils, nose tip,
+nostrils, mouth corners and lips, chin, jaw); `direction` is `forward|backward|frameForward|frameBackward` from the current
 time (or `start` / `end` in seconds). Like `track.analyze` it runs in the background in the app
 (progress in `track.status` under `mask`, `track.stop` cancels and keeps what was tracked) and is
-one undo step; `"wait": true` blocks.
+one undo step; `"wait": true` blocks. `track.extractFaceMeasurements {layer?}` (Extract & Copy Face
+Measurements) keys a Face Measurements effect from the Face Track Points (head position, scale and
+orientation X/Y/Z, eye openness, eyebrow distance from eye, mouth openness, width and offset), puts
+those keys on the keyframe clipboard and returns them per frame plus a tab-separated `clipboard`
+text.
 
 `mask.interpolate` (Window ▸ Mask Interpolation ▸ Apply) adds in-between Mask Path keys between
 each pair of selected Mask Path keys (or the existing keys at `times`, in seconds), giving both
@@ -138,7 +145,11 @@ the app re-analyses automatically (headless: run `warp.analyze` again). Settings
 properties under the instance's groups, e.g. `effects/#1/stabilization/result` (0 Smooth Motion,
 1 No Motion), `effects/#1/stabilization/smoothness`, `effects/#1/stabilization/method`,
 `effects/#1/borders/framing` (0 Stabilize Only … 3 Stabilize, Synthesize Edges),
-`effects/#1/borders/autoScale/maximumScale`, `effects/#1/advanced/showTrackPoints`.
+`effects/#1/borders/autoScale/maximumScale`, `effects/#1/advanced/showTrackPoints`. Method 3
+(Subspace Warp, the default) bends a mesh per frame to the smoothed feature trajectories, so
+scenes with depth stabilise where one perspective transform cannot;
+`effects/#1/advanced/rollingShutterRipple` (0 Automatic Reduction, 1 Enhanced Reduction: a finer,
+softer mesh) tunes it.
 
 ```sh
 effectcraft-cli run shaky.ecproj track.warpStabilizer '{"layer":"#1","wait":true}' \
@@ -150,13 +161,17 @@ effectcraft-cli exec warp.status '{"layer":"#1"}' shaky.ecproj --json
 ### 3D Camera Tracker
 
 Animation ▸ Track Camera and the Tracker panel's button (`track.camera {layer?, shotType?:
-fixed|variable|specify, aov?, solveMethod?: auto|typical|flat|tripod, detailed?, wait?}`) apply
+fixed|variable|specify, aov?, solveMethod?: auto|typical|flat|tripod, detailed?, lensDistortion?,
+undistort?, wait?}`) apply
 Effect ▸ Perspective ▸ 3D Camera Tracker (or reuse the layer's) and analyse it in the background:
 "Analyzing in background (step 1 of 2)" tracks features, "Solving camera" solves.
 `camera.analyze {layer?, effect?, wait?}` re-runs it, `camera.cancel` stops it without writing
 anything. `camera.solveStatus` reports `progress`/`banner`, `analyzed`, `solved`, `methodUsed`,
 `averageError` (pixels), `focalLength`, `horizontalAngleOfView`, the current frame's `camera`
-(position, orientation, zoom) and `groundPlane`. `camera.points {time?}` lists the solved points
+(position, orientation, zoom), `groundPlane` and the solved `lensDistortion` (`k1`, `k2`, normalised by
+half the frame diagonal) when Solve Lens Distortion (`effects/#1/advanced/lensDistortion`, also on
+with Detailed Analysis) is on; `effects/#1/advanced/undistort` renders the footage undistorted so
+the solved camera's 3D layers line up. `camera.points {time?}` lists the solved points
 visible now with their `id`, `comp` position, `depth`, `world` position and `error`.
 
 Select points with `camera.selectPoints {points, add?, toggle?}` (the viewer's click / Shift-click /
diff --git a/docs/architecture.md b/docs/architecture.md
index 7c341818..03c472a4 100644
--- a/docs/architecture.md
+++ b/docs/architecture.md
@@ -229,6 +229,13 @@ tracks them with pyramidal Lucas–Kanade and fits the chosen model (position 
 RANSAC; the motion moves the Mask Path's vertices and tangents and is keyed per frame. **Mask
 Interpolation** (`mask.interpolate`) uses `effectcraft-path`'s smart interpolation (arc-length
 vertex insertion, shape-context matching, rigid in-betweens) to key in-between shapes.
+The two **Face Tracking** methods (`effectcraft_track::face`) fit a face inside the mask instead:
+a skin colour model from the first frame, the face outline along rays from the skin component's
+centre (an area-moment ellipse gives centre, size and roll), facial features as the non-skin
+components inside it, and an active-shape point distribution model trained on synthetic face
+shapes (our own generator, no external weights) that fills and checks the landmarks. The outline
+keys the Mask Path; Detailed Features keys a Face Track Points effect, and
+`track.extractFaceMeasurements` derives a keyed Face Measurements effect (and copies its keys).
 
 **Warp Stabilizer** (`effects::warp_stab`, `engine::warp`): `warp.analyze` renders the layer's
 input to the effect (`Renderer::layer_input`: source, masks and the effects above it) for every
@@ -237,7 +244,10 @@ translation / similarity / homography fits) as JSON in the effect's hidden Analy
 with a key of the layer's source, In/Out, start, stretch and Time Remap. `Session::edit` clears
 analyses whose key no longer matches and queues them; the desktop app re-analyses them in the
 background. Rendering derives the stabilization plan (smoothed camera path or No Motion, framing,
-auto-scale; cached per analysis and settings) and warps each frame; Synthesize Edges fills the
+auto-scale; cached per analysis and settings) and warps each frame. Subspace Warp also stores long
+feature trajectories in the analysis, smooths them in a low-rank subspace (Liu et al. 2011) and
+warps each frame with a content-preserving mesh (Liu et al. 2009) fitted to the smoothed
+positions, rendered as a per-pixel mesh lookup; Synthesize Edges fills the
 borders from neighbouring frames read with `EffectHost::self_at`. The effect lives in the
 effects crate, hence the `effects → track` edge.
 
@@ -268,7 +278,10 @@ the essential matrix (normalised 8-point in RANSAC) or a planar homography decom
 resection and triangulation; a sparse Levenberg–Marquardt bundle adjustment (Schur complement over
 the points, Huber loss) with the focal length fixed, shared, or per frame; a log-spaced and
 golden-section focal search for Fixed Angle of View / Variable Zoom; a rotation-only model for tripod
-pans, chosen by Auto Detect when it explains the tracks as well. Tracks and solve are stored as JSON
+pans, chosen by Auto Detect when it explains the tracks as well; with Solve Lens Distortion (or
+Detailed Analysis) radial distortion `k1, k2` is adjusted jointly in a final bundle adjustment on
+the raw tracks and the solve repeated on tracks undistorted with that estimate; Undistort Footage
+renders the input through the inverse model. Tracks and solve are stored as JSON
 in hidden effect parameters with keys: changing the layer's frames clears both, changing Shot Type,
 Angle of View, Solve Method or deleting points only re-solves. The solve's canonical frame maps to
 comp space so the first frame's camera is the default comp camera (or so a chosen ground plane is
diff --git a/docs/effects.md b/docs/effects.md
index 84735286..fd29541b 100644
--- a/docs/effects.md
+++ b/docs/effects.md
@@ -3,11 +3,11 @@
 <!-- Generated from the effect registry (crates/effects/src/catalog.rs); do not edit by hand.
      Regenerate: UPDATE_DOCS=1 cargo test -p effectcraft-effects --lib effects_doc_is_current -->
 
-EffectCraft ships 304 effects with 2446 parameters in total, grouped into the same categories as the Effects & Presets panel. Every effect is our own implementation, written from public behaviour descriptions and standard image-processing literature. Parameter names, order, popup options, units, defaults and ranges follow the reference application so that projects, expressions and muscle memory carry over.
+EffectCraft ships 306 effects with 2488 parameters in total, grouped into the same categories as the Effects & Presets panel. Every effect is our own implementation, written from public behaviour descriptions and standard image-processing literature. Parameter names, order, popup options, units, defaults and ranges follow the reference application so that projects, expressions and muscle memory carry over.
 
 - **GPU**: 59 effects also run on the GPU compositor with identical results.
-- **32**: 304 effects process 32-bit float (HDR, overbright) pixels without clamping.
-- **Status**: 262 implemented in full, 42 partial (what is missing is listed).
+- **32**: 306 effects process 32-bit float (HDR, overbright) pixels without clamping.
+- **Status**: 264 implemented in full, 42 partial (what is missing is listed).
 
 Effects are addressed by id (`ec.<category>.<name>`) or display name in commands, scripts and the control channel, and every parameter by its id (`effects/#1/blurriness`).
 
@@ -302,7 +302,7 @@ Effects are addressed by id (`ec.<category>.<name>`) or display name in commands
 
 | Effect | Id | Params | GPU | 32 | Status |
 |---|---|---:|:-:|:-:|---|
-| 3D Camera Tracker | `ec.perspective.cameratracker` | 10 |  | 32 | Implemented |
+| 3D Camera Tracker | `ec.perspective.cameratracker` | 12 |  | 32 | Implemented |
 | 3D Glasses | `ec.perspective.3dglasses` | 8 |  | 32 | Implemented |
 | Bevel Alpha | `ec.perspective.bevelalpha` | 4 |  | 32 | Implemented |
 | Bevel Edges | `ec.perspective.beveledges` | 4 |  | 32 | Implemented |
@@ -416,6 +416,8 @@ Effects are addressed by id (`ec.<category>.<name>`) or display name in commands
 | CC Overbrights | `ec.utility.ccoverbrights` | 2 |  | 32 | Implemented |
 | Cineon Converter | `ec.utility.cineon` | 7 |  | 32 | Implemented |
 | Color Profile Converter | `ec.utility.colorprofileconverter` | 7 |  | 32 | Implemented |
+| Face Measurements | `ec.utility.facemeasurements` | 14 |  | 32 | Implemented |
+| Face Track Points | `ec.utility.facetrackpoints` | 26 |  | 32 | Implemented |
 | Grow Bounds | `ec.utility.growbounds` | 1 |  | 32 | Implemented |
 | HDR Compander | `ec.utility.hdrcompander` | 3 |  | 32 | Implemented |
 | HDR Highlight Compression | `ec.utility.hdrcompression` | 1 |  | 32 | Implemented |
diff --git a/docs/parity.md b/docs/parity.md
index 3974630f..877ed7ec 100644
--- a/docs/parity.md
+++ b/docs/parity.md
@@ -17,12 +17,12 @@ How close EffectCraft is to After Effects 2026, feature by feature, and how much
 
 By area: Layers 98%, Output 95%, Compositions 95%, Automation 96%, Paint 95%, Text 95%, Import 96%,
 Animation 94%, Masks 94%, Preview 94%, Interface 93%, Shapes 93%, 3D 91%, Audio 90%, Project 90%,
-Tracking 88%, Effects 85%, Web 65%.
+Tracking 96%, Effects 85%, Web 65%.
 
 What is left, in priority order: the web app's depth (threads, storage, audio, non-blocking
 renders); stroke taper/wave and multi-segment dashes; variable mask feather points; camera iris/bokeh and focus-link commands; Render Queue field render, crop/resize and templates;
-approximated effects (Subspace Warp, Key Cleaner) and Liquify's viewer brush; text in imported
-PDF/AI files and Illustrator procset EPS; more codecs; 59 preferences not yet wired; ScriptUI; face tracking; and the "better than After Effects" items (a plug-in API, branching history; GPU particles landed in M12.7).
+approximated effects (Key Cleaner) and Liquify's viewer brush; text in imported
+PDF/AI files and Illustrator procset EPS; more codecs; and the remaining "better than After Effects" items (ScriptUI, the plug-in API and branching history landed in M13.1; face tracking and Subspace Warp in M13.3; GPU particles in M12.7).
 
 The sections below are the original audit (morning of 2 October, ≈ 64%) and its updates, kept for
 history.
@@ -122,7 +122,7 @@ the disk cache in the formats wave; the 3D Camera Tracker in M12.6.)
 | Project | ≈ 68% | 5.5 | auto-save, folder moves, OCIO displays beyond the built-in tone map (Color Engine with OCIO/ACES working spaces, HDR compand/tone mapping, Rec. 2100 PQ/HLG output, Feet + Frames, display colour management, Simulate Output and the locked viewer landed in M7.7; proxies and Interpret Footage fields / pixel aspect / alpha guess landed: PRJ-8, PRJ-3) |
 | Masks & roto | 74% | 5.0 | Roto Brush's learned (3.0) segmentation model (variable-width mask feather points with the Mask Feather tool landed in M13.5; mask tracking and Mask Interpolation landed in M6.6; Roto Brush & Refine Edge with graph-cut segmentation, flow propagation, edge matting, decontamination and Freeze in M6.7) |
 | Preview | 66% | 4.0 | GPU bokeh depth of field, wireframes and Advanced 3D compositing (Classic 3D runs and adjustment layers composite on the GPU since M12.7; persistent disk cache with the blue cache bar landed; region of interest, snapshots, exposure and Fast Previews landed in M0.13) |
-| Tracking | ≈ 85% | 2.5 | face tracking, Subspace Warp's mesh warp, lens distortion in the camera solve (Rolling Shutter Repair landed in M9.11; point tracker, mask tracking, Warp Stabilizer and the 3D Camera Tracker landed in M6.x / M12.5 / M12.6) |
+| Tracking | ≈ 96% | 0.3 | face tracking is a classical (skin model + feature components + shape model) fitter, not a learned detector: profile views and occluded faces are weak; Rolling Shutter Ripple is approximated by Subspace Warp's mesh density (face tracking (Outline Only / Detailed Features with Face Track Points and Extract & Copy Face Measurements), Subspace Warp's content-preserving mesh warp on subspace-smoothed trajectories, and radial lens distortion (k1, k2) in the camera bundle adjustment with Undistort Footage landed in M13.3; Rolling Shutter Repair in M9.11; point tracker, mask tracking, Warp Stabilizer and the 3D Camera Tracker in M6.x / M12.5 / M12.6) |
 | Paint | 0% | 4.0 | Brush, Clone Stamp, Eraser |
 
 ## Effects still missing