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很難像以前那樣佔用時間。

希望這有助於。祝你好運:)