LucasKanade オプティカルフロー検出＠ MATLAB

この記事はメモです。
中のコードは以下のものをほとんど流用したので自分で書いたとはあまり言い張れないなぁ。
https://jp.mathworks.com/matlabcentral/fileexchange/48744-lucas-kanade-tutorial-example-1/content/LucasKanadeExample1/html/LKExample1.html

関数名：OpticalFlow
入力：前画像，後画像
出力：前画像から後画像へのオプティカルフローを図としてアウトプット

参照元からの変更分

関数への変更
Conv()からfilter2()関数の使用への変更

コード

%% Optical Flow
% This is not my original code
% https://jp.mathworks.com/matlabcentral/fileexchange/48744-lucas-kanade-tutorial-example-1/content/LucasKanadeExample1/html/LKExample1.html
% reviced conv() use to filter2

function OpticalFlow(img1,img2)

im1t = im2double(img1);
im1 = imresize(im1t, 0.5);
im2t = im2double(img2);
im2 = imresize(im2t, 0.5);

ww = 45;
w = round(ww/2);

% Lucas Kanade Here
% for each point, calculate I_x, I_y, I_t
fx=[-1,0,1;-1,0,1;-1,0,1];			% フィルタ行列の作成
fy = [-1,-1,-1;0,0,0;1,1,1];
Ix_m=filter2(fx,im1,'same');		% フィルタリングX方向微分
Iy_m=filter2(fy,im1,'same');		% フィルタリングY方向微分
It_m = im2 - im1;

% Ix_m = conv2(im1,[-1 1; -1 1], 'valid'); % partial on x
% Iy_m = conv2(im1, [-1 -1; 1 1], 'valid'); % partial on y
% It_m = conv2(im1, ones(2), 'valid') + conv2(im2, -ones(2), 'valid'); % partial on t
u = zeros(size(im1));
v = zeros(size(im2));

% within window ww * ww
for i = w+1:size(Ix_m,1)-w
   for j = w+1:size(Ix_m,2)-w
      Ix = Ix_m(i-w:i+w, j-w:j+w);
      Iy = Iy_m(i-w:i+w, j-w:j+w);
      It = It_m(i-w:i+w, j-w:j+w);

      Ix = Ix(:);
      Iy = Iy(:);
      b = -It(:); % get b here

      A = [Ix Iy]; % get A here
      nu = pinv(A)*b; % get velocity here

      u(i,j)=nu(1);
      v(i,j)=nu(2);
   end;
end;

% downsize u and v
u_deci = u(1:10:end, 1:10:end);
v_deci = v(1:10:end, 1:10:end);
% get coordinate for u and v in the original frame
[m, n] = size(im1t);
[X,Y] = meshgrid(1:n, 1:m);
X_deci = X(1:20:end, 1:20:end);
Y_deci = Y(1:20:end, 1:20:end);
figure();
imshow(img2);
hold on;
% draw the velocity vectors
quiver(X_deci, Y_deci, u_deci,v_deci, 'y')

end