% Convert the picture to black and white.
image = rgb2gray(imread('img.jpg'));

% Downscale the image to simplify.
image = im2double(imresize(image, 0.5));

% Compute the SVD of the image.
[U, S, V] = svd(image);

% The image is now represented as a matrix, 
sigmas = diag(S);
figure('position', [0, 0, 1000, 1000]); imshow(imread('img.jpg')); title('Original Image');
figure('position', [0, 0, 1000, 1000]); plot(sigmas); title('Pixel Intensity (y) vs Singular Value Count (x)');
figure('position', [0, 0, 1000, 1000]); plot(cumsum(sigmas) / sum(sigmas)); title('Image Representation % (y) vs Singular Value Count (x)');

% Show the original.
figure('position', [0, 0, 1000, 1000]); subplot(4, 2, 1), imshow(image), title('Original Image');

% Different ranks of the image.
ranks = [200, 100, 50, 30, 20, 10, 3];
for i = 1:length(ranks)
    % Only use the largest singular values.
    apx_sigmas = sigmas; apx_sigmas(ranks(i):end) = 0;

    % Make the matrix.
    ns = length(sigmas);
    apx_S = S; apx_S(1:ns, 1:ns) = diag(apx_sigmas);
    apx_image = U * apx_S * V';

    % Plot approximation
    subplot(4, 2, i + 1), imshow(apx_image), title(sprintf('Rank %d Image', ranks(i)));
end