%creates a noisy mxn pixel image of an elliptic band 
%uses diagonalization to create low-rank approximations of the image

close all
m = 150;
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(m,n);
for i=1:m    %set image = elliptic band of width ~w
    for j=1:n
        r = sqrt(1.5*(i-m/2)^2+(j-n/2)^2);
        if r > n/4 && r < n/4+w
            X(i,j) = 1;          %X(i,j) = 1 inside band, and 0 elsewhere
        end
    end
end

Y = X;
for i=1:m    %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

[U,Sig,V] = svd(Y);    %singular value decomposition
Z = cell(10,1);
Sig_dom = zeros(m,n);
for i=1:10
    Sig_dom(i,i) = Sig(i,i);         %D consists of i largest eigenvalues of Y
    Z{i} = U*Sig_dom*transpose(V);   %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
