% Two-dimensional deconvolution based on Tikhonov regularization.
% The routine DB1_data_comp.m should be computed before this file.
%
% Jennifer Mueller, Samuli Siltanen and Sanna Tyrvainen, October 2012
% Load precomputed results
load img_data psf noiselevel X mn p0
% Regularization parameter for Tikhonov regularization
delta = .009;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Demonstration of reconstruction using Tikhonov regularization and
% conjugate gradient method for the solution of the normal equations
%
% (A^T A + delta*I)x = A^T mn.
%
% The positive constant delta is the regularization parameter.
% We use Jacobi preconditioning, or replace matrix (A^T A + delta*I) by
%
% H := inv(P)*(A^T A + delta*I)
%
% and vector (A^T mn) by b = inv(P)*(A^T mn). Here P is a diagonal matrix
% of the simple form const*I, where const is the diagonal value of the matrix
% (A^T A + delta*I), which is p0^2+delta by construction.
const = p0^2+delta;
b = (1/const)*conv2(mn,psf,'same');
% Solve the minimization problem using conjugate gradient method.
% See Kelley: "Iterative Methods for Optimization", SIAM 1999, page 7.
K = 40; % maximum number of iterations
recn = zeros(size(X)); % initial iterate is the backprojected data
rho = zeros(K,1); % initialize parameters
% Compute residual
Arecn = conv2(recn,psf,'same');
ATArecn = conv2(Arecn,psf,'same');
Hrecn = (1/const)*(ATArecn + delta*recn);
r = b-Hrecn;
rho(1) = r(:).'*r(:);
% Start iteration
for kkk = 1:(K-1)
if kkk==1
p = r;
else
beta = rho(kkk)/rho(kkk-1);
p = r + beta*p;
end
Ap = conv2(p,psf,'same');
ATAp = conv2(Ap,psf,'same');
Hp = (1/const)*(ATAp + delta*p);
w = Hp;
a = rho(kkk)/(p(:).'*w(:));
recn = recn + a*p;
r = r - a*w;
rho(kkk+1) = r(:).'*r(:);
if mod(kkk,20)==0
disp([kkk K])
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Take a look at the reconstruction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2)
clf
subplot(1,3,1)
imagesc(X)%,[0 1])
colormap gray
title('Original image')
axis square
axis off
subplot(1,3,2)
imagesc(mn)
colormap gray
title('Blurry measurement')
axis square
axis off
subplot(1,3,3)
imagesc(recn)
colormap gray
title(['Reconstruction'])
axis square
axis off