Thursday, June 3, 2010

Tracking using Optical Flow



Optical flow with Hierarchical Lucas Kanade (using pyramids)

results are very strange!














Horn’s optical flow algorithm. Results are unclear.




Code:

function [u,v,cert] = HierarchicalLK(im1, im2, numLevels, windowSize, iterations, display)
%HIERARCHICALLK Hierarchical Lucas Kanade (using pyramids)
% [u,v]=HierarchicalLK(im1, im2, numLevels, windowSize, iterations, display)

if (size(im1,1)~=size(im2,1)) | (size(im1,2)~=size(im2,2))
error('images are not same size');
end;

if (size(im1,3) ~= 1) | (size(im2, 3) ~= 1)
error('input should be gray level images');
end;


% check image sizes and crop if not divisible
if (rem(size(im1,1), 2^(numLevels - 1)) ~= 0)
warning('image will be cropped in height, size of output will be smaller than input!');
im1 = im1(1:(size(im1,1) - rem(size(im1,1), 2^(numLevels - 1))), :);
im2 = im2(1:(size(im1,1) - rem(size(im1,1), 2^(numLevels - 1))), :);
end;

if (rem(size(im1,2), 2^(numLevels - 1)) ~= 0)
warning('image will be cropped in width, size of output will be smaller than input!');
im1 = im1(:, 1:(size(im1,2) - rem(size(im1,2), 2^(numLevels - 1))));
im2 = im2(:, 1:(size(im1,2) - rem(size(im1,2), 2^(numLevels - 1))));
end;

%Build Pyramids
pyramid1 = im1;
pyramid2 = im2;

for i=2:numLevels
im1 = reduce(im1);
im2 = reduce(im2);
pyramid1(1:size(im1,1), 1:size(im1,2), i) = im1;
pyramid2(1:size(im2,1), 1:size(im2,2), i) = im2;
end;

% base level computation
disp('Computing Level 1');
baseIm1 = pyramid1(1:(size(pyramid1,1)/(2^(numLevels-1))), 1:(size(pyramid1,2)/(2^(numLevels-1))), numLevels);
baseIm2 = pyramid2(1:(size(pyramid2,1)/(2^(numLevels-1))), 1:(size(pyramid2,2)/(2^(numLevels-1))), numLevels);
[u,v] = LucasKanade(baseIm1, baseIm2, windowSize);

for r = 1:iterations
[u, v] = LucasKanadeRefined(u, v, baseIm1, baseIm2);
end

%propagating flow 2 higher levels
for i = 2:numLevels
disp(['Computing Level ', num2str(i)]);
uEx = 2 * imresize(u,size(u)*2); % use appropriate expand function (gaussian, bilinear, cubic, etc).
vEx = 2 * imresize(v,size(v)*2);

curIm1 = pyramid1(1:(size(pyramid1,1)/(2^(numLevels - i))), 1:(size(pyramid1,2)/(2^(numLevels - i))), (numLevels - i)+1);
curIm2 = pyramid2(1:(size(pyramid2,1)/(2^(numLevels - i))), 1:(size(pyramid2,2)/(2^(numLevels - i))), (numLevels - i)+1);

[u, v] = LucasKanadeRefined(uEx, vEx, curIm1, curIm2);

for r = 1:iterations
[u, v, cert] = LucasKanadeRefined(u, v, curIm1, curIm2);
end
end;

if (display==1)
figure, quiver(reduce((reduce(medfilt2(flipud(u),[5 5])))), -reduce((reduce(medfilt2(flipud(v),[5 5])))), 0), axis equal
end


Code:

%Implementation based on Horn’s optical flow algorithm.
function [u, v] = LucasKanade(im1, im2, windowSize);
%LucasKanade lucas kanade algorithm, without pyramids (only 1 level);

%REVISION: NaN vals are replaced by zeros

[fx, fy, ft] = ComputeDerivatives(im1, im2);

u = zeros(size(im1));
v = zeros(size(im2));

halfWindow = floor(windowSize/2);
for i = halfWindow+1:size(fx,1)-halfWindow
for j = halfWindow+1:size(fx,2)-halfWindow
curFx = fx(i-halfWindow:i+halfWindow, j-halfWindow:j+halfWindow);
curFy = fy(i-halfWindow:i+halfWindow, j-halfWindow:j+halfWindow);
curFt = ft(i-halfWindow:i+halfWindow, j-halfWindow:j+halfWindow);

curFx = curFx';
curFy = curFy';
curFt = curFt';


curFx = curFx(:);
curFy = curFy(:);
curFt = -curFt(:);

A = [curFx curFy];

U = pinv(A'*A)*A'*curFt;

u(i,j)=U(1);
v(i,j)=U(2);
end;
end;

u(isnan(u))=0;
v(isnan(v))=0;

%u=u(2:size(u,1), 2:size(u,2));
%v=v(2:size(v,1), 2:size(v,2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [fx, fy, ft] = ComputeDerivatives(im1, im2);
%ComputeDerivatives Compute horizontal, vertical and time derivative
% between two gray-level images.

if (size(im1,1) ~= size(im2,1)) | (size(im1,2) ~= size(im2,2))
error('input images are not the same size');
end;

if (size(im1,3)~=1) | (size(im2,3)~=1)
error('method only works for gray-level images');
end;


fx = conv2(double(im1),double(0.25* [-1 1; -1 1])) + conv2(double(im2), double(0.25*[-1 1; -1 1]));
fy = conv2(double(im1), double(0.25*[-1 -1; 1 1])) + conv2(double(im2), double(0.25*[-1 -1; 1 1]));
ft = conv2(double(im1), double(0.25*ones(2))) + conv2(double(im2),double( -0.25*ones(2)));

% make same size as input
fx=fx(1:size(fx,1)-1, 1:size(fx,2)-1);
fy=fy(1:size(fy,1)-1, 1:size(fy,2)-1);
ft=ft(1:size(ft,1)-1, 1:size(ft,2)-1);



TO DO: More experiments




5 Frame difference optical flow using hierarchical lucas kanade on stabilized image sequence.
[u,v,w] = HierarchicalLK(ai1,ai2,3,1,1,1);
results show noise.




2 consecutive frames, stablized vs. Original
stabilized has more inconsistent movement! (Results are inconclusive)

1 comment:

  1. you haven't posted the implementation of this function LucasKanadeRefined can you help me ?

    ReplyDelete