clear all, close all, clc

% Use RPCA to remove noise, and then perform the eigenface example

load allFaces.mat


% Apply RPCA to the first 10 people
updatedFaces = zeros(m*n, sum(nfaces(1:10)));
for k=1:10
    baseind = 0;
    if (k>1)
        baseind = sum(nfaces(1:k-1));
    end
    ind = (1:nfaces(k)) + baseind;
    Y = faces(:,ind);
    [L,SS] = RPCA(Y);
    updatedFaces(:,ind) = L;
end

% We use the first 10 people for training data
trainingFaces = updatedFaces;
%trainingFaces = faces(:,1:sum(nfaces(1:10)));

avgFace = mean(trainingFaces,2);  % size n*m by 1;

% Compute eigenfaces on mean-subtracted training data
X = trainingFaces-avgFace*ones(1,size(trainingFaces,2));
[U,S,V] = svd(X,'econ');

figure(1)
imagesc(reshape(avgFace,n,m)), colormap(gray) % Plot avg face
figure(2)
imagesc(reshape(U(:,1),n,m)), colormap(gray)  % Plot first eigenface

%% Now show eigenface reconstruction of image that was omitted from test set

baseind = sum(nfaces(1:10));
ind = (1:nfaces(11)) + baseind;
Y = faces(:,ind);
[L,SS] = RPCA(Y);

%testFace = faces(:,1+sum(nfaces(1:10))); % First face of person 11
testFace = L(:,1);

ff = 3;
figure(ff)
imagesc(reshape(testFace,n,m)), colormap(gray)

testFaceMS = testFace - avgFace;
ff = 4;
for r=[30 60 100 200 300] 
    reconFace = avgFace + (U(:,1:r)*(U(:,1:r)'*testFaceMS));
    figure(ff)
    imagesc(reshape(reconFace,n,m)), colormap(gray)
    norm(reconFace-testFace)/norm(testFace)
    ff = ff + 1;
end
