%creates a noisy nxn pixel image of a circular band 
%uses SVD to create low-rank approximations of the image
%SVD is a generalization of diagonalization
%copyright D. Aristoff :)

close all
n = 100;     %size of image
w = 6;      %approximate width of circular band
p = 0.5;     %noise parameter
X = zeros(n,n);

for i=1:n    %set image = circular band of width w
    for j=1:n
        r = sqrt((i-n/2)^2+(j-n/2)^2);
        if r > n/4 && r < n/4+w
            X(i,j) = 1;
        end
    end
end

for i=1:n    %introduce noise into image
    for j=1:n
        if rand() < p
            X(i,j) = mod(X(i,j)+rand(),1);   %noisy value between 0 and 1
        end
    end
end

[U,S,V] = svd(X);    %X = USV is the singular value decomposition of X
sigmas = diag(S);    %vector of singular values, ordered greatest to least

D = zeros(n,n);
Y = cell(4,1);
for i=1:10
    D(i,i) = sigmas(i);       %D consists of i largest singular values
    Y{i} = U*D*transpose(V);  %Y{i} is a rank-i approximation of X
end


imagesc(X)
axis square
title('original image')
colorbar

in = input('process image? ','s');
if in == 'yes'
figure;
subplot(2,3,1)
imagesc(X)
axis square
title('original image')
subplot(2,3,2)
imagesc(Y{1})
axis square
title('rank 1 approximation')
subplot(2,3,3)
imagesc(Y{2})
axis square
title('rank 2 approximation')
subplot(2,3,4)
imagesc(Y{3})
axis square
title('rank 3 approximation')
subplot(2,3,5)
imagesc(Y{4})
axis square
title('rank 4 approximation')
subplot(2,3,6)
imagesc(Y{5})
axis square
title('rank 5 approximation')
end