%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% estimate PSF through blind deconvolution

cal_dir		= 'F:\mixer\lpt\calibration\21-06-30';
cal_pat		= '%s/CAM%d-%02d.mat';
cam_id		= 2;
cal_rng		= 11;%1 : 11;

x0			= [0;0]; % grid 
lambda		= 2;	% regularisation parameter
ksize		= 17;	% kernel size
n_grid		= 2;	% number of grid squares

%% load planar calibration
cal_file	= sprintf(cal_pat, cal_dir, cam_id, cal_rng);
fprintf('load %s\n', cal_file);
cal			= load(cal_file);

%% create a blurry image
PSF_ref		= fspecial('gaussian', ksize, ksize * 0.2);
[xm,ym]		= ndgrid(cal.x_axis, cal.y_axis);
im_sharp_ref= xor(mod(xm, abs(cal.grid_dx)*2) >= abs(cal.grid_dx), ...
				  mod(ym, abs(cal.grid_dy)*2) >= abs(cal.grid_dy));
im_full		= convn(im_sharp_ref, PSF_ref, 'same');
% create additive Gaussian noise
im_noise	= randn(size(xm)) * 0.01;
im_full     = im_full + im_noise;

%im_full		= cal.dewarped_image;

%% extract region surrounding grid point
dx			= cal.grid_dx;
dy			= cal.grid_dy;
x_axis		= cal.x_axis;
y_axis		= cal.y_axis;

x_ok		= abs((x_axis - x0(1))/dx) <= n_grid/2;
y_ok		= abs((y_axis - x0(2))/dy) <= n_grid/2;
im			= im_full(x_ok, y_ok);
x_trunc		= x_axis(x_ok);
y_trunc		= y_axis(y_ok);

%% estimate PSF
[PSF,im_sharp]=psf_estimate_local(x_trunc, y_trunc, im, dx, dy, ksize, lambda);
im_conv		= convn(im_sharp, PSF, 'same');
im_res		= im - im_conv;

%% display source images
figure(1);
clf;
j0			= 800;

layout		= tiledlayout(2,2);
nexttile;
imagesc(x_trunc, y_trunc, im.');
axis equal tight xy;
title('Reference');

nexttile
imagesc(x_trunc, y_trunc, im_sharp.');
axis equal tight xy;
title('Sharp');

nexttile
imagesc(x_trunc, y_trunc, im_conv.');
axis equal tight xy;
title('Convolved');

nexttile
imagesc(x_trunc, y_trunc, im_res.');
axis equal tight xy;
title('Residual');

%% compare PSF
figure(2);
clf;

layout		= tiledlayout(1,3);
nexttile;
imagesc(x_trunc, y_trunc, PSF_ref.');
axis equal tight;
title('Reference PSF');

nexttile;
imagesc(x_trunc, y_trunc, PSF.');
axis equal tight;
title('Recovered PSF');

nexttile;
% surface plot of PSF
S_ref = surf(PSF_ref);
hold on;
S     = surf(PSF);
set(S_ref, 'EdgeColor', 'r', 'FaceColor', 'none');
title('PSF Comparison');
daspect([1 1 0.5/ksize^2]);


