l'image de fond

fr       en

Optical flow

Lucas-Kanade optical flow

Method

To solve the system, we assume that for each pixel, all its neighbourhave got the same velocity. For a windows of size m*m, we can then write:

I x 1 V x + I y 1 V y = - I t 1 I x 2 V x + I y 2 V y = - I t 2 . . . I x n V x + I y n V y = - I t n

Where n is the number of pixel into the window. We have a over-determined system. So

[ I x 1 I y 1 I x 2 I y 2 . . . . . . I x n I y n ] ⋅ [ V x V y ] = [ - I t 1 - I t 2 . . . - I t n ]

It is a system easily solve with the least squares methods. With the previous equation we have:

A V ⟶ = - B

Then:

A T A V ⟶ = A T ( - B )

and then:

V ⟶ = ( A T A ) - 1 A T ( - B )

This equation can be solve with:

[ V x V y ] = [ ∑ I x i 2 ∑ I x i I y i ∑ I x i I y i ∑ I y i 2 ] - 1 [ - ∑ I x i I t i - ∑ I y i I t i ]

Program

clear all;
close all;


[filename, pathname] = uigetfile({'*.jpg;*.tif;*.png;*.gif;*.bmp','All Image Files';...
          '*.*','All Files' },'mytitle',...
          'C:\Work\myfile.jpg');

x = imread(filename);
x=double(x);
M1=double(x);
theta = 0.03 * pi/2;
n=256;
% original coordinates
[Y,X] = meshgrid(1:n,1:n);
% rotated coordinates
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
X1 = mod(X1-1,n)+1;
Y1 = mod(Y1-1,n)+1;
% interpolation
M2 = interp2(Y,X,M1,Y1,X1);
M2(isnan(M2)) = 0;

[l c]=size(M1);
l=l-1;
c=c-1;
%derivé temporel
dt = M1 - M2;
dt = dt(1:end-1,1:end-1);
% radient dx dy
dx = diff(M1,1,1);
dx = dx(:,1:end-1);
dy = diff(M1,1,2);
dy = dy(1:end-1,:);

[l c]=size(dt);
n=11;  
M = floor(l/n);
N = floor(c/n);

dm=floor(M/2);
dn=floor(N/2);
v=zeros(n,n,2);
for i = 1 : n
    for j = 1 : n
   
        mask_dx = dx((i*M)-(M-1):(i*M),(j*N)-(N-1):(j*N));
        mask_dx=mask_dx(1:N*M)';
       
        mask_dy = dy((i*M)-(M-1):(i*M),(j*N)-(N-1):(j*N));
        mask_dy=mask_dy(1:N*M)';
       
        mask_dt = dt((i*M)-(M-1):(i*M),(j*N)-(N-1):(j*N));
        mask_dt=mask_dt(1:N*M)';
       
        A=[mask_dx, mask_dy];
        b=mask_dt;
         
        V=inv(A'*A)*A'*(-b);
        v(i,j,1)=V(1);
        v(i,j,2)=V(2);
    end
end
v(isnan(v))=0;

figure(  'Name','Optical Flow: Lucas Kanade',...
            'NumberTitle','off',...
            'color',[0.3137 0.3137 0.5098]);
imshow(uint8(M2))
hold on
quiver([M-dm : M : n*M-dm]'*ones(1,n),ones(1,n)'*[N-dn : N : n*N-dn],v(:,:,1),v(:,:,2),'r')

 

Example

With a rotation,

Copyright © 2010-2014, tous droits réservés, contact : operationpixel@free.fr