%MTH 410/510 Numerical Optimization & DA
%sample code for image blurring and Picard condition/plot

load clown; %load the image  
%to use your image file myimage.jpg uncomment the lines below
%X = imread('myimage.jpg');
%X = double(X);
%X = X(:,:,1);

%get the matrix size
[m,n] = size(X);

%define the blur operator: 1-dim diffusion applied columnwise
b = zeros(m,1);
L=0.45; %specify a value in (0,0.5], higher values increase the diffusion 
b(1) = 1-2*L; b(2) = L;
B = toeplitz(b); %tridiagonal symmetric matrix
B1 = B^50; %50 time steps higher values increase blurring and cond number 
%check condition number 
cond(B1)
%visualize the original image and the blurred image
figure()
subplot(1,2,1)
imagesc(X); colormap(gray)
subplot(1,2,2)
BX = B1*X;
imagesc(BX); colormap(gray)

%perform SVD of the blurring operator and check Picard condition
[U,S,V] = svd(B1);
figure()

hold on
for i=1:n
    plot(log10(abs(U'*BX(:,i))),'go')
    plot(log10(abs(U'*BX(:,i))./diag(S)),'bo')
end
plot(log10(diag(S)),'ro')
%generate a noise vector 
e = randn(m,1);
%plot the corresponding noise impact on solution  
hold on; plot(log10(abs(U'*e)./diag(S)))

