%% This demonstrates basic use of MATLAB in image viewing and manipulation 
% and also of the SVD in image compression. 

% Requires JPEG file MsPotatoHead.JPG. 

% Load the JPEG image into an array with dimensions s1xs2xs3. 
MPH = imread('MsPotatoHead.JPG');
[s1,s2,s3] = size(MPH);

% View the image in MATLAB. 
figure(1); imagesc(MPH); axis equal; axis off;
title(['Color image, size ' num2str(s1) 'x' num2str(s2) 'x' num2str(s3)...
    ', storage ' num2str(s1*s2*s3)]);

fprintf('\n Size of full-color array: %d %d %d',s1,s2,s3); 
fprintf('\n Storage = %d\n',s1*s2*s3);

fprintf('\n Strike any key to continue.\n');pause; 

% The image needs to be rotated. Do this by permuting the first and 
% second indexes of the array and "flipping" the matrix left-to-right. 
Z = fliplr(permute(MPH,[2 1 3]));

figure(1); imagesc(Z); axis equal; axis off;
title(['Color image, size ' num2str(s1) 'x' num2str(s2) 'x' num2str(s3)...
    ', storage ' num2str(s1*s2*s3)]);

fprintf('\n Strike any key to continue.\n');pause; 

% Now create a black-and-white s1xs2 image by adding the components 
% along the 3rd index.
W = sum(Z,3);
[s1,s2] = size(W);
fprintf('\n Size of black-and-white array: %d %d ',s1,s2); 
fprintf('\n Storage = %d\n',s1*s2);

% View the B&W image. 
figure(1); imagesc(W); colormap('gray'); axis equal; axis off;
title(['B&W image, size ' num2str(s1) 'x' num2str(s2) ', storage ' num2str(s1*s2)]);

fprintf('\n Strike any key to continue.\n');pause; 

% Compute the SVD of W. This will take some time. 
fprintf('\n Computing SVD ... \n');
tic;[U,S,V] = svd(W);
fprintf('\n SVD compute time = %g\n',toc);

% Explore low-rank approximations of the image. 

fprintf('\n Exploring low-rank approximations of the image.\n');
fprintf('\n Full matrix storage = %d\n',size(W,1)*size(W,2));

rank = 1;

while rank, 
 rank = input('\n Type image-approximation rank or zero to quit.\n');
 if rank ~= 0;
    figure(1);figure(2); 
    imagesc(U(:,1:rank)*S(1:rank,1:rank)*V(:,1:rank)');
    colormap('gray');axis equal; axis off;
    storage = (size(U,1) + size(V,1) + 1)*rank;
    title(['Rank ',num2str(rank),' compressed image, storage ' num2str(storage)]);
    fprintf('\n Storage = %d\n',storage);
 else
     break;
 end
end
