% Example computations related to two-dimensional deconvolution.
% This file commits inverse crime.
%
% Samuli Siltanen February 2010

% Noise level
sigma = .02; 

% Construct example signal, a two-dimensional image. It is just an NxN
% matrix with pixel values zero apart from a square, where the values are
% one.
N = 32;
x = zeros(N,N);
x(round(N/4):round(N/2),round(N/4):round(N/2)) = 1;

% Construction of the measurement matrix A. We utilize here the routine
% oblur.m by Professor Per Christian Hansen, a Danish expert on 
% computational inverse problems. Here the number PSFradius gives the
% radius of the PSF mask; for example, PSFradius=3 leads to 5x5 PSF mask,
% and PSFradius=5 leads to 11x11 mask.
PSFradius = 3;
A = oblur(N,PSFradius);

% Construct ideal and noisy measurements m and mn. The noise level sigma is
% defined above and gives the standard deviation of the normally
% distributed additive noise. Note the command x(:) that reshapes the
% NxN matrix x into a vertical vector of length N^2. Also, note that the
% command m=reshape(m,N,N) reshapes a vertical vector m of length N^2 into
% a matrix of size NxN.
m  = A*x(:);
m  = reshape(m,N,N);
mn = m+sigma*randn(size(m));

% Demonstration of unsuccessful naive inversion (with inverse crime)
% Reconstruct from ideal and noisy data
Ainv = inv(A);
rec  = Ainv*m(:);
rec  = reshape(rec,N,N);
recn = Ainv*mn(:);
recn = reshape(recn,N,N);


% Take a look at the signal and matrix A
figure(1)
clf
subplot(1,2,1)
imagesc(x)
set(gca,'xtick',[1 N/2 N])
set(gca,'ytick',[1 N/2 N])
colormap gray
title('Signal x')
axis square
subplot(1,2,2)
spy(A)
title('Nonzero elements of matrix A')

% Take a look at both measurements
figure(2)
clf
subplot(1,3,1)
imagesc(x)
colormap gray
axis off
title('Signal x')
axis square
subplot(1,3,2)
imagesc(m)
colormap gray
axis off
title('Ideal measurement')
axis square
subplot(1,3,3)
imagesc(mn)
colormap gray
axis off
title('Noisy measurement')
axis square

% Take a look at both reconstructions
figure(3)
clf
subplot(1,3,1)
imagesc(x)
colormap gray
axis off
title('Signal x')
axis square
subplot(1,3,2)
imagesc(rec)
colormap gray
axis off
title('Naive reconstruction (ideal data)')
axis square
subplot(1,3,3)
imagesc(recn)
colormap gray
axis off
title('Naive reconstruction (noisy data)')
axis square
