Search code examples
algorithmmatlabimage-processingvideo-processingopticalflow

efficient matlab implementation for Lukas-Kanade step

I got an assignment in a video processing course - to implement the Lucas-Kanade algorithm. Since we have to do it in the pyramidal model, I first build a pyramid for each of the 2 input images, and then for each level I perform a number of LK iterations. in each step (iteration), the following code runs (note: the images are zero-padded so I can handle the image edges easily):

function [du,dv]= LucasKanadeStep(I1,I2,WindowSize)
It = I2-I1;
[Ix, Iy] = imgradientxy(I2);
Ixx = imfilter(Ix.*Ix, ones(5));
Iyy = imfilter(Iy.*Iy, ones(5));
Ixy = imfilter(Ix.*Iy, ones(5));
Ixt = imfilter(Ix.*It, ones(5));
Iyt = imfilter(Iy.*It, ones(5));
half_win = floor(WindowSize/2);
du = zeros(size(It));
dv = zeros(size(It));
A = zeros(2);
b = zeros(2,1);
%iterate only on the relevant parts of the images
for i = 1+half_win : size(It,1)-half_win  
    for j = 1+half_win : size(It,2)-half_win
          A(1,1) = Ixx(i,j);
          A(2,2) = Iyy(i,j);
          A(1,2) = Ixy(i,j);
          A(2,1) = Ixy(i,j);
          b(1,1) = -Ixt(i,j);
          b(2,1) = -Iyt(i,j);
          U = pinv(A)*b;
          du(i,j) = U(1);      
          dv(i,j) = U(2);
      end
  end
end

mathematically what I'm doing is calculating for every pixel (i,j) the following optical flow: enter image description here

as you can see, in the code I am calculating this for each pixel, which takes quite a long time (the whole processing for 2 images - including building 3 levels pyramids and 3 LK steps like the one above on each level - takes about 25 seconds (!) on a remote connection to my university servers).

My question: Is there a way to calculate this single LK step without the nested for loops? it must be more efficient because the next step of the assignment is to stabilize a short video using this algorithm.. thanks.

Solution

  • Eventually I was able to find a much more efficient solution to this problem. It is based on the formula shown in the question. The last 3 lines are what makes the difference - we get a loop-free code that works way faster. There were negligible differences from the looped version (~10^-18 or less in terms of absolute difference between the result matrices, ignoring the padding zone).

    Here is the code:

    function [du,dv]= LucasKanadeStep(I1,I2,WindowSize)

    half_win = floor(WindowSize/2);
    % pad frames with mirror reflections of itself
    I1 = padarray(I1, [half_win half_win], 'symmetric');
    I2 = padarray(I2, [half_win half_win], 'symmetric');

    % create derivatives (time and space)
    It = I2-I1;
    [Ix, Iy] = imgradientxy(I2, 'prewitt');

    % calculate dP = (du, dv) according to the formula
    Ixx = imfilter(Ix.*Ix, ones(WindowSize));
    Iyy = imfilter(Iy.*Iy, ones(WindowSize));
    Ixy = imfilter(Ix.*Iy, ones(WindowSize));
    Ixt = imfilter(Ix.*It, ones(WindowSize));
    Iyt = imfilter(Iy.*It, ones(WindowSize));

    % calculate the whole du,dv matrices AT ONCE!
    invdet = (Ixx.*Iyy - Ixy.*Ixy).^-1;
    du = invdet.*(-Iyy.*Ixt + Ixy.*Iyt);
    dv = invdet.*(Ixy.*Ixt - Ixx.*Iyt);

end