883 lines
35 KiB
Rust
883 lines
35 KiB
Rust
//! Video stabilization (the Warp Stabilizer's analysis and stabilization steps), clean-room from
|
||
//! the published approach of 2D stabilizers:
|
||
//!
|
||
//! 1. **Analysis** ([`Analyzer`]): Shi–Tomasi features are tracked between consecutive frames
|
||
//! with pyramidal Lucas–Kanade ([`crate::klt`]), and the dominant (background) motion between
|
||
//! each pair of frames is fitted robustly with RANSAC ([`crate::fit`]) for three models at
|
||
//! once: translation, similarity and homography. The result ([`WarpAnalysis`]) is small and
|
||
//! serialisable; it is stored in the effect instance.
|
||
//! 2. **Stabilization** ([`plan`]): the camera path is smoothed and per-frame corrective
|
||
//! transforms are derived. For *Smooth Motion* each frame is warped to the Gaussian-weighted
|
||
//! average of its neighbours' viewpoints — the motion-smoothing scheme of Matsushita et al.,
|
||
//! "Full-Frame Video Stabilization with Motion Inpainting" (PAMI 2006) — so the low-frequency
|
||
//! camera move stays and the jitter goes; *No Motion* maps every frame onto one reference
|
||
//! frame. Framing then crops to the region valid on every frame (a centred rectangle inside
|
||
//! every warped frame outline) and optionally scales it up to fill the frame; when that would
|
||
//! exceed Maximum Scale, the correction is relaxed towards identity on the offending frames —
|
||
//! trading smoothness for less crop, the constraint idea of Grundmann et al., "Auto-Directed
|
||
//! 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::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)]
|
||
#[serde(rename_all = "camelCase", default)]
|
||
pub struct FrameMotion {
|
||
/// Translation `[tx, ty]`.
|
||
pub t: [f64; 2],
|
||
/// Similarity `[a, b, tx, ty]` (`x' = a x − b y + tx`, `y' = b x + a y + ty`).
|
||
pub s: [f64; 4],
|
||
/// Homography, row-major without `h33 = 1`.
|
||
pub h: [f64; 8],
|
||
/// Correspondences that agreed with the homography / tracked in total.
|
||
pub inliers: u32,
|
||
pub features: u32,
|
||
/// A subset of this frame's tracked features (layer pixels), for Show Track Points.
|
||
#[serde(skip_serializing_if = "Vec::is_empty")]
|
||
pub points: Vec<[f32; 2]>,
|
||
}
|
||
|
||
impl Default for FrameMotion {
|
||
fn default() -> Self {
|
||
FrameMotion { t: [0.0; 2], s: [1.0, 0.0, 0.0, 0.0], h: [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0], inliers: 0, features: 0, points: vec![] }
|
||
}
|
||
}
|
||
|
||
impl FrameMotion {
|
||
/// The motion under a stabilization method.
|
||
pub fn motion(&self, method: Method) -> Homography {
|
||
match method {
|
||
Method::Position => Homography::translation(self.t),
|
||
Method::Similarity => {
|
||
let s = self.s;
|
||
Homography([[s[0], -s[1], s[2]], [s[1], s[0], s[3]], [0.0, 0.0, 1.0]])
|
||
}
|
||
Method::Perspective | Method::SubspaceWarp => {
|
||
let h = self.h;
|
||
Homography([[h[0], h[1], h[2]], [h[3], h[4], h[5]], [h[6], h[7], 1.0]])
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
/// The stored analysis of a clip.
|
||
#[derive(Clone, Debug, Default, PartialEq, Serialize, Deserialize)]
|
||
#[serde(rename_all = "camelCase", default)]
|
||
pub struct WarpAnalysis {
|
||
pub version: u32,
|
||
/// Layer time (seconds) of frame 0, and the frame duration.
|
||
pub start: f64,
|
||
pub frame_duration: f64,
|
||
/// Layer size (pixels).
|
||
pub size: [f64; 2],
|
||
/// Detailed Analysis was on.
|
||
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 {
|
||
pub fn is_empty(&self) -> bool {
|
||
self.frames.is_empty()
|
||
}
|
||
/// Frame index at a layer time (held at the ends).
|
||
pub fn frame_at(&self, layer_time: f64) -> Option<usize> {
|
||
if self.frames.is_empty() || self.frame_duration <= 0.0 {
|
||
return None;
|
||
}
|
||
let i = ((layer_time - self.start) / self.frame_duration).round();
|
||
Some(i.clamp(0.0, (self.frames.len() - 1) as f64) as usize)
|
||
}
|
||
pub fn to_json(&self) -> String {
|
||
serde_json::to_string(self).unwrap_or_default()
|
||
}
|
||
pub fn from_json(s: &str) -> Option<WarpAnalysis> {
|
||
if s.trim().is_empty() {
|
||
return None;
|
||
}
|
||
serde_json::from_str(s).ok()
|
||
}
|
||
}
|
||
|
||
/// Analysis options.
|
||
#[derive(Clone, Copy, Debug, Default)]
|
||
pub struct AnalyzeOpts {
|
||
/// Detailed Analysis: twice the resolution and more features.
|
||
pub detailed: bool,
|
||
}
|
||
|
||
/// Incremental analysis: push the clip's frames in order.
|
||
pub struct Analyzer {
|
||
opts: AnalyzeOpts,
|
||
prev: Option<GrayPyramid>,
|
||
frames: Vec<FrameMotion>,
|
||
size: [f64; 2],
|
||
trajectories: TrackAnalyzer,
|
||
}
|
||
|
||
fn round_pts(p: &[[f64; 2]], n: usize) -> Vec<[f32; 2]> {
|
||
let step = (p.len() / n.max(1)).max(1);
|
||
p.iter().step_by(step).take(n).map(|q| [((q[0] * 10.0).round() / 10.0) as f32, ((q[1] * 10.0).round() / 10.0) as f32]).collect()
|
||
}
|
||
|
||
impl Analyzer {
|
||
pub fn new(size: [f64; 2], opts: AnalyzeOpts) -> Analyzer {
|
||
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 {
|
||
self.frames.len()
|
||
}
|
||
pub fn is_empty(&self) -> bool {
|
||
self.frames.is_empty()
|
||
}
|
||
|
||
/// Analyse the next frame.
|
||
pub fn push(&mut self, frame: &Frame) {
|
||
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()
|
||
&& prev.height() == pyr.height()
|
||
&& prev.offset == pyr.offset
|
||
{
|
||
let md = (prev.width().max(prev.height()) as f64 / 40.0).clamp(4.0, 16.0);
|
||
let feats = good_features(prev, &CornerOpts { max_features: max_feat, min_distance: md, quality: 0.01, window: 2, border: 4 }, None);
|
||
let q = track(prev, &pyr, &feats, None, &LkOpts { radius: 5, fb_max: 0.75, ..Default::default() });
|
||
let (src, dst): (Vec<[f64; 2]>, Vec<[f64; 2]>) = feats.iter().zip(&q).filter_map(|(p, q)| q.map(|q| (prev.to_layer(*p), pyr.to_layer(q)))).unzip();
|
||
fm.features = src.len() as u32;
|
||
let thr = 1.0 * factor as f64;
|
||
let seed = self.frames.len() as u64 + 1;
|
||
if let Some(f) = ransac(Model::Translation, &src, &dst, thr * 1.5, 200, seed) {
|
||
let m = f.h.0;
|
||
fm.t = [m[0][2], m[1][2]];
|
||
}
|
||
if let Some(f) = ransac(Model::Similarity, &src, &dst, thr * 1.5, 300, seed) {
|
||
let m = f.h.0;
|
||
fm.s = [m[0][0], m[1][0], m[0][2], m[1][2]];
|
||
}
|
||
if let Some(f) = ransac(Model::Homography, &src, &dst, thr, 500, seed) {
|
||
let m = f.h.normalized().0;
|
||
// Reject wild projective fits (few features): keep the similarity then.
|
||
let sane = m[2][0].abs() < 0.01 && m[2][1].abs() < 0.01;
|
||
fm.h = if sane {
|
||
[m[0][0], m[0][1], m[0][2], m[1][0], m[1][1], m[1][2], m[2][0], m[2][1]]
|
||
} else {
|
||
[fm.s[0], -fm.s[1], fm.s[2], fm.s[1], fm.s[0], fm.s[3], 0.0, 0.0]
|
||
};
|
||
fm.inliers = f.inlier_count() as u32;
|
||
let inl: Vec<[f64; 2]> = dst.iter().zip(&f.inliers).filter(|(_, b)| **b).map(|(p, _)| *p).collect();
|
||
fm.points = round_pts(&inl, 48);
|
||
}
|
||
}
|
||
self.frames.push(fm);
|
||
self.prev = Some(pyr);
|
||
}
|
||
|
||
pub fn finish(self, start: f64, frame_duration: f64) -> WarpAnalysis {
|
||
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 }
|
||
}
|
||
}
|
||
|
||
// ---------------------------------------------------------------- stabilization
|
||
|
||
/// Stabilization ▸ Result.
|
||
#[derive(Clone, Copy, Debug, Default, PartialEq, Eq, Hash)]
|
||
pub enum StabResult {
|
||
#[default]
|
||
SmoothMotion,
|
||
NoMotion,
|
||
}
|
||
|
||
/// Stabilization ▸ Method.
|
||
#[derive(Clone, Copy, Debug, Default, PartialEq, Eq, Hash)]
|
||
pub enum Method {
|
||
Position,
|
||
/// Position, Scale, Rotation.
|
||
Similarity,
|
||
Perspective,
|
||
/// 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 {
|
||
StabilizeOnly,
|
||
StabilizeCrop,
|
||
#[default]
|
||
StabilizeCropAutoScale,
|
||
SynthesizeEdges,
|
||
}
|
||
|
||
/// Effect parameters that shape the stabilization.
|
||
#[derive(Clone, Copy, Debug, PartialEq)]
|
||
pub struct StabSettings {
|
||
pub result: StabResult,
|
||
/// Smoothness, percent (50 = default).
|
||
pub smoothness: f64,
|
||
pub method: Method,
|
||
pub preserve_scale: bool,
|
||
pub framing: Framing,
|
||
/// Maximum Scale, percent.
|
||
pub max_scale: f64,
|
||
/// Action Safe Margin, percent of the frame on each side.
|
||
pub action_safe: f64,
|
||
/// Additional Scale, percent.
|
||
pub additional_scale: f64,
|
||
/// Crop Less ↔ Smooth More (−100…100).
|
||
pub crop_less_smooth_more: f64,
|
||
/// Frames per second (for Smoothness).
|
||
pub fps: f64,
|
||
pub ripple: Ripple,
|
||
}
|
||
|
||
impl Default for StabSettings {
|
||
fn default() -> Self {
|
||
StabSettings {
|
||
result: StabResult::SmoothMotion,
|
||
smoothness: 50.0,
|
||
method: Method::SubspaceWarp,
|
||
preserve_scale: false,
|
||
framing: Framing::StabilizeCropAutoScale,
|
||
max_scale: 150.0,
|
||
action_safe: 0.0,
|
||
additional_scale: 100.0,
|
||
crop_less_smooth_more: 0.0,
|
||
fps: 30.0,
|
||
ripple: Ripple::Automatic,
|
||
}
|
||
}
|
||
}
|
||
|
||
impl StabSettings {
|
||
/// Gaussian smoothing sigma in frames.
|
||
pub fn sigma(&self) -> f64 {
|
||
let base = self.smoothness.max(0.0) / 100.0 * self.fps.max(1.0) * 0.6;
|
||
base * 2f64.powf(self.crop_less_smooth_more.clamp(-100.0, 100.0) / 100.0)
|
||
}
|
||
}
|
||
|
||
/// 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. 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]>,
|
||
/// Auto-scale factor (1 = none).
|
||
pub auto_scale: f64,
|
||
/// Fraction of the frame (per side length) valid on every frame before scaling.
|
||
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 {
|
||
let (a, b) = (q[i], q[(i + 1) % 4]);
|
||
let c = (b[0] - a[0]) * (p[1] - a[1]) - (b[1] - a[1]) * (p[0] - a[0]);
|
||
if c.abs() < 1e-12 {
|
||
continue;
|
||
}
|
||
if sign == 0.0 {
|
||
sign = c.signum();
|
||
} else if c.signum() != sign {
|
||
return false;
|
||
}
|
||
}
|
||
true
|
||
}
|
||
|
||
/// Largest `s ∈ [0, 1]` such that the centred `s·W × s·H` rectangle lies inside the frame
|
||
/// outline mapped by `b`.
|
||
pub fn valid_fraction(b: &Homography, size: [f64; 2]) -> f64 {
|
||
let [w, h] = size;
|
||
let quad = [b.apply([0.0, 0.0]), b.apply([w, 0.0]), b.apply([w, h]), b.apply([0.0, h])];
|
||
// A fold or a point at infinity: nothing is valid.
|
||
let m = b.normalized().0;
|
||
for p in [[0.0, 0.0], [w, 0.0], [w, h], [0.0, h]] {
|
||
if m[2][0] * p[0] + m[2][1] * p[1] + m[2][2] <= 1e-9 {
|
||
return 0.0;
|
||
}
|
||
}
|
||
let c = [w / 2.0, h / 2.0];
|
||
let fits = |s: f64| {
|
||
let (hw, hh) = (s * w / 2.0, s * h / 2.0);
|
||
[[c[0] - hw, c[1] - hh], [c[0] + hw, c[1] - hh], [c[0] + hw, c[1] + hh], [c[0] - hw, c[1] + hh]].iter().all(|p| convex_contains(&quad, *p))
|
||
};
|
||
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..30 {
|
||
let mid = 0.5 * (lo + hi);
|
||
if fits(mid) {
|
||
lo = mid;
|
||
} else {
|
||
hi = mid;
|
||
}
|
||
}
|
||
lo
|
||
}
|
||
|
||
fn remove_scale(b: &Homography, c: [f64; 2]) -> Homography {
|
||
// Local scale at the frame centre from the Jacobian of the map.
|
||
let e = 1.0;
|
||
let p0 = b.apply(c);
|
||
let px = b.apply([c[0] + e, c[1]]);
|
||
let py = b.apply([c[0], c[1] + e]);
|
||
let det = ((px[0] - p0[0]) * (py[1] - p0[1]) - (px[1] - p0[1]) * (py[0] - p0[0])).abs();
|
||
let s = det.sqrt();
|
||
if s < 1e-9 {
|
||
return *b;
|
||
}
|
||
Homography::scale_about(1.0 / s, p0).then_after(b)
|
||
}
|
||
|
||
/// Corrective transforms (source → stabilized, before framing) for every frame.
|
||
pub fn corrections(a: &WarpAnalysis, s: &StabSettings) -> Vec<Homography> {
|
||
let n = a.frames.len();
|
||
if n == 0 {
|
||
return vec![];
|
||
}
|
||
let method = s.method;
|
||
let inter: Vec<Homography> = a.frames.iter().map(|f| f.motion(method)).collect();
|
||
let inv: Vec<Homography> = inter.iter().map(|h| h.inverse().unwrap_or(Homography::IDENTITY)).collect();
|
||
let c = [a.size[0] / 2.0, a.size[1] / 2.0];
|
||
let mut out = Vec::with_capacity(n);
|
||
match s.result {
|
||
StabResult::NoMotion => {
|
||
// Cumulative path C_t (frame 0 → frame t); every frame maps onto the middle frame.
|
||
let mut cum = Vec::with_capacity(n);
|
||
let mut cur = Homography::IDENTITY;
|
||
for (k, h) in inter.iter().enumerate() {
|
||
if k > 0 {
|
||
cur = h.then_after(&cur);
|
||
}
|
||
cum.push(cur);
|
||
}
|
||
let r = n / 2;
|
||
for c_t in &cum {
|
||
let ci = c_t.inverse().unwrap_or(Homography::IDENTITY);
|
||
out.push(cum[r].then_after(&ci));
|
||
}
|
||
}
|
||
StabResult::SmoothMotion => {
|
||
let sigma = s.sigma();
|
||
let rad = ((3.0 * sigma).ceil() as usize).min(n.saturating_sub(1));
|
||
for t in 0..n {
|
||
if sigma < 0.3 || rad == 0 {
|
||
out.push(Homography::IDENTITY);
|
||
continue;
|
||
}
|
||
let mut acc = [[0.0; 3]; 3];
|
||
let mut wsum = 0.0;
|
||
let mut add = |h: &Homography, d: f64| {
|
||
let w = (-(d * d) / (2.0 * sigma * sigma)).exp();
|
||
let m = h.normalized().0;
|
||
for i in 0..3 {
|
||
for j in 0..3 {
|
||
acc[i][j] += w * m[i][j];
|
||
}
|
||
}
|
||
wsum += w;
|
||
};
|
||
add(&Homography::IDENTITY, 0.0);
|
||
// Forward: T_t^{i+1} = inter[i+1] ∘ T_t^i.
|
||
let mut tf = Homography::IDENTITY;
|
||
for i in t + 1..=(t + rad).min(n - 1) {
|
||
tf = inter[i].then_after(&tf);
|
||
add(&tf, (i - t) as f64);
|
||
}
|
||
// Backward: T_t^{i-1} = inter[i]⁻¹ ∘ T_t^i.
|
||
let mut tb = Homography::IDENTITY;
|
||
for i in (t.saturating_sub(rad)..t).rev() {
|
||
tb = inv[i + 1].then_after(&tb);
|
||
add(&tb, (t - i) as f64);
|
||
}
|
||
let m = acc.map(|r| r.map(|v| v / wsum));
|
||
out.push(Homography(m).normalized());
|
||
}
|
||
}
|
||
}
|
||
if s.preserve_scale && method != Method::Position {
|
||
out = out.iter().map(|b| remove_scale(b, c)).collect();
|
||
}
|
||
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 {
|
||
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: &[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);
|
||
let need_cover = 1.0 - 2.0 * margin;
|
||
let max_scale = (s.max_scale / 100.0).max(1.0);
|
||
if s.framing == Framing::StabilizeCropAutoScale && vmin * max_scale < need_cover && !b.is_empty() {
|
||
// Relax the correction on frames that need more than Maximum Scale.
|
||
let target = need_cover / max_scale;
|
||
let n = b.len();
|
||
let lam: Vec<f64> = (0..n)
|
||
.into_par_iter()
|
||
.map(|t| {
|
||
if fr[t] >= target {
|
||
return 1.0;
|
||
}
|
||
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);
|
||
let minf: Vec<f64> = (0..n).map(|t| lam[t.saturating_sub(r)..=(t + r).min(n - 1)].iter().copied().fold(1.0, f64::min)).collect();
|
||
let sig = r as f64 / 2.0;
|
||
let smooth: Vec<f64> = (0..n)
|
||
.map(|t| {
|
||
let (mut acc, mut ws) = (0.0, 0.0);
|
||
for i in t.saturating_sub(r)..=(t + r).min(n - 1) {
|
||
let d = i as f64 - t as f64;
|
||
let w = (-(d * d) / (2.0 * sig * sig)).exp();
|
||
acc += w * minf[i];
|
||
ws += w;
|
||
}
|
||
(acc / ws).min(lam[t])
|
||
})
|
||
.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);
|
||
}
|
||
let rect = |k: f64| [c[0] - k * size[0] / 2.0, c[1] - k * size[1] / 2.0, c[0] + k * size[0] / 2.0, c[1] + k * size[1] / 2.0];
|
||
let (auto, crop) = match s.framing {
|
||
Framing::StabilizeOnly | Framing::SynthesizeEdges => (1.0, None),
|
||
Framing::StabilizeCrop => (1.0, Some(rect(vmin * add))),
|
||
Framing::StabilizeCropAutoScale => {
|
||
let k = if vmin > 1e-6 { (need_cover / vmin).clamp(1.0, max_scale) } else { max_scale };
|
||
let vis = vmin * k * add;
|
||
(k, if vis >= 1.0 - 1e-9 { None } else { Some(rect(vis)) })
|
||
}
|
||
};
|
||
let z = Homography::scale_about(auto, c);
|
||
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]>) {
|
||
let pos = |k: usize| {
|
||
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![], tracks: vec![] };
|
||
let path: Vec<[f64; 2]> = (0..n).map(pos).collect();
|
||
for k in 0..n {
|
||
let mut f = FrameMotion::default();
|
||
if k > 0 {
|
||
// Scene points move opposite to the camera.
|
||
let d = [path[k - 1][0] - path[k][0], path[k - 1][1] - path[k][1]];
|
||
f.t = d;
|
||
f.s = [1.0, 0.0, d[0], d[1]];
|
||
f.h = [1.0, 0.0, d[0], 0.0, 1.0, d[1], 0.0, 0.0];
|
||
}
|
||
a.frames.push(f);
|
||
}
|
||
(a, path)
|
||
}
|
||
|
||
fn jitter(p: &[[f64; 2]]) -> f64 {
|
||
// Mean magnitude of the second difference (acceleration).
|
||
let mut s = 0.0;
|
||
for k in 1..p.len() - 1 {
|
||
s += (p[k + 1][0] - 2.0 * p[k][0] + p[k - 1][0]).hypot(p[k + 1][1] - 2.0 * p[k][1] + p[k - 1][1]);
|
||
}
|
||
s / (p.len() - 2) as f64
|
||
}
|
||
|
||
#[test]
|
||
fn smoothing_removes_jitter_and_no_motion_locks() {
|
||
let (a, path) = shaky(90);
|
||
// A scene point at (320, 180) on frame 0 appears at (320, 180) - path[k] on frame k.
|
||
let scene = |k: usize| [320.0 - path[k][0] + path[0][0], 180.0 - path[k][1] + path[0][1]];
|
||
let before: Vec<[f64; 2]> = (0..90).map(scene).collect();
|
||
for method in [Method::Position, Method::Similarity, Method::Perspective] {
|
||
let s = StabSettings { method, framing: Framing::StabilizeOnly, ..Default::default() };
|
||
let p = plan(&a, &s);
|
||
let after: Vec<[f64; 2]> = (0..90).map(|k| p.warps[k].apply(scene(k))).collect();
|
||
let (j0, j1) = (jitter(&before), jitter(&after));
|
||
assert!(j1 < 0.2 * j0, "{method:?}: jitter {j0} → {j1}");
|
||
}
|
||
let s = StabSettings { result: StabResult::NoMotion, framing: Framing::StabilizeOnly, ..Default::default() };
|
||
let p = plan(&a, &s);
|
||
let fixed: Vec<[f64; 2]> = (0..90).map(|k| p.warps[k].apply(scene(k))).collect();
|
||
for q in &fixed {
|
||
assert!((q[0] - fixed[45][0]).abs() < 1e-6 && (q[1] - fixed[45][1]).abs() < 1e-6);
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn framing_crops_and_scales() {
|
||
let (a, _) = shaky(60);
|
||
let p = plan(&a, &StabSettings { framing: Framing::StabilizeCrop, ..Default::default() });
|
||
let crop = p.crop.unwrap();
|
||
assert!(p.auto_scale == 1.0 && crop[0] > 0.0 && crop[2] < 640.0);
|
||
// Every frame covers the crop rectangle.
|
||
for w in &p.warps {
|
||
let inv = w.inverse().unwrap();
|
||
for q in [[crop[0], crop[1]], [crop[2], crop[1]], [crop[2], crop[3]], [crop[0], crop[3]]] {
|
||
let s = inv.apply(q);
|
||
assert!(s[0] >= -1e-6 && s[1] >= -1e-6 && s[0] <= 640.0 + 1e-6 && s[1] <= 360.0 + 1e-6, "{s:?}");
|
||
}
|
||
}
|
||
let p = plan(&a, &StabSettings::default());
|
||
assert!(p.auto_scale > 1.0 && p.crop.is_none(), "{} {:?}", p.auto_scale, p.crop);
|
||
for w in &p.warps {
|
||
let inv = w.inverse().unwrap();
|
||
for q in [[0.0, 0.0], [640.0, 0.0], [640.0, 360.0], [0.0, 360.0]] {
|
||
let s = inv.apply(q);
|
||
assert!(s[0] >= -1e-3 && s[1] >= -1e-3 && s[0] <= 640.001 && s[1] <= 360.001, "{s:?}");
|
||
}
|
||
}
|
||
// A tight Maximum Scale relaxes the correction instead of showing borders.
|
||
let p = plan(&a, &StabSettings { max_scale: 101.0, ..Default::default() });
|
||
assert!(p.auto_scale <= 1.01 + 1e-9);
|
||
assert!(p.valid_fraction * p.auto_scale >= 0.999 || p.crop.is_some());
|
||
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);
|
||
}
|
||
}
|