Lukas-Kanade步骤的高效matlab实现

Question

我在视频处理课程中获得了一个任务 - 实现了Lucas-Kanade算法。由于我们必须在金字塔模型中做到这一点，因此我首先为每个输入图像构建一个金字塔，然后为每个级别执行一系列LK迭代。在每一个步骤（迭代），下面的代码运行（注意：图片是零填充，所以我可以很容易地处理图像边缘）：

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 

数学我在做什么是计算为每个像素 （I，J）以下光流：

，你可以看到，在代码我计算这对每一个像素，这需要相当长的时间（整个处理2个图像 - 包括建筑物3倍的水平金字塔和3 LK步骤，每个级别都需要大约25秒（！）远程连接到我的大学服务器）。

我的问题：有没有办法计算这个单一的LK步骤没有嵌套for循环？它必须更高效，因为下一步的任务是使用这种算法来稳定短视频。谢谢。

Answer 1

最终我能够找到一个更有效的解决方案来解决这个问题。 它基于问题中显示的公式。最后3行代表了不同之处 - 我们得到了一个无循环的代码，其运行速度更快。与循环版本（根据结果矩阵之间的绝对差异，忽略填充区域）差异可以忽略不计（〜10^-18或更小）。

下面是代码：

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 

Answer 2

我在我的系统上运行了你的代码并进行了分析。这是我得到的。

正如你所看到的反转矩阵（PINV）正在采取的大部分时间。我猜你可以尝试并引导你的代码，但我不知道该怎么做。但我知道一个提高计算时间的技巧。您必须利用矩阵A的最小方差。也就是说，只有在A的​​最小方差大于某个阈值时才计算逆。这将提高速度，因为您不会对所有像素的矩阵进行反转。

您可以通过将您的代码修改为下面显示的代码来实现此目的。

function [du,dv]= LucasKanadeStep(I1,I2,WindowSize) 
It = double(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); 
     % +++++++++++++++++++++++++++++++++++++++++++++++++++ 
     % Code I added , threshold better be outside the loop. 
     lambda = eig(A); 
     threshold = 0.2 
     if (min(lambda)> threshold) 
      U = A\B; 
      du(i,j) = U(1); 
      dv(i,j) = U(2); 
     end 
     % end of addendum 
     % +++++++++++++++++++++++++++++++++++++++++++++++++++ 


%   U = pinv(A)*B; 
%   du(i,j) = U(1);  
%   dv(i,j) = U(2); 
     end 
    end 
end 

我已将阈值设置为0.2。你可以尝试一下。通过使用特征值技巧，我能够将计算时间从37秒缩短到10秒（如下所示）。使用特征，pinv很难像以前那样占用时间。

希望这有助于。祝你好运:)