%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Eigenfaces.m
%
% This file is my version of the eigenfaces approach to 
% face recognition.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all, close all
%-------------------------------------------------------------------------
% Load and plot (if desired) the faces. Mean subtract the data since we are
% computing a principle component analysis (PCA).
%-------------------------------------------------------------------------
Dim       = 112;
N_training = 64;
Images=zeros([Dim*Dim N_training]);
disp('1. First we load in and view the faces. Hit any key to continue.'), pause
figure(1)
for I=1:N_training
   name=sprintf('./faces/face%dn.pcx',I);
   [tempIm,map]=imread(name);
   imagesc(tempIm), colormap(gray) 
   title(['Face number ', num2str(I)])
   pause(0.2)
   Images(:,I)=tempIm(:);
end
% Mean subtract the data.
A = (1/N_training)*sum(Images,2);
Data = zeros(size(Images));
for i = 1:N_training
  Data(:,i) = Images(:,i)-A;
end

%-------------------------------------------------------------------------
% The eigenfaces are then the left singular vectors Data. The total power
% in each eigenface over the entire data set is determined by the size of
% the singular values. So we compute the SVD of the data matrix Data.
%-------------------------------------------------------------------------
disp('2. Next, we compute the skinny SVD of the Data. The left singular vectors')
disp('   are the PCA vectors and the singular values are proportional to the')
disp('   square root of the PCA spectral values. Hit any key to continue.'), pause
[U,S,V] = svd(Data,0);
power = diag(S).^2/N_training;
total_power = sum(power);
disp('3. Now, we compress the data by taking a subset of the left singular vectors')
disp('   determined by a choice for the percentage of total power wanted in the basis.')
alpha = input('Percentage of total power = ');
figure(2)
  plot(power), hold on, plot(alpha*ones(size(power)),'r'), hold off
  legend('PCA spectrum','truncation level')
r = 0;
ipower = 0;
while ipower <= alpha & r < N_training
    r = r + 1;
    ipower = ipower + 100*power(r)/total_power;
end 
fprintf('Number of singular vectors/values saved = %d\n',r)
% Eigenfaces/Compressed basis.
Ur = U(:,1:r);

%-------------------------------------------------------------------------
% Note what needs to be stored in order to reconstruct the faces over the
% PCA basis is the matrix Ur and the reconstruction coefficients Ur'*Data.
% New faces can be tested against the library by looking at ||Ur'*face||.
% If this norm is large enough, presumably the face corresponds to one in
% the library. We now plot the eigenfaces and the projections of the faces
% onto the PCA basis. The reduction is storage is as follows: without 
% compression (nxp) where n=112^2 and p=64; with compression Ur (nxr) and
% Ur'*Data (rxp). If r<<p, the storage might be significant. 
%-------------------------------------------------------------------------
disp('4. Now we view the eigenfaces/compressed basis and the projection of')
disp('   the data (original faces) onto this basis. Hit any key to continue.'), pause
figure(3)
for I=1:r
   imagesc(reshape(Ur(:,I),112,112)), colormap(gray) 
   title(['Eigenface, number ', num2str(I)])
   pause(0.2)
end
% Compressed faces.
DataComp = Ur*(Ur'*Data);
figure(4)
for I=1:N_training
   imagesc(reshape(DataComp(:,I)+A,112,112)), colormap(gray) 
   title(['Compressed face, number ', num2str(I)])
   pause(0.2)
end

%-------------------------------------------------------------------------
% Now load a test set of faces that is not included in the training set, 
% mean subtract the data (using the mean computed above), and project the 
% new faces onto reduced PCA bases. Then have a look at how well the faces
% are reconstructed by the bases.
%-------------------------------------------------------------------------
disp('5. Finally, we load in a test set of new images, view them, and view')
disp('   their projections onto the compressed basis. Hit any key to continue.'), pause
N_test = 14;
TestData=zeros([Dim*Dim N_test]);
for I=1:N_test
   name=sprintf('./faces/face%dtn.pcx',I);
   [tempIm,map]=imread(name);
   TestData(:,I)=double(tempIm(:))-A;
end
Projections = U*(U'*TestData);
for i = 1:N_test
  Projections(:,i) = Projections(:,i)+A;
end
% Plot the reconstructions.
figure(5);
for i = 1:N_test
  Ei = reshape(Projections(:,i),Dim,Dim);
  imagesc(Ei), colormap(gray)
  pause(0.2)
end