%% Charge une image JPEG
clear;
load myimage.mat
M = double(myimage(:,:,1));

% Decomposition en valeurs singulieres
[U,S,V] = svd(M);

figure(1); clf;
imagesc(M);
axis equal
axis tight
title('original')
print -dpng svd_original.png

%% Genere et affiche les approximations bas-rang successives
p = rank(S);
r = [1 2 4 8 16 32 64 p];
for i = 1:8
    % Approximation avec une matrice de rang r
    idx = find(S);
    SS = zeros(size(S));
    SS(idx(1:r(i))) = S(idx(1:r(i)));

    MM = U*SS*V';

    figure(1); clf;
    imagesc(MM)
    axis equal
    axis tight
    title(['r = ' num2str(r(i))])
    print('-dpng',['svd_rank_' num2str(r(i)) '.png'])
end