Lucas-Kanade Optical Flow

令I(x,y,t)表示在t時間時座標(x,y)的圖像值,假設dt時間內該目標移動了(dx,dy)

讓右邊經過一階泰勒展開

消去I(x,y,t)可得以下等式

同除dt並令u=(dx/dt),v=(dy/dt),(u,v)即為像素值在該時間的速度,稱之為光流

上式可改寫為

假設只討論小區域內的(u,v)光流分析

改寫上式為

目標是解 min(|Au-b|^2),所以可以得到u矩陣即為所求

clear;
im1=single((rgb2gray(imread('data\46.jpg'))));
im2=single((rgb2gray(imread('data\47.jpg'))));
numLevels=3; %作幾層金字塔
window=9;    %小區域計算光流的區域呎吋
iterations=1;
alpha = 0.001;

hw = floor(window/2);   
t0=clock;
pyramid1 = im1;
pyramid2 = im2;
for i=2:numLevels  %1最大3最小
    im1 = impyramid(im1, 'reduce');
    im2 = impyramid(im2, 'reduce');
    pyramid1(1:size(im1,1), 1:size(im1,2), i) = im1;
    pyramid2(1:size(im2,1), 1:size(im2,2), i) = im2;
end;
for p = 1:3
       im1 = pyramid1(1:(size(pyramid1,1)/(2^(numLevels - p))), 1:(size(pyramid1,2)/(2^(numLevels - p))), (numLevels - p)+1);
       im2 = pyramid2(1:(size(pyramid2,1)/(2^(numLevels - p))), 1:(size(pyramid2,2)/(2^(numLevels - p))), (numLevels - p)+1);
  if p==1
    u=zeros(size(im1));
    v=zeros(size(im1));
    else  
    %resizing
    u = 2 * imresize(u,size(u)*2,'bilinear');   
    v = 2 * imresize(v,size(v)*2,'bilinear');
    end
   for r = 1:iterations
    u=round(u);
    v=round(v);
        for i = 1+hw:size(im1,1)-hw
            for j = 1+hw:size(im2,2)-hw
        %作小區域的圖像patch1,patch2
        patch1 = im1(i-hw:i+hw, j-hw:j+hw);
      
            %moved patch 
            lr = i-hw+v(i,j);
            hr = i+hw+v(i,j);
            lc = j-hw+u(i,j);
            hc = j+hw+u(i,j);
           
                  if (lr < 1)||(hr > size(im1,1))||(lc < 1)||(hc > size(im1,2))  
                  %Regularized least square processing
                  else
                  patch2 = im2(lr:hr, lc:hc);
      
                  fx = conv2(patch1, 0.25* [-1 1; -1 1]) + conv2(patch2, 0.25*[-1 1; -1 1]);
                  fy = conv2(patch1, 0.25* [-1 -1; 1 1]) + conv2(patch2, 0.25*[-1 -1; 1 1]);
                  ft = conv2(patch1, 0.25*ones(2)) + conv2(patch2, -0.25*ones(2));

      
                  Fx = fx(2:window-1,2:window-1)';
                  Fy = fy(2:window-1,2:window-1)';
                  Ft = ft(2:window-1,2:window-1)';

                  A = [Fx(:) Fy(:)];      
                  G=A'*A;
              
                  G(1,1)=G(1,1)+alpha; G(2,2)=G(2,2)+alpha;
                  %U = G反矩陣*A'*(-Ft)
                  U=1/(G(1,1)*G(2,2)-G(1,2)*G(2,1))*[G(2,2) -G(1,2);-G(2,1) G(1,1)]*A'*-Ft(:);
                  u(i,j)=u(i,j)+U(1); v(i,j)=v(i,j)+U(2);
                  end
            end
        end
   end
end
x=135;
y=170;
figure(1)
im1=((rgb2gray(imread('data\46.jpg'))));
imshow(im1);  
 hold on
 plot(x,y,'r');
figure(2)
im1=((rgb2gray(imread('data\47.jpg'))));
imshow(im1);  
 hold on
 plot(x+u(x,y),y+v(x,y),'g');
 u(x,y)
 v(x,y)