%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% create synthetic particle images 
addpath(genpath([pwd '/../../core']));

%% input parameters
% source data from DNS
src_dir			= 'L:/shadow/reference/AR8.00/';
src_pat         = '%s/%d.mat';
% scale and offset 
src_scale_length= 50 / 1.364;		% mm / code unit
%src_scale_time  = 25e-3 / 0.0424;	% seconds / code unit
src_idx			= 1;
src_subset		= 1 : 50000;		% subset of particles to view
src_offset		= [pi;pi;pi];		% offset in code units domain

dst_dir         = 'L:/shadow/images/AR8.00_3C/';
dst_pat			= '%s/S%04d_F%05d.cam%d.tif';
mat_pat         = '%s/S%04d_F%05d.mat';

calib_file		= [dst_dir '/setup.mat'];
n_particles     = 100;
particle_radius = 0.05; % mm
particle_length = 8*particle_radius; % mm
particle_type   = 'rod';

n_min_cams		= 3;

mkdir(dst_dir);

%% load setup
fprintf('loading camera setup from %s\n', calib_file);
setup			= importdata(calib_file);
camera			= setup.camera;
n_cameras		= setup.n_cameras;
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;

im_corners		= [i_lim([1 1 2 2]); j_lim([1 2 2 1])];

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);

%% report size of measurable volume
re_bb_lim		= setup.re_bb_lim;
fprintf('reconstruction limits are:\n');
fprintf('  x = %+6.1f - %+6.1f mm\n', re_bb_lim(:,1));
fprintf('  y = %+6.1f - %+6.1f mm\n', re_bb_lim(:,2));
fprintf('  z = %+6.1f - %+6.1f mm\n', re_bb_lim(:,3));

% extended limits for simulation domain
re_sim_lim		= diag([1.1 1.1 1.1])*re_bb_lim.';

%% load initial condition from DNS
ic_file			= sprintf('%s/ic.mat', src_dir);
ic				= load(ic_file);
n_t				= length(ic.t);
n_particles		= length(src_subset);

%% define viewing frustum of each camera
M_ary			= zeros(4, 4, n_cameras);
U_ary			= zeros(3, 4, n_cameras);
xc_ary			= zeros(3, n_cameras);
for ic = 1 : n_cameras
	[M_ary(:,:,ic),U_ary(:,:,ic),xc_ary(:,ic)] = frustum(camera(ic));
end

%% iterate over timesteps
t_start			= tic;
for n = 1 : n_t
	%% load particle distribution
	t_frame     = tic;
	ds_file		= sprintf(src_pat, src_dir, n);
	ds			= load(ds_file);
	
	%% convert to lab coordinate system
	re_pos		= (ds.x(:,src_subset) - src_offset) * src_scale_length;
	re_radius   = ones(1, n_particles) * particle_radius;
	re_length	= squeeze(ds.R(:,1,src_subset)) * particle_length;

	%% identify subset of particles within measurement volume
	% check for visibility in each camera's viewing frustum
	b_visible		= false(n_cameras, n_particles);
	for ic = 1 : n_cameras
		pndist		= M_ary(:,:,ic) * [re_pos; ones(1,n_particles)];
		ok			= all(pndist>=0,1);
		b_visible(ic,:)=ok;
	end
	% check also if within limits of depth of field
	b_in_limits		= all(re_pos >= re_sim_lim(:,1) & re_pos <= re_sim_lim(:,2), 1);
	b_visibleN		= sum(b_visible,1)>=n_min_cams & b_in_limits;
	n_visible		= nnz(b_visibleN);
	fprintf('%d particles visible in at least %d cameras\n', n_visible, n_min_cams);

	% select only visible particles
	prt_id		= find(b_visibleN);
	re_pos		= re_pos(:,b_visibleN);
	re_radius   = re_radius(:,b_visibleN);
	re_length   = re_length(:,b_visibleN);

	%% create a random distribution of rods
	particles3d = struct('re_pos', num2cell(re_pos, 1), ...
    	're_radius', num2cell(re_radius, 1), ...
		're_length', num2cell(re_length, 1), ...
		'id', num2cell(prt_id), ...
    	'type', particle_type);

	%% create synthetic image for each camera
	img_ary         = cell(1, n_cameras);
	p2d_ary         = cell(1, n_cameras);
	for c = 1 : n_cameras
    	fprintf('projecting %d particles in camera %d\n', n_visible, c);
    	[p2d_ary{c}, img_ary{c}] = shad_project(setup.camera(c), particles3d); 
	end
	
	%% save images and reference dataset
	dst_file  = sprintf(mat_pat, dst_dir, src_idx, n);
	fprintf('   save %s\n', dst_file);
	save(dst_file, 'particles3d');

	for c = 1 : n_cameras
		dst_file = sprintf(dst_pat, dst_dir, src_idx, n, c-1);
		fprintf('   save %s\n', dst_file);
		imwrite(img_ary{c}.', dst_file);
	end
	fprintf('frame %d took %0.0f seconds\n', n, toc(t_frame));
	fprintf('%2.0f h %2.0f m %2.0f s elapsed\n', hms(toc(t_start)))

	%% display projected images
	continue;
	for c = 1 : n_cameras
    	figure(100+c);
    	clf;
    	imagesc(i_axis, j_axis, img_ary{c}.');
		hold on;
		im_pos = [p2d_ary{c}.im_pos];
		plot(im_pos(1,:), im_pos(2,:), 'r.')

		axis equal tight xy;
		xlim(sort(i_axis([1 end])));
		ylim(sort(j_axis([1 end])));
    	colormap gray;
	end
end


%% plot particles
re_endpoints = cat(1, re_pos+re_length/2, re_pos-re_length/2);
figure(1);
clf;

% plot rods
plot3(re_endpoints([1 4], :), re_endpoints([2 5], :), re_endpoints([3 6], :), ...
    'r-', 'linewidth', 2);
hold on;

% plot viewing frustums
cmap			= hsv(n_cameras);
for ic = 1 : n_cameras
	xc			= camera(ic).calib.xc;
	plot3(xc(1), xc(2), xc(3), 'ko');
	re_corners	= back_project_1c(camera(ic).calib, [0 0 -1 re_bb_lim(2,3)], im_corners, true);
	P			= patch(re_corners(1,:), re_corners(2,:), re_corners(3,:), cmap(ic,:));
	set(P, 'FaceAlpha', 0.5, 'edgecolor', cmap(ic,:));
	hold on;
	for k = 1 : 4
		plot3([xc(1) re_corners(1,k)], [xc(2) re_corners(2,k)], [xc(3) re_corners(3,k)], '-', 'color', cmap(ic,:));
	end
	% point sources
	otf			= camera(ic).otf;
	plot3(otf.re_source(1,:), otf.re_source(2,:), otf.re_source(3,:), '.', 'color', cmap(ic,:));
end

axis equal tight xy;
axis(setup.re_bb_lim(:)*2);
camup([0 1 0]);
view(0,0);