% Eigenfaces sample code. 
% Written by: Alexander Stoytchev
%
% The face images were cropped from the PDF version of the paper
% "Beyond Eigenfaces: Probabilistic Matching for Face Recognition"
% by Moghaddam, Wahid, and Pentland. 1998
%
% For a gentle introduction to the math behind eigenfaces see
% Dana H. Ballard (1999). ``An Introduction to Natural Computation 
% (Complex Adaptive Systems)'', Chapter 4, pp 70-94, MIT Press. 


%% NOTE: All face images must have the same size

N = 214; % height of each face image
M = 144; % width of each face image

K = 9;  % number of faces

% Alocate some variable that are needed by the code
Face=cell(K);
dFace=cell(K);
X=cell(K);
eigenFace=cell(K);
reconstructedFace = cell(K);

for i=1:K
   fname= [ '.\faces\F', (i+48), '_a.bmp' ];
   Face{i}=imread(fname);

   %Also, convert to double and scale to [0,1]
   dFace{i}=double(Face{i})/255.0;
end


% Just for visualization; show the first 4 faces
figure
RES=[];
for i=1:K
   RES= [RES, Face{i}];
end
imshow(RES)
imwrite(RES, 'original_faces.jpg');
title('Original Face Images');


%%Calculate Average Face
avgFace = zeros(N, M);
for i=1:K
   avgFace = avgFace + dFace{i};
end
avgFace = avgFace ./ K;

%Show the average Face
figure
imshow(uint8((avgFace)*255))
imwrite(uint8((avgFace)*255), 'average_face.jpg');
title('Average Face');


%% Normalize and convert to vectors row-wise
%% Each of the Xs is 6400x1
for i=1:K
  X{i}=reshape((dFace{i}-avgFace)', 1,N*M)';
end


% Just for visualization; show the normalized Faces
% They have negative values and need to be rescaled.
Normalized=[];
for i=1:K
  Normalized = [Normalized, uint8((dFace{i}-avgFace+1)*255/2.0)];
end
figure
imshow(Normalized)
imwrite(uint8(Normalized), 'normalized_faces.jpg');
title('Normalized Faces');


%%Construct the matrix A ((N*M)x4)  -> A = [ X{1}, X{2}, ... , X{K}]
A=[];
for i=1:K
  A = [A, X{i}];
end
A = A ./K;


%And the matrix R ( K x K = K x(N*M) * (N*M) x K)
R=A' * A;


%Find the Eigenvectors and Eigenvalues of R
[V,D]=eig(R) %% -- returns them in incresing order
%[V,D]=eigs(R) %% - return only top 6 in decreasing order 


%Use them to find the eigenvectors of A*A' (which is (N*M) x (N*M))
U = A *V;

%u1=A*V(:,1);
%u2=A*V(:,2);
%...
%uK=A*V(:,K);

%Convert the eigenvectors back to images (N x M) for visualization
for i=1:K
  eigenFace{i} = reshape(U(:, i), M, N)';
end

%eigenFace{1}=reshape(u1, M, N)';
%eigenFace{2}=reshape(u2, M, N)';
%...
%eigenFace{K}=reshape(uk, M, N)';

figure
EigenImages=[];
for i=1:K
   EigenImages= [EigenImages, uint8((eigenFace{i}+1)*255.0/2.0)];
end
imshow(EigenImages);
imwrite(EigenImages, 'eigen_images.jpg');
title('Eigenfaces')

%% Do a morph between two eigenfaces
eigMorph=cell(5);
imgMorph=cell(5);
RES = [];
for i=1:5
    eigMorph{i} = (5-i)*0.25* U(:, 1) + (i-1)*0.25* U(:, 2); 
    imgMorph{i} = reshape(eigMorph{i}, M, N)';
    RES = [RES, imgMorph{i}];
end
figure 
imshow(uint8((RES+1)*255/2))
title('Morphing in between the 1st and the 2nd eigenfaces');


%
% Reconstruct the original images using the eigenvectors 
%
% For each normalized face calculate its projection into face space
for i=1:K
   omega{i} = U' * X{i};
end

% Uncomment if you want to Reconstruct using only 2 eigenfaces
%D(3,3)=0.000000000001;

%% Reconstruct the face images using the eigenfaces and the face space coordinates
for i=1:K
   temp = U * (omega{i} ./ diag(D)); % divide by the eigenvalues
   reconFace{i} = reshape(temp,  M, N)';
   reconFace{i} = uint8((reconFace{i} + avgFace)*255);
end


figure
original=[];
reconstructed =[];
for i=1:K
   original    =  [original, Face{i}];
   reconstructed = [reconstructed, reconFace{i}];
end;
imshow( [original; reconstructed]);
imwrite([original; reconstructed], 'original_vs_reconstructed.jpg');
title('Top Row: Original Faces;  bottom Row: Reconstructed Faces');


figure
difference    = uint8(abs(double(original) - double(reconstructed)));
imshow( difference);
title('Reconstruction Error');


