Lucas–Kanade step by step: iterative subpixel motion estimation

Following the slope of the error

We look for the displacement $\mathbf{u} = (u, v)$ that moves a window of pixels $\mathbf{x}_i$ in frame $I_0$ onto the same content in frame $I_1$ , by minimizing the sum of squared differences

E(\mathbf{u}) = \sum_i \left[ I_1(\mathbf{x}_i + \mathbf{u}) - I_0(\mathbf{x}_i) \right]^2 = \sum_i e_i^2 .

A full search evaluates $E$ at every shift on a grid, which is expensive and only as accurate as the grid. Lucas and Kanade instead start from an initial guess $\mathbf{u}_0$ , often zero or the motion of the previous frame, and improve it step by step using the image gradient. Shifts between pixels are no problem, since $I_1$ can be interpolated anywhere.

Linearizing the image

Close to the current estimate $\mathbf{u}$ , the shifted image is approximated by its first-order Taylor expansion, with the image gradient $\mathbf{J}_1 = \nabla I_1 = (\partial I_1 / \partial x, \partial I_1 / \partial y)$ :

E(\mathbf{u} + \Delta\mathbf{u}) = \sum_i \left[ I_1(\mathbf{x}_i + \mathbf{u} + \Delta\mathbf{u}) - I_0(\mathbf{x}_i) \right]^2 \approx \sum_i \left[ \mathbf{J}_1(\mathbf{x}_i + \mathbf{u})\, \Delta\mathbf{u} + e_i \right]^2

with the current errors $e_i = I_1(\mathbf{x}_i + \mathbf{u}) - I_0(\mathbf{x}_i)$ . This is linear in the unknown update $\Delta\mathbf{u}$ . For $\mathbf{u} = 0$ , $e_i$ is the temporal derivative $I_t$ , and setting each term to zero gives the optical flow constraint $I_x u + I_y v + I_t = 0$ : one equation per pixel for two unknowns. A single pixel can only tell the motion along its gradient, which is why Lucas–Kanade assumes the same motion for the whole window.

A 2 × 2 least-squares problem

Stacking the gradients of all $N$ window pixels as the rows of a matrix $A$ and the negated errors as a vector $\mathbf{b}$ turns the approximation into $\| A\, \Delta\mathbf{u} - \mathbf{b} \|^2$ , a least-squares problem like line fitting:

A = \begin{bmatrix} I_x(\mathbf{x}_1 + \mathbf{u}) & I_y(\mathbf{x}_1 + \mathbf{u}) \\ \vdots & \vdots \\ I_x(\mathbf{x}_N + \mathbf{u}) & I_y(\mathbf{x}_N + \mathbf{u}) \end{bmatrix}, \qquad \mathbf{b} = - \begin{bmatrix} e_1 \\ \vdots \\ e_N \end{bmatrix}

Its minimum solves the normal equations

A^\top A\, \Delta\mathbf{u} = A^\top \mathbf{b}, \qquad \begin{bmatrix} \sum I_x^2 & \sum I_x I_y \\ \sum I_x I_y & \sum I_y^2 \end{bmatrix} \begin{bmatrix} \Delta u \\ \Delta v \end{bmatrix} = - \begin{bmatrix} \sum I_x e_i \\ \sum I_y e_i \end{bmatrix} .

$A^\top A$ is the second moment matrix $M$ of the autocorrelation surface and the Harris detector, and the dashed ellipse on the error surface is the contour of the linearized error through the current estimate. Its center is where the linear model expects the minimum, at $\mathbf{u} + \Delta\mathbf{u}$ .

Iterating

Because the image is not linear, the model's minimum is only an approximation. Lucas–Kanade moves there, or a fraction $\alpha$ of the way, linearizes again and repeats,

\mathbf{u}_{k+1} = \mathbf{u}_k + \alpha\, \Delta\mathbf{u}_k ,

until the update is smaller than a tolerance, here $0.001$ px. This is the Gauss–Newton method. Near the solution, full steps ( $\alpha = 1$ ) make the distance to the truth drop steeply on the logarithmic plot. A smaller $\alpha$ is more cautious and shrinks the distance by roughly a constant factor per step, while $\alpha > 1$ overshoots and zigzags across the minimum, and too large a step never settles.

The gradient $\mathbf{J}_1(\mathbf{x}_i + \mathbf{u})$ has to be interpolated and $A^\top A$ rebuilt in every step. Since $I_1(\mathbf{x}_i + \mathbf{u}^*) = I_0(\mathbf{x}_i)$ at the solution, the gradient $\mathbf{J}_0(\mathbf{x}_i) = \nabla I_0(\mathbf{x}_i)$ is a good substitute near it, and it can be computed once, together with $A^\top A$ and its inverse.

When it fails

Aperture problem: on a straight edge all gradients point the same way, $A^\top A$ has one eigenvalue near zero, and only the motion across the edge is determined. In a flat region both eigenvalues vanish. The windows that track well are the ones a corner detector would pick.
Large motion: the linearization only holds over roughly the size of the image structures. From a start farther away, the steps can point the wrong way and end in a different local minimum. Smoothing both frames widens the basin of convergence at the cost of detail, which is the idea behind estimating the motion coarse to fine on an image pyramid.
Changed appearance: the method assumes that a point keeps its brightness (brightness constancy). Noise, lighting changes and occlusion leave errors that no shift can explain.

Try this

Press Play from u₀ with the defaults: after 7 steps the estimate is within 0.02 px of the true shift (2.6, −1.8), far below the grid of a full search. The small rest comes from interpolation: $I_1$ was made by resampling $I_0$ , which the smoothing does not undo exactly. With $\sigma = 0$ , the estimate ends at (2.600, −1.800).
Click on the error surface at about (−4, 4): the path follows the slope down to the minimum. Start in a corner, e.g. (−7, −7) or (7, 7): the steps run off to a wrong minimum outside the surface shown. Raise the pre-smoothing $\sigma$ to 3 and try the same starts again.
Set the step size $\alpha$ to 0.3: every step covers only part of the way, and the distance falls along a nearly straight line on the log plot. At 1.3 the path zigzags across the minimum and needs about 25 steps; at 1.5 it no longer settles.
Set the window radius to 4 and select a window on a straight edge, such as the right wall of the house below its window: $\lambda_2 \ll \lambda_1$ , the dashed ellipse becomes long and thin, and the estimate stops somewhere along the dark valley. Only the motion across the edge is right. Click the flat sky: $A^\top A$ is singular and there is no step at all.
Set $\sigma$ to 0 and start at (−4, 4): without smoothing, the surface has flat stretches with tiny gradients, and the estimate crawls along and runs out of iterations.
Switch the gradient to $\nabla I_0(\mathbf{x}_i)$ : from $\mathbf{u}_0 = 0$ it converges as before. From (−4, 4) it gets lost, since $\nabla I_0$ only describes $I_1$ near the solution.
Add noise $\sigma_n = 0.05$ : the error image stays grainy and $E$ cannot reach zero, but the estimate still ends within a few hundredths of a pixel, since the window averages the noise over all its pixels.