%% Tutorial 3. Face Recognition Analysis using Singular Value Decomposition.
%
% In this script, we will demonstrate how to perform face recognition
% using Principal Component Analysis. Our dataset consists of headshots
% for current CBE graduate students, which will be analyzed to form a
% "training" set containing the most prominent basis facial features for
% the group.
%
% Singular Value Decomposition can then be utilized in the reconstruction
% of a different set of images, called the eigenfaces, made up of
% "proportions" of the features contained in the training set.
%
% For more information about this technique, please refer to the following
% wikipedia article
% http://en.wikipedia.org/wiki/Eigenface
%% 1. Assemble all student pictures into a training set.
%
% The first step in our analysis will be to collect all the photographs and
% store them into an array. Data for color pictures is typically represented
% by specifying intensity of three different channels (red, green and
% blue) and storing them into a matrix composed of three columns. The size
% of the matrix will depend on the number of pixels that the original
% picture contains.

clc
clear
echo on

% Perform a query of all jpeg files in the current directory and determine
% how many images there are in it.
ifiles = dir('*.jpeg');
nimages = length(ifiles);

% All of the pictures are 256 x 256 pixels. Create a matrix x to store the
% data.
x = zeros(256,256);

for i = 1:nimages
    
    % In this loop, we first import the data using the imread command and then
    % turn every picture to black and white. Using this approach it becomes
    % possible to add a linear combination of the red, green and blue intensities
    % of the original data, which reduces our matrix dimension to a single
    % column.
    img = double(imread(ifiles(i).name));
    xn = 0.2126*img(:,:,1) + 0.7152*img(:,:,2) + 0.0722*img(:,:,3);
    
    % It is then possible to store the black and white image data into
    % another array.
    idata(i,:,:) = double(xn);
    
    % Finally, we use our original matrix x to store a running sum of
    % images used later in the script.
    x = x + double(xn);
    echo off
end

% Display a sample of one of the images that we captured.
figure(1);
imagesc(xn);
title('Sample headshot for student')
axis 'off'

%% 2. Compute an average face for a CBE grad student.
% A typical approach used in analyzing a training set is to compute an
% average of the image data to standardize all of the individuals in the
% group. In this part of the script we display what the average CBE
% graduate student would look like.

echo on
% Compute an average of the vector x and then display the results in a figure.
x = x/nimages;
figure(2);
imagesc(x);
title('Average face for a CBE graduate student')
axis 'square'
colormap('gray');
axis 'off'

% In this loop we prepare a matrix that contains color intensity in the
% rows and number of pictures in the columns.
xdata = [];
for i = 1:nimages
    x = idata(i,:,:);
    xdata(:,i) = x(:);
    echo off
end

%% 3. Reconstruct an face from the training set.
% The first approach that we will use for face recognition will employ the
% backslash operator. The idea will be to set up a linear system where
% we calculate what are the best fitted prediction variables that
% reconstruct a given student's face from the calculated mean.

% First, select an image from the training set.
echo on
k = 10;
ki = [1:(k-1),(k+1):nimages];

% Use the backslash operator to approximate the student's face.
A = xdata(:,ki);
b = xdata(:,k);
z = A\b;

% Finally, generate a plot that compares the student with its calculated
% approximation.
z1 = reshape(idata(k,:,:),[256,256]);
z2 = reshape(A*z,[256,256]);

figure(3)
subplot(121)
imagesc(z1);
axis 'square'
title('Original face')

subplot(122)
imagesc(z2);
colormap('gray');
axis 'square'
title('Face recognition using the backslash operator')
%% 4. Eigenface analysis using SVD.

% In this part we recalculate the average graduate student employing the xdata
% matrix obtained from part 2. We then display the recalculated average
% student.
[~,m] = size(xdata);
xmean = sum(xdata')'/m;
figure(4)
imagesc(reshape(xmean,[256,256]));
colormap 'gray'
axis 'square'
axis 'off'
title('Average face for a CBE graduate student')

% Calculate a matrix of deviations. This matrix specifies how far away
% each image is from the mean.
[n,m] = size(xdata);
xdev = zeros(n,m);
for j = 1:m
    xdev(:,j) = xdata(:,j) - xmean;
    echo off
end

% Perform singular value decomposition on the matrix of deviations.
echo on
[u,s,v] = svd(xdev,0);
s = diag(s);

% Construct a subplot containing the average face as well as the first 11
% eigenfaces, which are obtained from the columns of the matrix U
% calculated using the SVD command.
figure(5)
subplot(3,4,1);
imagesc(reshape(xmean,[256,256]));
axis 'square'
axis 'off'
colormap('gray');
title('Average CBE student')

i = 1:11;
for j = i
    subplot(3,4,1+j);
    imagesc(reshape(u(:,j),[256,256]));
    axis 'square'
    axis 'off'
    echo off
    colormap('gray');
    title(['eigenface ',num2str(j)])
end
%% 5. Principal component analysis of a single face.
%
% In the last section of the script, we will test how well our face
% recognition approach works. The basic idea behind this method is that
% starting with the average face we can map any individual by computing the
% dot product between columns in U and the matrix of deviations.

% Select a student from the group at random.
k = round(1 + (nimages-1)*rand);

% Construct a subplot displaying the student's face and its projection.
figure(6)
z1 = reshape(idata(k,:,:),[256,256]);
subplot(121)
imagesc(z1);
axis 'square' 'off'
title('Original student face')

% Calculate a projection utilizing principal component analysis. Display it
% in a figure.
x = u(:,i)'*xdev(:,k);
subplot(122)
imagesc(reshape(xmean + u(:,i)*x,[256,256]));
colormap('gray');
axis 'square' 'off'
title('Student Approximation')
