%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% create a set of artificial particle images based on RGB LED illumination

addpath(genpath([pwd '/../core/']));
addpath(([pwd '/../../misc/']));

%% input parameters
% path to camera model
calib_dir		= './calib/';
calib_pat		= '%s/linear3.mat';

% displacement field
DeltaX			= 0.8; % ~ 0.1 mm/px
U_mean			= [DeltaX; 0; 0]; % displacement, mm

% volume limits for particles
re_x_lim		= [-200 200+134*2];
re_y_lim		= [-150 150];
re_z_lim		= [-25 25];
% seeding concentration, particles / mm^3
rho_p			= 0.001 * 12e6 / (300*400*re_z_lim*[-1;1]);
% mean noise level, counts
im_noise_level	= 2;			

%% load setup
cal_file		= sprintf(calib_pat, calib_dir);
fprintf('loading camera setup from %s\n', cal_file);
setup			= importdata(cal_file);
camera			= setup.camera;
n_cameras		= setup.n_cameras;
otf				= setup.otf;

i_axis			= camera(1).i_axis;
j_axis			= camera(1).j_axis;
di				= i_axis(2)-i_axis(1);
dj				= j_axis(2)-j_axis(1);
n_i				= length(i_axis);
n_j				= length(j_axis);
i_lim			= [min(i_axis) max(i_axis)];
j_lim			= [min(j_axis) max(j_axis)];

fprintf('# cameras:    %d\n', n_cameras);
fprintf('resolution:   %d x %d\n', n_i, n_j);

%% create particle volume
fprintf('creating synthetic particle volume:\n');
re_lim			= [re_x_lim; re_y_lim; re_z_lim];
FOV				= re_lim * [-1; 1];
V_meas			= prod(FOV);
n_part			= round(rho_p * V_meas);
re_part_t1		= diag(FOV)*rand(3,n_part) + re_lim(:,1);
re_part_t0		= re_part_t1 - U_mean;
re_part_t2		= re_part_t1 + U_mean;
re_part			= {re_part_t0, re_part_t1, re_part_t2};
% particle image intensity
I_part			= 255 * ones(1, n_part);
fprintf('# particles:   %d\n', n_part);

%% create particle images
fprintf('creating synthetic particle images ... \n');
img_ary			= cell(1, n_cameras);
img_ary_rgb		= cell(1, n_cameras);
img_ary_noisy	= cell(1, n_cameras);

for ic = 1 : n_cameras
	%% check if particle is visible in this camera
	fprintf('camera %d:\n', ic);
	
	%% accumulate images over three exposures
	img_rgb		= zeros(n_i, n_j, 3);
	for rgb = 1 : 3
		im_part		= pinhole_model(re_part{rgb}, camera(ic).calib, true);
		b_visible	= im_part(1,:)>=i_lim(1) & im_part(1,:) <= i_lim(2) ...
					& im_part(2,:)>=j_lim(1) & im_part(2,:) <= j_lim(2);
		n_visible	= nnz(b_visible);
		fprintf('   %d particles visible t=%d\n', n_visible, rgb);
		% project only visible 
		[~,~,~,img_rgb(:,:,rgb)] = rgb_project(camera(ic), otf(ic), ...
			re_part{rgb}(:,b_visible), ...
			I_part(b_visible), rgb);
	end
	img				= sum(img_rgb, 3);
	
	%% add noise
	if im_noise_level > 0
		fprintf('   adding noise...\n');
		noise		= poissrnd_fast(im_noise_level, n_i, n_j);
		img_noisy	= img + noise;
	else
		img_noisy	= img;
	end
	
	%% store in array
	img_ary{ic}		= img;
	img_ary_rgb{ic} = img_rgb;
	img_ary_noisy{ic}=img_noisy;
end

%% demosaic images
figure(1);
clf;
layout			= tiledlayout(2,2);
ax = [];

for ic = 1 : n_cameras
	nexttile;
	img				= img_ary{ic};
	
	img_demosaic	= demosaic(uint8(img), 'bggr');
	imshow(permute(img_demosaic, [2 1 3]));
	axis equal tight;
	title(sprintf('Camera %d', ic));
end

nexttile;
imagesc(img.');
axis equal tight;
title(sprintf('Camera %d raw', ic));

%% plot particle distribution
figure(10);
clf;
plot3(re_part_t0(1,:), re_part_t0(2,:), re_part_t0(3,:), 'r.');
hold on;
plot3(re_part_t1(1,:), re_part_t1(2,:), re_part_t1(3,:), 'g.');
plot3(re_part_t2(1,:), re_part_t2(2,:), re_part_t2(3,:), 'b.');

for ic = 1 : n_cameras
	xc=camera(ic).calib.xc;
	plot3(xc(1), xc(2), xc(3), 'ko');
end

axis equal tight xy;
%xlim(re_x_lim);
%ylim(re_y_lim);
%zlim(re_z_lim);
camup([0 1 0]);

xlabel('$x$ [mm]', 'interpreter', 'latex', 'fontsize', font_size);
ylabel('$y$ [mm]', 'interpreter', 'latex', 'fontsize', font_size);
zlabel('$z$ [mm]', 'interpreter', 'latex', 'fontsize', font_size);
set(gca,'fontsize',font_size,'TickLabelInterpreter','latex');




