Implementing Lukas and Kanade’s Optical Flow

By Mikel Rodriguez

Introduction:

Optical flow is a method used for estimating motion of objects across a series of frames. The method is based on an assumption which states that points on the same object location (therefore the corresponding pixel values) have constant brightness over time.

Optical flow can be said to have two components, normal flow and parallel flow. Normal flow can be computed directly without any further constraints. However, Parallel Flow cannot be computed this way. Several methods have been proposed to compute Parallel Flow; in this project a method proposed by Lucas and Kanade is implemented.

Results:

From the original image sequence, two frames were selected (60,61). A Gaussian pyramid of 3 levels was built for each of these frames:

Frame 60 Level 3 Level 2 Level 1

Frame 61 Level 3 Level 2 Level 1

For each level of the Gaussian pyramid, the optical flow was computed for each pair of frames, the following is the resulting plot of the optical flow vectors.

Level 3:

Level 2:

Level 1:

Source Code

%%*****************************************

%% Mikel Rodriguez {mikel at cs dot ucf dot edu} *

%% Lucas and Kanade's Optical flow *

%%*****************************************

function optical_flow();

clear;

gaus_sigma = 1;

%%Image Variables:

%%I precompute the priramyds using my code from the previous assignment

ORIGINAL_IMAGE_1=imread('flower_frame_0060_LEVEL0.jpg');

ORIGINAL_IMAGE_2=imread('flower_frame_0061_LEVEL0.jpg');

[height,width]=size(ORIGINAL_IMAGE_1);

ORIGINAL_IMAGE_1=im2double(ORIGINAL_IMAGE_1);

ORIGINAL_IMAGE_2=im2double(ORIGINAL_IMAGE_2);

IMAGE_1_SMOOTHED=zeros(height,width);

IMAGE_2_SMOOTHED=zeros(height,width);

%%Derivate Variables:

Dx_1=zeros(height,width);

Dy_1=zeros(height,width);

Dx_2=zeros(height,width);

Dy_2=zeros(height,width);

Ix=zeros(height,width);

Iy=zeros(height,width);

It=zeros(height,width);

%%Optical flow variables

neighborhood_size=5;

A=zeros(2,2);

B=zeros(2,1);

output1=zeros(height,width);

output2=zeros(height,width);

%%Kernel Variables:

Kernel_Size = 6*gaus_sigma+1;

k = (Kernel_Size-1)/2;

gaus_kernel_x=zeros(Kernel_Size,Kernel_Size);

gaus_kernel_y=zeros(Kernel_Size,Kernel_Size);

kernel=zeros(Kernel_Size,Kernel_Size);

%Make a kernel for partial derivatve

%Of of gaussian with respect to x (for computing Dx)

for i=1:Kernel_Size

for j=1:Kernel_Size

gaus_kernel_x(i,j) = -( (j-k-1)/( 2* pi * gaus_sigma^3 ) ) * exp ( - ( (i-k-1)^2 + (j-k-1)^2 )/ (2*gaus_sigma^2) );

end

end

%Make a kernel for partial derivatve

%Of of gaussian with respect to y (for computing Dy)

for i=1:Kernel_Size

for j=1:Kernel_Size

gaus_kernel_y(i,j) = -( (i-k-1)/( 2* pi * gaus_sigma^3 ) ) * exp ( - ( (i-k-1)^2 + (j-k-1)^2 )/ (2*gaus_sigma^2) );

end

end

%%Compute x and y derivates for both images:

Dx_1 = filter2(gaus_kernel_x,ORIGINAL_IMAGE_1);

Dy_1 = filter2(gaus_kernel_y,ORIGINAL_IMAGE_1);

Dx_2 = filter2(gaus_kernel_x,ORIGINAL_IMAGE_2);

Dy_2 = filter2(gaus_kernel_y,ORIGINAL_IMAGE_2);

Ix = (Dx_1 + Dx_2) / 2;

Iy = (Dy_1 + Dy_2) / 2;

%%Build a gaussian kernel to smooth images for computing It

for i=1:Kernel_Size

for j=1:Kernel_Size

kernel(i,j) = (1/(2*pi*(gaus_sigma^2))) * exp (-((i-k-1)^2 + (j-k-1)^2)/(2*gaus_sigma^2));

end

end

IMAGE_1_SMOOTHED = filter2(kernel,ORIGINAL_IMAGE_1);

IMAGE_2_SMOOTHED = filter2(kernel,ORIGINAL_IMAGE_2);

It = IMAGE_2_SMOOTHED - IMAGE_1_SMOOTHED;

for i=(1+floor(neighborhood_size/2)):(height-floor(neighborhood_size/2))

for j=(1+floor(neighborhood_size/2)):(width-floor(neighborhood_size/2))

A=zeros(2,2);

B=zeros(2,1);

for m=i-floor(neighborhood_size/2):i+floor(neighborhood_size/2)

for n=j-floor(neighborhood_size/2):j+floor(neighborhood_size/2)

A(1,1)=A(1,1) + Ix(m,n)*Ix(m,n);

A(1,2)=A(1,2) + Ix(m,n)*Iy(m,n);

A(2,1)=A(2,1) + Ix(m,n)*Iy(m,n);

A(2,2)=A(2,2) + Iy(m,n)*Iy(m,n);

B(1,1)=B(1,1) + Ix(m,n)*It(m,n);

B(2,1)=B(2,1) + Iy(m,n)*It(m,n);

end

end

Ainv=pinv(A); %%Pseudo inverse

result=Ainv*(-B);

output1(i,j)=result(1,1);

output2(i,j)=result(2,1);

end

end

output1=flipud(output1); % Flip matrix in up/down direction

output2=flipud(output2); % Same

quiver(output1, output2); %plot optical flow vectors as arrows