% HOW TO USE THIS PROGRAM:
% This program allows you to study any jpg-image, for example, an image of
% your face. The image does not have to be square. 
%
% Down load the image you want to study and put in the same directory as this
% file. Make sure Matlab works in this directory. You can check this
% by typing ls at the Matlab prompt. You should then 
% see both this file and your image file. If not, you can change directory
% by using the command cd. 
%
% For exampel, on a PC, type
%    cd d:
% to go to the D-drive.
%
% Once everything is in the same directory type
%    imrank
% to run this program. Follow the instructions in the command window 
% carefully!
%
%
% By Kristian Sandberg November 19, 2001.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5



clear;
close all;



% ================= READ AND CONVERT IMAGE ========================
%

fname=input('Give file name within single quotes: ');
colorflag=input('Enter 1 for a color image, 0 otherwise: ');

I=imread(fname);          % Read the image.
if colorflag == 1
   I=rgb2gray(I);         % Convert the RGB-format (color) to gray scale format.
end
I=double(I);              % Convert numbers into double precision (real)
			  % numbers.

figure(1)
imshow(mat2gray(I))       % Display the image. mat2gray is not necessary but
			  % sometimes useful. This command re-scales the image
			  % so that all values fall within 0 and 1. 
title(['Gray scale version of ' fname])

%
% ==================================================================




% ============== SINGULAR VALUE DECOMPOSITION OF IMAGE =============
%


disp('==========================================================')
disp('We will now study the singular value decompostion of the image.')
disp(' ')
disp('Press any key to compute the singular value decomposition.')
pause
disp('Please wait...')

[U S V]=svd(I,0);

disp('The singular value decomposition has been computed.')
disp('The output contains three matrices U, S and V.')
whos
disp(' ')
disp('Press any key to continue')
pause


disp('==========================================================')
disp('We will now look at the singular values.')
disp('The singular values are given along the diagonal of S.')
disp('Notice the rapid decay!')

figure(2)
plot(diag(S))
title('The singular values of the image')
ylabel('Magnitude of singular values')

disp(' ')
disp('Press any key to continue.')
pause

disp('==========================================================')
disp('The columns of U contain an orthogonal basis for the ') 
disp('column space of the image.')
disp('The columns of V contain an orthogonal basis for the ')
disp('row space of the image.')
disp(' ')
disp('Press any key to continue.')
pause


%
% ==================================================================




% ============== STUDY RANK REDUCED VERSIONS OF THE IMAGE ==========
%



disp('==========================================================')
disp('Let us look at a rank one approximation of the image.')


Ssp=sparse(S);

[M,N]=size(U);
Utemp=zeros(M,N);
Utemp(:,1)=U(:,1);

[M,N]=size(V);
Vtemp=zeros(M,N);
Vtemp(:,1)=V(:,1);

Irank1=Utemp*Ssp*Vtemp';
figure(3)
imshow(mat2gray(Irank1))
title('A rank one approximation of the image')

disp('Note that all columns are just multiples of a single column vector!')

disp(' ')
disp('Press any key to continue.')
pause


disp('==========================================================')
disp('Let us look at a rank two approximation of the image.')


[M,N]=size(U);
Utemp=zeros(M,N);
Utemp(:,1:2)=U(:,1:2);

[M,N]=size(V);
Vtemp=zeros(M,N);
Vtemp(:,1:2)=V(:,1:2);

Irank2=Utemp*Ssp*Vtemp';
figure(4)
imshow(mat2gray(Irank2))
title('A rank two approximation of the image')

disp('All columns are linear combination of just two column vectors.')


disp(' ')
disp('Press any key to continue.')
pause


disp('==========================================================')
disp('Let us look at a rank four approximation of the image.')


[M,N]=size(U);
Utemp=zeros(M,N);
Utemp(:,1:4)=U(:,1:4);

[M,N]=size(V);
Vtemp=zeros(M,N);
Vtemp(:,1:4)=V(:,1:4);

Irank4=Utemp*Ssp*Vtemp';
figure(5)
imshow(mat2gray(Irank4))
title('A rank four approximation of the image')

disp('All columns are linear combination of just four column vectors.')
disp('Despite using only four basis vectors, you should be able ')
disp('to see some structure in your image.')


disp(' ')
disp('Press any key to continue.')
pause



disp('==========================================================')
disp('Let us look at a rank eight approximation of the image.')


[M,N]=size(U);
Utemp=zeros(M,N);
Utemp(:,1:8)=U(:,1:8);

[M,N]=size(V);
Vtemp=zeros(M,N);
Vtemp(:,1:8)=V(:,1:8);

Irank8=Utemp*Ssp*Vtemp';
figure(6)
imshow(mat2gray(Irank8))
title('A rank eight approximation of the image')

disp('All columns are linear combination of eight column vectors.')



disp(' ')
disp('Press any key to continue.')
pause
disp('==========================================================')
disp('Now choose your own rank!')
Nrank=input('Enter rank you want to study: ')


Ssp=sparse(S);

[M,N]=size(U);
Utemp=zeros(M,N);
Utemp(:,1:Nrank)=U(:,1:Nrank);

[M,N]=size(V);
Vtemp=zeros(M,N);
Vtemp(:,1:Nrank)=V(:,1:Nrank);

Irank=Utemp*Ssp*Vtemp';
figure(7)
imshow(mat2gray(Irank))
title(['A rank eight ', num2str(Nrank), ' approximation of the image.'])




disp(' ')
disp('Press any key to continue.')
pause
disp('==========================================================')
disp('Now try to explain everything you have seen!')
disp('In particular, how does the rapid decay of the singular ')
disp('values relate to the fact that even a relative low rank ')
disp('can reproduce a decent approximation of the image?')