%%

db = 'konstanz';
load(db);

%%
NN  = size(F,1);        % Overal number of face images.
fps = NN/N;             % Number of faces per subject.

% Display faces
exprNum = 2;            % The number of the particular face expression.

nrows = round(sqrt(N));
ncols = ceil(N/nrows);
figure(1); clf; set(gcf, 'Name', 'Faces');
for ii=1:N
    subplot(nrows,ncols,ii);
    tmp = (ii-1)*fps + exprNum;
    imagesc( reshape(F(tmp,:),m,n) );
    colormap gray; axis equal tight off;
    title(num2str(ii));
end

%%
[U S V] = svd(F, 'econ');
PC = U*S;                                % the principle components
figure(2); clf; semilogy(diag(S), '.');

%%
nrows = 2;
ncols = 3;
figure(3); clf; set(gcf, 'Name', 'EigenFaces');
for ii=1:nrows*ncols
    subplot(nrows,ncols,ii);
    imagesc( reshape(V(:,ii),m,n) );
    colormap gray; axis equal tight off
    title(['\sigma_' num2str(ii) ' = ' num2str(S(ii,ii), '%1.2e') ]);
end

%%
nDim = 20;                   % dimensionality of the eigenfaces subspace
V1 = V(:,1:nDim);            % reduced basis of eigenfaces

fn = randint(1,1,[1, NN]);            % Select face at random
f   = F(fn,:)';
% YOUR CODE HERE:
pc1 = PC(fn,1:nDim);                 % project it on the first nDim principle axes
f1 = V1*pc1';                        % reconstructed face
rmse = norm(f-f1)/norm(f);           % rms reprojection error

% Plot the original and reconstructed images
figure(4); clf;
subplot(1,2,1); imagesc(reshape(f,m,n)); title('Original');
axis equal tight off
subplot(1,2,2); imagesc(reshape(f1, m,n));
title(['||\epsilon_{projection}|| = ' num2str(rmse*100,'%2.2f') ' %']);
colormap gray; axis equal tight off

%%
PC_mean = zeros(N, nDim);       % average feature vectors for each subject
for ii=1:N
    ix1 = fps*(ii-1)+1; ix2 = ix1+fps-1;
    PC_mean(ii,:) = mean(PC(ix1:ix2,1:nDim));
end
% Plot faces reconstructed from the feature vectors
nrows = round(sqrt(N));
ncols = ceil(N/nrows);
figure(1); set(gcf, 'Name', 'Face averages');
for ii=1:N
    subplot(nrows,ncols,ii);
    imagesc(reshape(V1*PC_mean(ii,:)',m,n));
    colormap gray; axis equal tight off
    title(num2str(ii));
end

%% Face Recognition
% After we have constructed the feature space, the algorithm for 
% face recognition is straightforward.

% For a (vectorized) face image f:

% * Compute its feature vector, i.e., its coordinates in the 
% (reduced) eigenface basis;
% * Compute the distances between the feature vectors 
% corresponding to f and to the subjects in our database;
% * Find the nearest subject.
% Task 4: Face recognition
% * Implement the face recognition algorithm described above for a 
% test face selected at random from our Face matrix;
% * In a single figure, display the test and the recognized faces;
% * Plot the inter-feature-vectors distances as a bar plot.
% Your plot should look like the one below.
% YOUR CODE HERE:


% Select a face at random and compute its feature vector:
fn = randint(1,1,[1, NN]);
f = F(fn,:);                        % test face
pc1 = V1'*f';                       % its principle component

% Compute distances in the feature space
Dist = PC_mean - ones(N,1)*pc1';
dist = sqrt(mean(Dist.^2, 2));
[min_dist, min_ix] = min(dist);

% Plot
figure(5); clf; set(gcf, 'Name', 'Face recognition')
subplot(2,2,1); imagesc( reshape(f,m,n) );
axis equal tight off; colormap gray;
title('Test face');
subplot(2,2,2); imagesc( reshape( F((min_ix-1)*fps+7,:), m,n) );
axis equal tight off; colormap gray;
title('Recognized as');
subplot(2,2,[3 4]); bar(dist, 'FaceColor', 'b');
xlabel('Subject'); ylabel('Distance [a.u.]');

%% Real-time face recognition
% In reality, the face recognition with this simple 
% approach is not very reliable. Notice the distances in the plot below.

camF = double(t);
camf = camF(:);

pc1 = V1'*camf;                       % its principle component
% Compute distances in the feature space
Dist = PC_mean - ones(N,1)*pc1';
dist = sqrt(mean(Dist.^2, 2));
[min_dist, min_ix] = min(dist);

% Plot
figure(5); clf; set(gcf, 'Name', 'Face recognition')
subplot(2,2,1); imagesc( camF );
axis equal tight off; colormap gray;
title('Camera image');
subplot(2,2,2); imagesc( reshape( F((min_ix-1)*fps+5,:), m,n) );
axis equal tight off; colormap gray;
title('Recognized as');
subplot(2,2,[3 4]); bar(dist, 'facecolor', 'b');
xlabel('Subject'); ylabel('Distance [a.u.]');
