12 Optical Flow and Motion

Chapter goal. This chapter develops optical flow and motion analysis from a mathematical point of view. The central theme is that motion is not directly observed; it is inferred from temporal changes in image intensity, geometry, and learned features. We will derive the optical-flow constraint equation, explain why it is ill-posed, and study local least-squares, global variational, multiscale, robust, and modern learned approaches.

Computer vision becomes especially difficult when the scene moves. A still image contains spatial information. A video contains spatial and temporal information, and the goal is often to estimate how visual content moves from one frame to the next.

The central object of this chapter is a motion field or optical flow field

\[ w(x,y)=(u(x,y),v(x,y)), \]

where $u$ is horizontal displacement and $v$ is vertical displacement. If $I_0$ and $I_1$ are two consecutive frames, optical flow tries to describe the correspondence

\[ I_0(x,y) \approx I_1(x+u(x,y),y+v(x,y)). \]

This looks simple, but it is mathematically subtle. A single pixel gives one brightness equation but two unknown motion components. This is the aperture problem. To obtain a useful estimate, one must add assumptions: local constant motion, smoothness, feature consistency, geometric constraints, probabilistic priors, or learned representations.

For MA Applied Mathematics and MS Statistics students, optical flow is a beautiful case study in:

differential calculus and Taylor approximation;
linear least squares and conditioning;
variational regularization and partial differential equations;
robust statistics and outlier modeling;
multiscale numerical algorithms;
modern deep learning with cost volumes and iterative refinement.

12.1 Learning objectives

After studying this chapter, you should be able to:

define sparse and dense optical flow fields;
derive the optical-flow constraint equation from brightness constancy;
explain the aperture problem geometrically and statistically;
formulate Lucas–Kanade optical flow as a weighted least-squares problem;
connect Lucas–Kanade conditioning to the structure tensor eigenvalues;
formulate Horn–Schunck flow as a variational regularization problem;
explain pyramids, warping, and coarse-to-fine estimation;
describe robust penalties, occlusion handling, and motion boundaries;
connect classical flow to KLT tracking, video stabilization, action analysis, and learned models such as FlowNet and RAFT;
implement small Python experiments for optical-flow estimation and motion-based learning.

12.2 11.1 Videos as functions of space and time

A grayscale video may be modeled as a function

\[ I:\Omega\times[0,T]\to\mathbb{R}, \qquad (x,y,t)\mapsto I(x,y,t), \]

where $\Omega\subset\mathbb{R}^2$ is the image domain. A color video is

\[ I:\Omega\times[0,T]\to\mathbb{R}^3. \]

A digital video is a sampled version:

\[ I[i,j,k]\approx I(x_i,y_j,t_k), \]

where $(i,j)$ indexes pixels and $k$ indexes frames.

A motion field is a time-dependent vector field

\[ w(x,y,t)=(u(x,y,t),v(x,y,t)). \]

For two frames separated by $\Delta t$, the displacement version is often written

\[ w(x,y)=\big(u(x,y),v(x,y)\big), \]

where $u$ and $v$ are measured in pixels per frame. If one wants velocity instead of displacement, then

\[ \dot x(x,y,t)=\frac{u(x,y,t)}{\Delta t}, \qquad \dot y(x,y,t)=\frac{v(x,y,t)}{\Delta t}. \]

12.2.1 Motion field versus optical flow

It is useful to distinguish two related ideas.

Object	Meaning
3D motion field	projected motion of actual 3D scene points
optical flow	apparent image motion inferred from intensity or feature changes

The two are not always the same. Specular highlights, shadows, occlusions, transparency, illumination changes, and repeated texture can produce apparent image motion that does not correspond to physical object motion.

In this chapter, we first study optical flow as a 2D image-based inference problem. Later chapters connect motion to camera geometry, stereo, 3D reconstruction, and structure from motion.

12.3 11.2 Brightness constancy

The classical assumption is brightness constancy:

\[ I(x(t),y(t),t)=\text{constant along a moving point}. \]

If a visible point moves through the image according to

\[ \frac{dx}{dt}=u(x,y,t), \qquad \frac{dy}{dt}=v(x,y,t), \]

then brightness constancy says

\[ \frac{d}{dt} I(x(t),y(t),t)=0. \]

By the chain rule,

\[ I_x\frac{dx}{dt}+I_y\frac{dy}{dt}+I_t=0. \]

Therefore,

\[ I_x u + I_y v + I_t = 0. \]

This is the optical-flow constraint equation.

12.3.1 Vector form

Let

\[ \nabla I=(I_x,I_y)^T, \qquad w=(u,v)^T. \]

Then the constraint becomes

\[ \nabla I^T w + I_t=0. \]

Equivalently,

\[ \nabla I^T w = -I_t. \]

At each pixel, this is one linear equation in two unknowns. Hence the problem is underdetermined unless additional information is used.

12.3.2 Small-displacement derivation

For two frames $I_0$ and $I_1$, brightness constancy is often written

\[ I_1(x+u,y+v)=I_0(x,y). \]

For small $(u,v)$, a first-order Taylor expansion gives

\[ I_1(x+u,y+v) \approx I_1(x,y)+I_x(x,y)u+I_y(x,y)v. \]

Thus

\[ I_xu+I_yv+I_1(x,y)-I_0(x,y)\approx 0. \]

\[ I_t\approx I_1-I_0, \]

we recover

\[ I_xu+I_yv+I_t\approx 0. \]

This derivation is local and linear. It works best when motion is small, the image is differentiable, and brightness is approximately preserved.

A statistical warning

Brightness constancy is a model, not a law of nature. In practice, a more realistic equation is

\[ I_xu+I_yv+I_t=\varepsilon, \]

where $\varepsilon$ includes illumination changes, discretization error, noise, occlusion, and model mismatch.

12.4 11.3 Interactive demo: brightness constancy and the flow constraint

The following demo shows a moving Gaussian blob. The same physical motion generates spatial derivatives, a temporal derivative, and a linear constraint on $(u,v)$.

12.4.1 Moving Gaussian and optical-flow constraint

Move the sliders to change the true horizontal and vertical displacement. The displayed constraint line is

\[ I_x u + I_y v + I_t = 0. \]

Horizontal displacement $u$: 1.0 Vertical displacement $v$: 0.5

12.5 11.4 The aperture problem

At one pixel, the equation

\[ \nabla I^T w=-I_t \]

restricts $w$ to a line in the $(u,v)$ plane, assuming $\nabla I\neq 0$.

Only the component of motion in the gradient direction is determined. Let

\[ n=\frac{\nabla I}{\|\nabla I\|_2}. \]

Then

\[ n^T w = -\frac{I_t}{\|\nabla I\|_2}. \]

The quantity

\[ w_n=\left(-\frac{I_t}{\|\nabla I\|_2}\right)n \]

is the normal flow. Any tangential component can be added without changing the brightness equation:

\[ w = w_n + \alpha t, \]

where $t$ is perpendicular to $n$.

This is the aperture problem: viewing a moving edge through a small aperture reveals motion perpendicular to the edge, but not motion parallel to the edge.

12.5.1 Connection with local image structure

The aperture problem is closely related to the structure tensor from Chapter 7. If gradients inside a window all point in nearly the same direction, then the local least-squares problem is ill-conditioned. If gradients point in two independent directions, then motion can be estimated more reliably.

Local pattern	Gradient directions	Flow estimate
flat region	very small gradients	impossible/noisy
straight edge	one dominant direction	aperture ambiguity
corner/texture	two directions	well-conditioned

12.6 11.5 Interactive demo: aperture problem geometry

12.6.1 One equation, two unknowns

Choose the spatial gradient direction and temporal derivative. The line shows all possible flow vectors satisfying the local constraint.

Gradient angle: 45 degrees Temporal derivative $I_t$: -1.0

12.7 11.6 Lucas–Kanade optical flow

The Lucas–Kanade method resolves the aperture problem by assuming that motion is approximately constant in a small neighborhood $W$ around a pixel.

For each pixel $p=(x,y)$ in the window, the constraint is

\[ I_x(p)u+I_y(p)v=-I_t(p). \]

Stack these equations over the window. Let

\[ A= \begin{bmatrix} I_x(p_1) & I_y(p_1)\\ I_x(p_2) & I_y(p_2)\\ \vdots & \vdots\\ I_x(p_m) & I_y(p_m) \end{bmatrix}, \qquad b=- \begin{bmatrix} I_t(p_1)\\ I_t(p_2)\\ \vdots\\ I_t(p_m) \end{bmatrix}. \]

Then the local constant flow estimate solves

\[ \widehat w = \arg\min_{w\in\mathbb{R}^2}\|Aw-b\|_2^2. \]

The normal equations are

\[ A^T A\widehat w=A^T b. \]

Thus

\[ \widehat w=(A^TA)^{-1}A^Tb, \]

provided $A^TA$ is invertible.

12.7.1 Weighted Lucas–Kanade

Usually nearby pixels should matter more than far pixels. Let $\omega(p_i)\ge 0$ be weights, often Gaussian. The weighted objective is

\[ E(w)=\sum_{p\in W}\omega(p)\left(I_x(p)u+I_y(p)v+I_t(p)\right)^2. \]

Define

\[ G= \sum_{p\in W}\omega(p) \begin{bmatrix} I_x(p)^2 & I_x(p)I_y(p)\\ I_x(p)I_y(p) & I_y(p)^2 \end{bmatrix}, \]

and

\[ h=-\sum_{p\in W}\omega(p)I_t(p) \begin{bmatrix} I_x(p)\\ I_y(p) \end{bmatrix}. \]

Then

\[ G\widehat w=h. \]

The matrix $G$ is exactly the local second-moment matrix or structure tensor.

12.7.2 Conditioning and eigenvalues

Let $\lambda_1\ge \lambda_2\ge 0$ be eigenvalues of $G$.

If both eigenvalues are small, the window is nearly flat.
If $\lambda_1$ is large and $\lambda_2$ is small, the window is edge-like and has aperture ambiguity.
If both eigenvalues are reasonably large, the window contains corner-like or textured information.

This is why KLT tracking often selects points with two large structure-tensor eigenvalues.

12.7.3 Statistical interpretation

Suppose the local equations satisfy

\[ b=Aw+\varepsilon, \qquad \varepsilon\sim N(0,\sigma^2 I). \]

Then Lucas–Kanade is the maximum likelihood estimator. Its covariance is

\[ \operatorname{Cov}(\widehat w)=\sigma^2(A^TA)^{-1}. \]

Thus ill-conditioned windows produce high-variance flow estimates. This gives a clean statistical meaning to the structure-tensor eigenvalues.

12.8 11.7 Interactive demo: Lucas–Kanade conditioning

12.8.1 Local least squares and gradient diversity

The slider controls how diverse the local gradient directions are. When gradients are nearly parallel, the least-squares estimate becomes unstable. When gradients point in multiple directions, the estimate is well-conditioned.

Gradient diversity: 0.60 Noise level: 0.05

12.9 11.8 From local windows to tracking: KLT

The Kanade–Lucas–Tomasi tracker tracks selected feature points across frames. The idea is:

select points with strong local structure tensor eigenvalues;
estimate small displacement using Lucas–Kanade;
iterate because real displacement is not exactly infinitesimal;
use image pyramids for larger motion.

Let $T(x)$ be a template patch in the first frame and let $I(x;p)$ be a warped patch in the second frame, depending on transformation parameters $p$. The alignment problem is

\[ \min_p \sum_{x\in W}\left[I(W(x;p))-T(x)\right]^2. \]

For pure translation, $p=(u,v)$. For affine tracking,

\[ W(x;p)=Ax+b, \]

where $A\in\mathbb{R}^{2\times 2}$ and $b\in\mathbb{R}^2$.

A Gauss–Newton step linearizes the residual and solves a least-squares problem. This makes KLT a direct application of nonlinear least squares.

12.10 11.9 Large motion: pyramids and warping

The derivation of the optical-flow constraint assumes small motion. If an object moves 20 pixels between frames, the Taylor approximation is poor. A standard remedy is a coarse-to-fine image pyramid.

Let

\[ I^{(0)},I^{(1)},\ldots,I^{(L)} \]

be images at increasingly coarse scales, where $I^{(L)}$ is the smallest image. The algorithm proceeds:

estimate flow at the coarsest scale;
upsample the flow to the next finer scale;
warp one image using the current estimate;
estimate a residual flow correction;
repeat until the original scale is reached.

Mathematically, if $w_0$ is the current flow estimate, then one linearizes around a warp:

\[ I_1(x+w_0(x)+\delta w(x))\approx I_0(x). \]

The unknown becomes the correction $\delta w$, not the full motion.

This idea appears in classical flow, image registration, stereo matching, and modern neural correspondence models.

12.11 11.10 Horn–Schunck optical flow

Lucas–Kanade assumes locally constant motion. Horn–Schunck instead estimates a dense flow field by minimizing a global energy:

\[ E(u,v)= \int_\Omega \left(I_xu+I_yv+I_t\right)^2 +\alpha^2\left(\|\nabla u\|_2^2+\|\nabla v\|_2^2\right) \,dx\,dy. \]

The first term enforces brightness constancy. The second term enforces smoothness. The parameter $\alpha>0$ controls the trade-off.

12.11.1 Euler–Lagrange equations

The Euler–Lagrange equations are

\[ I_x(I_xu+I_yv+I_t)-\alpha^2\Delta u=0, \]

and

\[ I_y(I_xu+I_yv+I_t)-\alpha^2\Delta v=0. \]

Here $\Delta$ is the Laplacian operator. These equations show that optical flow estimation is a PDE-regularized inverse problem.

12.11.2 Discrete interpretation

On a pixel grid, the regularization term becomes

\[ \sum_{i,j}\left[(u_{i+1,j}-u_{i,j})^2+(u_{i,j+1}-u_{i,j})^2 +(v_{i+1,j}-v_{i,j})^2+(v_{i,j+1}-v_{i,j})^2\right]. \]

This penalizes neighboring flow vectors that differ strongly. It encourages smooth motion fields, but it can oversmooth motion boundaries.

12.11.3 Bayesian interpretation

Horn–Schunck can be seen as MAP estimation. The data model is

\[ I_xu+I_yv+I_t\approx 0, \]

and the prior is a Gaussian smoothness prior on the vector field:

\[ p(u,v)\propto \exp\left[-\frac{\alpha^2}{2}\int_\Omega \|\nabla u\|^2+\|\nabla v\|^2\,dx\,dy\right]. \]

Thus $\alpha$ controls prior strength.

12.12 11.11 Interactive demo: smoothness regularization

12.12.1 Data fitting versus smooth flow

This demo shows a one-dimensional toy version of flow smoothing. The observed local flow estimates are noisy. The regularized estimate solves

\[ \min_u \sum_i (u_i-z_i)^2 + \alpha^2 \sum_i (u_{i+1}-u_i)^2. \]

Smoothness $\alpha$: 1.00

12.13 11.12 Robust optical flow

The squared error is sensitive to outliers. But optical flow has many outliers: occluded pixels, shadows, reflections, repeated texture, motion boundaries, and violations of brightness constancy.

A robust energy replaces the square by a robust penalty:

\[ E(u,v)= \int_\Omega \rho_d(I_xu+I_yv+I_t) + \alpha\rho_s(\nabla u,\nabla v) \,dx\,dy. \]

Examples include:

\[ \rho(r)=\sqrt{r^2+\epsilon^2}, \]

\[ \rho(r)=|r|. \]

A common robust smoothness term is total variation:

\[ \int_\Omega \|\nabla u\|+\|\nabla v\|\,dx\,dy. \]

This encourages piecewise-smooth flow while allowing sharp motion boundaries.

12.13.1 Occlusion

A pixel in frame 0 may not be visible in frame 1. Then no corresponding brightness value exists. Classical optical flow often uses forward-backward consistency:

\[ w_{0\to 1}(x)+w_{1\to 0}(x+w_{0\to 1}(x))\approx 0. \]

Large inconsistency indicates possible occlusion or mismatch.

12.14 11.13 Motion boundaries and segmentation

Optical flow often has discontinuities at object boundaries. If a foreground object moves differently from the background, then

\[ w_{\text{foreground}}(x)\neq w_{\text{background}}(x). \]

The gradient of the flow field can reveal motion boundaries:

\[ \|\nabla w\|^2 = \|\nabla u\|^2+\|\nabla v\|^2. \]

This connects optical flow to segmentation. Instead of segmenting by color alone, one can segment by motion. In video, motion is a powerful cue for separating objects from background.

12.15 11.14 Interactive demo: motion boundaries

12.15.1 Piecewise motion and flow gradients

The demo creates a simple two-region motion field. The flow gradient highlights the boundary where motion changes.

Foreground horizontal motion: 2.0 Background horizontal motion: -0.5

12.16 11.15 Optical flow as image registration

Optical flow is also a form of image registration. Let $\Phi:\Omega\to\Omega$ be a deformation map

\[ \Phi(x)=x+w(x). \]

Then one seeks

\[ I_1(\Phi(x))\approx I_0(x). \]

A common registration objective is

\[ \min_\Phi \int_\Omega \left[I_1(\Phi(x))-I_0(x)\right]^2\,dx +\lambda R(\Phi), \]

where $R$ penalizes non-smooth or nonphysical deformations.

For small deformations,

\[ \Phi(x)=x+w(x). \]

For larger transformations, one may parameterize $\Phi$ as affine, projective, spline-based, or diffeomorphic.

This viewpoint is important in medical imaging, remote sensing, video stabilization, and 3D reconstruction.

12.17 11.16 Motion induced by camera movement

Suppose a static 3D scene is observed by a moving camera. A 3D point $X=(X,Y,Z)^T$ projects to normalized image coordinates

\[ x=\frac{X}{Z}, \qquad y=\frac{Y}{Z}. \]

If the camera translates or rotates, the image point moves. The resulting optical flow depends on scene depth. Points closer to the camera often move more rapidly in the image than far points.

This explains why motion is useful for 3D perception: image motion contains depth information. Later chapters will connect this to epipolar geometry, stereo vision, and structure from motion.

12.17.1 Focus of expansion

For forward camera motion, flow vectors often radiate outward from a point called the focus of expansion. This point corresponds to the direction of camera translation.

12.18 11.17 Interactive demo: radial flow and depth

12.18.1 Camera motion intuition

This toy demo shows radial flow. Larger depth values reduce image motion; closer points move more.

Forward motion strength: 1.0 Depth variation: 1.0

12.19 11.18 Python example: synthetic translation and flow constraints

The following example creates two synthetic frames and checks the optical-flow constraint numerically.

import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import shift, sobel

n = 96
y, x = np.mgrid[-1:1:complex(n), -1:1:complex(n)]
I0 = np.exp(-((x + 0.2)**2 + (y - 0.1)**2) / 0.08)

true_u_pixels = 3.0
true_v_pixels = -2.0
I1 = shift(I0, shift=(true_v_pixels, true_u_pixels), order=1, mode="nearest")

Ix = sobel((I0 + I1) / 2, axis=1) / 8.0
Iy = sobel((I0 + I1) / 2, axis=0) / 8.0
It = I1 - I0

residual = Ix * true_u_pixels + Iy * true_v_pixels + It

fig, axes = plt.subplots(1, 4, figsize=(12, 3))
for ax, img, title in zip(
    axes,
    [I0, I1, It, residual],
    ["frame 0", "frame 1", "temporal derivative", "flow-constraint residual"]
):
    ax.imshow(img, cmap="gray")
    ax.set_title(title)
    ax.axis("off")
plt.tight_layout()
plt.show()

print("Mean absolute residual:", np.mean(np.abs(residual)))

The residual is not exactly zero because the derivation uses linearization and the images are discrete. Smaller motion and smoother images usually reduce the residual.

12.20 11.19 Python example: Lucas–Kanade in one window

import numpy as np
from scipy.ndimage import gaussian_filter, sobel, shift

rng = np.random.default_rng(7243)

n = 80
y, x = np.mgrid[-1:1:complex(n), -1:1:complex(n)]
I0 = np.exp(-((x + 0.25)**2 + (y + 0.05)**2) / 0.06)
I0 += 0.4 * np.exp(-((x - 0.25)**2 + (y - 0.15)**2) / 0.03)
I0 = gaussian_filter(I0, 1.0)

true_flow = np.array([1.5, -1.0])
I1 = shift(I0, shift=(true_flow[1], true_flow[0]), order=1, mode="nearest")
I1 += 0.01 * rng.normal(size=I1.shape)

Im = (I0 + I1) / 2
Ix = sobel(Im, axis=1) / 8.0
Iy = sobel(Im, axis=0) / 8.0
It = I1 - I0

cy, cx = n // 2, n // 2
r = 9
sl = np.s_[cy-r:cy+r+1, cx-r:cx+r+1]

A = np.column_stack([Ix[sl].ravel(), Iy[sl].ravel()])
b = -It[sl].ravel()

# Optional Gaussian window weights
wy, wx = np.mgrid[-r:r+1, -r:r+1]
w = np.exp(-(wx**2 + wy**2) / (2 * (r/2)**2)).ravel()
Aw = A * np.sqrt(w[:, None])
bw = b * np.sqrt(w)

flow_hat, *_ = np.linalg.lstsq(Aw, bw, rcond=None)
G = Aw.T @ Aw
evals = np.linalg.eigvalsh(G)

print("True flow (u, v):", true_flow)
print("Estimated flow:", flow_hat)
print("Structure tensor eigenvalues:", evals)
print("Condition number:", evals[-1] / max(evals[0], 1e-12))

12.21 11.20 Python example: Horn–Schunck iteration

The classical Horn–Schunck update uses local averages of the current flow estimates. A simple implementation is useful pedagogically.

import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import convolve, sobel, shift, gaussian_filter

n = 96
y, x = np.mgrid[-1:1:complex(n), -1:1:complex(n)]
I0 = np.exp(-((x + 0.2)**2 + y**2) / 0.08)
I0 += 0.6 * np.exp(-((x - 0.3)**2 + (y + 0.2)**2) / 0.04)
I0 = gaussian_filter(I0, 1.0)
I1 = shift(I0, shift=(-1.0, 2.0), order=1, mode="nearest")

Ix = sobel((I0 + I1) / 2, axis=1) / 8.0
Iy = sobel((I0 + I1) / 2, axis=0) / 8.0
It = I1 - I0

alpha = 1.0
num_iter = 200
u = np.zeros_like(I0)
v = np.zeros_like(I0)

avg_kernel = np.array([[1/12, 1/6, 1/12],
                       [1/6,  0.0, 1/6],
                       [1/12, 1/6, 1/12]])

for _ in range(num_iter):
    u_bar = convolve(u, avg_kernel, mode="nearest")
    v_bar = convolve(v, avg_kernel, mode="nearest")
    residual = Ix * u_bar + Iy * v_bar + It
    denom = alpha**2 + Ix**2 + Iy**2
    u = u_bar - Ix * residual / denom
    v = v_bar - Iy * residual / denom

step = 8
plt.figure(figsize=(5, 5))
plt.imshow(I0, cmap="gray")
plt.quiver(
    np.arange(0, n, step), np.arange(0, n, step),
    u[::step, ::step], -v[::step, ::step],
    color="yellow"
)
plt.title("Horn--Schunck flow field")
plt.axis("off")
plt.show()

print("Mean estimated u:", np.mean(u))
print("Mean estimated v:", np.mean(v))

12.22 11.21 Python example: OpenCV optical flow

For practical projects, students may use OpenCV’s implementations. This example uses Farneback dense optical flow.

import cv2
import numpy as np
import matplotlib.pyplot as plt

# I0 and I1 should be grayscale images with values in [0, 255].
I0_u8 = (255 * I0 / I0.max()).astype(np.uint8)
I1_u8 = (255 * I1 / I1.max()).astype(np.uint8)

flow = cv2.calcOpticalFlowFarneback(
    I0_u8, I1_u8,
    None,
    pyr_scale=0.5,
    levels=3,
    winsize=21,
    iterations=3,
    poly_n=5,
    poly_sigma=1.2,
    flags=0
)

u = flow[..., 0]
v = flow[..., 1]
mag = np.sqrt(u**2 + v**2)

plt.figure(figsize=(5, 4))
plt.imshow(mag, cmap="magma")
plt.colorbar(label="flow magnitude")
plt.title("Dense optical-flow magnitude")
plt.axis("off")
plt.show()

The goal of this book is not to treat OpenCV as a black box. Instead, students should understand what assumptions and approximations each implementation uses.

12.23 11.22 AI component: optical flow in deep learning

Modern optical-flow methods learn correspondences from data. They typically combine several ideas:

feature extraction: encode images into learned feature maps;
matching or correlation: compare features across frames;
cost volume: store matching evidence for candidate displacements;
regularization/refinement: enforce spatial coherence;
iterative updates: repeatedly improve the flow estimate.

A learned optical-flow model can be written schematically as

\[ \widehat w = f_\theta(I_0,I_1), \]

where $f_\theta$ is trained from supervised flow labels or self-supervised photometric objectives.

12.23.1 Supervised loss

If ground-truth flow $w^*$ is available, a common loss is endpoint error:

\[ \mathcal{L}_{\text{EPE}}(\theta) = \frac{1}{|\Omega|}\sum_{x\in\Omega} \left\|f_\theta(I_0,I_1)(x)-w^*(x)\right\|_2. \]

12.23.2 Self-supervised photometric loss

Without ground-truth flow, one can warp $I_1$ back to $I_0$:

\[ \widehat I_0(x)=I_1(x+\widehat w(x)). \]

Then minimize

\[ \mathcal{L}_{\text{photo}} = \sum_x \rho\left(I_1(x+\widehat w(x))-I_0(x)\right). \]

This is the same brightness-constancy idea, now embedded inside a differentiable learning pipeline.

12.23.3 Cost volume

A cost volume compares features $F_0$ and $F_1$:

\[ C(x,d)=\langle F_0(x),F_1(x+d)\rangle, \]

where $d$ ranges over candidate displacements. The cost volume makes matching explicit and is central in many learned correspondence architectures.

12.24 11.23 TensorFlow example: differentiable photometric loss

import tensorflow as tf


def bilinear_sample(image, x, y):
    """Minimal bilinear sampler for one image, pedagogical only.

    image: tensor of shape [H, W, C]
    x, y: floating coordinate tensors of shape [H, W]
    """
    H = tf.shape(image)[0]
    W = tf.shape(image)[1]

    x0 = tf.cast(tf.floor(x), tf.int32)
    x1 = x0 + 1
    y0 = tf.cast(tf.floor(y), tf.int32)
    y1 = y0 + 1

    x0c = tf.clip_by_value(x0, 0, W - 1)
    x1c = tf.clip_by_value(x1, 0, W - 1)
    y0c = tf.clip_by_value(y0, 0, H - 1)
    y1c = tf.clip_by_value(y1, 0, H - 1)

    Ia = tf.gather_nd(image, tf.stack([y0c, x0c], axis=-1))
    Ib = tf.gather_nd(image, tf.stack([y1c, x0c], axis=-1))
    Ic = tf.gather_nd(image, tf.stack([y0c, x1c], axis=-1))
    Id = tf.gather_nd(image, tf.stack([y1c, x1c], axis=-1))

    x0f = tf.cast(x0, tf.float32)
    x1f = tf.cast(x1, tf.float32)
    y0f = tf.cast(y0, tf.float32)
    y1f = tf.cast(y1, tf.float32)

    wa = (x1f - x) * (y1f - y)
    wb = (x1f - x) * (y - y0f)
    wc = (x - x0f) * (y1f - y)
    wd = (x - x0f) * (y - y0f)

    return (
        wa[..., None] * Ia +
        wb[..., None] * Ib +
        wc[..., None] * Ic +
        wd[..., None] * Id
    )


def photometric_loss(I0, I1, flow):
    """I0, I1: [H, W, C]; flow: [H, W, 2] with components (u, v)."""
    H = tf.shape(I0)[0]
    W = tf.shape(I0)[1]
    yy, xx = tf.meshgrid(
        tf.range(H, dtype=tf.float32),
        tf.range(W, dtype=tf.float32),
        indexing="ij"
    )
    u = flow[..., 0]
    v = flow[..., 1]
    I1_warped = bilinear_sample(I1, xx + u, yy + v)
    return tf.reduce_mean(tf.abs(I1_warped - I0))

This example shows why differentiable warping is central in learned optical flow, video prediction, spatial transformers, and self-supervised depth estimation.

12.25 11.24 Evaluation metrics

Given estimated flow $\widehat w(x)$ and ground truth $w^*(x)$, the most common metric is endpoint error:

\[ \operatorname{EPE}(x)=\|\widehat w(x)-w^*(x)\|_2. \]

The average endpoint error is

\[ \operatorname{AEPE} = \frac{1}{|\Omega|}\sum_{x\in\Omega}\|\widehat w(x)-w^*(x)\|_2. \]

Another metric is angular error. Treat the flow as a 3D vector $(u,v,1)$ and compare directions:

\[ \operatorname{AE}(x) = \arccos \frac{\widehat u u^*+\widehat v v^*+1} {\sqrt{\widehat u^2+\widehat v^2+1}\sqrt{(u^*)^2+(v^*)^2+1}}. \]

Evaluation should also consider:

errors near motion boundaries;
errors in occluded regions;
small-motion versus large-motion cases;
textureless regions;
computational cost;
robustness to illumination changes.

12.26 11.25 Common failure modes

Failure mode	Mathematical cause	Possible remedy
flat region	$\nabla I\approx 0$	feature priors, regularization, larger context
straight edge	rank-deficient structure tensor	corners, texture, global smoothness
large motion	Taylor expansion invalid	pyramids, matching, learned cost volume
occlusion	no true correspondence	occlusion masks, forward-backward checks
illumination change	brightness constancy violated	gradient constancy, learned features
motion boundary	smoothness prior incorrect	robust/edge-aware regularization
repetitive pattern	ambiguous matching	larger context, semantic features
nonrigid motion	simple parametric model too weak	dense deformable models

12.27 11.26 Summary

Optical flow is a foundational problem in video analysis. The basic equation

\[ I_xu+I_yv+I_t=0 \]

is simple but underdetermined. The main mathematical issue is not writing the equation; it is making the inverse problem well-posed.

Lucas–Kanade adds a local constant-motion assumption and solves a least-squares problem. Horn–Schunck adds a global smoothness prior and solves a variational problem. Robust methods handle outliers and motion boundaries. Pyramids and warping handle larger displacement. Modern learned methods add feature learning, cost volumes, and iterative refinement, but they still preserve the classical ideas of correspondence, regularization, and motion consistency.

The conceptual message is:

\[ \text{motion estimation} = \text{brightness or feature matching} + \text{regularization} + \text{optimization or learning}. \]

12.28 AI-assisted study prompts

Use AI tools as a mathematical tutor, not as a replacement for derivation.

Ask an AI system to derive the optical-flow constraint equation. Then identify where differentiability and brightness constancy enter.
Ask for an explanation of the aperture problem. Check whether the answer distinguishes normal flow from full flow.
Ask the AI to compare Lucas–Kanade and Horn–Schunck in terms of objective functions.
Ask for a pseudocode implementation of coarse-to-fine optical flow. Then label which steps are warping, linearization, and residual correction.
Ask the AI to explain why learned optical-flow models still use cost volumes or correlations. Relate the answer to classical matching.
Ask the AI to generate a failure case for optical flow and classify it as a statistical, geometric, or optimization failure.

12.29 Exercises

12.29.1 Conceptual exercises

Explain the difference between optical flow and physical 3D motion.
Why does brightness constancy fail for shadows and specular highlights?
What is the aperture problem? Give an example using a moving straight edge.
Why are corners better than flat regions for tracking?
Explain why motion boundaries create difficulty for smoothness regularization.
Compare sparse tracking and dense optical flow.

12.29.2 Mathematical exercises

Starting from $I(x(t),y(t),t)=c$, derive $I_xu+I_yv+I_t=0$.
Suppose $\nabla I=(3,4)^T$ and $I_t=-10$. Find the normal flow.
For a Lucas–Kanade window with design matrix $A$, prove that the least-squares solution satisfies $A^TA\widehat w=A^Tb$.
Show that if all gradient vectors in a window are scalar multiples of one vector, then $A^TA$ is rank one.
Derive the Euler–Lagrange equation for the Horn–Schunck $u$ variable.
Consider the one-dimensional regularized problem \[ \min_u \sum_i (u_i-z_i)^2+\alpha^2\sum_i(u_{i+1}-u_i)^2. \] Write the corresponding linear system.
Show how endpoint error changes under a global scaling of flow vectors.

12.29.3 Python exercises

Generate two frames containing a shifted Gaussian blob and estimate the shift using Lucas–Kanade.
Compare Lucas–Kanade estimates for a corner-like image patch and an edge-like image patch.
Implement Horn–Schunck flow and test different values of $\alpha$.
Use OpenCV Farneback flow on a short video clip and visualize magnitude and angle.
Create a synthetic piecewise-motion scene and show that smooth Horn–Schunck regularization blurs the motion boundary.
Implement forward-backward consistency for occlusion detection.

12.29.4 AI/modeling exercises

Design a self-supervised optical-flow loss using photometric error, smoothness, and occlusion masking.
Explain how a cost volume $C(x,d)$ relates to nearest-neighbor matching in feature space.
Compare a CNN-based optical-flow model with a transformer-based video correspondence model.
Ask an AI system to write a RAFT-style pseudocode description. Then rewrite the answer using only mathematical notation.
Propose a final project that uses optical flow for medical imaging, traffic analysis, sports analytics, or microscopy.

12.30 Notes and references

Classical optical flow begins with variational and differential methods such as Horn–Schunck and Lucas–Kanade (Horn and Schunck 1981; Lucas and Kanade 1981). KLT tracking is closely related to local least-squares image alignment and the conditioning ideas in the structure tensor (Shi and Tomasi 1994). A useful modern survey perspective appears in Szeliski’s computer vision text (Szeliski 2022). For learned optical flow, FlowNet introduced end-to-end CNN-based flow estimation (Dosovitskiy et al. 2015), and RAFT introduced a strong recurrent all-pairs field transformation architecture (Teed and Deng 2020).

Dosovitskiy, Alexey, Philipp Fischer, Eddy Ilg, et al. 2015. “FlowNet: Learning Optical Flow with Convolutional Networks.” Proceedings of the IEEE International Conference on Computer Vision, 2758–66.

Horn, Berthold K. P., and Brian G. Schunck. 1981. “Determining Optical Flow.” Artificial Intelligence 17 (1–3): 185–203.

Lucas, Bruce D., and Takeo Kanade. 1981. “An Iterative Image Registration Technique with an Application to Stereo Vision.” Proceedings of the 7th International Joint Conference on Artificial Intelligence, 674–79.

Shi, Jianbo, and Carlo Tomasi. 1994. “Good Features to Track.” Proceedings of IEEE Conference on Computer Vision and Pattern Recognition, 593–600.

Szeliski, Richard. 2022. Computer Vision: Algorithms and Applications. 2nd ed. Springer. https://szeliski.org/Book/.

Teed, Zachary, and Jia Deng. 2020. “RAFT: Recurrent All-Pairs Field Transforms for Optical Flow.” European Conference on Computer Vision, 402–19.

# Optical Flow and Motion {#sec-optical-flow-motion} > **Chapter goal.** This chapter develops optical flow and motion analysis from a mathematical point of view. The central theme is that motion is not directly observed; it is inferred from temporal changes in image intensity, geometry, and learned features. We will derive the optical-flow constraint equation, explain why it is ill-posed, and study local least-squares, global variational, multiscale, robust, and modern learned approaches. Computer vision becomes especially difficult when the scene moves. A still image contains spatial information. A video contains spatial and temporal information, and the goal is often to estimate how visual content moves from one frame to the next. The central object of this chapter is a **motion field** or **optical flow field** $$ w(x,y)=(u(x,y),v(x,y)), $$ where $u$ is horizontal displacement and $v$ is vertical displacement. If $I_0$ and $I_1$ are two consecutive frames, optical flow tries to describe the correspondence $$ I_0(x,y) \approx I_1(x+u(x,y),y+v(x,y)). $$ This looks simple, but it is mathematically subtle. A single pixel gives one brightness equation but two unknown motion components. This is the **aperture problem**. To obtain a useful estimate, one must add assumptions: local constant motion, smoothness, feature consistency, geometric constraints, probabilistic priors, or learned representations. For MA Applied Mathematics and MS Statistics students, optical flow is a beautiful case study in: - differential calculus and Taylor approximation; - linear least squares and conditioning; - variational regularization and partial differential equations; - robust statistics and outlier modeling; - multiscale numerical algorithms; - modern deep learning with cost volumes and iterative refinement. ## Learning objectives After studying this chapter, you should be able to: 1. define sparse and dense optical flow fields; 2. derive the optical-flow constraint equation from brightness constancy; 3. explain the aperture problem geometrically and statistically; 4. formulate Lucas--Kanade optical flow as a weighted least-squares problem; 5. connect Lucas--Kanade conditioning to the structure tensor eigenvalues; 6. formulate Horn--Schunck flow as a variational regularization problem; 7. explain pyramids, warping, and coarse-to-fine estimation; 8. describe robust penalties, occlusion handling, and motion boundaries; 9. connect classical flow to KLT tracking, video stabilization, action analysis, and learned models such as FlowNet and RAFT; 10. implement small Python experiments for optical-flow estimation and motion-based learning. ## 11.1 Videos as functions of space and time A grayscale video may be modeled as a function $$ I:\Omega\times[0,T]\to\mathbb{R}, \qquad (x,y,t)\mapsto I(x,y,t), $$ where $\Omega\subset\mathbb{R}^2$ is the image domain. A color video is $$ I:\Omega\times[0,T]\to\mathbb{R}^3. $$ A digital video is a sampled version: $$ I[i,j,k]\approx I(x_i,y_j,t_k), $$ where $(i,j)$ indexes pixels and $k$ indexes frames. A motion field is a time-dependent vector field $$ w(x,y,t)=(u(x,y,t),v(x,y,t)). $$ For two frames separated by $\Delta t$, the displacement version is often written $$ w(x,y)=\big(u(x,y),v(x,y)\big), $$ where $u$ and $v$ are measured in pixels per frame. If one wants velocity instead of displacement, then $$ \dot x(x,y,t)=\frac{u(x,y,t)}{\Delta t}, \qquad \dot y(x,y,t)=\frac{v(x,y,t)}{\Delta t}. $$ ### Motion field versus optical flow It is useful to distinguish two related ideas. | Object | Meaning | |---|---| | **3D motion field** | projected motion of actual 3D scene points | | **optical flow** | apparent image motion inferred from intensity or feature changes | The two are not always the same. Specular highlights, shadows, occlusions, transparency, illumination changes, and repeated texture can produce apparent image motion that does not correspond to physical object motion. In this chapter, we first study optical flow as a 2D image-based inference problem. Later chapters connect motion to camera geometry, stereo, 3D reconstruction, and structure from motion. ## 11.2 Brightness constancy The classical assumption is **brightness constancy**: $$ I(x(t),y(t),t)=\text{constant along a moving point}. $$ If a visible point moves through the image according to $$ \frac{dx}{dt}=u(x,y,t), \qquad \frac{dy}{dt}=v(x,y,t), $$ then brightness constancy says $$ \frac{d}{dt} I(x(t),y(t),t)=0. $$ By the chain rule, $$ I_x\frac{dx}{dt}+I_y\frac{dy}{dt}+I_t=0. $$ Therefore, $$ I_x u + I_y v + I_t = 0. $$ This is the **optical-flow constraint equation**. ### Vector form Let $$ \nabla I=(I_x,I_y)^T, \qquad w=(u,v)^T. $$ Then the constraint becomes $$ \nabla I^T w + I_t=0. $$ Equivalently, $$ \nabla I^T w = -I_t. $$ At each pixel, this is one linear equation in two unknowns. Hence the problem is underdetermined unless additional information is used. ### Small-displacement derivation For two frames $I_0$ and $I_1$, brightness constancy is often written $$ I_1(x+u,y+v)=I_0(x,y). $$ For small $(u,v)$, a first-order Taylor expansion gives $$ I_1(x+u,y+v) \approx I_1(x,y)+I_x(x,y)u+I_y(x,y)v. $$ Thus $$ I_xu+I_yv+I_1(x,y)-I_0(x,y)\approx 0. $$ If $$ I_t\approx I_1-I_0, $$ we recover $$ I_xu+I_yv+I_t\approx 0. $$ This derivation is local and linear. It works best when motion is small, the image is differentiable, and brightness is approximately preserved. ::: {.callout-note} ## A statistical warning Brightness constancy is a model, not a law of nature. In practice, a more realistic equation is $$ I_xu+I_yv+I_t=\varepsilon, $$ where $\varepsilon$ includes illumination changes, discretization error, noise, occlusion, and model mismatch. ::: ## 11.3 Interactive demo: brightness constancy and the flow constraint The following demo shows a moving Gaussian blob. The same physical motion generates spatial derivatives, a temporal derivative, and a linear constraint on $(u,v)$. ::: {.interactive-card} ### Moving Gaussian and optical-flow constraint Move the sliders to change the true horizontal and vertical displacement. The displayed constraint line is $$ I_x u + I_y v + I_t = 0. $$ <div class="controls"> <label>Horizontal displacement $u$: <span id="ch11-flow-u-label">1.0</span></label> <input id="ch11-flow-u" type="range" min="-3" max="3" step="0.25" value="1.0" /> <label>Vertical displacement $v$: <span id="ch11-flow-v-label">0.5</span></label> <input id="ch11-flow-v" type="range" min="-3" max="3" step="0.25" value="0.5" /> </div> <div id="ch11-brightness-plot" class="plotly-demo"></div> ::: <script> (function(){ function gauss(x,y,cx,cy,s){ return Math.exp(-((x-cx)*(x-cx)+(y-cy)*(y-cy))/(2*s*s)); } function grid(n){ const xs=[], ys=[]; for(let i=0;i<n;i++) xs.push(-5+10*i/(n-1)); for(let j=0;j<n;j++) ys.push(-5+10*j/(n-1)); return {xs,ys}; } function draw(){ const u=parseFloat(document.getElementById('ch11-flow-u').value); const v=parseFloat(document.getElementById('ch11-flow-v').value); document.getElementById('ch11-flow-u-label').innerText=u.toFixed(2); document.getElementById('ch11-flow-v-label').innerText=v.toFixed(2); const {xs,ys}=grid(45); const z0=[], z1=[], diff=[]; for(let j=0;j<ys.length;j++){ const row0=[], row1=[], rowd=[]; for(let i=0;i<xs.length;i++){ const a=gauss(xs[i],ys[j],-1.0,-0.5,1.4); const b=gauss(xs[i],ys[j],-1.0+u,-0.5+v,1.4); row0.push(a); row1.push(b); rowd.push(b-a); } z0.push(row0); z1.push(row1); diff.push(rowd); } const x0=0.0, y0=0.0, cx=-1.0, cy=-0.5, s=1.4; const I0=gauss(x0,y0,cx,cy,s); const Ix=-(x0-cx)/(s*s)*I0; const Iy=-(y0-cy)/(s*s)*I0; const It=gauss(x0,y0,cx+u,cy+v,s)-I0; const vv=[]; const uu=[]; for(let a=-4;a<=4.01;a+=0.1){ if(Math.abs(Ix)>1e-6){ vv.push(a); uu.push((-It-Iy*a)/Ix); } } const traces=[ {z:z0, x:xs, y:ys, type:'contour', contours:{coloring:'heatmap'}, showscale:false, name:'frame 0'}, {z:z1, x:xs, y:ys, type:'contour', contours:{coloring:'lines'}, line:{width:2}, showscale:false, name:'frame 1'}, {x:[x0], y:[y0], mode:'markers', marker:{size:10}, name:'sample pixel'}, {x:[0,u], y:[0,v], mode:'lines+markers', name:'true flow vector'}, {x:uu, y:vv, mode:'lines', name:'constraint line in (u,v) space', xaxis:'x2', yaxis:'y2'} ]; const layout={ margin:{t:30,l:40,r:20,b:40}, grid:{rows:1,columns:2,pattern:'independent'}, xaxis:{title:'x'}, yaxis:{title:'y', scaleanchor:'x'}, xaxis2:{title:'candidate u'}, yaxis2:{title:'candidate v'}, title:'Frame contours and one local optical-flow constraint' }; Plotly.react('ch11-brightness-plot',traces,layout,{responsive:true}); } ['ch11-flow-u','ch11-flow-v'].forEach(id=>document.getElementById(id).addEventListener('input',draw)); draw(); })(); </script> ## 11.4 The aperture problem At one pixel, the equation $$ \nabla I^T w=-I_t $$ restricts $w$ to a line in the $(u,v)$ plane, assuming $\nabla I\neq 0$. Only the component of motion in the gradient direction is determined. Let $$ n=\frac{\nabla I}{\|\nabla I\|_2}. $$ Then $$ n^T w = -\frac{I_t}{\|\nabla I\|_2}. $$ The quantity $$ w_n=\left(-\frac{I_t}{\|\nabla I\|_2}\right)n $$ is the **normal flow**. Any tangential component can be added without changing the brightness equation: $$ w = w_n + \alpha t, $$ where $t$ is perpendicular to $n$. This is the aperture problem: viewing a moving edge through a small aperture reveals motion perpendicular to the edge, but not motion parallel to the edge. ### Connection with local image structure The aperture problem is closely related to the structure tensor from Chapter 7. If gradients inside a window all point in nearly the same direction, then the local least-squares problem is ill-conditioned. If gradients point in two independent directions, then motion can be estimated more reliably. | Local pattern | Gradient directions | Flow estimate | |---|---|---| | flat region | very small gradients | impossible/noisy | | straight edge | one dominant direction | aperture ambiguity | | corner/texture | two directions | well-conditioned | ## 11.5 Interactive demo: aperture problem geometry ::: {.interactive-card} ### One equation, two unknowns Choose the spatial gradient direction and temporal derivative. The line shows all possible flow vectors satisfying the local constraint. <div class="controls"> <label>Gradient angle: <span id="ch11-ap-angle-label">45</span> degrees</label> <input id="ch11-ap-angle" type="range" min="0" max="180" step="5" value="45" /> <label>Temporal derivative $I_t$: <span id="ch11-ap-it-label">-1.0</span></label> <input id="ch11-ap-it" type="range" min="-3" max="3" step="0.25" value="-1.0" /> </div> <div id="ch11-aperture-plot" class="plotly-demo"></div> ::: <script> (function(){ function draw(){ const angle=parseFloat(document.getElementById('ch11-ap-angle').value); const It=parseFloat(document.getElementById('ch11-ap-it').value); document.getElementById('ch11-ap-angle-label').innerText=angle.toFixed(0); document.getElementById('ch11-ap-it-label').innerText=It.toFixed(2); const th=angle*Math.PI/180; const Ix=Math.cos(th), Iy=Math.sin(th); const xs=[], ys=[]; for(let a=-4;a<=4.001;a+=0.05){ if(Math.abs(Iy)>0.05){ xs.push(a); ys.push((-It-Ix*a)/Iy); } else { xs.push(-It/Ix); ys.push(a); } } const normalScale=-It; const traces=[ {x:xs,y:ys,mode:'lines',name:'constraint line'}, {x:[0,Ix],y:[0,Iy],mode:'lines+markers',name:'gradient direction'}, {x:[0,normalScale*Ix],y:[0,normalScale*Iy],mode:'lines+markers',name:'normal flow'} ]; const layout={ margin:{t:30,l:40,r:20,b:40}, title:'Aperture problem: all flows on a line fit one pixel', xaxis:{title:'u', range:[-4,4], zeroline:true}, yaxis:{title:'v', range:[-4,4], scaleanchor:'x', zeroline:true} }; Plotly.react('ch11-aperture-plot',traces,layout,{responsive:true}); } ['ch11-ap-angle','ch11-ap-it'].forEach(id=>document.getElementById(id).addEventListener('input',draw)); draw(); })(); </script> ## 11.6 Lucas--Kanade optical flow The Lucas--Kanade method resolves the aperture problem by assuming that motion is approximately constant in a small neighborhood $W$ around a pixel. For each pixel $p=(x,y)$ in the window, the constraint is $$ I_x(p)u+I_y(p)v=-I_t(p). $$ Stack these equations over the window. Let $$ A= \begin{bmatrix} I_x(p_1) & I_y(p_1)\\ I_x(p_2) & I_y(p_2)\\ \vdots & \vdots\\ I_x(p_m) & I_y(p_m) \end{bmatrix}, \qquad b=- \begin{bmatrix} I_t(p_1)\\ I_t(p_2)\\ \vdots\\ I_t(p_m) \end{bmatrix}. $$ Then the local constant flow estimate solves $$ \widehat w = \arg\min_{w\in\mathbb{R}^2}\|Aw-b\|_2^2. $$ The normal equations are $$ A^T A\widehat w=A^T b. $$ Thus $$ \widehat w=(A^TA)^{-1}A^Tb, $$ provided $A^TA$ is invertible. ### Weighted Lucas--Kanade Usually nearby pixels should matter more than far pixels. Let $\omega(p_i)\ge 0$ be weights, often Gaussian. The weighted objective is $$ E(w)=\sum_{p\in W}\omega(p)\left(I_x(p)u+I_y(p)v+I_t(p)\right)^2. $$ Define $$ G= \sum_{p\in W}\omega(p) \begin{bmatrix} I_x(p)^2 & I_x(p)I_y(p)\\ I_x(p)I_y(p) & I_y(p)^2 \end{bmatrix}, $$ and $$ h=-\sum_{p\in W}\omega(p)I_t(p) \begin{bmatrix} I_x(p)\\ I_y(p) \end{bmatrix}. $$ Then $$ G\widehat w=h. $$ The matrix $G$ is exactly the local second-moment matrix or structure tensor. ### Conditioning and eigenvalues Let $\lambda_1\ge \lambda_2\ge 0$ be eigenvalues of $G$. - If both eigenvalues are small, the window is nearly flat. - If $\lambda_1$ is large and $\lambda_2$ is small, the window is edge-like and has aperture ambiguity. - If both eigenvalues are reasonably large, the window contains corner-like or textured information. This is why KLT tracking often selects points with two large structure-tensor eigenvalues. ### Statistical interpretation Suppose the local equations satisfy $$ b=Aw+\varepsilon, \qquad \varepsilon\sim N(0,\sigma^2 I). $$ Then Lucas--Kanade is the maximum likelihood estimator. Its covariance is $$ \operatorname{Cov}(\widehat w)=\sigma^2(A^TA)^{-1}. $$ Thus ill-conditioned windows produce high-variance flow estimates. This gives a clean statistical meaning to the structure-tensor eigenvalues. ## 11.7 Interactive demo: Lucas--Kanade conditioning ::: {.interactive-card} ### Local least squares and gradient diversity The slider controls how diverse the local gradient directions are. When gradients are nearly parallel, the least-squares estimate becomes unstable. When gradients point in multiple directions, the estimate is well-conditioned. <div class="controls"> <label>Gradient diversity: <span id="ch11-lk-div-label">0.60</span></label> <input id="ch11-lk-div" type="range" min="0.02" max="1.50" step="0.02" value="0.60" /> <label>Noise level: <span id="ch11-lk-noise-label">0.05</span></label> <input id="ch11-lk-noise" type="range" min="0" max="0.5" step="0.01" value="0.05" /> </div> <div id="ch11-lk-plot" class="plotly-demo"></div> ::: <script> (function(){ function randn(seed){ const x=Math.sin(seed*999.13)*10000; return x-Math.floor(x); } function draw(){ const div=parseFloat(document.getElementById('ch11-lk-div').value); const noise=parseFloat(document.getElementById('ch11-lk-noise').value); document.getElementById('ch11-lk-div-label').innerText=div.toFixed(2); document.getElementById('ch11-lk-noise-label').innerText=noise.toFixed(2); const trueW=[1.2,-0.6]; const gx=[], gy=[], b=[]; for(let i=0;i<40;i++){ const base=0.15; const ang=base + div*(randn(i+1)-0.5); const mag=0.4+0.8*randn(i+13); const ix=mag*Math.cos(ang), iy=mag*Math.sin(ang); const eps=noise*(randn(i+31)-0.5)*2; gx.push(ix); gy.push(iy); b.push(ix*trueW[0]+iy*trueW[1]+eps); } let g11=0,g12=0,g22=0,h1=0,h2=0; for(let i=0;i<gx.length;i++){ g11+=gx[i]*gx[i]; g12+=gx[i]*gy[i]; g22+=gy[i]*gy[i]; h1+=gx[i]*b[i]; h2+=gy[i]*b[i]; } const det=g11*g22-g12*g12; let wh=[NaN,NaN]; if(Math.abs(det)>1e-8){ wh=[(g22*h1-g12*h2)/det,(-g12*h1+g11*h2)/det]; } const tr=g11+g22; const disc=Math.sqrt(Math.max(0,(g11-g22)*(g11-g22)+4*g12*g12)); const lam1=(tr+disc)/2, lam2=(tr-disc)/2; const theta=[]; const ell=[]; for(let a=0;a<2*Math.PI;a+=0.02){ const x=Math.cos(a), y=Math.sin(a); const val=g11*x*x+2*g12*x*y+g22*y*y; theta.push(a); ell.push(val); } const traces=[ {x:gx,y:gy,mode:'markers',name:'local gradients'}, {x:[0,trueW[0]],y:[0,trueW[1]],mode:'lines+markers',name:'true flow',xaxis:'x2',yaxis:'y2'}, {x:[0,wh[0]],y:[0,wh[1]],mode:'lines+markers',name:'estimated flow',xaxis:'x2',yaxis:'y2'}, {x:theta,y:ell,type:'scatterpolar',mode:'lines',name:'quadratic form',subplot:'polar'} ]; const layout={ margin:{t:45,l:40,r:30,b:40}, title:`Lucas--Kanade conditioning: λ1=${lam1.toFixed(2)}, λ2=${lam2.toFixed(2)}`, grid:{rows:1,columns:3,pattern:'independent'}, xaxis:{title:'Ix'}, yaxis:{title:'Iy', scaleanchor:'x'}, xaxis2:{title:'u', range:[-3,3]}, yaxis2:{title:'v', range:[-3,3], scaleanchor:'x2'}, polar:{domain:{x:[0.70,1.0],y:[0,1]}, radialaxis:{visible:true}} }; Plotly.react('ch11-lk-plot',traces,layout,{responsive:true}); } ['ch11-lk-div','ch11-lk-noise'].forEach(id=>document.getElementById(id).addEventListener('input',draw)); draw(); })(); </script> ## 11.8 From local windows to tracking: KLT The Kanade--Lucas--Tomasi tracker tracks selected feature points across frames. The idea is: 1. select points with strong local structure tensor eigenvalues; 2. estimate small displacement using Lucas--Kanade; 3. iterate because real displacement is not exactly infinitesimal; 4. use image pyramids for larger motion. Let $T(x)$ be a template patch in the first frame and let $I(x;p)$ be a warped patch in the second frame, depending on transformation parameters $p$. The alignment problem is $$ \min_p \sum_{x\in W}\left[I(W(x;p))-T(x)\right]^2. $$ For pure translation, $p=(u,v)$. For affine tracking, $$ W(x;p)=Ax+b, $$ where $A\in\mathbb{R}^{2\times 2}$ and $b\in\mathbb{R}^2$. A Gauss--Newton step linearizes the residual and solves a least-squares problem. This makes KLT a direct application of nonlinear least squares. ## 11.9 Large motion: pyramids and warping The derivation of the optical-flow constraint assumes small motion. If an object moves 20 pixels between frames, the Taylor approximation is poor. A standard remedy is a **coarse-to-fine image pyramid**. Let $$ I^{(0)},I^{(1)},\ldots,I^{(L)} $$ be images at increasingly coarse scales, where $I^{(L)}$ is the smallest image. The algorithm proceeds: 1. estimate flow at the coarsest scale; 2. upsample the flow to the next finer scale; 3. warp one image using the current estimate; 4. estimate a residual flow correction; 5. repeat until the original scale is reached. Mathematically, if $w_0$ is the current flow estimate, then one linearizes around a warp: $$ I_1(x+w_0(x)+\delta w(x))\approx I_0(x). $$ The unknown becomes the correction $\delta w$, not the full motion. This idea appears in classical flow, image registration, stereo matching, and modern neural correspondence models. ## 11.10 Horn--Schunck optical flow Lucas--Kanade assumes locally constant motion. Horn--Schunck instead estimates a dense flow field by minimizing a global energy: $$ E(u,v)= \int_\Omega \left(I_xu+I_yv+I_t\right)^2 +\alpha^2\left(\|\nabla u\|_2^2+\|\nabla v\|_2^2\right) \,dx\,dy. $$ The first term enforces brightness constancy. The second term enforces smoothness. The parameter $\alpha>0$ controls the trade-off. ### Euler--Lagrange equations The Euler--Lagrange equations are $$ I_x(I_xu+I_yv+I_t)-\alpha^2\Delta u=0, $$ and $$ I_y(I_xu+I_yv+I_t)-\alpha^2\Delta v=0. $$ Here $\Delta$ is the Laplacian operator. These equations show that optical flow estimation is a PDE-regularized inverse problem. ### Discrete interpretation On a pixel grid, the regularization term becomes $$ \sum_{i,j}\left[(u_{i+1,j}-u_{i,j})^2+(u_{i,j+1}-u_{i,j})^2 +(v_{i+1,j}-v_{i,j})^2+(v_{i,j+1}-v_{i,j})^2\right]. $$ This penalizes neighboring flow vectors that differ strongly. It encourages smooth motion fields, but it can oversmooth motion boundaries. ### Bayesian interpretation Horn--Schunck can be seen as MAP estimation. The data model is $$ I_xu+I_yv+I_t\approx 0, $$ and the prior is a Gaussian smoothness prior on the vector field: $$ p(u,v)\propto \exp\left[-\frac{\alpha^2}{2}\int_\Omega \|\nabla u\|^2+\|\nabla v\|^2\,dx\,dy\right]. $$ Thus $\alpha$ controls prior strength. ## 11.11 Interactive demo: smoothness regularization ::: {.interactive-card} ### Data fitting versus smooth flow This demo shows a one-dimensional toy version of flow smoothing. The observed local flow estimates are noisy. The regularized estimate solves $$ \min_u \sum_i (u_i-z_i)^2 + \alpha^2 \sum_i (u_{i+1}-u_i)^2. $$ <div class="controls"> <label>Smoothness $\alpha$: <span id="ch11-hs-alpha-label">1.00</span></label> <input id="ch11-hs-alpha" type="range" min="0" max="5" step="0.1" value="1.0" /> </div> <div id="ch11-hs-plot" class="plotly-demo"></div> ::: <script> (function(){ const n=80; const x=[]; const truth=[]; const z=[]; for(let i=0;i<n;i++){ x.push(i); const t=i<38 ? 1.0 : -0.8; truth.push(t); z.push(t+0.45*Math.sin(0.73*i)+0.22*Math.sin(2.17*i)); } function smooth(alpha){ let u=z.slice(); const tau=0.12; for(let iter=0;iter<400;iter++){ const old=u.slice(); for(let i=1;i<n-1;i++){ const lap=old[i-1]-2*old[i]+old[i+1]; u[i]=old[i]-tau*((old[i]-z[i])-alpha*alpha*lap); } } return u; } function draw(){ const alpha=parseFloat(document.getElementById('ch11-hs-alpha').value); document.getElementById('ch11-hs-alpha-label').innerText=alpha.toFixed(2); const u=smooth(alpha); const traces=[ {x:x,y:truth,mode:'lines',name:'true piecewise flow'}, {x:x,y:z,mode:'markers',name:'noisy local estimates'}, {x:x,y:u,mode:'lines',name:'regularized estimate'} ]; const layout={ margin:{t:35,l:45,r:20,b:40}, title:'Horn--Schunck idea in 1D: smoothness can help but may blur boundaries', xaxis:{title:'pixel index'}, yaxis:{title:'flow component'} }; Plotly.react('ch11-hs-plot',traces,layout,{responsive:true}); } document.getElementById('ch11-hs-alpha').addEventListener('input',draw); draw(); })(); </script> ## 11.12 Robust optical flow The squared error is sensitive to outliers. But optical flow has many outliers: occluded pixels, shadows, reflections, repeated texture, motion boundaries, and violations of brightness constancy. A robust energy replaces the square by a robust penalty: $$ E(u,v)= \int_\Omega \rho_d(I_xu+I_yv+I_t) + \alpha\rho_s(\nabla u,\nabla v) \,dx\,dy. $$ Examples include: $$ \rho(r)=\sqrt{r^2+\epsilon^2}, $$ or $$ \rho(r)=|r|. $$ A common robust smoothness term is total variation: $$ \int_\Omega \|\nabla u\|+\|\nabla v\|\,dx\,dy. $$ This encourages piecewise-smooth flow while allowing sharp motion boundaries. ### Occlusion A pixel in frame 0 may not be visible in frame 1. Then no corresponding brightness value exists. Classical optical flow often uses forward-backward consistency: $$ w_{0\to 1}(x)+w_{1\to 0}(x+w_{0\to 1}(x))\approx 0. $$ Large inconsistency indicates possible occlusion or mismatch. ## 11.13 Motion boundaries and segmentation Optical flow often has discontinuities at object boundaries. If a foreground object moves differently from the background, then $$ w_{\text{foreground}}(x)\neq w_{\text{background}}(x). $$ The gradient of the flow field can reveal motion boundaries: $$ \|\nabla w\|^2 = \|\nabla u\|^2+\|\nabla v\|^2. $$ This connects optical flow to segmentation. Instead of segmenting by color alone, one can segment by motion. In video, motion is a powerful cue for separating objects from background. ## 11.14 Interactive demo: motion boundaries ::: {.interactive-card} ### Piecewise motion and flow gradients The demo creates a simple two-region motion field. The flow gradient highlights the boundary where motion changes. <div class="controls"> <label>Foreground horizontal motion: <span id="ch11-mb-uf-label">2.0</span></label> <input id="ch11-mb-uf" type="range" min="-3" max="3" step="0.25" value="2.0" /> <label>Background horizontal motion: <span id="ch11-mb-ub-label">-0.5</span></label> <input id="ch11-mb-ub" type="range" min="-3" max="3" step="0.25" value="-0.5" /> </div> <div id="ch11-motion-boundary-plot" class="plotly-demo"></div> ::: <script> (function(){ function draw(){ const uf=parseFloat(document.getElementById('ch11-mb-uf').value); const ub=parseFloat(document.getElementById('ch11-mb-ub').value); document.getElementById('ch11-mb-uf-label').innerText=uf.toFixed(2); document.getElementById('ch11-mb-ub-label').innerText=ub.toFixed(2); const n=26; const xs=[], ys=[], X=[], Y=[], U=[], V=[]; for(let i=0;i<n;i++){ xs.push(-3+6*i/(n-1)); ys.push(-3+6*i/(n-1)); } const mag=[]; for(let j=0;j<n;j++){ const row=[]; for(let i=0;i<n;i++){ const x=xs[i], y=ys[j]; const fg=(x*x+y*y)<2.0; row.push(fg?Math.abs(uf-ub):0); if(i%3===0 && j%3===0){ X.push(x); Y.push(y); U.push(fg?uf:ub); V.push(fg?0.6:-0.1); } } mag.push(row); } const quiverX=[], quiverY=[]; for(let k=0;k<X.length;k++){ quiverX.push(X[k],X[k]+0.25*U[k],null); quiverY.push(Y[k],Y[k]+0.25*V[k],null); } const traces=[ {z:mag,x:xs,y:ys,type:'heatmap',showscale:true,name:'flow-gradient magnitude'}, {x:quiverX,y:quiverY,mode:'lines',name:'flow vectors'} ]; const layout={ margin:{t:35,l:45,r:20,b:40}, title:'Motion boundary from discontinuity in flow', xaxis:{title:'x'}, yaxis:{title:'y', scaleanchor:'x'} }; Plotly.react('ch11-motion-boundary-plot',traces,layout,{responsive:true}); } ['ch11-mb-uf','ch11-mb-ub'].forEach(id=>document.getElementById(id).addEventListener('input',draw)); draw(); })(); </script> ## 11.15 Optical flow as image registration Optical flow is also a form of image registration. Let $\Phi:\Omega\to\Omega$ be a deformation map $$ \Phi(x)=x+w(x). $$ Then one seeks $$ I_1(\Phi(x))\approx I_0(x). $$ A common registration objective is $$ \min_\Phi \int_\Omega \left[I_1(\Phi(x))-I_0(x)\right]^2\,dx +\lambda R(\Phi), $$ where $R$ penalizes non-smooth or nonphysical deformations. For small deformations, $$ \Phi(x)=x+w(x). $$ For larger transformations, one may parameterize $\Phi$ as affine, projective, spline-based, or diffeomorphic. This viewpoint is important in medical imaging, remote sensing, video stabilization, and 3D reconstruction. ## 11.16 Motion induced by camera movement Suppose a static 3D scene is observed by a moving camera. A 3D point $X=(X,Y,Z)^T$ projects to normalized image coordinates $$ x=\frac{X}{Z}, \qquad y=\frac{Y}{Z}. $$ If the camera translates or rotates, the image point moves. The resulting optical flow depends on scene depth. Points closer to the camera often move more rapidly in the image than far points. This explains why motion is useful for 3D perception: image motion contains depth information. Later chapters will connect this to epipolar geometry, stereo vision, and structure from motion. ### Focus of expansion For forward camera motion, flow vectors often radiate outward from a point called the **focus of expansion**. This point corresponds to the direction of camera translation. ## 11.17 Interactive demo: radial flow and depth ::: {.interactive-card} ### Camera motion intuition This toy demo shows radial flow. Larger depth values reduce image motion; closer points move more. <div class="controls"> <label>Forward motion strength: <span id="ch11-radial-s-label">1.0</span></label> <input id="ch11-radial-s" type="range" min="0" max="3" step="0.1" value="1.0" /> <label>Depth variation: <span id="ch11-radial-d-label">1.0</span></label> <input id="ch11-radial-d" type="range" min="0" max="2" step="0.1" value="1.0" /> </div> <div id="ch11-radial-flow-plot" class="plotly-demo"></div> ::: <script> (function(){ function draw(){ const s=parseFloat(document.getElementById('ch11-radial-s').value); const d=parseFloat(document.getElementById('ch11-radial-d').value); document.getElementById('ch11-radial-s-label').innerText=s.toFixed(2); document.getElementById('ch11-radial-d-label').innerText=d.toFixed(2); const x=[], y=[], qx=[], qy=[]; for(let i=-4;i<=4;i++){ for(let j=-4;j<=4;j++){ if(i===0 && j===0) continue; const xx=i/2, yy=j/2; const depth=2.0+d*(0.5+0.5*Math.sin(1.2*xx+0.7*yy)); const u=s*xx/depth, v=s*yy/depth; qx.push(xx,xx+0.4*u,null); qy.push(yy,yy+0.4*v,null); x.push(xx); y.push(yy); } } const traces=[ {x:x,y:y,mode:'markers',marker:{size:5},name:'image points'}, {x:qx,y:qy,mode:'lines',name:'flow'} ]; const layout={ margin:{t:35,l:45,r:20,b:40}, title:'Toy radial flow: nearby points move more for the same camera motion', xaxis:{title:'x', range:[-2.5,2.5]}, yaxis:{title:'y', range:[-2.5,2.5], scaleanchor:'x'} }; Plotly.react('ch11-radial-flow-plot',traces,layout,{responsive:true}); } ['ch11-radial-s','ch11-radial-d'].forEach(id=>document.getElementById(id).addEventListener('input',draw)); draw(); })(); </script> ## 11.18 Python example: synthetic translation and flow constraints The following example creates two synthetic frames and checks the optical-flow constraint numerically. ```{python} #| label: py-ch11-synthetic-translation #| eval: false import numpy as np import matplotlib.pyplot as plt from scipy.ndimage import shift, sobel n = 96 y, x = np.mgrid[-1:1:complex(n), -1:1:complex(n)] I0 = np.exp(-((x + 0.2)**2 + (y - 0.1)**2) / 0.08) true_u_pixels = 3.0 true_v_pixels = -2.0 I1 = shift(I0, shift=(true_v_pixels, true_u_pixels), order=1, mode="nearest") Ix = sobel((I0 + I1) / 2, axis=1) / 8.0 Iy = sobel((I0 + I1) / 2, axis=0) / 8.0 It = I1 - I0 residual = Ix * true_u_pixels + Iy * true_v_pixels + It fig, axes = plt.subplots(1, 4, figsize=(12, 3)) for ax, img, title in zip( axes, [I0, I1, It, residual], ["frame 0", "frame 1", "temporal derivative", "flow-constraint residual"] ): ax.imshow(img, cmap="gray") ax.set_title(title) ax.axis("off") plt.tight_layout() plt.show() print("Mean absolute residual:", np.mean(np.abs(residual))) ``` The residual is not exactly zero because the derivation uses linearization and the images are discrete. Smaller motion and smoother images usually reduce the residual. ## 11.19 Python example: Lucas--Kanade in one window ```{python} #| label: py-ch11-lucas-kanade-window #| eval: false import numpy as np from scipy.ndimage import gaussian_filter, sobel, shift rng = np.random.default_rng(7243) n = 80 y, x = np.mgrid[-1:1:complex(n), -1:1:complex(n)] I0 = np.exp(-((x + 0.25)**2 + (y + 0.05)**2) / 0.06) I0 += 0.4 * np.exp(-((x - 0.25)**2 + (y - 0.15)**2) / 0.03) I0 = gaussian_filter(I0, 1.0) true_flow = np.array([1.5, -1.0]) I1 = shift(I0, shift=(true_flow[1], true_flow[0]), order=1, mode="nearest") I1 += 0.01 * rng.normal(size=I1.shape) Im = (I0 + I1) / 2 Ix = sobel(Im, axis=1) / 8.0 Iy = sobel(Im, axis=0) / 8.0 It = I1 - I0 cy, cx = n // 2, n // 2 r = 9 sl = np.s_[cy-r:cy+r+1, cx-r:cx+r+1] A = np.column_stack([Ix[sl].ravel(), Iy[sl].ravel()]) b = -It[sl].ravel() # Optional Gaussian window weights wy, wx = np.mgrid[-r:r+1, -r:r+1] w = np.exp(-(wx**2 + wy**2) / (2 * (r/2)**2)).ravel() Aw = A * np.sqrt(w[:, None]) bw = b * np.sqrt(w) flow_hat, *_ = np.linalg.lstsq(Aw, bw, rcond=None) G = Aw.T @ Aw evals = np.linalg.eigvalsh(G) print("True flow (u, v):", true_flow) print("Estimated flow:", flow_hat) print("Structure tensor eigenvalues:", evals) print("Condition number:", evals[-1] / max(evals[0], 1e-12)) ``` ## 11.20 Python example: Horn--Schunck iteration The classical Horn--Schunck update uses local averages of the current flow estimates. A simple implementation is useful pedagogically. ```{python} #| label: py-ch11-horn-schunck #| eval: false import numpy as np import matplotlib.pyplot as plt from scipy.ndimage import convolve, sobel, shift, gaussian_filter n = 96 y, x = np.mgrid[-1:1:complex(n), -1:1:complex(n)] I0 = np.exp(-((x + 0.2)**2 + y**2) / 0.08) I0 += 0.6 * np.exp(-((x - 0.3)**2 + (y + 0.2)**2) / 0.04) I0 = gaussian_filter(I0, 1.0) I1 = shift(I0, shift=(-1.0, 2.0), order=1, mode="nearest") Ix = sobel((I0 + I1) / 2, axis=1) / 8.0 Iy = sobel((I0 + I1) / 2, axis=0) / 8.0 It = I1 - I0 alpha = 1.0 num_iter = 200 u = np.zeros_like(I0) v = np.zeros_like(I0) avg_kernel = np.array([[1/12, 1/6, 1/12], [1/6, 0.0, 1/6], [1/12, 1/6, 1/12]]) for _ in range(num_iter): u_bar = convolve(u, avg_kernel, mode="nearest") v_bar = convolve(v, avg_kernel, mode="nearest") residual = Ix * u_bar + Iy * v_bar + It denom = alpha**2 + Ix**2 + Iy**2 u = u_bar - Ix * residual / denom v = v_bar - Iy * residual / denom step = 8 plt.figure(figsize=(5, 5)) plt.imshow(I0, cmap="gray") plt.quiver( np.arange(0, n, step), np.arange(0, n, step), u[::step, ::step], -v[::step, ::step], color="yellow" ) plt.title("Horn--Schunck flow field") plt.axis("off") plt.show() print("Mean estimated u:", np.mean(u)) print("Mean estimated v:", np.mean(v)) ``` ## 11.21 Python example: OpenCV optical flow For practical projects, students may use OpenCV's implementations. This example uses Farneback dense optical flow. ```{python} #| label: py-ch11-opencv-farneback #| eval: false import cv2 import numpy as np import matplotlib.pyplot as plt # I0 and I1 should be grayscale images with values in [0, 255]. I0_u8 = (255 * I0 / I0.max()).astype(np.uint8) I1_u8 = (255 * I1 / I1.max()).astype(np.uint8) flow = cv2.calcOpticalFlowFarneback( I0_u8, I1_u8, None, pyr_scale=0.5, levels=3, winsize=21, iterations=3, poly_n=5, poly_sigma=1.2, flags=0 ) u = flow[..., 0] v = flow[..., 1] mag = np.sqrt(u**2 + v**2) plt.figure(figsize=(5, 4)) plt.imshow(mag, cmap="magma") plt.colorbar(label="flow magnitude") plt.title("Dense optical-flow magnitude") plt.axis("off") plt.show() ``` The goal of this book is not to treat OpenCV as a black box. Instead, students should understand what assumptions and approximations each implementation uses. ## 11.22 AI component: optical flow in deep learning Modern optical-flow methods learn correspondences from data. They typically combine several ideas: 1. **feature extraction:** encode images into learned feature maps; 2. **matching or correlation:** compare features across frames; 3. **cost volume:** store matching evidence for candidate displacements; 4. **regularization/refinement:** enforce spatial coherence; 5. **iterative updates:** repeatedly improve the flow estimate. A learned optical-flow model can be written schematically as $$ \widehat w = f_\theta(I_0,I_1), $$ where $f_\theta$ is trained from supervised flow labels or self-supervised photometric objectives. ### Supervised loss If ground-truth flow $w^*$ is available, a common loss is endpoint error: $$ \mathcal{L}_{\text{EPE}}(\theta) = \frac{1}{|\Omega|}\sum_{x\in\Omega} \left\|f_\theta(I_0,I_1)(x)-w^*(x)\right\|_2. $$ ### Self-supervised photometric loss Without ground-truth flow, one can warp $I_1$ back to $I_0$: $$ \widehat I_0(x)=I_1(x+\widehat w(x)). $$ Then minimize $$ \mathcal{L}_{\text{photo}} = \sum_x \rho\left(I_1(x+\widehat w(x))-I_0(x)\right). $$ This is the same brightness-constancy idea, now embedded inside a differentiable learning pipeline. ### Cost volume A cost volume compares features $F_0$ and $F_1$: $$ C(x,d)=\langle F_0(x),F_1(x+d)\rangle, $$ where $d$ ranges over candidate displacements. The cost volume makes matching explicit and is central in many learned correspondence architectures. ### RAFT-style iterative refinement RAFT-style models build a correlation volume and iteratively update flow estimates using a recurrent update operator. The schematic form is $$ w^{(k+1)}=w^{(k)}+\Delta_\theta\left(w^{(k)}, C, h^{(k)}\right), $$ where $C$ is a correlation representation and $h^{(k)}$ is a hidden state. The important mathematical point is that learned methods still reflect classical structure: matching, regularization, warping, and iterative correction. ## 11.23 TensorFlow example: differentiable photometric loss ```{python} #| label: py-ch11-tensorflow-photometric-loss #| eval: false import tensorflow as tf def bilinear_sample(image, x, y): """Minimal bilinear sampler for one image, pedagogical only. image: tensor of shape [H, W, C] x, y: floating coordinate tensors of shape [H, W] """ H = tf.shape(image)[0] W = tf.shape(image)[1] x0 = tf.cast(tf.floor(x), tf.int32) x1 = x0 + 1 y0 = tf.cast(tf.floor(y), tf.int32) y1 = y0 + 1 x0c = tf.clip_by_value(x0, 0, W - 1) x1c = tf.clip_by_value(x1, 0, W - 1) y0c = tf.clip_by_value(y0, 0, H - 1) y1c = tf.clip_by_value(y1, 0, H - 1) Ia = tf.gather_nd(image, tf.stack([y0c, x0c], axis=-1)) Ib = tf.gather_nd(image, tf.stack([y1c, x0c], axis=-1)) Ic = tf.gather_nd(image, tf.stack([y0c, x1c], axis=-1)) Id = tf.gather_nd(image, tf.stack([y1c, x1c], axis=-1)) x0f = tf.cast(x0, tf.float32) x1f = tf.cast(x1, tf.float32) y0f = tf.cast(y0, tf.float32) y1f = tf.cast(y1, tf.float32) wa = (x1f - x) * (y1f - y) wb = (x1f - x) * (y - y0f) wc = (x - x0f) * (y1f - y) wd = (x - x0f) * (y - y0f) return ( wa[..., None] * Ia + wb[..., None] * Ib + wc[..., None] * Ic + wd[..., None] * Id ) def photometric_loss(I0, I1, flow): """I0, I1: [H, W, C]; flow: [H, W, 2] with components (u, v).""" H = tf.shape(I0)[0] W = tf.shape(I0)[1] yy, xx = tf.meshgrid( tf.range(H, dtype=tf.float32), tf.range(W, dtype=tf.float32), indexing="ij" ) u = flow[..., 0] v = flow[..., 1] I1_warped = bilinear_sample(I1, xx + u, yy + v) return tf.reduce_mean(tf.abs(I1_warped - I0)) ``` This example shows why differentiable warping is central in learned optical flow, video prediction, spatial transformers, and self-supervised depth estimation. ## 11.24 Evaluation metrics Given estimated flow $\widehat w(x)$ and ground truth $w^*(x)$, the most common metric is **endpoint error**: $$ \operatorname{EPE}(x)=\|\widehat w(x)-w^*(x)\|_2. $$ The average endpoint error is $$ \operatorname{AEPE} = \frac{1}{|\Omega|}\sum_{x\in\Omega}\|\widehat w(x)-w^*(x)\|_2. $$ Another metric is angular error. Treat the flow as a 3D vector $(u,v,1)$ and compare directions: $$ \operatorname{AE}(x) = \arccos \frac{\widehat u u^*+\widehat v v^*+1} {\sqrt{\widehat u^2+\widehat v^2+1}\sqrt{(u^*)^2+(v^*)^2+1}}. $$ Evaluation should also consider: - errors near motion boundaries; - errors in occluded regions; - small-motion versus large-motion cases; - textureless regions; - computational cost; - robustness to illumination changes. ## 11.25 Common failure modes | Failure mode | Mathematical cause | Possible remedy | |---|---|---| | flat region | $\nabla I\approx 0$ | feature priors, regularization, larger context | | straight edge | rank-deficient structure tensor | corners, texture, global smoothness | | large motion | Taylor expansion invalid | pyramids, matching, learned cost volume | | occlusion | no true correspondence | occlusion masks, forward-backward checks | | illumination change | brightness constancy violated | gradient constancy, learned features | | motion boundary | smoothness prior incorrect | robust/edge-aware regularization | | repetitive pattern | ambiguous matching | larger context, semantic features | | nonrigid motion | simple parametric model too weak | dense deformable models | ## 11.26 Summary Optical flow is a foundational problem in video analysis. The basic equation $$ I_xu+I_yv+I_t=0 $$ is simple but underdetermined. The main mathematical issue is not writing the equation; it is making the inverse problem well-posed. Lucas--Kanade adds a local constant-motion assumption and solves a least-squares problem. Horn--Schunck adds a global smoothness prior and solves a variational problem. Robust methods handle outliers and motion boundaries. Pyramids and warping handle larger displacement. Modern learned methods add feature learning, cost volumes, and iterative refinement, but they still preserve the classical ideas of correspondence, regularization, and motion consistency. The conceptual message is: $$ \text{motion estimation} = \text{brightness or feature matching} + \text{regularization} + \text{optimization or learning}. $$ ## AI-assisted study prompts Use AI tools as a mathematical tutor, not as a replacement for derivation. 1. Ask an AI system to derive the optical-flow constraint equation. Then identify where differentiability and brightness constancy enter. 2. Ask for an explanation of the aperture problem. Check whether the answer distinguishes normal flow from full flow. 3. Ask the AI to compare Lucas--Kanade and Horn--Schunck in terms of objective functions. 4. Ask for a pseudocode implementation of coarse-to-fine optical flow. Then label which steps are warping, linearization, and residual correction. 5. Ask the AI to explain why learned optical-flow models still use cost volumes or correlations. Relate the answer to classical matching. 6. Ask the AI to generate a failure case for optical flow and classify it as a statistical, geometric, or optimization failure. ## Exercises ### Conceptual exercises 1. Explain the difference between optical flow and physical 3D motion. 2. Why does brightness constancy fail for shadows and specular highlights? 3. What is the aperture problem? Give an example using a moving straight edge. 4. Why are corners better than flat regions for tracking? 5. Explain why motion boundaries create difficulty for smoothness regularization. 6. Compare sparse tracking and dense optical flow. ### Mathematical exercises 1. Starting from $I(x(t),y(t),t)=c$, derive $I_xu+I_yv+I_t=0$. 2. Suppose $\nabla I=(3,4)^T$ and $I_t=-10$. Find the normal flow. 3. For a Lucas--Kanade window with design matrix $A$, prove that the least-squares solution satisfies $A^TA\widehat w=A^Tb$. 4. Show that if all gradient vectors in a window are scalar multiples of one vector, then $A^TA$ is rank one. 5. Derive the Euler--Lagrange equation for the Horn--Schunck $u$ variable. 6. Consider the one-dimensional regularized problem $$ \min_u \sum_i (u_i-z_i)^2+\alpha^2\sum_i(u_{i+1}-u_i)^2. $$ Write the corresponding linear system. 7. Show how endpoint error changes under a global scaling of flow vectors. ### Python exercises 1. Generate two frames containing a shifted Gaussian blob and estimate the shift using Lucas--Kanade. 2. Compare Lucas--Kanade estimates for a corner-like image patch and an edge-like image patch. 3. Implement Horn--Schunck flow and test different values of $\alpha$. 4. Use OpenCV Farneback flow on a short video clip and visualize magnitude and angle. 5. Create a synthetic piecewise-motion scene and show that smooth Horn--Schunck regularization blurs the motion boundary. 6. Implement forward-backward consistency for occlusion detection. ### AI/modeling exercises 1. Design a self-supervised optical-flow loss using photometric error, smoothness, and occlusion masking. 2. Explain how a cost volume $C(x,d)$ relates to nearest-neighbor matching in feature space. 3. Compare a CNN-based optical-flow model with a transformer-based video correspondence model. 4. Ask an AI system to write a RAFT-style pseudocode description. Then rewrite the answer using only mathematical notation. 5. Propose a final project that uses optical flow for medical imaging, traffic analysis, sports analytics, or microscopy. ## Notes and references Classical optical flow begins with variational and differential methods such as Horn--Schunck and Lucas--Kanade [@hornSchunck1981; @lucasKanade1981]. KLT tracking is closely related to local least-squares image alignment and the conditioning ideas in the structure tensor [@shiTomasi1994]. A useful modern survey perspective appears in Szeliski's computer vision text [@szeliski2022computer]. For learned optical flow, FlowNet introduced end-to-end CNN-based flow estimation [@flownet2015], and RAFT introduced a strong recurrent all-pairs field transformation architecture [@raft2020].