Optical Flow Computation¶

Important: Please read the installation page for details about how to install the toolboxes. $\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\enscond}[2]{\lbrace #1, #2 \rbrace}$ $\newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\umax}[1]{\underset{#1}{\max}\;}$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\uargmin}[1]{\underset{#1}{argmin}\;}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\diag}[1]{{diag}\left( #1 \right)}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\qarrq}{ \quad \Longrightarrow \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\CC}{\mathbb{C}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\Zz}{\mathcal{Z}}$ $\newcommand{\Ww}{\mathcal{W}}$ $\newcommand{\Vv}{\mathcal{V}}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\NN}{\mathcal{N}}$ $\newcommand{\Hh}{\mathcal{H}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\Ee}{\mathcal{E}}$ $\newcommand{\Cc}{\mathcal{C}}$ $\newcommand{\Gg}{\mathcal{G}}$ $\newcommand{\Ss}{\mathcal{S}}$ $\newcommand{\Pp}{\mathcal{P}}$ $\newcommand{\Ff}{\mathcal{F}}$ $\newcommand{\Xx}{\mathcal{X}}$ $\newcommand{\Mm}{\mathcal{M}}$ $\newcommand{\Ii}{\mathcal{I}}$ $\newcommand{\Dd}{\mathcal{D}}$ $\newcommand{\Ll}{\mathcal{L}}$ $\newcommand{\Tt}{\mathcal{T}}$ $\newcommand{\si}{\sigma}$ $\newcommand{\al}{\alpha}$ $\newcommand{\la}{\lambda}$ $\newcommand{\ga}{\gamma}$ $\newcommand{\Ga}{\Gamma}$ $\newcommand{\La}{\Lambda}$ $\newcommand{\si}{\sigma}$ $\newcommand{\Si}{\Sigma}$ $\newcommand{\be}{\beta}$ $\newcommand{\de}{\delta}$ $\newcommand{\De}{\Delta}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$

This numerical tour explores the computation of optical flow between two images. It is at the heart of video coding.

In [2]:

addpath('toolbox_signal')
addpath('toolbox_general')
addpath('solutions/multidim_5_opticalflow')

Loading Warped Images¶

To evaluate the performance of optical flow computation, we compute a pair of image obtained by a smooth warping (small deformation), here a simple rotation.

Original frame #1.

In [3]:

n = 256;
name = 'lena';
M1 = rescale( load_image(name,n) );

The second image |M2| is obtaind by rotating the first one.

angle of rotation

In [4]:

theta = .03 * pi/2;

original coordinates

In [5]:

[Y,X] = meshgrid(1:n,1:n);

rotated coordinates

In [6]:

X1 = (X-n/2)*cos(theta) + (Y-n/2)*sin(theta) + n/2;
Y1 =-(X-n/2)*sin(theta) + (Y-n/2)*cos(theta) + n/2;

boundary handling

In [7]:

X1 = mod(X1-1,n)+1;
Y1 = mod(Y1-1,n)+1;

interpolation

In [8]:

M2 = interp2(Y,X,M1,Y1,X1);
M2(isnan(M2)) = 0;

Display the two images.

In [9]:

clf;
imageplot(M1, 'Frame #1', 1,2,1);
imageplot(M2, 'Frame #2', 1,2,2);

Optical Flow Computation with Regularization¶

A first approach to optical flow computation is to solve a ill posed problem corresponding to the optical flow equation constraint (consistency of gray level intensity when moving along the flow).

Compute the derivatives in time and space.

In [10]:

global D;
Dt = M1-M2;
D = grad(M1);

Display them.

In [11]:

clf;
imageplot(Dt, 'd/dt', 1,3,1);
imageplot(D(:,:,1), 'd/dx', 1,3,2);
imageplot(D(:,:,2), 'd/dy', 1,3,3);

The optical flow constraint asks for consistency of gray levels when moving along the flow |v=[v1,v2]|. This is written as a linear equation

|Dt + v1.D1 + v2.D2=0|

This equation does not constrain enough the flow (one equation for two unknown). One thus needs to add other constraints, and this is achieved by performing a Sobolev regularization, as first proposed by Horn and Schunck in the paper:

Horn, B.K.P., and Schunck, B.G., Determining Optical Flow, AI(17), No. 1-3, August 1981, pp. 185-203

This corresponds to a quadratic regularization with a Sobolev prior:

|min_{v} norm(Dt + v1.D1 + v2.D2)^2 + lambdanorm(grad(v1))^2 + lambdanorm(grad(v2))^2|

Its solution is computed by solving a linear system resolution, which sets to zero the gradient of the functional. It can be computed using a gradient descent, or, better, a conjugate gradient descent. We first detail the gradient descent, and shows that is not very efficient.

Regularization strength.

In [12]:

global lambda;
lambda = .1;

Gradient step size.

In [13]:

tau = .2;

Initialization.

In [14]:

v = zeros(n,n,2);

Compute the gradient of the functional. First compute |Dt+v1D1+v2D2|

In [15]:

U = Dt + sum(v.*D,3);

Then compute the Laplacian |L(:,:,k)| of each channel |v(:,:,k)| of the vector field

In [16]:

L = cat(3, div(grad(v(:,:,1))), div(grad(v(:,:,2))));

And gather everything together to build the gradient of the functional.

In [17]:

G = D.*repmat(U, [1 1 2])  - lambda * L;

Perform the descent.

In [18]:

v = v - tau*G;

Exercise 1

Perform the gradient descent of the energy, and display the decay of the energy.

In [19]:

exo1()

In [20]:

%% Insert your code here.

A much faster algorithm is the conjugate gradient. Several variant are implemented within matlab, and can be used by implementing a callback function.

Set up parameters for the CG algorithm (tolerance and maximum number of iterations.

In [21]:

tol = 1e-5;
maxit = 200;

Right hand side of the linear system.

In [22]:

b = -D.*cat(3,Dt,Dt);

Resolution by conjugate gradient.

In [23]:

[v,flag,relres,it,resvec] = cgs(@callback_optical_flow,b(:),tol,maxit);
v = reshape(v, [n n 2]);

Display the flow as a color image and as arrows.

In [24]:

clf;
imageplot(v, '', 1,2,1);
subplot(1,2,2);
w = 12; m = ceil(n/w);
t = w/2 + ((0:m-1)*w);
[V,U] = meshgrid(t,t);
hold on;
imageplot(M1);
quiver(t,t,v(1:w:n,1:w:n,2), v(1:w:n,1:w:n,1));
axis('ij');

Compute the image warped along the flow.

compute the grid, translated along the flow

In [25]:

[Y,X] = meshgrid(1:n,1:n);
X = clamp(X+v(:,:,1),1,n);
Y = clamp(Y+v(:,:,2),1,n);

compute the first fame, translated along the flow

In [26]:

Ms = interp2( 1:n,1:n, M1, Y,X );

One can compare the residual with and without the flow

residual without flow

In [27]:

R0 = M2-M1;

residual along the flow

In [28]:

R = Ms-M1;

ensure same dynamic range (just for display)

In [29]:

v = max( [max(abs(R0(:))) max(abs(R(:)))] );
R(1)=v; R(2)=-v; R0(1)=v; R0(2)=-v;

display

In [30]:

clf;
imageplot(R0, 'Residual without flow', 1,2,1);
imageplot(R, 'Residual with flow', 1,2,2);