%Script file for the eigenface routine.  I haven't
%been able to get permission to distribute the
%actual face data used in the example, so I will 
%try to get some more data to post here.  This
%script can be modified to accommodate your own
%data, however.


%Load the images into Matlab:
X=imread('face01.tif');
Y(:,1)=double(X(:));
X=imread('face02.tif');
Y(:,2)=double(X(:));
X=imread('face03.tif');
Y(:,3)=double(X(:));
X=imread('face04.tif');
Y(:,4)=double(X(:));
X=imread('face05.tif');
Y(:,5)=double(X(:));
X=imread('face06.tif');
Y(:,6)=double(X(:));
X=imread('face07.tif');
Y(:,7)=double(X(:));
X=imread('face08.tif');
Y(:,8)=double(X(:));
X=imread('face09.tif');
Y(:,9)=double(X(:));
X=imread('face10.tif');
Y(:,10)=double(X(:));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Compute the mean face and subtract
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

meanpic=mean(Y')';
H=Y-repmat(meanpic,1,10);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Plotting routines
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(1)
colormap(gray)
for i=1:10
subplot(2,5,i)
imagesc(vectopic(Y(:,i))');
axis square
axis off
end


figure(2)
colormap(gray)
for i=1:6
subplot(3,2,i)
imagesc(vectopic(H(:,i))');
axis square
axis off
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%  The heart of the code is in the next 4 lines!
%%
%%  Think of the 10 64 x 64 images as
%%    4,096 points in 10 dimensional space,
%%    and compute the covariance matrix.
%%  If we let X=USV^T be the SVD of X,
%%    then the vectors in U (4096 dimensions)
%%    are the eigenfaces, and the vectors in V
%%    (10 dimensions) can be used to construct
%%    U.  The following computations first
%%    construct V, then compute U=XVS^(-1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

C=(1/4096)*H'*H;
[V,S,V]=svd(C);

inS=diag(1./sqrt(diag(S(1:9,1:9))));

%In the next computation, 64=sqrt(4096)
U=(1/64)*H*V(:,1:9)*inS;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Next round of commands puts up the
%%  Eigenfaces and mean face.
%%
%%  "vectopic" is not a built in Matlab
%%  command- its purpose is to take the
%   4096-vector and output a 64 x 64 matrix
%   for viewing.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Am=vectopic(meanpic);
A1=vectopic(U(:,1));
A2=vectopic(U(:,2));
A3=vectopic(U(:,3));
A4=vectopic(U(:,4));
A5=vectopic(U(:,5));
A6=vectopic(U(:,6));
A7=vectopic(U(:,7));

figure(3)
colormap(gray);
subplot(2,4,1);
imagesc(Am');
axis square
axis off
subplot(2,4,2);
imagesc(A1');
axis square
axis off
subplot(2,4,3);
imagesc(A2');
axis square
axis off
subplot(2,4,4);
imagesc(A3');
axis square
axis off
subplot(2,4,5);
imagesc(A4');
axis square
axis off
subplot(2,4,6);
imagesc(A5');
axis square
axis off
subplot(2,4,7);
imagesc(A6');
axis square
axis off
subplot(2,4,8);
imagesc(A7');
axis square
axis off


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Reconstruction Section.  The matrix Coef
%%   holds the coefficients of each image in
%%   terms of the eigenface basis.  Recon holds
%%   the 7-term approximations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Coef=U(:,1:7)'*H;
Recon=U(:,1:7)*Coef;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Graphics section:  Change the vector to
%%    an image, and add the mean back in.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:10
R{i}=vectopic(Recon(:,i));
R{i}=R{i}+Am;
Z{i}=vectopic(Y(:,i));
end

figure(4)
colormap(gray)

for i=1:4
subplot(2,4,i)
imagesc(Z{i}');
axis square
axis off
subplot(2,4,i+4)
imagesc(R{i}')
axis square
axis off
end

figure(5)
colormap(gray)
for i=1:4
subplot(2,4,i)
imagesc(Z{i+4}');
axis square
axis off
subplot(2,4,i+4);
imagesc(R{i+4}');
axis square
axis off
end

figure(6)
colormap(gray)
for i=1:2
subplot(2,2,i)
imagesc(Z{i+8}');
axis square
axis off
subplot(2,2,i+2);
imagesc(R{i+8}');
axis square
axis off
end