algorithm - Lukas-Kanade 步骤的高效 matlab 实现

Question

我在视频处理类(class)中接到了一项任务 - 实现 Lucas-Kanade 算法。由于我们必须在金字塔模型中进行，我首先为 2 个输入图像中的每一个构建一个金字塔，然后为每个级别执行多次 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 秒(!))。

我的问题:有没有办法在没有嵌套 for 循环的情况下计算这个单个 LK 步骤？它必须更有效，因为作业的下一步是使用此算法稳定短视频。谢谢。

Answer 1

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

如您所见，反转矩阵 (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 秒(如下所示)。使用eigen，pinv几乎不像以前那样占用时间。

希望这对您有所帮助。祝你好运 :)