Horn–Schunck Method

The Horn–Schunck method computes a dense optical flow field, one motion vector per pixel, by minimizing a single energy over the whole image. The energy combines the brightness constancy constraint at every pixel with a penalty on how quickly the flow changes from pixel to pixel. Berthold Horn and Brian Schunck introduced it in 1981. It founded the variational approach to flow, from which many of the most accurate pre-deep-learning methods descend.

Problem

Brightness constancy fixes only the motion component along the image gradient (the aperture problem), and nothing in uniform regions. Where the Lucas–Kanade method assumes constant flow in a small window, Horn and Schunck assume that the flow varies smoothly almost everywhere, as it does for opaque objects moving rigidly or deforming, and impose this over the entire image.

Inputs and Outputs

Inputs: two consecutive frames; a smoothness weight $\alpha$; a number of iterations; optionally an initial flow field.

Outputs: the flow components $u$ and $v$ at every pixel, with no per-pixel confidence.

Intuition

Think of the flow field as an elastic sheet. At each pixel, the image gradient pulls the flow toward its constraint line, more strongly where the gradient is large, while the smoothness term ties each pixel to its neighbours. In texture the data pull dominates. Where the image is flat, or has a single edge orientation, some directions feel no pull and the sheet takes its shape from the surroundings, so motion measured at corners and in texture spreads along edges and into uniform regions.

Algorithm

1. Estimate the derivatives $I_x$, $I_y$, $I_t$ at every pixel from the two frames.
2. Initialize $u$ and $v$, usually to zero.
3. Repeat for a fixed number of iterations:
   1. Compute the local averages $\bar{u}$ and $\bar{v}$ of the current flow at every pixel.
   2. Move every pixel from the local average toward its constraint line, using the formula below.
4. Return $u$ and $v$.

Mathematical Formulation

Energy

Horn and Schunck write brightness as $E$; this article uses $I$. The method minimizes

$$\mathcal{E}(u, v) = \iint \Big[ \big( I_x u + I_y v + I_t \big)^2 + \alpha^2 \big( \lVert \nabla u \rVert^2 + \lVert \nabla v \rVert^2 \big) \Big] \, \mathrm{d}x \, \mathrm{d}y .$$

The first term, the data term, is the squared error in the linearized brightness constancy equation. The second, the smoothness term, is the squared gradient magnitude of each flow component, weighted by $\alpha^2$. The paper mentions the sum of squared Laplacians of $u$ and $v$ as an alternative, the measure stated in the 1980 MIT AI memo version.

Euler–Lagrange equations

By the calculus of variations, a minimizer satisfies

$$\begin{aligned} I_x^2\, u + I_x I_y\, v &= \alpha^2 \nabla^2 u - I_x I_t, \\ I_x I_y\, u + I_y^2\, v &= \alpha^2 \nabla^2 v - I_y I_t, \end{aligned}$$

where $\nabla^2$ is the Laplacian. At image borders, the natural boundary condition is a zero derivative of the flow normal to the border.

Discretization

The Laplacian is approximated by the difference between a local average and the value at the pixel, $\nabla^2 u \approx \kappa (\bar{u} - u)$, with

$$\bar{u}_{i,j} = \tfrac{1}{6}\big( u_{i-1,j} + u_{i,j+1} + u_{i+1,j} + u_{i,j-1} \big) + \tfrac{1}{12}\big( u_{i-1,j-1} + u_{i-1,j+1} + u_{i+1,j+1} + u_{i+1,j-1} \big)$$

and $\kappa = 3$ for these weights and a grid spacing of one pixel; $\bar{v}$ is defined the same way. Substituting into the Euler–Lagrange equations gives a $2 \times 2$ linear system at each pixel. Strictly, the weight becomes $\kappa \alpha^2$; the paper’s equations carry $\alpha^2$ alone (the coefficient determinant it gives is $\alpha^2(\alpha^2 + I_x^2 + I_y^2)$), which amounts to absorbing $\kappa$ into the weight. Rearranged, the system reads

$$\begin{aligned} (\alpha^2 + I_x^2 + I_y^2)(u - \bar{u}) &= -I_x \big( I_x \bar{u} + I_y \bar{v} + I_t \big), \\ (\alpha^2 + I_x^2 + I_y^2)(v - \bar{v}) &= -I_y \big( I_x \bar{u} + I_y \bar{v} + I_t \big). \end{aligned}$$

Geometrically, the solution lies on the line through the local average $(\bar{u}, \bar{v})$ perpendicular to the constraint line, at a distance proportional to how badly the average violates the constraint.

Iterative solution

The equations for all pixels form a very large, sparse linear system. Rather than solving it directly, Horn and Schunck, citing iterative methods such as Gauss–Seidel, compute new estimates from the averages of the previous ones:

$$\begin{aligned} u^{n+1} &= \bar{u}^{\,n} - \frac{I_x \big( I_x \bar{u}^{\,n} + I_y \bar{v}^{\,n} + I_t \big)}{\alpha^2 + I_x^2 + I_y^2}, \\ v^{n+1} &= \bar{v}^{\,n} - \frac{I_y \big( I_x \bar{u}^{\,n} + I_y \bar{v}^{\,n} + I_t \big)}{\alpha^2 + I_x^2 + I_y^2}. \end{aligned}$$

As the paper notes, the new value at a pixel does not depend directly on its own previous value; as written, with all averages taken from the previous iterate, this is a Jacobi-type iteration. Where the average needs points outside the image, velocities are copied from adjacent points further in.

Derivative estimates

Horn and Schunck required $I_x$, $I_y$, and $I_t$ to refer to the same point in space and time, and estimated all three at the centre of a $2 \times 2 \times 2$ cube of pixels spanning two frames, each as the average of four first differences along the cube’s parallel edges.

Parameters

• Smoothness weight $\alpha$. It matters mainly where the gradient is small, where $\alpha^2$ dominates the denominator and keeps the update near the local average, suppressing noise in the derivatives. Horn and Schunck suggested setting $\alpha^2$ roughly equal to the expected noise in the estimate of $I_x^2 + I_y^2$, so it must change with the intensity scale. Larger values give smoother fields and more bleeding across motion boundaries.
• Number of iterations. Filling in uniform regions from their border is analogous to heat diffusion. Horn and Schunck advised using more iterations than the number of pixels across the largest region to be filled in. Because diffusion distance grows only with the square root of the number of steps, full convergence can take far more, as the example below shows.
• Initialization. For video, the paper also proposes a single iteration per new frame, starting from the previous flow.
• Pre-smoothing. Smoothing the images before differentiation reduces noise and aliasing.

Complexity

Each iteration costs $O(N)$ for $N$ pixels (two $3 \times 3$ convolutions and pointwise arithmetic), and memory is $O(N)$. The total is $O(N T)$ for $T$ iterations, and $T$ grows with the size of the regions to be filled in, so plain iteration is slow on large images. Successive over-relaxation and multigrid solvers need far fewer sweeps.

Implementation

A textured square with a uniform hole moves by $(1, 0.5)$ pixels over a static textured background. The code follows the paper’s derivative estimates, local average, and update, and reports the mean endpoint error on the object, the background, a band around the object’s boundary, and the hole.

import numpy as np
from scipy import ndimage

N = 128
y, x = np.mgrid[0:N, 0:N].astype(float)

def texture(seed, pad=20):
    t = np.random.default_rng(seed).standard_normal((N + 2 * pad, N + 2 * pad))
    t = ndimage.gaussian_filter(t, sigma=3)
    return t / t.std()

background, object_texture = texture(1)[20:-20, 20:-20], texture(2)
true_w = np.array([1.0, 0.5])  # object motion (u, v) in pixels per frame

def frame(t):
    """A textured square with a uniform hole, moving over a static textured background."""
    dx, dy = t * true_w
    moved = ndimage.shift(object_texture, (dy, dx), order=3)[20:-20, 20:-20]
    r = np.maximum(np.abs(x - 64 - dx), np.abs(y - 64 - dy))  # Chebyshev distance to centre
    image = np.where(r < 30, moved, background)
    return np.where(r < 10, 0.0, image)

def derivatives(E0, E1):
    """E_x, E_y, E_t as averages of four first differences over a 2x2x2 cube."""
    def cube(E, k):  # k: 2x2 weights placed on pixels (i, j) .. (i+1, j+1)
        return ndimage.correlate(E, np.pad(k, ((1, 0), (1, 0))), mode="nearest")
    kx = np.array([[-1, 1], [-1, 1]]) / 4
    ky = np.array([[-1, -1], [1, 1]]) / 4
    kt = np.full((2, 2), 1 / 4)
    return cube(E0, kx) + cube(E1, kx), cube(E0, ky) + cube(E1, ky), cube(E1, kt) - cube(E0, kt)

# Local average: weight 1/6 for the four edge neighbours, 1/12 for the four diagonal ones.
AVERAGE = np.array([[1, 2, 1], [2, 0, 2], [1, 2, 1]]) / 12

def horn_schunck(E0, E1, alpha, iterations):
    Ex, Ey, Et = derivatives(E0, E1)
    u, v = np.zeros_like(E0), np.zeros_like(E0)
    for _ in range(iterations):
        u_bar = ndimage.correlate(u, AVERAGE, mode="nearest")
        v_bar = ndimage.correlate(v, AVERAGE, mode="nearest")
        step = (Ex * u_bar + Ey * v_bar + Et) / (alpha**2 + Ex**2 + Ey**2)
        u, v = u_bar - Ex * step, v_bar - Ey * step
    return u, v

I0, I1 = frame(0), frame(1)
r = np.maximum(np.abs(x - 64), np.abs(y - 64))
gt_u, gt_v = np.where(r < 30, true_w[0], 0.0), np.where(r < 30, true_w[1], 0.0)
regions = {
    "object": (r > 13) & (r < 26),
    "background": (r > 34) & (r < 56),
    "boundary": np.abs(r - 30) <= 3,
    "hole": r < 7,
}

print("alpha  iters   mean endpoint error by region               flow in hole")
for alpha, iterations in [(0.5, 10), (0.5, 100), (0.5, 1000), (0.1, 1000), (2.0, 1000)]:
    u, v = horn_schunck(I0, I1, alpha, iterations)
    epe = np.hypot(u - gt_u, v - gt_v)
    errors = "  ".join(f"{name} {epe[m].mean():.2f}" for name, m in regions.items())
    hole = regions["hole"]
    print(f"{alpha:5.1f} {iterations:6d}   {errors}   ({u[hole].mean():.2f}, {v[hole].mean():.2f})")

Output:

alpha  iters   mean endpoint error by region               flow in hole
  0.5     10   object 0.58  background 0.00  boundary 0.40  hole 1.08   (0.04, 0.02)
  0.5    100   object 0.06  background 0.01  boundary 0.39  hole 0.40   (0.65, 0.35)
  0.5   1000   object 0.03  background 0.02  boundary 0.39  hole 0.10   (0.99, 0.52)
  0.1   1000   object 0.02  background 0.01  boundary 0.36  hole 0.14   (1.01, 0.48)
  2.0   1000   object 0.13  background 0.08  boundary 0.47  hole 0.06   (0.97, 0.52)

On the textured object, the error falls from 0.58 to 0.03 pixels. The uniform hole, 20 pixels across, has no local information, yet after 1,000 iterations it carries the object’s motion, filled in from its border; after 100 the filling-in is only partly done. Near the outline the error stays near 0.4 pixels whatever the settings, because the smoothness term blends the two motions, and a larger $\alpha$ spreads the object’s motion further into the static background.

Properties and Behavior

• Dense output and filling in. Every pixel gets a flow vector. In uniform regions the iteration reduces to repeated averaging, and the result solves Laplace’s equation with the surrounding flow as boundary values; along a single edge orientation, the normal component comes from the data and the tangential component from the neighbourhood.
• Unique solution. The energy is a convex quadratic in $(u, v)$, so the discrete problem has a single minimizer unless every image gradient is parallel or zero, in which case the flow along the common edge direction is undetermined.
• Behaviour of the iteration. From zero flow, the first iteration gives vectors along the gradient direction, a shrunken normal flow; later iterations bring in the tangential component from neighbours. Horn and Schunck reported robustness to coarse quantization and additive noise in synthetic tests.

Limitations

• Blurred motion boundaries. The quadratic smoothness term penalizes large flow differences heavily, so the solution varies smoothly across occluding boundaries instead of jumping. Horn and Schunck anticipated that occluding edges would cause difficulties.
• Sensitivity to outliers. With a quadratic data term, pixels that violate brightness constancy (occlusions, highlights, illumination changes) pull the flow of their neighbourhood, and the smoothness term spreads the error.
• Small motions only. The data term is linearized around zero flow, so the original method is accurate only for displacements of about a pixel or less, depending on the image texture. Larger motions need coarse-to-fine estimation with warping.
• Parameter sensitivity. The best $\alpha$ depends on image content, noise, and intensity scaling.
• No confidence measure. Propagated values look the same as measured ones.

Variants

The Horn–Schunck energy is the template for variational optical flow, a data term plus a regularizer minimized over the whole image, and with Lucas–Kanade one of the two classical starting points of optical flow estimation. Its descendants change one of the terms or the minimization.

Robust penalties. Black and Anandan replaced the quadratic penalties in both terms with robust functions such as the Lorentzian, whose influence decreases for large residuals, so pixels violating brightness constancy are treated as outliers and the flow can change abruptly at motion boundaries. The energy is not convex; they minimized it with graduated non-convexity, a continuation method.

Image-driven and anisotropic smoothness. Other regularizers reduce smoothing across image edges, on the assumption that motion boundaries usually coincide with intensity edges. The oriented smoothness constraint of Nagel and Enkelmann is an early example.

Warping and better data terms. Brox, Bruhn, Papenberg, and Weickert kept the data term without linearization, added a gradient constancy assumption, and used a convex robust penalty close to the absolute value in both terms. They showed that the outer loop of their nested fixed-point minimization amounts to warping, which, embedded in a coarse-to-fine scheme, handles much larger displacements.

TV-L1. Zach, Pock, and Bischof used an absolute-value (L1) data term with a total variation regularizer, the sum of $\lVert \nabla u \rVert$ and $\lVert \nabla v \rVert$ rather than their squares. Total variation allows sharp discontinuities and the L1 term resists outliers. Their duality-based scheme alternates pointwise thresholding and denoising steps and ran in real time on a GPU.

Combined local–global methods. Bruhn, Weickert, and Schnörr replaced the pointwise data term with the window-integrated data term of Lucas–Kanade, combining the noise robustness of the local method with the dense output of the global one.