% Kronecker blur.
% TSVD deblur

% 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 = 5;%input(' The standard deviation of the noise = ');
b = Ac*xtrue*Ar' + noise*randn(m,n);
figure(2), imagesc(b), axis image, colormap(gray), colorbar

% Compute naive solution and plot it.
[Uc,Sc,Vc]=svd(Ac);
[Ur,Sr,Vr]=svd(Ar);
alpha = input(' Truncation level = ');
Sc_trunc = (Sc>alpha).*Sc;
Sr_trunc = (Sr>alpha).*Sr;
xtsvd = ( ( Vc * ( pinv(Sc_trunc) * (Uc'*b) ) * Vr ) * pinv(Sr_trunc) ) * Ur';

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

D = full(diag(kron(sparse(Sc),sparse(Sr))));
figure(4) 
  semilogy(sort(D,1,'descend'))
  hold on, semilogy(alpha*ones(size(D)),'r'), hold off
  title('Singular Values')