function [im_part, pimg_ary, pimg_map, img, pimg_jac] = rgb_project(camera, otf, re_part, I_part, channel)
	% [im_part, pimg_ary, pimg_map, img, pimg_jac] = rgb_project(camera, otf, re_part, I_part, channel)
	% create synthetic reconstruction of the particle distribution
	% using 
	
	% output arguments
	b_img					= nargout > 3;
	b_jacobian				= nargout > 4;
	
	%% extract model parameters
	n_part					= size(re_part, 2);
	n_support				= otf.n_support;
	i_axis					= camera.i_axis;
	j_axis					= camera.j_axis;
	bayer					= camera.bayer;
	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)];
	
	%% convert to gpuArray
	%re_part					= gpuArray(re_part);
	%I_part					= gpuArray(I_part);
	
	%% geometric projection of particles using camera model
	% this gives us the image space position of the particle im_part
	% and the jacobian of this quantity with respect to object space
	% coordinate
	[im_part, im_jac_re]	= pinhole_model(re_part, camera.calib, true);
	
	%% obtain individual particle images
	% pimg_ary, pimg_row and pimg_col are n_support^2 x n_part arrays
	% representing the particle image in value,row,col format
	% pimg_jac is [n_support^2 2 n_part] is jacobian with respect to image
	% space coordinate
	
	if strcmp(otf(1).type, 'diffraction')
		[pimg_ary, pimg_jac_im, pimg_row, pimg_col]	= otf_evaluate_diffraction(camera, otf, re_part, im_part, channel);
	elseif strcmp(otf(1).type, 'diffraction_interpolant')
		[pimg_ary, pimg_jac_im, pimg_row, pimg_col]	= otf_evaluate_diffraction_interpolant(camera, otf, re_part, im_part, channel);
	elseif strcmp(otf(1).type, 'diffraction_bicubic')
		[pimg_ary, pimg_jac_im, pimg_row, pimg_col]	= otf_evaluate_diffraction_bicubic(camera, otf, re_part, im_part, channel);
	else
		error('unsupported model %s', otf.type);
	end	
	% define map with single array
	pimg_map			= cat(3, pimg_row, pimg_col);
	
	%% modulate image array by particle brightness
	pimg_ary			= pimg_ary .* reshape(I_part, [1 n_part]);
	pimg_jac_im			= pimg_jac_im .* reshape(I_part, [1 1 n_part]);
	
	%% accumulate sum of particle images
	% this is an optional output
	if b_img
		ok				= pimg_row >= 1 & pimg_row <= n_i ...
						& pimg_col >= 1 & pimg_col <= n_j;
		img				= accumarray({pimg_row(ok), pimg_col(ok)}, pimg_ary(ok), [n_i n_j]);
	end
	
	%% calculate jacobian
	% this is an optional output
	% jacobian of particle image with respect to object space position
	% this line does
	%	pimg_jac_re(:,:,p) = pimg_jac_im(:,:,p) * im_jac_re(:,:,p)
	% which gives the derivative of the OTF with respect to object space
	% coordinate
	if b_jacobian
		pimg_jac_re			= mtimesx(pimg_jac_im, im_jac_re);
		pimg_jac_re1		= reshape(pimg_jac_re(:,1,:), [n_support^2 n_part]);
		pimg_jac_re2		= reshape(pimg_jac_re(:,2,:), [n_support^2 n_part]);
		pimg_jac_re3		= reshape(pimg_jac_re(:,3,:), [n_support^2 n_part]);
		pimg_jac			= cat(3, pimg_jac_re1, pimg_jac_re2, pimg_jac_re3);
	end
end