Contents

% Jana Kosecka, George Mason University, Spring 2005.
% Homework 2, CS 223b, computer vision

% Purely translational flow model based on Kanade-Tomasi tracker
% (Chapter 4 of MASKS)
clear; close all;

prefilt = [0.223755 0.552490 0.223755];
derivfilt = [-0.453014 0 0.45301];

% alternatively you can try to use derivatives of gaussians for
% computing the image derivatives here are some examples of
% of them:
% gaussderiv = [0.0047 0.4905 0 -0.4905 -0.0047 ];
% gauss = [0.0006    0.1262    0.7465    0.1262    0.0006];
% gaussderiv = [0.5000  0   -0.5000];
% gauss = [0.0177 0.9647 0.0177];


% blur serves as weighted sum - weighting the middle pixel most
% you can achieve integration (sumation over a window) by
% convolution with a weight filter

% blur  = [1 6 15 20 15 6 1];
blur    = [1 1 1 1 1 1 1 1 1 1];
blur    = blur / sum(blur);

% must match the filter size for computation of temporal derivatives
nframes = 3;

% Get image sequence
[dimy,dimx]=size(readpgm('yos.10.pnm'));
seq(1).im = zeros(dimy,dimx);
for i=1:nframes
   filename=sprintf('yos.%d.pnm',9+i);
%    disp(['Reading image ',filename]);
   seq(i).im=readpgm(filename);
end

% view motion sequence
% if( 1 )
%         clear M;
%         for i =1 : nframes
%                 %imshow( seq(i).im,256 );
%                 imagesc( seq(i).im);
%                 colormap('gray');
%                 M(:,i) = getframe;
%         end
%         movie(M,-10,30);
% end

Spatial and temporal derivative

f       = zeros(dimy,dimx);

% prefilter in the temporal diection
for i = 1 : nframes
        f = f + prefilt(i)*seq(i).im;
end

% compute derivatives in x and y, prefilter in the other dimension
fx      = conv2( conv2( f, prefilt', 'same' ), derivfilt, 'same' );
fy      = conv2( conv2( f, prefilt, 'same' ), derivfilt', 'same' );

ft      = zeros(dimy,dimx);
% compute derivative in time
for i = 1 : nframes
        ft = ft + derivfilt(i)*seq(i).im;
end

% and prefilter in the other directions
ft      = conv2( conv2( ft, prefilt', 'same' ), prefilt, 'same' );
blur    = [1 6 15 20 15 6 1];
blur    = blur / sum(blur);
fx2     = conv2( conv2( fx .* fx, blur', 'same' ), blur, 'same'  );
fy2     = conv2( conv2( fy .* fy, blur', 'same'  ), blur, 'same'  );
fxy     = conv2( conv2( fx .* fy, blur', 'same'  ), blur, 'same'  );
fxt     = conv2( conv2( fx .* ft, blur', 'same'  ), blur, 'same'  );
fyt     = conv2( conv2( fy .* ft, blur', 'same'  ), blur, 'same'  );

% COMPUTE FLOW (at each s-th pixel)
s       = 4;

[ydim,xdim] = size( fx );
Ux      = zeros( ydim/s, ceil(xdim/s) );
Uy      = zeros( ydim/s, ceil(xdim/s) );
cx      = 1;


%%%%  your code comes here
x = zeros( ydim/s, ceil(xdim/s) );
y = zeros( ydim/s, ceil(xdim/s) );

for i=1:size(Ux, 1)
    % Start and end indices for the ROWS of the window
    i_start = (i-1)*s + 1;
    i_end = i*s;
    for j=1:size(Ux, 2)
        % Start and end indices for the COLUMNS of the window
        j_start = (j-1)*s + 1;
        j_end = j*s;
        % The matrix G (as in part a of the homework)
        G = [sum(sum(fx2(i_start:i_end, j_start:j_end))), ...
             sum(sum(fxy(i_start:i_end, j_start:j_end))); ...
             sum(sum(fxy(i_start:i_end, j_start:j_end))), ...
             sum(sum(fy2(i_start:i_end, j_start:j_end)))];
        % The vector b (as in part a of the homework)
        b = [sum(sum(fxt(i_start:i_end, j_start:j_end))); ...
             sum(sum(fyt(i_start:i_end, j_start:j_end)))];
        % Solving for u and decomposing into x and y components
        u = G\(-b);
        Ux(i,j) = u(1);
        Uy(i,j) = u(2);
        % Locations for the arrows (center of the window)
        x(i,j) = j_end - s/2;
        y(i,j) = i_end - s/2;
    end
end

%%% SHOW FLOW FIELD using quiver
imagesc(seq(2).im);
colormap('gray');
hold on;
quiver(x, y, Ux, Uy, 'r');
set(gcf, 'color', 'w');