%IMAGE COMPRESSION USING SVD

clear all, close all, clc

% the following reads a bit map image and converts it to numeral format
a=double(imread('lena128.bmp'));
sza=size(a);
nx=sza(1);
ny=sza(2);

% This is the storage space required for the original image (in computer
% words):
disp([num2str(nx*ny) ' numbers represent the original image'])
colormap(gray(256));

disp('See the original image:')
image(a)
pause

figure(2)
colormap(gray(256));
am=mean(a);
disp('See the average image:')
image(ones(ny,1)*am)
pause

a1=double(a)-ones(ny,1)*am;     %subtract the mean
disp('See the image with zero average:')
image(a1)
pause

[u,s,v]=svd(a1); % SVD of the image row vectors

figure(2)
colormap(gray(256));
for k=1:ny/2-1
% the number of eigenvectors used in compression/reconstruction:
   Neigenvectors=k                                  
   a2=ones(ny,1)*am+u(:,1:k)*s(1:k,1:k)*v(:,1:k)'; % reconstructed image
   image(a2);     % display reconstructed compressed image
% nx words to store the mean vector, k*ny words to store matrix U
% k words to store diagonal of S, k*nx words to store matrix V:
disp([num2str(nx+ny*k+k+nx*k) ' numbers are needed to obtain this image.'])
disp(['The average difference between the original and compressed'])
disp(['images is ' num2str(mean(mean(abs(a2-double(a))))) '. Storage saving is ' num2str(100-100.*(nx+ny*k+k+nx*k)/nx/ny) '%'])
   pause
end