%% Setup

% Deconvolution will normally require that your data be in double format.
I_coins = double(imread('coins.png')); 
I_coins = I_coins./255;

%% Blurring / Convolution

% Make a convolution kernel that we know about. Sometimes you can derive
% this through some specific measurements using your imaging system.
% Different imaging conditions can lead to different sorts of point spread
% functions.
blur_kernel_gauss = fspecial('gaussian',[7,7],2);
blur_kernel_aniso = fspecial('motion',9,30); % Kernel

% Display the kernels
figure; colormap gray;
tiledlayout(1,2);
nexttile;
imagesc(blur_kernel_gauss); title('Gaussian blur'); axis equal tight;
nexttile;
imagesc(blur_kernel_aniso); title('Anisotropic motion blur'); axis equal tight;

% Apply.
coins_blurred_gauss = imfilter(I_coins,blur_kernel_gauss);
coins_blurred_aniso = imfilter(I_coins,blur_kernel_aniso);

% Display the filtered/blurred images.
figure; colormap gray;
imagesc(I_coins); title('Original image'); axis equal tight off;
figure; colormap gray;
imagesc(coins_blurred_gauss); title('Gaussian blurred'); axis equal tight off;
figure; colormap gray;
imagesc(coins_blurred_aniso); title('Anisotropic motion blurred'); axis equal tight off;

% Add noise for later.
n_mean = 0;
n_var = 1e-3; 
cbg_noise = imnoise(coins_blurred_gauss,'gaussian',n_mean,n_var);
cba_noise = imnoise(coins_blurred_aniso,'gaussian',n_mean,n_var);

figure; imagesc(cbg_noise); title('Gaussian blur + noise'); axis equal tight off; colormap gray;
figure; imagesc(cba_noise); title('Motion blur + noise'); axis equal tight off; colormap gray;


%
%% deconvwnr - Wiener deconvolution

% Wiener deconvolution is similar to Wiener filtering, and benefits from an
% estimate of the noise:signal ratio (NSR). Here, we already know that the
% noise variance (n_var in the section above) is 0.001 .

n_var_est = 0.001; % Vary this to see the effects of inaccurate noise estimates.
nsr_cbg = n_var_est ./ var(cbg_noise(:));
nsr_cba = n_var_est ./ var(cba_noise(:));

% Deconvolve each of the blurred images using the matched PSF...
coins_unblurred_g_w_g = deconvwnr(cbg_noise,blur_kernel_gauss,nsr_cbg);
coins_unblurred_a_w_a = deconvwnr(cba_noise,blur_kernel_aniso,nsr_cba);

figure; colormap gray;
tiledlayout(2,2)
nexttile; imagesc(cbg_noise); title('Gaussian blur + noise'); axis equal tight off; 
nexttile; imagesc(cba_noise); title('Motion blur + noise'); axis equal tight off; 
nexttile; imagesc(coins_unblurred_g_w_g); title('Matched Gaussian deblur'); axis equal tight off; 
nexttile; imagesc(coins_unblurred_a_w_a); title('Matched motion deblur'); axis equal tight off; 

% ... and the wrong PSF.
coins_unblurred_g_w_a = deconvwnr(cbg_noise,blur_kernel_aniso,nsr_cbg);
coins_unblurred_a_w_g = deconvwnr(cba_noise,blur_kernel_gauss,nsr_cba);

figure; colormap gray;
tiledlayout(2,2)
nexttile; imagesc(cbg_noise); title('Gaussian blur + noise'); axis equal tight off; 
nexttile; imagesc(cba_noise); title('Motion blur + noise'); axis equal tight off; 
nexttile; imagesc(coins_unblurred_g_w_a); title('Unmatched Gaussian deblur'); axis equal tight off; 
nexttile; imagesc(coins_unblurred_a_w_g); title('Unmatched motion deblur'); axis equal tight off; 



%
%% deconvreg - Regularized deconvolution

% You can use regularized methods to reduce the deconvolution process's
% sensitivity to noise.

n_var_est = 0.001;
n_pwr_est = n_var_est .* numel(cbg_noise);

% deconvreg can take either a noise power estimate or not.
coins_unblurred_reg_gaussian = deconvreg(cbg_noise,blur_kernel_gauss,n_pwr_est);
coins_unblurred_reg_gaussian_no_est = deconvreg(cbg_noise,blur_kernel_gauss);

figure; colormap gray;
imagesc(cbg_noise); title('Gaussian + noise'); axis equal tight off; 

figure; colormap gray;
imagesc(coins_unblurred_reg_gaussian); title('deconvreg w/ est'); axis equal tight off; 

figure; colormap gray;
imagesc(coins_unblurred_reg_gaussian_no_est); title('deconvreg no est'); axis equal tight off; 


%
%% deconvblind - Blind deconvolution

% If you know nothing about your image's point spread function, you can use
% blind deconvolution to attempt to recover both the deblurred image and an
% estimate of the point spread function.

% Blind methods are particularly ill-posed, and suffer greatly due to
% noise. Here we'll use the less noisy image. Notice we're using the
% anisotropic blurred image, but the Gaussian blur kernel/point spread
% function. A Gaussian is usually a good first estimate of the PSF, but you
% should always try to have a good estimate to start with.
I_blur = coins_blurred_aniso;

[I_deblur,psf_recover] = deconvblind(I_blur,blur_kernel_gauss,100);

figure;
tiledlayout(1,3);
nexttile; imagesc(I_coins); axis equal tight off; title('Original');
nexttile; imagesc(I_blur); axis equal tight off; title('Blurred');
% Notice that the recovered image has high-value outliers where noise is
% amplified.
nexttile; imagesc(I_deblur); axis equal tight off; title('Recovered image');


figure;
tiledlayout(1,3);
nexttile; imagesc(blur_kernel_aniso); axis equal tight off; title('True PSF');
nexttile; imagesc(blur_kernel_gauss); axis equal tight off; title('Estimated PSF');
nexttile; imagesc(psf_recover); axis equal tight off; title('Recovered PSF');

% You can also use the estimate of the PSF as a seed for another run of the
% blind deconvolution algorithm. 




