% Principal Component Analysis of Human Faces
% Cesar Marcondes, cesar@cs.ucla.edu

function [] = myEigenfaces (work_dir)
% 'D:\\matlabR12\\work\\stats\\face'

% Variables: 
% \Gamma - a particular image
% \GammaSet - the set of all images

% obtaining the training set directory
cd(work_dir);
bmp_dir = dir('face*.bmp');

[n m] = size(bmp_dir);
K = 150;    % separation of training to testing samples
training_set = 150;
testing_set = n;

% load the set of n-dimensional vectors (one per image)
tic;
progressbar
for i=1:1:n
    progressbar(i/n)
    cur_image = bmp_dir(i).name;    % obtain the filename by position
    data = (imread(cur_image));     % read the image as matrix data
    dimension = size(data);         % read the dimensionality of data
    Gamma = reshape(data,dimension(1)*dimension(2),1); % one sample n-dimensional vector
    GammaSet(:,i) = Gamma;          % insert vector (65536) in the column i
end;
disp( ['dataset loaded in ' num2str(toc) 's.'] );

% call the subrotine that calculates the average and variance matrix
tic;
[avgface, PhiSet] = myPhiSet(GammaSet(:,1:training_set));
disp( ['average and covariance matrix generated in ' num2str(toc) 's.'] );

% call the subrotine myPCA that uses the matrix trick from the 
% Paper "Eigenfaces for Recognition" 
tic;
[eigenfaces, eigenvalues] = myPCA(PhiSet, dimension, K);
disp( ['obtaining eigenfaces through PCA in ' num2str(toc) 's.'] );

% Display/Save the mean face
tic; figure;
mean_face_matrix = reshape(avgface,dimension); % reshape to 256x256 figure
% imwrite(mean_face_matrix,'avgface.bmp','bmp');
imshow(255*mean_face_matrix,[0 255]); title('Average Face');
disp( ['display average face in ' num2str(toc) 's.'] );

eigenfaces2 = eigenfaces*eigenvalues;
% display K-eigenfaces (K=20)

% reshape from dimension 65356 back to 255x255 matrix - eigenface
% the matrix could contain negative values and floating numbers
% therefore, in order to visualize, normalize it to [0,255]

tic; figure;
for q=1:1:20
    % set all numbers to positive and as a minimum 0
    % divide by the range to normalize (0-1) and multiply by 255

    eigenfaces2(:,q) = eigenfaces2(:,q) - min(eigenfaces2(:,q));
    eigenfaces2(:,q) = 255*(eigenfaces2(:,q)/range(eigenfaces2(:,q)));
    Eigenface_matrix = reshape(eigenfaces2(:,q),dimension);
    subplot(5,4,q); imshow(Eigenface_matrix, [0 255]); title(['eig ',num2str(q),'']);
end;
disp( ['display the first 20 eigenfaces in ' num2str(toc) 's.'] );

[TestPhiSet] = TestingPhiSet(GammaSet(:,training_set+1:testing_set),avgface);

% Projection Phase and Error Calculation Varying the Range of Eigenfaces
tic;
progressbar
for q=1:1:150
    
    % reconstruction according to PCA equation
    % from Paper "Low-dimensional procedures for the character. of human faces"
    % weights - projection of weights onto the eigenfaces (linear combinations)
    progressbar(q/training_set)
    weights = (TestPhiSet'*eigenfaces(:,1:q))';
    project_images = eigenfaces(:,1:q)*weights;

    % calculate the difference per pixel value between projected faces and original faces
    % Since both are based on the same average face it turns out that I can compare
    % the project_image with the variance vector of the original image from the mean
    
    error_pixel_images = (abs(project_images - TestPhiSet));
    
    % from all 27 image vectors of 65536 pixels we obtain the average error per pixel
    % doing mean(error_pixel_images), we repeat the process by summing together
    % the average error per pixel for every image and dividing by number of images
    
    error(:,q) = sum((mean(error_pixel_images)))/size(error_pixel_images,2);
end;
disp( ['calculate the reconstruction error in ' num2str(toc) 's.'] );

% Plot the reconstruction error per pixel over the number of eigen-faces K
tic; figure; 
plot(255*error, '-');
disp( ['plotting reconstruction error in ' num2str(toc) 's.'] );

% Using 20 eigenfaces reconstruct the remaining 27 testing faces
% Project the TestPhiSet (testing images) onto the eigenfaces subspace
tic; figure; 
weights = (TestPhiSet'*eigenfaces(:,1:20))';
project_images = eigenfaces(:,1:20)*weights;

for q=1:1:(testing_set-training_set)
    reconstructed_images(:,q) = avgface + project_images(:,q);
    reconstructed_images(:,q) = reconstructed_images(:,q) - min(reconstructed_images(:,q));
    reconstructed_images(:,q) = (reconstructed_images(:,q)/range(reconstructed_images(:,q)));
    reconstructed_face_matrix = reshape(reconstructed_images(:,q),dimension);
    subplot(5,6,q);
    % just in case we want to save the reconstructed faces
    % filename = ['reconst_' num2str(q) '_img.jpg'];
    % imwrite(reconstructed_face_matrix,filename,'jpg');
    imshow(255*reconstructed_face_matrix, [0 255]);
end;
disp( ['reconstruct 28 testfaces using 20-eigenfaces in ' num2str(toc) 's.'] );