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