Toy program to compress greyscale images using a singular value decomposition

svd.m

clear;

RGB = imread('escher.jpg');
I = double(rgb2gray(RGB));

[m, n] = size(I);
[V, S, W] = svd(I);

% Achieves a 2.04 compression ratio while still being visible
slimit = mean(mean(S)) + 4 * std(std(S));

% The number of non-zero singular values is equal to the rank of I.
% This way we avoid calling rank(I) which is expensive.
maxRank = min(m, n);
r = maxRank;

for i=1:maxRank
    if (S(i, i) == 0)
        r = i;
        break
    end
end

% Compute the compressed matrix
R = zeros(m, n);
used = 0;

for i=1:r
    svalue = S(i, i);
    
    if svalue < slimit
        used = i;
        fprintf('Reached singular value s_%d = %f, l = %f, r = %d\n', used, svalue, slimit, r);
        break
    end
    
    R = R + S(i, i) * V(:,i) * W(:,i)';
end

originalCount = m * n;
valueCount = used * (m + n);
compRatio = originalCount / valueCount;
error = norm(I - R);

fprintf('compression ratio = %f (original requires %d, compressed has %d); error = %.2f\n', compRatio, originalCount, valueCount, error);

displayRange = [0 255];

f1 = figure(1);
f1.Name = 'Original';
imshow(I, displayRange);

f2 = figure(2);
f2.Name = 'Compressed';
imshow(R, displayRange);

Original:

Relativity by M. C. Escher

Compressed (4 * std):

error = 692.92, compression ratio = 2.038.

Compressed (12 * std):

error = 2054.70, compression ratio = 6.502.