Goal

Estimate a dense, piecewise-smooth optical flow field $(u(s), v(s))$ at every pixel site $s$, or recover one or more parametric (affine) motion models from a two-frame image sequence, in the presence of multiple competing motions. Input: two consecutive grayscale frames $I(\mathbf{x}, t)$ and $I(\mathbf{x}, t{+}1)$, where $\mathbf{x} = (x, y)$. Output: either (a) a dense flow field with flow in pixels per frame, or (b) a set of affine parameter vectors $\mathbf{a}$ together with a per-pixel outlier indicator. Both outputs are obtained by replacing the quadratic (L2) penalty in the data conservation and spatial smoothness terms with a redescending robust M-estimator, solved via Iteratively Reweighted Least Squares (IRLS) within a Graduated Non-Convexity (GNC) continuation schedule.

Algorithm

Let $I_x$, $I_y$, $I_t$ denote the partial spatial and temporal image derivatives at site $s$. Let $u_s$, $v_s$ denote the horizontal and vertical flow at $s$. Let $N(s)$ denote the 4-connected spatial neighbourhood of $s$. Let $\lambda > 0$ weight the smoothness term relative to the data term. Let $\rho(\cdot, \sigma)$ denote a robust penalty function with scale parameter $\sigma$, and $\psi(x, \sigma) = \partial \rho / \partial x$ its influence function.

Two penalty functions are used throughout.

Piecewise-smooth (dense) variant

The robust data energy penalises violation of the brightness-constancy constraint at each pixel.

The robust smoothness energy penalises large flow differences between neighbouring sites.

Minimisation proceeds by SOR within a GNC continuation schedule. The SOR update for $u_s$ at iteration $n$ is

$u_s^{n+1} = u_s^n - \frac{\omega}{T(u_s)} \cdot \frac{\partial E}{\partial u_s},$

where $\omega \in (0, 2)$ is the overrelaxation parameter and $T(u_s)$ is an upper bound on $\partial^2 E / \partial u_s^2$. The partial derivative expands as

$\frac{\partial E}{\partial u_s} = I_x\,\psi_D(I_x u_s + I_y v_s + I_t,\; \sigma_D) + \lambda \sum_{n \in N(s)} \psi_S(u_s - u_n,\; \sigma_S).$

The optimal overrelaxation parameter for an $n \times n$ image is

$\beta_{\max} = \cos\!\left(\frac{\pi}{n+1}\right), \qquad \omega_{\mathrm{opt}} = \frac{2}{1 + \sqrt{1 - \beta_{\max}^2}}.$

To avoid local minima from the non-convex $\rho$, the scale $\sigma$ is initialised large enough that the objective is convex. For the Lorentzian, convexity holds when $\rho_L'' > 0$ for all $x$, which requires $\sigma_{\mathrm{init}} \geq r_{\max} / \sqrt{2}$, where $r_{\max}$ is the maximum expected residual magnitude. $\sigma$ is then reduced linearly toward the target scale over 6 continuation stages. A 3-level Gaussian image pyramid is used; at each level, 20 SOR iterations are performed per continuation stage and the flow estimate is upsampled as initialisation for the finer level. Reference parameter ranges from §6: $\sigma_D$ from $18/\sqrt{2}$ to $5/\sqrt{2}$; $\sigma_S$ from $3/\sqrt{2}$ to $0.03/\sqrt{2}$; $\lambda_D = 5$, $\lambda_S = 1$.

Parametric (affine) variant

The 6-parameter affine flow model at pixel $\mathbf{x} = (x, y)$ is

$\mathbf{u}(\mathbf{x}; \mathbf{a}) = \begin{pmatrix} a_0 + a_1 x + a_2 y \\ a_3 + a_4 x + a_5 y \end{pmatrix}.$

The robust parametric energy over a region $R$ is

$E_p(\mathbf{a}) = \sum_{\mathbf{x} \in R} \rho\!\left((\nabla I)^T \mathbf{u}(\mathbf{x}; \mathbf{a}) + I_t,\; \sigma\right).$

Minimisation proceeds by SOR over the six parameters $a_0, \ldots, a_5$ with the GNC $\sigma$-schedule. The SOR update for each $a_i$ is

$a_i^{n+1} = a_i^n - \frac{\omega}{T_{a_i}} \cdot \frac{\partial E_p}{\partial a_i},$

where $T_{a_i} = \sum_{\mathbf{x}} \{x^2, y^2, 1, \ldots\} \cdot \max_x \rho''(x, \sigma)$ bounds the second derivative. The overrelaxation parameter is fixed at $\omega = 1.995$ and the GNC schedule reduces $\sigma$ geometrically: $\sigma_{i+1} = 0.95\,\sigma_i$.

After recovering the dominant motion parameters $\mathbf{a}^*$, pixels where $|r_\mathbf{x}| = |(\nabla I)^T \mathbf{u}(\mathbf{x}; \mathbf{a}^*) + I_t| \geq \tau$ are collected as outliers. The residual set is re-submitted to the same robust regression to recover a second motion; the process repeats until no consistent motion remains in the residuals.

Implementation

The per-pixel IRLS update for the dense variant in Rust, implementing the data-plus-smoothness gradient $\partial E / \partial u_s$ and the SOR step. The snippet covers one horizontal-component update; the vertical component $v_s$ follows the same pattern with $I_y$ in place of $I_x$.

/// Lorentzian IRLS weight: w(x, σ) = 2σ² / (2σ² + x²).
#[inline]
fn w_lorentzian(x: f32, sigma: f32) -> f32 {
    let s2 = 2.0 * sigma * sigma;
    s2 / (s2 + x * x)
}

/// One SOR update for u_s (horizontal flow at site s).
fn sor_update_u(
    r: f32,            // brightness-constancy residual I_x·u_s + I_y·v_s + I_t
    ix: f32,           // spatial gradient I_x at s
    neighbours: &[f32],// (u_n) for n ∈ N(s); length 4 on a 4-connected grid
    u_s: f32,
    sigma_d: f32,
    sigma_s: f32,
    lambda: f32,
    omega: f32,
    t_upper: f32,      // upper bound on ∂²E/∂u_s²
) -> f32 {
    // ψ_D(r, σ_D) · I_x, with ψ = w · x.
    let grad_data = ix * (w_lorentzian(r, sigma_d) * r);
    // λ · Σ_{n∈N(s)} ψ_S(u_s − u_n, σ_S).
    let grad_smooth: f32 = lambda * neighbours.iter().map(|&u_n| {
        let diff = u_s - u_n;
        w_lorentzian(diff, sigma_s) * diff
    }).sum::<f32>();
    u_s - (omega / t_upper) * (grad_data + grad_smooth)
}

The vectorised inner loop in Python:

import numpy as np

def w_lorentzian(x, sigma):
    s2 = 2.0 * sigma ** 2
    return s2 / (s2 + x ** 2)

def sor_update_u(r, ix, u, u_nbrs, sigma_d, sigma_s, lam, omega, t_upper):
    grad_data = ix * (w_lorentzian(r, sigma_d) * r)
    diffs = u[None] - u_nbrs                              # (4, H, W)
    grad_smooth = lam * (w_lorentzian(diffs, sigma_s) * diffs).sum(axis=0)
    return u - (omega / t_upper) * (grad_data + grad_smooth)

Remarks

• $O(K \cdot M \cdot N)$ work per pyramid level for the dense variant, where $K$ is the number of GNC continuation stages, $M$ is the SOR iteration count per stage, and $N$ is the number of pixels. Reference values are $K = 6$, $M = 20$ for the dense variant, and $M = 30$ per pyramid level for the parametric variant.
• Both the data term and the smoothness term must use robust penalties. Robustifying the smoothness term alone — the approach taken by classical line-process methods — leaves the data term quadratic and can increase flow error relative to the fully-quadratic baseline.
• The Lorentzian outlier threshold is $\tau_L = \sqrt{2}\,\sigma$, derived from where $\rho_L'' = 0$. The Geman-McClure threshold is $\tau_{GM} = \sigma / \sqrt{3}$. Sites where the final residual magnitude exceeds $\tau$ are classified as motion outliers; this yields an explicit outlier map as a by-product of estimation.
• Local minima in the non-convex energy are avoided via the GNC $\sigma$-schedule: $\sigma$ is initialised at $\sigma_{\mathrm{init}} \geq r_{\max} / \sqrt{2}$ (Lorentzian) or $\sigma_{\mathrm{init}} \geq r_{\max} \cdot \sqrt{3}$ (Geman-McClure) so the objective is globally convex, then lowered until the true non-convex penalty is recovered. A geometric schedule $\sigma_{i+1} = 0.95\,\sigma_i$ is used for the parametric variant; a linear ramp is preferred for the dense variant.
• The Gaussian pyramid handles displacements larger than one pixel per frame. Motions exceeding the coarse-level resolution cause hard failure; transparency that does not yield a dominant motion in a local region and global illumination changes are not handled by the brightness-constancy data term.
• The redescending M-estimator is equivalent to an analog (continuous-valued) outlier process: an explicit line-process layer on the flow field is therefore unnecessary when a redescending $\rho$ is used.

References

1. M. J. Black and P. Anandan. The Robust Estimation of Multiple Motions: Parametric and Piecewise-Smooth Flow Fields. Computer Vision and Image Understanding, 63(1)
   –104, 1996. doi
   .1006/cviu.1996.0006
2. B. K. P. Horn and B. G. Schunck. Determining Optical Flow. Artificial Intelligence, 17
   –203, 1981. hdl
   .1/6337
3. B. D. Lucas and T. Kanade. An Iterative Image Registration Technique with an Application to Stereo Vision. Proc. IJCAI, 1981. pdf