% Benjamin Garrett
% CIS 750 - Video Mining Seminar
% May 2003
% Final Project - Pfinder implementation
% Based on algorithms described in:
% Wren et al., Pfinder: Real-Time Tracking of the Human Body. 
% IEEE Transactions on Pattern Analysis and Machine Intelligence
% July 1997, vol 19, no 7, pp. 780-785

clear; close all;

%===================================
% PARAMETERS
%===================================
IMAGES_LOADED = 0; % to skip processing of images
ALPHA = 0.01; % adaptive learning constant
BETA = 0.01; % percentage change to detect movement
K_MODELS = 3; % number models in the mixture
INITIAL_FRAMES = 10; % number of frames to read before looking for changes
NORMALIZATION = 0; % to include or not the section containing normalized values
START_FRAME = 70; % start processing with this frame
STOP_FRAME = 100; % stop processing with this frame
FIRST_MOVEMENT = 90; %this is first frame containing some movement
TWO_PI = 6.283185307179586476925286766559; % equal to 2*pi

% Read in video stream to be analyzed, count frames and get frames size
if (IMAGES_LOADED==0)
    
	movieObj = aviread('security1.avi');
  	noOfFrameImages = size(movieObj,2);
    numrows = size(movieObj(1).cdata, 1)
	numcols = size(movieObj(1).cdata, 2)
	
	% Convert images to double precision for calculations
    % Cut off the first START_FRAME images to save some time
	d_movieObj=[];
	for i = START_FRAME:STOP_FRAME
        d_movieObj(i-(START_FRAME-1)).cdata = double(movieObj(i).cdata);
	end;
  	noOfFrameImages = size(d_movieObj,2);
       
end;

% CCIR 601 Recommendation Standard 
% YUV Mask for converting RGB images to YUV images

% |Y| |  0.299  0.587  0.114 | |R|
% |U|=| -0.169 -0.331  0.500 |x|G|
% |V| |  0.500 -0.419 -0.081 | |B|

% Y = 0.299R + 0.587G + 0.114B 
% Cr = R - Y 
% Cb = B - Y 

%=================================
% INITIALIZATION STAGE
% This part will execute until a 'significant' change has been found in the scene
%=================================

% Convert all images from RGB to YUV and U*V* (luminance normalization) color spaces
if (IMAGES_LOADED==0)

for i = 1:noOfFrameImages
    im.y(:,:,i) = (0.299 * d_movieObj(i).cdata(:,:,1) + 0.587 * d_movieObj(i).cdata(:,:,2) + 0.114 * d_movieObj(i).cdata(:,:,3));                
    im.u(:,:,i) = (d_movieObj(i).cdata(:,:,1) - im.y(:,:,i));
    im.v(:,:,i) = (d_movieObj(i).cdata(:,:,3) - im.y(:,:,i));
%     if (NORMALIZATION==1)
%         if (not(im.y(:,:,i) == 0))
%             im.u_st(:,:,i) = im.u(:,:,i) ./ im.y(:,:,i);
%         end;
%         if (not(im.y(:,:,i) == 0))
%             im.v_st(:,:,i) = im.v(:,:,i) ./ im.y(:,:,i);
%         end;    
%     end;
end;     
  
% Go through the first 20 images and build mean statistics
im.my = mean(im.y, 3);
im.mu = mean(im.u, 3);
im.mv = mean(im.v, 3);
% if (NORMALIZATION==1)
%     im.mu_st = mean(im.u_st, 3);
% 	im.mv_st = mean(im.v_st, 3);
% end;

% Get covariance matrix for first set of frames and use adaptive learning
im.k0=[1 1 1;
       1 1 1;
       1 1 1];
for k=4:INITIAL_FRAMES
    for i=1:numrows
        for j=1:numcols
            w=[im.y(i,j,k)       im.u(i,j,k)      im.v(i,j,k);
               im.y(i,j,(k-1))   im.u(i,j,(k-1))  im.v(i,j,(k-1));
               im.y(i,j,(k-2))   im.u(i,j,(k-2))  im.v(i,j,(k-2))];
            im.k(:,:,k)=cov(w);
            im.k0=(ALPHA.*im.k(:,:,k))+((1-ALPHA).*im.k0);
        end;
    end;
end;

% if (NORMALIZATION==1)
%     w_nor=[im.u_st(:,:,(INITIAL_FRAMES-1))    im.v_st(:,:,(INITIAL_FRAMES-1));
%            im.u_st(:,:,(INITIAL_FRAMES-2))    im.v_st(:,:,(INITIAL_FRAMES-2));
%            im.u_st(:,:,(INITIAL_FRAMES-3))    im.v_st(:,:,(INITIAL_FRAMES-3))];
%     im.k_nor=cov(w_nor)
% end;

end; % end of if(IMAGES_LOADED==0)
     % after an initial run this part will be skipped
     
% For each pixel, compute mean and covariances
% Then check their Mahalanobis distance (dk) in the color space
% NOTE: I've defined here FIRST_MOVEMENT to be where I know a change occurs
%       in the scene.
for k=INITIAL_FRAMES:noOFrameImages
    for i=1:numrows
        for j=1:numcols
            % calculate the mean for each pixel, and for each spectral component
            im.my(i,j,k)=mean([im.y(i,j,k) im.y(i,j,(k-1)) im.y(i,j,(k-2))]);
            im.mu(i,j,k)=mean([im.u(i,j,k) im.u(i,j,(k-1)) im.u(i,j,(k-2))]);
            im.mv(i,j,k)=mean([im.v(i,j,k) im.v(i,j,(k-1)) im.v(i,j,(k-2))]);
            
            % calculate the covariance matrix for each pixel
            w=[im.y(i,j,k)     im.u(i,j,k)     im.v(i,j,k);
               im.y(i,j,(k-1)) im.u(i,j,(k-1)) im.v(i,j,(k-1));
               im.y(i,j,(k-2)) im.u(i,j,(k-2)) im.v(i,j,(k-2))];
            im.k(:,:,k)=cov(w);
            
            % start calculating distance
            p1=(-0.5).*([im.y(i,j,k) im.u(i,j,k) im.v(i,j,k)] - [im.my(i,j,k) im.mu(i,j,k) im.mv(i,j,k)]);
            p2=1./(im.k(:,:,k));
            p3=p1*p2;
            p4=p3*(([im.y(i,j,k) im.u(i,j,k) im.v(i,j,k)] - [im.my(i,j,k) im.mu(i,j,k) im.mv(i,j,k)])');
            p5=(0.5).*(log(abs(im.k(:,:,k))));
            p6=p4-p5;
            p7=(0.5).*(log(TWO_PI));
            dk=p6-p7;            
        end;
    end;
end;

whos
