% Example computations related to two-dimensional deconvolution.
% This file attempts naive inversion of noisy data, avoiding inverse crime
% at least to some extent. Namely, the measurement data is simulated at a
% double resolution and then downsampled back to original resolution.
% Properly avoiding inverse crime would involve using grids that are not
% subgrids of each other at all, but this approach will do for now.
%
% Samuli Siltanen February 2010

% Noise level
sigma = .02; 

% Construct the same simple example signal than was used in the routine
% conv2D_01_iversecrime.m: a two-dimensional image with a white square 
% on a black background. This time we need to construct the image on two 
% resolutions: the coarse resolution image x is an NxN matrix, and
% the fine resolution image xf is a (2N)x(2N) matrix. 
% The construction of images is done in a bit involved way:
% first, create a set of x- and y-coordinates of the centerpoints of pixels
% in the two images of different resolutions. We think of both images being
% defined on the unit square 0<x<1 and 0<y<1. After defining the
% coordinates of centerpoints of pixels at both resolutions, the white
% square can be simply defined as inequalities concerning the coordinates.
% Since both coordinate systems are defined inside the unit square, the
% images are comparable, as shown in Figure 1.

N       = 16; % It is best to take N to be a power of 2. Try 16 and 32.
dx      = 1/N;
t       = dx/2+[0:(N-1)]*dx;
[t1,t2] = meshgrid(t);
x       = zeros(size(t1));
ind     = (t1>1/4) & (t1<1/2) & (t2>1/4) & (t2<1/2);
x(ind)  = 1;

Nf        = 2*N;
dxf       = 1/Nf;
tf        = dxf/2+[0:(Nf-1)]*dxf;
[t1f,t2f] = meshgrid(tf);
xf        = zeros(size(t1f));
indf      = (t1f>1/4) & (t1f<1/2) & (t2f>1/4) & (t2f<1/2);
xf(indf)  = 1;


% Construction of the measurement matrix A at both resolutions. 
% We utilize here the routine oblur.m by Professor Per Christian Hansen, 
% a renowned 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. In our process of avoiding inverse crime, we use roughly
% twice the size of PSF on the fine grid. 
PSFradius = 3;
%A = 1/(N^2)*oblur(N,PSFradius);
A = oblur(N,PSFradius);
PSFradius_f = 2*PSFradius-1;
%Af = 1/(Nf^2)*oblur(2*N,PSFradius_f);
Af = oblur(2*N,PSFradius_f);

%figure(6);
%imagesc(A);
%colormap gray;

% Construct ideal measurements o both resolutions. 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);
mf = Af*xf(:);
mf = reshape(mf,Nf,Nf);

% Downsample the fine-resolution measurement mf (of size (2N)x(2N))
% to the coarse resolution and call the result m2 (size NxN). 
% Downsampling here means taking the average of four neighboring pixels
% while reducing resolution.
m2 = mf(:);
m2 = reshape(m2,2,length(m2)/2);
m2 = sum(m2);
m2 = reshape(m2,N,2*N);
m2 = m2.';
m2 = m2(:);
m2 = reshape(m2,2,length(m2)/2);
m2 = sum(m2);
m2 = reshape(m2,N,N);
m2 = m2.';
m2 = m2/4;

% Calculate the relative sup error percentage between m and m2
relerr = norm(m2-m)/norm(m);
relerr = round(relerr*100);

% Construct noisy measurement m2n. The noise level sigma is
% defined above and gives the standard deviation of the normally
% distributed additive noise.
m2n = m2+sigma*randn(size(m));

% Demonstration of unsuccessful naive inversion (without inverse crime)
% Reconstruct from ideal and noisy data
Ainv = inv(A);
rec  = Ainv*m2(:);
rec  = reshape(rec,N,N);
if isempty(find(isnan(rec)));
    relerr_ideal = norm(x-rec)/norm(x);
else 
    relerr_ideal = NaN;
end
relerr_ideal = round(relerr_ideal*100);
recn = Ainv*m2n(:);
recn = reshape(recn,N,N);
if isempty(find(isnan(recn)));
    relerr_noisy = norm(x-recn)/norm(x);
else
    relerr_noisy = NaN;
end
relerr_noisy = round(relerr_noisy*100);



% Take a look at the signals at both resolutions
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 on coarse resolution')
axis square
set(gca,'xtick',[1 size(x,2)/2 size(x,2)])
set(gca,'ytick',[1 size(x,1)/2 size(x,1)])
subplot(1,2,2)
imagesc(xf)
set(gca,'xtick',[1 N/2 N])
set(gca,'ytick',[1 N/2 N])
colormap gray
title('Signal x on fine resolution')
axis square
set(gca,'xtick',[1 size(xf,2)/2 size(xf,2)])
set(gca,'ytick',[1 size(xf,1)/2 size(xf,1)])

% Take a look at the matrices A and Af
figure(2)
clf
subplot(1,2,1)
spy(A)
title('Nonzero elements of matrix A')
set(gca,'xtick',[1 size(A,2)/2 size(A,2)])
set(gca,'ytick',[1 size(A,1)/2 size(A,1)])
subplot(1,2,2)
spy(Af)
title('Nonzero elements of matrix Af')
set(gca,'xtick',[1 size(Af,2)/2 size(Af,2)])
set(gca,'ytick',[1 size(Af,1)/2 size(Af,1)])


% Take a look at the blurred signals at both resolutions
figure(3)
clf
subplot(1,2,1)
imagesc(m)
set(gca,'xtick',[1 N/2 N])
set(gca,'ytick',[1 N/2 N])
colormap gray
title('Blurred signal m on coarse resolution')
axis square
set(gca,'xtick',[1 size(x,2)/2 size(x,2)])
set(gca,'ytick',[1 size(x,1)/2 size(x,1)])
subplot(1,2,2)
imagesc(mf)
set(gca,'xtick',[1 N/2 N])
set(gca,'ytick',[1 N/2 N])
colormap gray
title('Blurred signal mf on fine resolution')
axis square
set(gca,'xtick',[1 size(xf,2)/2 size(xf,2)])
set(gca,'ytick',[1 size(xf,1)/2 size(xf,1)])


% Compare blurred signals m and m2 at coarse resolution
figure(4)
clf
imagesc([m,m2,abs(m-m2)])
colormap gray
title(['Blurred image m, blurred image m2, and their difference. Relative error',num2str(relerr),'%'])
axis equal
axis off



% Take a look at both reconstructions
figure(5)
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), error ',num2str(relerr_ideal),'%'])
axis square
subplot(1,3,3)
imagesc(recn)
colormap gray
axis off
title(['Naive reconstruction (noisy data), error ',num2str(relerr_noisy),'%'])
axis square


%%%%%%%%%
% SVD reconstruction / Lauri

[U,D,V] = svds(A,100); % Compute the first 100 singular values
Ns = 30;         % Number of singular values used
[row,col]        = size(D.');
Dplus            = sparse(row,col);
svals            = diag(D);
Dplus(1:Ns,1:Ns) = diag(1./svals(1:Ns));
recsvd              = V*Dplus*U.'*m2n(:);
recsvd             = reshape(recsvd,N,N);

relerr_svd = norm(x-recsvd)/norm(x)

figure(6);
plot(svals);

figure(7);
imagesc(recsvd);
axis square
colormap gray;

figure;
plot(V(:,1));
figure;
plot(V(:,2));
figure;
plot(V(:,3));
figure;
plot(V(:,4));
figure;
plot(V(:,5));
