%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% synthetic image generation for RGB images from JHU turbulence database
% point clouds
function rgb_sig_spill(sig_file)
addpath(genpath([pwd '../core/']));
addpath(([pwd '/core/misc/']));

%if nargin < 1
sig_file='/Volumes/home/test/SIG_C4_f20_N4.0_interpolant.mat';
%end

%% input parameters
%clc;
load(sig_file);
fprintf('------------------------------------------------------------------\n');
fprintf('calibration:      %s\n', calib_file);
fprintf('src:              %s\n', src_dir);
fprintf('src pat:          %s\n', src_pat);
fprintf('src range (time): %d-%d\n', src_rng_t([1 end]));
fprintf('src range (z):    %d-%d\n', src_rng_z([1 end]));
fprintf('dst:              %s\n', dst_dir);
fprintf('dst pat:          %s\n', dst_pat);
fprintf('image noise:      %0.1f counts\n', im_noise_level);
fprintf('# particles:      <%d\n', n_subsample);
fprintf('------------------------------------------------------------------\n');

font_size		= 24;
%% load setup
fprintf('loading camera setup from %s\n', calib_file);
setup			= importdata(calib_file);
camera			= setup.camera;
n_cameras		= setup.n_cameras;
otf				= setup.otf;
delta_t			= setup.delta_t;

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)];
t				= (-1 : 1) * delta_t;

fprintf('# cameras:    %d\n', n_cameras);
fprintf('resolution:   %d x %d\n', n_i, n_j);
fprintf('delta t:      %0.3f ms\n', delta_t*1e3);

%% load particle volume
[iz,it]			= ndgrid(src_rng_z, src_rng_t);
n_src			= numel(iz);
for s = 1 : n_src
	%% load data
	src_name	= sprintf(src_pat, src_dir, iz(s), it(s));
	t_start		= tic;
	if exist(src_name, 'file') == 0
		warning('%s not found, skipping', src_name);
		continue;
	end
	fs			= load(src_name);
	fprintf('[z%03d,t%03d] %s\n', iz(s), it(s), src_name);
	fprintf('   z = %0.3f - %0.3f, t = %0.3f\n', fs.z_lim, fs.t);
	
	%% subsample data
	if fs.n_part > n_subsample
		warning('subsampling...\n');
		i_select	= 1 : n_subsample;% randperm(fs.n_part, n_subsample);
		fs.n_part	= n_subsample;
		fs.x		= fs.x(:, i_select);
		fs.u		= fs.u(:, i_select);
	end
	re_part			= fs.x;
	u_part			= fs.u;
	n_part			= fs.n_part;
	
	%% calculate particle displacement
	re_part_t0		= re_part - delta_t * u_part;
	re_part_t1		= re_part;
	re_part_t2		= re_part + delta_t * u_part;
	re_part_t		= {re_part_t0, re_part_t1, re_part_t2};
	I_part			= im_mean_bright * ones(1, n_part);
	
	%% sanity checking
	u_mean			= mean(fs.u,2);
	re_dx_mean		= u_mean * delta_t;
	fprintf('   dx mean: %0.3f %0.3f %0.3f mm\n', re_dx_mean);

	%% 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

        % this is probably the place to add colour spill
		img_rgb		= zeros(n_i, n_j, 3);
		for rgb = 1 : 3
            clean_img = zeros(n_i, n_j); % image without colour spill
			im_part		= pinhole_model(re_part_t{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 (%0.3f x 10^3 part/px)\n', n_visible, n_visible/(n_i*n_j)*1e3);
			% project only visible 
			[~,~,~,clean_img(:,:)] = rgb_project(camera(:,ic), otf(:,ic), ...
				re_part_t{rgb}(:,b_visible), ...
				I_part(b_visible), rgb);

            for spill_channel = 1:3
                img_rgb(:,:, spill_channel) = img_rgb(:,:, spill_channel) + (clean_img(:,:) * otf(spill_channel).spillage(rgb, spill_channel));
            end
            

		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
		% quantisation noise
		img_noisy		= uint8(img_noisy);

		%% store in array
		img_ary{ic}		= img;
		img_ary_rgb{ic} = img_rgb;
		img_ary_noisy{ic}=img_noisy;
		
		%% save images
		dst_name		= sprintf(dst_pat, dst_dir, iz(s), it(s), ic);
		fprintf('   save %s\n', dst_name);
		imwrite(img_noisy, dst_name);
	end
	
	%continue;
	%% demosaic images
	figure(1);
	clf;
	layout			= tiledlayout(1, 3);
	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));
end



