%%
clear
clc

%%
% We make an image from some text:
%
% tobj = text(0, 0, 'MATH 3511', 'fontsize', 90);
% ex = get(tobj, 'extent');
% axis([ex(1) ex(1)+ex(3) ex(2) ex(2)+ex(4)]);
% axis off
% saveas(gcf, 'math3511.png')

%%
% Reload image as a matrix.
% 
% The image has three ``layers'' for red, green, and blue.
% We convert it to a single matrix indicating shades of gray from black (0) to white (255).
% Then we have to tell MATLAB to represent the entries as floating-point
% values rather than integers. 
A = imread('math3511.png');
A = double(rgb2gray(A));
imagesc(A)  
colormap gray
[m, n] = size(A)

%%
% Next we show that the singular values decrease exponentially,
% until they are about sigma_1 * macheps. For all numerical purposes, 
% this determines the rank of the matrix.
[U, S, V] = svd(A);
sigma = diag(S);
semilogy(sigma, '.')
title('singular values')
xlabel('i') 
ylabel('\sigma_i')
grid on
r = find(sigma/sigma(1) > 10*eps, 1, 'last')

%%
% The rapid decrease suggests that we can get fairly good low-rank
% approximations. 
for i = 1:4
    subplot(2, 2, i)
    k = 3*i;
    Ak = U(:,1:k)*S(1:k,1:k)*V(:,1:k)';
    imshow(Ak, [0 255])
    title(sprintf('rank = %d', k))
end


%%
% Consider how little data is needed to reconstruct these images. For rank
% 8, for instance, we have 8 left and right singular vectors plus 8
% singular values, for a compression ratio of better than 25:1.  
compression = 12*(m+n+1) / (m*n)