%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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		= 1 : 11;
n_support	= 160;

n			= 1;

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

%% create synthetic checkerboard matrix
[xm,ym]		= ndgrid(cal.x_axis, cal.y_axis);
% create synthetic checkerboard
im_check_ref= xor(mod(xm, abs(cal.grid_dx)*2) >= abs(cal.grid_dx), ...
				  mod(ym, abs(cal.grid_dy)*2) >= abs(cal.grid_dy));

% 
% label regions on the checkerboard
[im_check_w,n_white]= bwlabel( im_check_ref,4);
[im_check_b,n_black]= bwlabel(~im_check_ref,4);
% create index for all tiles
im_check_b  = im_check_b + n_white;
im_check_b(im_check_ref) = 0;
im_check_ind= im_check_w+ im_check_b;
n_tiles		= n_white + n_black;
% average value over checkerboard tile
im_bw_val   = accumarray(im_check_ind(:), im_dewarp(:), [n_tiles 1], @mean);
im_sharp    = im_bw_val(im_check_ind);

%% estimate PSF 
im_sharp_ft = fft2(im_sharp);
im_deconv_ft= fft2(im_deconv);
im_psf_ft   = im_deconv_ft ./ im_sharp_ft;
im_psf_fg   = ifft2(im_psf_ft);
im_psf_fg   = im_psf_fg(1:n_support, 1:n_support);

%% deconvolve
fprintf('deconvolving...\n');
% PSF0		= zeros(n_support);
% i0			= ceil(n_support/2);
% j0			= i0;
% PSF0(i0,j0) = 1;
%PSF0		= fspecial('gaussian', n_support);
PSF0		= im_psf_fg;
[im_deconv, PSF] = deconvblind(im_dewarp, PSF0);

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

layout		= tiledlayout(1,2);
nexttile;
imagesc(cal.x_axis, cal.y_axis, im_dewarp.');
hold on;
plot(cal.x_axis([1 end]), cal.y_axis(j0)*[1 1], 'r--')

axis equal tight xy;


title('original');

nexttile
% plot line profile
plot(cal.x_axis, im_dewarp(:,j0));
hold on;
plot(cal.x_axis, im_sharp(:,j0), 'r--');

pbaspect([4/3 1 1]);

%% display source and deconvolved images
figure(3);
clf;
layout		= tiledlayout(2,2);
nexttile;
imagesc(cal.x_axis, cal.y_axis, im_dewarp.');
axis equal tight xy;
title('original');
caxis([0 2^15]);

nexttile
imagesc(cal.x_axis, cal.y_axis, im_sharp.');
axis equal tight xy;
title('sharp');
caxis([0 2^15]);

nexttile
imagesc(cal.x_axis, cal.y_axis, im_deconv.');
axis equal tight xy;
title('deconvolved');
%caxis([0 2^15]);

%% display PSF
figure(2);
clf;
imagesc(PSF.');
axis equal tight ij;
title('PSF');