Lukas-Kanade 步骤的高效 matlab 实现

Question

我在视频处理课程中接到了一个任务——实现 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几乎不像以前那样占用时间。

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

Answer 2

最终我找到了解决这个问题的更有效的方法。它基于问题中显示的公式。最后 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

Lukas-Kanade 步骤的高效 matlab 实现

efficient matlab implementation for Lukas-Kanade step

algorithm

matlab

image-processing

video-processing

opticalflow