%NNPCA:Nearest neighbor PCA-based face presentation and recognition
%This demo presents a implemenation of PCA-base face recognition;
%The face database used here is ORL database, for detail information about
%it, one can refer to 

%F. Samaria and A. Harter 
%"Parameterisation of a stochastic model for human face identification"
%2nd IEEE Workshop on Applications of Computer Vision
%December 1994, Sarasota (Florida).

%Also, the eigenfaces and the face represenation using PCA are illustrated.
%To implement this demo, please change the directory according to your
%computer

%Nov.1st,2010


clc
clear all;
allsamples=[];%all training images
alltests=[];
tic
%% input all face sampes
for i=1:40
    for j=1:5
        a=imread(strcat('D:\work\face_PCA\ORL\s',num2str(i),'\',num2str(i*10+j),'.pgm'));
        b=a(1:112*92); % size (b)=1*N,N=10304, first colum, then row.
        b=double(b);
%       size(allsamples)
        allsamples=[allsamples; b];  % allsamples is a matrix of size M*N, every row of allsample is an training image, M=200
    end
end

samplemean=mean(allsamples); % Compute the mean image of all samples; Note, If A is a matrix, mean(A) treats the columns of A as vectors, returning a row vector of mean values
figure(1)
imshow(reshape(samplemean,112,92),[])
title('the mean face')
for i=1:200 
    xmean(i,:)=allsamples(i,:)-samplemean; % xmean is a M*N matrix, of which each row is the sample image minus mean image;
end;

%% generate eigenvectors and eignevalues;
sigma=xmean*xmean';   % M*M matrix;
[v d]=eig(sigma); % produces matrices of eigenvalues (d) and eigenvectors (v) of matrix sigma, so that sigma*v = v*d. 
d1=diag(d);% puts d on the main diagonal
[d2 index]=sort(d1); % sorts the elements along different dimensions of d1, and arranges those elements in ascending order.
cols=size(v,2);% #columns of eigenvectors;
dsort=[];
for i=1:cols
    vsort(:,i) = v(:, index(cols-i+1) ); % vsort size is M*col (note:cols generally equals M). Every column is a eigenvectors in decending order.
    dsort(i)   = d1( index(cols-i+1) );  % dsort 's elements are eigenvalues in decending order. 
end  % ascending order
%preserving 90 percents energy. 
dsum = sum(dsort);
dsum_extract = 0;
p = 0;
while( dsum_extract/dsum < 0.9)
    p = p + 1;
    dsum_extract = sum(dsort(1:p));
end
i=1;
%% Training
% (traning step)
while (i<=p && dsort(i)>0.9)
    base(:,i) = dsort(i)^(-1/2) * xmean' * vsort(:,i);   % base is an N*p matrix, divided by dsort(i)^(-1/2) for normalization
    i = i + 1;
end

% add by wolfsky  projection
allcoor = allsamples * base;
%accu = 0;
%%%%Eigenfaces
%Eigenfaces (the principal components of face space provide a low-dimensional representation of any face
figure(2)
title('eigenfaces')
for f=1:1:6
    eigface=reshape(base(:,f),112,92);
    subplot(2,3,f)
     imshow(eigface,[]);
     title(['eigenface',int2str(f)])
end

%%PCA for face representation
%original face
figure(3)
title('face representation using PCA')
face_original=allsamples(1,:);
subplot(1,2,1)
imshow(reshape(face_original,112,92),[])
title('original face');
c=allcoor(1,:);
recon=zeros(112*92,1);

for r=1:1:10        %PCA approximates a mean-subtracted face, z = x ¨C m, as a weighted sum of the first k eigenfaces
    temp=c(r)*base(:,r);
    recon=recon+temp;
end
subplot(1,2,2)
imshow(reshape(recon+samplemean',112,92),[])
title('Approximating a face using the first 10 eigenfaces')



%% testing step face recognition
accu=0;
for i=1:40
        a=imread(strcat('D:\work\face_PCA\ORL\s',num2str(i),'\',num2str(i*10+j),'.pgm'));
        b=a(1:10304);
        b=double(b);
        tcoor= b * base; %calculate the coordinates
        for k=1:200 
            mdist(k)=norm(tcoor-allcoor(k,:));
        end;
        %third order nearest neighbor
        [dist,index2]=sort(mdist);
        class1=floor(index2(1)/5)+1;
        class2=floor(index2(2)/5)+1;
        class3=floor(index2(3)/5)+1;
        if class1~=class2 && class2~=class3
            class=class1;
        elseif class1==class2
            class=class1;
        elseif class2==class3
            class=class2;
        end;
        if class==i
            accu=accu+1;
        end;
%     end;
end;
accuracy=accu/40 %accuracy
toc
