% 2D imaging example: using the Truncated Singular Value Decomposition
% to find an approximate solution to an image of the moon

% Load the image of the moon and plot it.
clear all, close all
xtrue = imread('moon.tif');
xtrue = double(xtrue);
figure(1), imagesc(xtrue), axis image, colormap(gray), colorbar

% Create the first type of blur. 
% Here Ac corresponds to the column blur 
% and Ar corresponds to the row blurr. 
[m,n] = size(xtrue);
c = zeros(m,1);
c(1:5) = [5:-1:1]'/25;
Ac = toeplitz(c);
c = zeros(n,1);
c(1:10) = [5:-.5:.5]'/50;
Ar = toeplitz(c);

% Create blurred noise data and plot it.
noise=input('Noise level of RHS = ');
if isempty(noise)
  noise=.01
end

disp('... implementing column blur and row blur transformations and adding nose ...') 

b = Ac*xtrue*Ar' + noise*randn(m,n);
figure(2), imagesc(b), axis image, colormap(gray), colorbar


% Compute naive solution and plot it.
disp('... computing svd of both the column blur and row blur matrices ...') 
[Uc,Sc,Vc]=svd(Ac);
[Ur,Sr,Vr]=svd(Ar);
diagSc=diag(Sc);
diagSr=diag(Sr);

figure(3)
subplot(1,2,1)
plot(diagSc)
title('spectrum of the column blur')
subplot(1,2,2)
plot(diagSr)
title('spectrum of the row blur')

alpha = input(' Truncation level = ');

tic
Sc_trunc = find(diagSc>alpha);
Sr_trunc = find(diagSr>alpha);
xtsvd = ( ( Vc(:,Sc_trunc) * diag(1./diagSc(Sc_trunc)) * (Uc(:,Sc_trunc)'*b) ) * Vr(:,Sr_trunc) ) * diag(1./diagSr(Sr_trunc))  * Ur(:,Sr_trunc)';
toc

figure(4) 
  imagesc(xtsvd)
  axis image
  colormap(gray), colorbar



  figure(5)
  
  subplot(1,3,1)
  plot(Ur(:,2))
  title('2nd')
  
  subplot(1,3,2)
  title('Singular vectors of the row blur matrix')
  plot(Ur(:,10))
  title('10th')
  
  subplot(1,3,3)
  plot(Ur(:,end))
  title('Last')