function [particles2d, img] = shad_project(camera, particles3d)
	% [particles2d, img] = shad_project(camera, particles3d)
	%
	% create synthetic shadowgraph images of rod- and sphereical particles
	% based on the camera model and particle coordinates specified
	%
	% Input:
	%	camera		a camera calibration structure as generated by shad_calib
	%	particles3d 1 x N structure array containing the following fields
	%		re_pos		3 x 1 object space position of particle
	%		re_radius	1 x 1 sphere radius (for spheres) or cylinder
	%					radius (for rods)
	%		re_length	3 x 1 particle length/orientation (for rods)
	%		

	%% sanity check
	assert(isstruct(particles3d) && isvector(particles3d), 'particles3d must be structure array');
	assert(isfield(particles3d, 're_pos'), 'particles3d must have re_pos field');
	assert(isfield(particles3d, 're_radius'), 'particles3d must have re_radius field');
	assert(isfield(particles3d, 're_length'), 'particles3d must have re_length field');
	n_particles		= numel(particles3d);

	%% define support
	n_support		= camera.otf.n_support;
	irel			= -floor(n_support/2) : floor(n_support/2);
	irel			= irel(:);
	jrel			= irel.';
	n_support		= [length(irel) length(jrel)];
	[irmat,jrmat]	= ndgrid(irel, jrel);
	n_support_px	= numel(irmat);

	i_axis			= camera.i_axis;
	j_axis			= camera.j_axis;
	n_i				= length(i_axis);
	n_j				= length(j_axis);
	
	%% otf
	re_source		= camera.otf.re_source;
	W_source		= camera.otf.W_source;
	n_source		= size(re_source, 2);
	re_source		= reshape(re_source, [3 1 n_source]);
	
	
	%% extract particle position and orientation into matrices
	particles3d		= reshape(particles3d, [1 n_particles]);
	re_pos			= cat(2, particles3d.re_pos);
	re_length		= cat(2, particles3d.re_length);
	re_radius		= cat(2, particles3d.re_radius);
	b_is_rod		= isfinite(re_length(1,:));

	%% project particle centers
	[im_pos, im_jac]= pinhole_model(re_pos, camera.calib, true);
	im_length		= pagemtimes(im_jac, reshape(re_length, [3 1 n_particles]));
	im_length		= squeeze(im_length);
	im_pos0			= round(im_pos);
    
    % TODO: implement calculation of radius in image space
    im_radius       = nan(1, n_particles);

	%% project particles locally
    particles2d     = cell(1, n_particles);
    img             = ones(n_i, n_j);
	for n = 1 : n_particles
		%% identify set of pixels near particle image
		% image space coordinates of nearby pixels
		im_near		= [irmat(:)+im_pos0(1,n), jrmat(:)+im_pos0(2,n)].';
		% corresponding pixels in image
		[i_ok,row]	= ismembertol(im_pos0(1,n)+irel, i_axis, 0.5, 'DataScale', 1);
		[j_ok,col]	= ismembertol(im_pos0(2,n)+jrel, j_axis, 0.5, 'DataScale', 1);
		% index and flag for mapping
		map         = (col-1)*n_i + row;
		map_ok      = i_ok & j_ok;
        map(~map_ok)= 0;
		
		%% create particle image
		if ~any(map_ok, 'all')
			im_part = ones([n_support_px 1]);
		else
			%% evaluate object space coordinates of focal plane
			% back project rays from image plane to focal plane
			re_focal_plane	= back_project_1c(camera.calib, camera.otf.M_focus, im_near, true);
	
        	%% check if ray is blocked
        	% form 
        	R           = particles3d(n).re_radius;
        	if strcmp(particles3d(n).type, 'sphere')
            	%% for spheres: evaluate closest distance from point to line
            	% n_source rays originate from
            	%	re_source (3 x 1 x n_source)
            	% in direction 
            	%   e_ray (3 x n_support_px x n_source)
            	% and are tested for their nearest distance to the particle at
            	%   re_pos(:,n)
            	e_ray		= re_focal_plane - re_source;
            	e_ray		= e_ray ./ sqrt(sum(e_ray.^2,1));
            	re_rel_pos  = re_source - re_pos(:,n);
				d_ray_sphere= sum(e_ray .* re_rel_pos, 1).^2;
				d_src_sphere= sum(re_rel_pos.^2, 1);
				r_ray_sphere= squeeze(sqrt(d_src_sphere - d_ray_sphere));
	
            	im_part     = sum(r_ray_sphere >= R, 2) / n_source;
			elseif strcmp(particles3d(n).type, 'rod')
            	%% for rods: evaluate closest distance between two lines
            	% we find the minimum distance between lines
            	%  xA = re_focal_plane(:,i) + A*e_ray(:,i)
            	%  xB = re_pos(:,n) + B*re_length(:,n)
            	% and identify where the closest point is along line B
            	% the shortest distance is along the direction normal to both
            	% line A and B,
            	e_ray		= re_focal_plane - re_source;
            	e_ray		= e_ray ./ sqrt(sum(e_ray.^2,1));
            	re_rel_pos  = re_focal_plane - re_pos(:,n);
            	e_rod       = repmat(re_length(:,n), [1 n_support_px n_source]);
            	
            	e_normal	= cross(e_ray, e_rod, 1);
            	e_normal    = e_normal ./ sqrt(sum(e_normal.^2,1));

				% check for points at ends of rod
				re_rel_end1 = re_source - (re_pos(:,n) + re_length(:,n)/2);
				re_rel_end2 = re_source - (re_pos(:,n) - re_length(:,n)/2);
				d_ray_end1  = sum(e_ray .* re_rel_end1, 1).^2;
				d_ray_end2  = sum(e_ray .* re_rel_end2, 1).^2;
				d_src_end1  = sum(re_rel_end1.^2, 1);
				d_src_end2  = sum(re_rel_end2.^2, 1);
				r_ray_end1  = squeeze(d_src_end1 - d_ray_end1);
				r_ray_end2  = squeeze(d_src_end2 - d_ray_end2);
				b_blocked   = r_ray_end1 < R^2 | r_ray_end2 < R^2;

            	% solve a set of equations
            	% xA0 + s*e_A + r*e_AB = xB0 + t*e_B
            	% xA0 - xB0 = - r*e_AB + t*e_B - s*e_A 
            	e_normal    = reshape(e_normal, [3 1 n_support_px n_source]);
            	e_ray       = reshape(e_ray   , [3 1 n_support_px n_source]);
            	e_rod       = reshape(e_rod   , [3 1 n_support_px n_source]);
            	RHS         = cat(2, -e_normal, e_rod, -e_ray);
            	LHS         = repmat(reshape(re_rel_pos, [3 1 n_support_px]), [1 1 1 n_source]);
            	%rts         = pagemldivide(RHS, LHS);
            	rts         = mmx('backslash', RHS, LHS);
            	d_ray_rod   = squeeze(rts(1,1,:,:));
            	d_along_rod = squeeze(rts(2,1,:,:));
            	b_blocked   = b_blocked ...
							| (abs(d_ray_rod) < R & abs(d_along_rod) < 1/2);

				
            	im_part     = sum(~b_blocked, 2) / n_source;
			else
				error('unrecognised particle type "%s"', particles3d(n).type);
			end

			%% accumulate multiplicatively in image
			img(map(map_ok)) = img(map(map_ok)) .* im_part(map_ok);
		end	

        %% add to structure
        p2d     = struct('map', map, ...
            'pimg', im_part, ...
            'im_pos', im_pos(:,n), ...
            'im_length', im_length(:,n), ...
            'im_radius', im_radius(n));
        particles2d{n} = p2d;

		%% debug
        continue;
		figure(1);
        clf;
        im_part     = reshape(im_part, n_support);
        imagesc(irel, jrel, im_part');
        %imagesc(i_axis, j_axis, img.');
        axis equal tight ij;
        colormap gray;
    end

    %% aggregate
    particles2d = cell2mat(particles2d);
end
