% 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.2.pnm'));
seq(1).im = zeros(dimy,dimx);
for i=1:nframes
   filename=sprintf('yos.%d.pnm',1+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



%%% SHOW FLOW FIELD using quiver
%%%







