Lucas–Kanade Method

The Lucas–Kanade method estimates how a small image window moved between two frames. It assumes that every pixel in the window shares the same displacement, turns the brightness constancy constraints of all those pixels into a small least-squares problem, and refines the answer by warping and repeating. It is the standard way to compute sparse optical flow, the tracking step of the KLT feature tracker, and the prototype of direct image alignment.

Bruce Lucas and Takeo Kanade published it in 1981 as an image registration technique that uses the spatial intensity gradient to direct the search for the best match instead of testing displacements one by one. Their application was stereo vision, where it refined object depths and camera parameters. Applied to small windows in consecutive video frames, the same step computes optical flow.

Problem

Given frames $I_0$ and $I_1$ and a point $\mathbf{x}_0$, find the displacement $\mathbf{w} = (u, v)$ such that the window around $\mathbf{x}_0$ in $I_0$ matches the window around $\mathbf{x}_0 + \mathbf{w}$ in $I_1$. One brightness constancy equation per pixel cannot fix two unknowns (the aperture problem), so Lucas–Kanade assumes the flow is constant over the window.

Inputs and Outputs

Inputs: two grayscale frames; the points at which flow is wanted (or every pixel, for dense output); a window size and weighting; an initial displacement, often zero or a coarser pyramid level’s estimate.

Outputs: a displacement $(u, v)$ per point, usually with a status flag or a quality measure such as the smaller eigenvalue of the gradient matrix or the final residual.

Intuition

Each pixel in the window contributes one constraint line of possible velocities. If the window contains edges at varied orientations, the lines cross near one point, the common motion, and the method takes the velocity closest to all of them in the least-squares sense.

The estimate is exact only if the image is linear in the displacement, which holds for motions much smaller than the texture scale, so the method warps the second image by the current estimate, measures the remaining motion, and repeats, much like Newton’s method for finding a root.

Algorithm

For each point:

1. Compute the gradients $I_x, I_y$ in the window and the weighted $2 \times 2$ gradient matrix $\mathbf{M}$.
2. If the smaller eigenvalue of $\mathbf{M}$ is below a threshold, report failure (flat window or single edge orientation).
3. Initialize $\mathbf{w}$ (zero, or a prediction).
4. Repeat until the update is small or a maximum number of iterations is reached:
   1. Warp: sample $I_1$ at $\mathbf{x} + \mathbf{w}$ for every window pixel $\mathbf{x}$, by interpolation.
   2. Compute the residual $I_t = I_1(\mathbf{x} + \mathbf{w}) - I_0(\mathbf{x})$.
   3. Solve $\mathbf{M}\, \Delta\mathbf{w} = -\sum g\, I_t \nabla I_0$ for an update and set $\mathbf{w} \leftarrow \mathbf{w} + \Delta\mathbf{w}$.
5. Optionally, reject points with a large final residual.

For larger motions, the procedure runs on an image pyramid, from coarse to fine.

Mathematical Formulation

Least squares over a window

The linearized brightness constancy constraint at a pixel is $I_x u + I_y v + I_t = 0$. Lucas–Kanade minimizes the weighted sum of squared residuals over a window $\Omega$ centred on the point:

$$E(u, v) = \sum_{\mathbf{x} \in \Omega} g(\mathbf{x}) \big( I_x u + I_y v + I_t \big)^2,$$

where $g$ is a uniform or Gaussian weighting (a Gaussian favours pixels near the tracked point). (Lucas and Kanade’s one-dimensional derivation instead weighted pixels by an estimate of how linear the image is there.) Setting the derivatives with respect to $u$ and $v$ to zero gives the normal equations

$$\underbrace{\begin{bmatrix} \sum g I_x^2 & \sum g I_x I_y \\ \sum g I_x I_y & \sum g I_y^2 \end{bmatrix}}_{\mathbf{M}} \begin{bmatrix} u \\ v \end{bmatrix} = - \begin{bmatrix} \sum g I_x I_t \\ \sum g I_y I_t \end{bmatrix},$$

so $\mathbf{w} = -\mathbf{M}^{-1} \mathbf{b}$ with $\mathbf{b} = \sum g\, I_t \nabla I$. The matrix $\mathbf{M}$ is the weighted structure tensor of the window. It is invertible exactly when the window contains gradients in at least two directions, and its eigenvalues separate flat regions, edges, and corners (aperture problem); in practice the smaller eigenvalue must also be large compared with the noise.

Iterative refinement

The linear solution is a single step of an iterative minimization of the actual alignment error,

$$E(\mathbf{w}) = \sum_{\mathbf{x} \in \Omega} g(\mathbf{x}) \big[ I_1(\mathbf{x} + \mathbf{w}) - I_0(\mathbf{x}) \big]^2 .$$

Linearizing $I_1$ around the current estimate, $I_1(\mathbf{x} + \mathbf{w} + \Delta\mathbf{w}) \approx I_1(\mathbf{x} + \mathbf{w}) + \nabla I_1(\mathbf{x} + \mathbf{w})^\top \Delta\mathbf{w}$, and minimizing over $\Delta\mathbf{w}$ gives

$$\Delta\mathbf{w} = -\mathbf{M}_{\mathbf{w}}^{-1} \sum_{\mathbf{x} \in \Omega} g\, \nabla I_1(\mathbf{x} + \mathbf{w}) \big[ I_1(\mathbf{x} + \mathbf{w}) - I_0(\mathbf{x}) \big],$$

where $\mathbf{M}_{\mathbf{w}}$ is built from the gradients of the warped image. This is a Gauss–Newton step; Lucas and Kanade called it “a type of Newton–Raphson iteration”. Using the gradients of $I_0$ instead, as in the algorithm above, keeps $\mathbf{M}$ fixed across iterations and is equivalent to first order (see Variants).

General warps

The same derivation applies to any differentiable warp $\mathbf{W}(\mathbf{x}; \mathbf{p})$, such as affine motion or a homography. The gradient is multiplied by the warp’s Jacobian $\partial \mathbf{W} / \partial \mathbf{p}$, and $\mathbf{M}$ becomes the Gauss–Newton approximation to the Hessian, one row and column per parameter. Lucas and Kanade already proposed a general linear transformation of the coordinates, plus contrast and brightness parameters for the intensities.

Parameters

• Window size. Larger windows contain more gradient orientations and average out noise, but more often straddle a motion boundary. OpenCV’s default is 21×21 pixels.
• Iterations and stopping threshold. A few to a few dozen iterations per level, stopping when $\lVert \Delta\mathbf{w} \rVert$ falls below about 0.01 pixels.
• Pyramid levels. Each level halves the motion, so three or four levels extend the range to tens of pixels (coarse-to-fine estimation). Lucas and Kanade already proposed smoothing and a coarse-to-fine strategy to extend the range of convergence.
• Eigenvalue threshold. The minimum acceptable smaller eigenvalue of $\mathbf{M}$, relative to the image’s gradient scale.
• Pre-smoothing. Blurring before differentiation reduces noise and widens the convergence range, at the cost of fine detail.

Complexity

For an $n \times n$ window, each iteration costs $O(n^2)$ per point for interpolation and sums; the $2 \times 2$ solve is constant time. Tracking $K$ points for $T$ iterations over $L$ pyramid levels costs $O(K T L n^2)$, plus $O(N)$ to build the pyramid of an $N$-pixel image. Dense, non-iterated Lucas–Kanade sums the five gradient products over windows with separable filters, costing $O(N)$ regardless of window size.

Implementation

Lucas–Kanade at four points of a synthetic image pair, with Gaussian weighting and iterative refinement. The left part of the image is random texture; the right part is vertical stripes. The gradients are taken from $I_0$ rather than from the warped $I_1$, so $\mathbf{M}$ is computed once per point; for a pure translation this is the inverse compositional algorithm described under Variants.

import numpy as np
from scipy import ndimage

rng = np.random.default_rng(1)
y, x = np.mgrid[0:160, 0:160].astype(float)

# Frame 0: smooth random texture on the left, vertical stripes on the right.
texture = ndimage.gaussian_filter(rng.standard_normal(x.shape), sigma=2.5)
texture /= texture.std()
I0 = np.where(x < 100, texture, np.sin(2 * np.pi * x / 16))

# Frame 1: frame 0 translated by (u, v) = (2.3, -1.4) pixels.
true_w = np.array([2.3, -1.4])
I1 = ndimage.shift(I0, (true_w[1], true_w[0]), order=3, mode="nearest")

Iy, Ix = np.gradient(I0)  # template gradients, computed once

def lucas_kanade(px, py, half=10, sigma=5.0, iterations=10, min_eig=1e-3):
    """Displacement of the window centred at (px, py) from I0 to I1."""
    win = (slice(py - half, py + half + 1), slice(px - half, px + half + 1))
    yy, xx = y[win], x[win]
    g = np.exp(-((xx - px) ** 2 + (yy - py) ** 2) / (2 * sigma**2))  # Gaussian weights
    gx, gy, template = Ix[win], Iy[win], I0[win]
    # Weighted structure tensor; it does not change between iterations.
    M = np.array([[np.sum(g * gx * gx), np.sum(g * gx * gy)],
                  [np.sum(g * gx * gy), np.sum(g * gy * gy)]])
    lam_min = np.linalg.eigvalsh(M)[0] / g.sum()
    if lam_min < min_eig:
        return None, None, lam_min              # ill-conditioned: aperture problem
    w, history = np.zeros(2), []
    for _ in range(iterations):
        # Warp I1 back by the current estimate and solve for the residual motion.
        warped = ndimage.map_coordinates(I1, [yy + w[1], xx + w[0]], order=3, mode="nearest")
        It = warped - template
        b = -np.array([np.sum(g * gx * It), np.sum(g * gy * It)])
        w = w + np.linalg.solve(M, b)
        history.append(w)
    return history[0], history[-1], lam_min

print(f"true (u, v) = {true_w}")
for px, py in [(40, 40), (70, 110), (50, 80), (130, 80)]:
    first, final, lam_min = lucas_kanade(px, py)
    if first is None:
        print(f"point ({px:3d}, {py:3d})  min eigenvalue {lam_min:.3f}  rejected")
    else:
        print(f"point ({px:3d}, {py:3d})  min eigenvalue {lam_min:.3f}  "
              f"one step {np.round(first, 2)}  10 iterations {np.round(final, 2)}")

Output:

true (u, v) = [ 2.3 -1.4]
point ( 40,  40)  min eigenvalue 0.058  one step [ 1.65 -1.43]  10 iterations [ 2.3 -1.4]
point ( 70, 110)  min eigenvalue 0.044  one step [ 2.08 -0.73]  10 iterations [ 2.3 -1.4]
point ( 50,  80)  min eigenvalue 0.023  one step [ 1.64 -1.16]  10 iterations [ 2.3 -1.4]
point (130,  80)  min eigenvalue 0.000  rejected

The single linear step is off by up to about 0.7 pixels; the iterated estimate is exact to the printed precision. The point on the stripes has only horizontal gradients, so its eigenvalue is zero and it is rejected. OpenCV’s calcOpticalFlowPyrLK provides a pyramidal implementation (feature tracking).

Properties and Behavior

• Local and sparse. Each estimate depends only on one window, so points are independent, parallel, and cheap; on well-textured windows with small motion, iterated estimates reach sub-pixel accuracy.
• Self-diagnosing. The eigenvalues of $\mathbf{M}$ indicate in advance whether a point can be tracked, and the final residual whether it was. Shi and Tomasi made this a selection rule, keeping points whose smaller eigenvalue exceeds a threshold, the basis of the KLT tracker (feature tracking). They also compared each feature with its first appearance under an affine model, to detect features that have become unreliable, for example through occlusion.
• Convergence. Gauss–Newton converges quickly within a basin set by the image’s dominant spatial frequencies; outside it, the iteration may lock onto a wrong alignment.

Limitations

• Small motion. Without a pyramid, displacements beyond a pixel or two on fine texture fail; with one, small fast-moving structures that vanish at coarse levels are still lost.
• Constant flow in the window. Windows spanning a motion boundary or occlusion, or undergoing rotation or scaling, produce a blend of motions.
• Brightness constancy. Lighting changes and specular highlights bias the least-squares fit, which a few outlier pixels can dominate.
• Aperture problem and textureless regions. Nothing is propagated into windows on edges or in flat areas, so dense Lucas–Kanade output is unreliable over much of a typical image.

Variants

Pyramidal Lucas–Kanade runs the iteration at every pyramid level, passing the scaled estimate to the next finer level; most feature trackers, including OpenCV’s, use it. Bergen, Anandan, Hanna, and Hingorani combined coarse-to-fine estimation with Gauss–Newton refinement of parametric motion models such as affine flow.

Baker and Matthews’ unifying framework classifies Lucas–Kanade algorithms by how the warp is updated. The original is the forwards additive algorithm, which recomputes the gradient and Hessian at every iteration. The forwards compositional algorithm composes the current warp with an incremental one. The inverse compositional algorithm estimates the incremental warp on the template and composes its inverse with the current warp, so the gradient, Jacobian, and Hessian depend only on the template and are precomputed; it requires warps that form a group. The inverse additive algorithm of Hager and Belhumeur also precomputes, for a narrower class of warps. Baker and Matthews show that the four are equivalent to first order in the parameter update, and that the inverse algorithms are substantially cheaper per iteration.

General warps (affine, homography) handle rotation, scaling, and shear of a tracked patch, or align planar regions and whole images for mosaicking and stabilization. Robust variants replace the squared residual with a robust penalty, solved by iteratively reweighted least squares, to down-weight occluded or specular pixels; others add gain and bias parameters for illumination.

Local–global combinations. Bruhn, Weickert, and Schnörr combined the local least-squares data term of Lucas–Kanade with the global smoothness term of the Horn–Schunck method, producing dense flow that is robust to noise yet fills in textureless regions. The two methods are the classical starting points of optical flow estimation.