%creates a noisy nxn pixel image of a circular band 
%uses SVD to create low-rank approximations of the image

close all
n = 100;     %size of image
w = 6;       %approximate width of circular band
p = .4;     %noise parameter

f = @(x) heaviside(x-n/4)*(1-heaviside(x-2*n/4));

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

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

[Q,D] = eig(Y);      %Y = QDQ^T = orthogonal diagonalization of Y
evalues = diag(D);   %the next few lines are for ordering the evalues of Y
[a,order] = sort(abs(diag(D)),'descend');  %in descending order 
Q = Q(:,order);      %and then reordering the eigenvectors appropriately
evalues = evalues(order);

Z = cell(10,1);
D_dom = zeros(n,n);
for i=1:10
    D_dom(i,i) = evalues(i);  %D_dom consists of i largest eigenvalues of Y
    Z{i} = Q*D_dom*transpose(Q);    %Z{i} is a "rank-i" approximation of Y
end

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

figure;
imagesc(Y)
axis square
title('noisy image')
colorbar

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