function [pimg_ary_cpu, pimg_jac_cpu, pimg_row_cpu, pimg_col_cpu] = otf_evaluate_diffraction_gpu(camera, otf, re_pos, im_pos, channel)
	% render points ij_pos in image space acccording to diffraction limited imaging
	% we also compute the jacobian of the particle image
	% with respect to the image space coordinates

	% OTF structure contains
	%  resolution  resolution for integrating profile [pixels]
	%  Z_0         sensor plane position[mm]
	%  f           focal length         [mm]
	%  D           aperture diameter    [mm]
	%  pitch       pixel ptich          [um]
	%  lambda      wavelength           [um]
	%  phi         dimensionless axial distance
	%  psi         dimensionless radial coordinate
	%  I           diffraction pattern I(phi,psi)
	%  gradI       gradient of I with respect to [phi,psi]


	%% sanity check
	n_part					= size(im_pos,2);
	assert(size(im_pos,1) == 2);
	
	% default colour channel to evaluate is first one (red)
	if nargin < 5
		channel					= 1;
	end
	
	%% convert to gpuArray if not already
	re_pos					= gpuArray(re_pos);
	im_pos					= gpuArray(im_pos);
	otf.phi					= gpuArray(otf.phi);
	otf.psi					= gpuArray(otf.psi);
	otf.I					= gpuArray(otf.I);

	%% support
	support						= otf.support;
	S							= size(support);
	n_support					= numel(support);
	assert(S(1)==S(2), 'Can only handle square support right now');
	b_fancy_support				= any(support ~= support(1)); % any point in multiplying?

	%% extract details of image grid
	i_axis						= camera.i_axis;
	j_axis						= camera.j_axis;
	delta_i						= i_axis(2) - i_axis(1);
	delta_j						= j_axis(2) - j_axis(1);
	n_i							= length(i_axis);
	n_j							= length(j_axis);
	Nint						= round(1 / otf.resolution);
	bayer						= camera.bayer(:,:,channel);
	assert(size(bayer,1) == 2 && size(bayer,2) == 2, 'bayer pattern was unexpected size');
	res							= 1 / Nint;
	assert(abs(delta_i)==1);
	assert(abs(delta_j)==1);

	%% model parameters
	f                       = otf.f;				% focal length
	D                       = otf.D;				% aperture
	wavelength              = otf.lambda(channel);	% wavelength
	Z_0                     = otf.Z_0;				% sensor plane position
	pitch                   = otf.pitch;			% pixel pitch

	%% evaluate focal distance for each particle
	P                       = camera.calib.P;
	z                       = abs(P(3,1:3) * re_pos + P(3,4));
	Z                       = (f*z)./(z-f);% thin lens equation, map to image plane
	phi_part                = 2*pi*(Z_0 - Z)/wavelength*(D/(2*f))^2;
                        
	%% define grid for integration of diffraction limited profile	
	% define the grid of neighbouring pixels
	irel					= (0 : S(1)-1) - floor(S(1)/2);
	jrel					= (0 : S(2)-1) - floor(S(2)/2);	
	irel					= gpuArray(irel);
	jrel					= gpuArray(jrel);
	[imat,jmat]				= ndgrid(irel, jrel);
	% and a fine-grained grid for integration
	irel_fine				= linspace(irel(1)-0.5+res/2, irel(end)+0.5-res/2, S(1)*Nint);
	jrel_fine				= linspace(jrel(1)-0.5+res/2, jrel(end)+0.5-res/2, S(2)*Nint);
	[imat_fine,jmat_fine]	= ndgrid(irel_fine, jrel_fine);
	n_fine					= numel(imat_fine);
	% for the pixel n with coordinate (i_n,j_n), we integrate over the
	% region 
	%	i_n - 0.5 < i < i_n + 0.5
	%  j_n - 0.5 < i < j_n + 0.5
	% or more succinctly where:
	%  |i - i_n| < 0.5 and |j -j_n| < 0.5
	%
	integration_kernel		= (abs(bsxfun(@minus, imat(:), imat_fine(:)'))<=0.5) ...
							.*(abs(bsxfun(@minus, jmat(:), jmat_fine(:)'))<=0.5);
	% normalise kernel
	% this ensures when we integrate over a unit field
	% we get a unit response for each pixel
	integration_kernel		= bsxfun(@rdivide, integration_kernel, sum(integration_kernel,2));
	% take advantage of sparsity structure
	integration_kernel		= sparse(integration_kernel);
	
	%% break down into smaller chunks to fit into GPU memory
	type					= 'double';
	pimg_ary_cpu			= zeros([n_support n_part], type);
	pimg_jac_cpu			= zeros([n_support 2 n_part], type);
	pimg_row_cpu			= zeros([n_support n_part], type);
	pimg_col_cpu			= zeros([n_support n_part], type);
	% TODO: automatically calculate this based on available GPU memory
	n_part_max				= 1e4;
	%disp(gpuDevice(1))
	
	for ip0 = 1 : n_part_max : n_part
		%% range of particles
		ip					= ip0 : min(n_part, ip0+n_part_max);
		n_part_chunk		= length(ip);
		
		%% evaluate OTF on image plane
		% this gives us a particle image measured in 2D image coordinates
		% see Lab Book 68, page 86-89 for the maths behind this 
		ij0_pos					= round(im_pos(:,ip));
		igrid_part				= bsxfun(@plus, ij0_pos(1,:), imat(:));
		jgrid_part				= bsxfun(@plus, ij0_pos(2,:), jmat(:));
		% y^ = (y - Y) / sigma
		%    = (y - Y0 - (Y-Y0) / sigma
		dij_pos					= im_pos(:,ip) - ij0_pos;
		ifine_part				= bsxfun(@minus, imat_fine(:), dij_pos(1,:));
		jfine_part				= bsxfun(@minus, jmat_fine(:), dij_pos(2,:));
		rfine_part              = sqrt(ifine_part.^2 + jfine_part.^2);
		psi_fine_part           = pi * (pitch/wavelength) * (D/f) * rfine_part;
		phi_fine_part           = repmat(phi_part(1,ip), [n_fine 1]);

		% evaluate particle image and its jacobian on finely spaced grid
		% result is of size n_fine x n_part
		pfine_part              = interpn(otf.phi, otf.psi, otf.I, phi_fine_part, psi_fine_part, 'linear', 0);
		pfine_part_jac_dr       = pi * (pitch/wavelength) * (D/f) * ...
								 interpn(otf.phi, otf.phi, otf.gradI(:,:,2), phi_fine_part, psi_fine_part, 'linear', 0);   
		% jacobian with respect to image space coordinates 
		pfine_part_jac_dx		= -pfine_part_jac_dr .* ifine_part ./ rfine_part;
		pfine_part_jac_dy		= -pfine_part_jac_dr .* jfine_part ./ rfine_part;
		% integrate image and its jacobian
		pimg_ary				= integration_kernel * pfine_part;
		pimg_jac_x				= integration_kernel * pfine_part_jac_dx;
		pimg_jac_y				= integration_kernel * pfine_part_jac_dy;

		%% create mapping function to image
		% determine the correspondence between values in otf_ary
		% and the row and column indices of pixels in the image volume
		% ii = row index
		% jj = column index
		pimg_row						= 1 + (igrid_part - i_axis(1)) / delta_i;
		pimg_col						= 1 + (jgrid_part - j_axis(1)) / delta_j;

		%% modulate by bayer filter for this channel
		pimg_idx						= sub2ind([2 2], mod(pimg_row,2)+1, mod(pimg_col,2)+1);
		pimg_bayer						= reshape(bayer(pimg_idx), [n_support n_part_chunk]);
		pimg_ary						= pimg_ary .* pimg_bayer;
		pimg_jac_x						= pimg_jac_x .* pimg_bayer;
		pimg_jac_y						= pimg_jac_y .* pimg_bayer;

		%% weight by support
		if b_fancy_support
			pimg_ary					= bsxfun(@times, support, pimg_ary);
			pimg_jac_x					= bsxfun(@times, support, pimg_jac_x);
			pimg_jac_y					= bsxfun(@times, support, pimg_jac_y);
		end

		%% reshape arrays
		% pimg_ary -> [S^2   n_part]
		% pimg_jac -> [S^2 2 n_part]
		% pimg_row -> [S^2   n_part]
		% pimg_col -> [S^2   n_part]
		pimg_ary						= reshape(pimg_ary  , [n_support   n_part_chunk]);
		pimg_row						= reshape(pimg_row  , [n_support   n_part_chunk]);
		pimg_col						= reshape(pimg_col  , [n_support   n_part_chunk]);
		pimg_jac_x						= reshape(pimg_jac_x, [n_support 1 n_part_chunk]);
		pimg_jac_y						= reshape(pimg_jac_y, [n_support 1 n_part_chunk]);
		%pimg_jac						= cat(2, pimg_jac_x, pimg_jac_y);
		
		%% gather to CPU
		% gather operation is implicit
		pimg_ary_cpu(:,ip)				= pimg_ary;
		pimg_row_cpu(:,ip)				= pimg_row;
		pimg_col_cpu(:,ip)				= pimg_col;
		pimg_jac_cpu(:,1,ip)			= pimg_jac_x;
		pimg_jac_cpu(:,2,ip)			= pimg_jac_y;
		%disp(gpuDevice(1))
	end
	%% debug plot
	%{
	figure(1);
	clf;
	k=22;
	imagesc(igrid_part(:,k), jgrid_part(:,k), otf_ary(:,:,k)');
	hold on;
	plot(ij_pos(1,k), ij_pos(2,k), 'kx', 'MarkerSize', 8);
	axis equal tight xy;
	colormap jet;
	%}
end