clear all;

NrOfSubjects=20;
NrOfImages=10;
nbEigenFacesToShow=2;


% L O A D   T R A I N I N G   D A T A 
ImageData = [];
for s=1:NrOfSubjects
    s
    for n=1:NrOfImages
        FileName=strcat('FaceDatabase/S',...
            num2str(s*100+n),'.bmp');
        dummy=imread(FileName);
        [sy,sx]=size(dummy);
        dummy=dummy(1:4:sy,1:4:sx);
        
        [sy,sx]=size(dummy);
        TrainingData(:,n,s)=double(reshape(dummy,sy*sx,1));
        ImageData = [ImageData,TrainingData(:,n,s)];
    end
end

% L O A D   T E S T   D A T A 
for s=1:NrOfSubjects
    FileName=strcat('FaceDatabase/T',...
        num2str(s),'.bmp');
    dummy=imread(FileName);
    [sy,sx]=size(dummy);
    dummy=dummy(1:4:sy,1:4:sx);
    
    [sy,sx]=size(dummy);
    TestData(:,s)=double(reshape(dummy,sy*sx,1));
end

% get the covariance matrix and the eigenvectors
K=cov(ImageData');
[V,D]=eig(K);
[dy,dx]=size(V);


% P L O T   T H E   F E A T U R E   V E C T O R S
figure(1);clf;hold on;

% here we could be using the PCAData matrix, but if the number of
% eigenfaces has been set to one, then the reduced data won't provide
% enough features (1 instead of 2) for plotting in 2D
PlotVector = [V(:,dx)';V(:,dx-1)']; % first two eigenvectors
PlotData = PlotVector*ImageData; % reduced data

% for each subject, plot each face
for s=1:NrOfSubjects
    % trick for giving different colors per class
    r=bitand(s,1)/1;
    g=bitand(s,2)/2;
    b=bitand(s,4)/4;
    d=bitand(s,8)/8;
    
    for n=1:NrOfImages
        i=(s-1)*NrOfImages+n;
        % avoid 7 and 15 to end up with [r g b] = [1 1 1] => white
        if r+g+b+d==4
            r=0;g=0;b=0;d=0; %black
        end
        
        h=plot(PlotData(1,i),PlotData(2,i),'ro');
        set(h,'Color',(.5+.5*d)*[r g b]);
        h=text(PlotData(1,i),PlotData(2,i),num2str(s));
        set(h,'Color',(.5+.5*d)*[r g b]);
    end
end


% S H O W   F I R S T   T W O   E I G E N F A C E S

% Get the first two eigenfaces
for i=1:nbEigenFacesToShow
    % reshape into a matrix of the size of the rescaled training images
    f=reshape(V(:,dx-i+1),sy,sx);
    
    % normalise and rescale
    maxEF=max(max(f));
    minEF=min(min(f));
    im = uint8(255*((f-minEF)/(maxEF-minEF)));
    
    %if Jean-Michel, write, else show
    if(isunix)
        % give it a name, for writing on disk
        fn=strcat('eigenfaces/e',...
            num2str(i),'.jpg')
        imwrite(im,fn,'JPEG'); %write on disk
    else
        figure(2);
        imshow(im); %show on screen
    end
end


% S T A R T   T H E   R E C O G N I T I O N
for nbEigenFaces=1:6
    T=[];
    for i=0:nbEigenFaces-1
        T=[T;V(:,dy-i)'];
    end
    
    %reduce
    PCAData = T * ImageData;
    PCADataTest = T * TestData;
    
    m=[];
    c=[];
    for s=1:NrOfSubjects
        %offset where the eigenfaces for the current subject start and stop
        first=1+(s-1)*NrOfImages;
        last=s*NrOfImages;
        c(:,:,s) = cov(PCAData(:,first:last));
        m(:,s) = mean(PCAData(:,first:last),2);
    end
    
    
    
    
    
    % E U C L I D E A N   D I S T A N C E   M E A S U R E    
    
    %calc euclidean dist. for each subject of test set
    euclideanFRCounter=0; % false recognition counter
    
    % For each (unknown/test) subject, we compute his euclidean distance
    % to the mean of each class
    
    for unknownSubjectIndex=1:NrOfSubjects   % for each unknown/test subject
        for trainingSubject=1:NrOfSubjects   % for each class
            vec = PCADataTest(:,unknownSubjectIndex) - m(:,trainingSubject);
            distances(trainingSubject) = norm(vec); %store all distances
        end
        
        % get the index of the closest (lowest euclidean distance) class
        [d,classIndex] = min(distances);
        eucDistance(unknownSubjectIndex)=classIndex;
        
        % we already know which class *should* be found, we verify the
        % recognition and increment the counter if false
        if(unknownSubjectIndex ~= classIndex)
            euclideanFRCounter = euclideanFRCounter + 1;
        end
    end
    
    euclideanFRCounter
    eucDistance
    
    
    
    % M E A S U R E   W I T H   P R O B A B I L I T I E S
    
    part02FRCounter=0;
    part03FRCounter=0;
    
    % For each (unknown/test) subject, we compute his mahalonobis distance
    % to the cluster of each class and then the probability p(O|s)
    
    for unknownSubjectIndex=1:NrOfSubjects
        unknownSubject=PCADataTest(:,unknownSubjectIndex);
        for subject=1:NrOfSubjects
            start = 1+(subject-1)*NrOfImages;
            stop = subject*NrOfImages;
            subjectData=PCAData(:,start:stop);
            classMean=mean(subjectData');
            classCov=cov(subjectData',1);
            %proba of observing a face O given subject s
            mahaDistance=(unknownSubject - classMean')' * inv(classCov) * (unknownSubject -classMean');
            prob(unknownSubjectIndex,subject) = exp(-mahaDistance);
        end
    end
    
    prob2=zeros(size(prob));
    
    for unknownSubjectIndex=1:NrOfSubjects
        prob2(unknownSubjectIndex,:) = prob(unknownSubjectIndex,:)/sum(prob(unknownSubjectIndex,:));
    end
    
    for unknownSubjectIndex=1:NrOfSubjects
        
        [probability03,index03] = max(prob2(:,unknownSubjectIndex));
        [probability02,index02] = max(prob(:,unknownSubjectIndex));
        
        %% UNCOMMENT TO USE A TRESHOLD PROBABILITY
        %   treshold=0.5;
        %   if(probability03 < treshold)
        %     unknownSubjectIndex
        %     index03=-1
        %   end
        
        % here we should multiply by the number of faces of that subject in
        % the training set divided by the number of images. but all subject are
        % equally represented in the training set, it is useless
        
        % count the false recognition for each method
        if(unknownSubjectIndex ~= index02)
            part02FRCounter = part02FRCounter + 1;
        end
        if(unknownSubjectIndex ~= index03)
            part03FRCounter = part03FRCounter + 1;
        end
        
        %store the recognised (highest probability) face for each method
        part02Results(unknownSubjectIndex)=index02;
        part03Results(unknownSubjectIndex)=index03;
    end
    
    
    
    part02Results
    part03Results
    
    part02FRCounter
    part03FRCounter
    
    errors(:,nbEigenFaces) = [euclideanFRCounter;part02FRCounter;part03FRCounter]
end

figure(3);
plot(errors');