%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ImageEigenvectors.m 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This script (by MPeszynska, for MTH 420/520) is based on the book 
%% "Data Driven Modeling & Scinetific Computation", by J.N. Kutz
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% The spectral decomposition and reconstruction could be done by SVD
%% instead of with eigs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% The script assumes you have various images p?.bmp and others.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% To use the script, you must have the images in your folder 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ImageEigenvectors
function pimage(A,imtitle)
    pcolor(flipud(reshape(A,120,80))),shading interp, colormap(gray);
    if nargin>1,title(imtitle);end
end
function bimage(p,btitle)  
    bar(p(2:end));ylim([min(p) max(p)]);
    if nargin>1,title(btitle);end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% First, read the images for the training set
B = zeros(20,120*80);
for j=1:20 
     pname=sprintf('p%d',j);
     A=imresize(double(rgb2gray(imread(pname,'bmp'))),[120 80]);
     pimage(A,sprintf('Training image %d',j));pause
     B(j,:)=reshape(A,1,120*80);
end
%%
%% Show the size of the image matrix
%%
fprintf('Images read: %d x %d\n',size(B));
%% 
%% construct the "average" image of each person
%%
fprintf('Showing average faces\n');
clooney  = sum(B(1:5,:))/5;pimage(clooney,'Clooney ave');pause
obama    = sum(B(6:10,:))/5;pimage(obama,'Obama ave');pause
thatcher = sum(B(11:15,:))/5;pimage(thatcher,'Thatcher ave');pause
damon    = sum(B(16:20,:))/5;pimage(damon,'Damon ave');pause
%%
%% Now extract (dimension reduction) information 
%% and reconstruct the images
%% 
%%
%% Construct the correlation matrix
%%
fprintf('Computing correlation matrix ...');
tic,
C = B'*B;
toc,
%% take only first 20 eigenmodes and do spectral dcomposition
fprintf('Performing spectral decomposition');
tic,
[V,D] = eigs(C,20,'lm'); 
toc,
%%
%% plot the eigenvalues from largest to smallest
%%
plot(diag(D),'*');title('Image eigenvalues');
pause;
%%
%% reconstruct average images with the averages
%% 
%% find projection (coefficients of expansion of the image in the V basis)
myimage = clooney;
proj = myimage*V; 
figure(1);bar(proj);ylim([min(proj) max(proj)]);pause
bimage(proj,sprintf('Clooney SIGNATURE. First eig=%g',proj(1)));
pause
%% reconstruction
recon = V*proj';
figure(2);pimage(myimage,'Clooney training ave');
figure(3);pimage(recon,'Clooney reconstructed');
pause;
%% repeat for other faces
myimage = obama;
proj = myimage*V; 
figure(1);bimage(proj,sprintf('Obama SIGNATURE. First eig=%g',proj(1)));
recon = V*proj';
figure(2);pimage(myimage,'Clooney training ave');
figure(3);pimage(recon,'Clooney reconstructed');
pause;
%%
%% repeat for other faces
myimage = thatcher;
proj = myimage*V; 
figure(1);bimage(proj,sprintf('Thatcher SIGNATURE. First eig=%g',proj(1)));
recon = V*proj';
figure(2);pimage(myimage,'Thatcher training ave');
figure(3);pimage(recon,'Thatcher reconstructed');
pause;
%%
myimage = damon;
proj = myimage*V; 
figure(1);bimage(proj,sprintf('Damon SIGNATURE. First eig=%g',proj(1)));
recon = V*proj';
figure(2);pimage(myimage,'Damon training ave');
figure(3);pimage(recon,'Damon reconstructed');
pause;
%%
clf(1);clf(2);clf(3);
%%
%% Now try image recognition for faces NOT in the training set
%%
%% read a new image
fprintf('Processing a new image/face not in the raining set');
streep=reshape(imresize(double(rgb2gray(imread('streep','jpg'))),[120 80]),1,120*80);
figure(1);pimage(streep,'Streep as Thatcher');
pause
%% check its coefficients in V basis
%% compare to others
figure(1);bimage(clooney*V,'Clooney');
figure(2);bimage(obama*V,'Obama');
figure(3);bimage(thatcher*V,'Thatcher');
figure(4);bimage(damon*V,'Damon');
figure(5);bimage(streep*V,'Streep as Thatcher');
pause
%% now compare reconstructions
figure(1);pimage(V*(clooney*V)','Clooney');
figure(2);pimage(V*(obama*V)','Obama');
figure(3);pimage(V*(thatcher*V)','Thatcher');
figure(4);pimage(V*(damon*V)','Damon');
figure(5);pimage(V*(streep*V)','Streep as Thatcher');
pause
%%
%% Now check yourself ! 
%% 
me = reshape(imresize(double(rgb2gray(imread('MP','jpg'))),[120 80]),1,120*80);
figure(1);pimage(me,'ME');
pause
%% check its coefficients in V basis
%% compare to others
figure(1);bimage(clooney*V,'Clooney');
figure(2);bimage(obama*V,'Obama');
figure(3);bimage(thatcher*V,'Thatcher');
figure(4);bimage(damon*V,'Damon');
figure(5);bimage(me*V,'ME');
pause
%% now compare reconstructions
figure(1);pimage(V*(clooney*V)','Clooney');
figure(2);pimage(V*(obama*V)','Obama');
figure(3);pimage(V*(thatcher*V)','Thatcher');
figure(4);pimage(V*(damon*V)','Damon');
figure(5);pimage(V*(me*V)','ME');
fprintf('THE IMAGE RECOGNITION IS POOR.\n');
fprintf(' W H Y ???? What can you do to improve ?\n');
pause
for j=1:5, clf(j);end
end