%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   matlab script to calculate the SVD of an image
%   This script first reads in an image (a detail of an etching by Durer of a
%   magic square) and plots it.
%
%      It then takes the SVD of the  image and calculates the
%     spectrum of singular values.  Note that only the first 50 or so of the
%     singular values are large (and really only the first 4 or so).  Because the
%     high singular values are negligible, we can reconstruct much of the image just using
%    the first 50 "Empirical orthogonal Functions" (which are just the first 50 Eigenvectors
%    in V.  This script builds up the reconstruction one EOF at a time.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
close all;
%%%%%%
%  load and display the original image
%%%%%%
load detail
[m,n]=size(X);
imagesc(X); colormap(map); axis image; axis off;
set(gca,'fontweight','bold','fontsize',[14]);
title(sprintf('%s, %dx%d image',caption(1,:),m,n));

pause;

%%%%%
%  now calculate the SVD of the image
%     and set VT=V transpose
%%%%%

[U,S,V]=svd(X,0);
VT=V';

%%%%%%
%  plot the spectrum of the Singular values sigma_1-sigma_m
%%%%%%

figure
semilogy(diag(S),'b-o');
set(gca,'fontsize',[16],'fontweight','bold');
title('Singular Values');
grid;

pause;

%%%%%
% Now calculate all of the coefficients of the principal components
%%%%%
C=U*S;  

%%%%%
%  and start reconstructing the image, one component of V at a
%  time up to 50 modes
%%%%%

figure
i=1;
Xi=C(:,1:i)*VT(1:i,:);
compressionPercent = (1-i*(m+n+1)/m/n)*100;
imagesc(Xi); colormap(map); axis image; axis off;
set(gca,'fontsize',[16],'fontweight','bold');
title(sprintf('EOF reconstruction with %d modes, Compression=%g%%',i,compressionPercent));
disp('Hit return for more modes')
pause;
for i=2:50
  Xi=C(:,1:i)*VT(1:i,:);
  compressionPercent = (1-i*(m+n+1)/m/n)*100;
  imagesc(Xi); colormap(map); axis image; axis off;
  set(gca,'fontsize',[16],'fontweight','bold')
  title(sprintf('EOF reconstruction with %d modes, Compression=%g%%',i,compressionPercent));
  pause;
  drawnow;
end
%%%%
%  and now just show for 100 modes and 200 modes
%%%%
figure
i=100;
Xi=C(:,1:i)*VT(1:i,:);
compressionPercent = (1-i*(m+n+1)/m/n)*100;
imagesc(Xi); colormap(map); axis image; axis off;
set(gca,'fontsize',[16],'fontweight','bold');
title(sprintf('EOF reconstruction with %d modes, Compression=%g%%',i,compressionPercent));
pause;

figure
i=200;
Xi=C(:,1:i)*VT(1:i,:);
compressionPercent = (1-i*(m+n+1)/m/n)*100;
imagesc(Xi); colormap(map); axis image; axis off;
set(gca,'fontsize',[16],'fontweight','bold');
title(sprintf('EOF reconstruction with %d modes, Compression=%g%%',i,compressionPercent));