function [pimg_ary, pimg_jac, pimg_row, pimg_col] = otf_evaluate_diffraction_interpolant(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
	%  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
	%  irel        i coordinate, relative to particle center, px
	%  jrel        i coordinate, relative to particle center, px
	%  I           particle image integrated on sensor I(irel, jrel, phi)
	%  gradI       gradient of I with respect to [irel,jrel]


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

	%% support
	otf							= otf(channel);
	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?

	%% convert to gpuArray if necessary
	if isa(im_pos,'gpuArray')
		otf.irel				= gpuArray(otf.irel);
		otf.jrel				= gpuArray(otf.jrel);
		otf.phi					= gpuArray(otf.phi);
		otf.I					= gpuArray(otf.I);
		otf.gradI				= gpuArray(otf.gradI);
	end

	%% 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);
	bayer						= camera.bayer(:,:,channel);
	assert(size(bayer,1) == 2 && size(bayer,2) == 2, 'bayer pattern was unexpected size');
	assert(abs(delta_i)==1);
	assert(abs(delta_j)==1);

	%% model parameters
	f                       = otf.f;				% focal length
	D                       = otf.D;				% aperture
	wavelength              = otf.lambda;			% 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);						
	[imat,jmat]				= ndgrid(irel, jrel);
	
	%% 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);
	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 - ij0_pos;
	imat_part				= bsxfun(@minus, imat(:), dij_pos(1,:));
	jmat_part				= bsxfun(@minus, jmat(:), dij_pos(2,:));
	phimat_part				= repmat(phi_part(:).', [n_support 1]);

	% evaluate particle image and its jacobian on finely spaced grid
	% result is of size n_fine x n_part
	pimg_ary				= interpn(otf.irel, otf.jrel, otf.phi, otf.I, ...
		imat_part, jmat_part, phimat_part, 'linear', 0);
	% jacobian with respect to image space coordinates
	pimg_jac_x				= -interpn(otf.irel, otf.jrel, otf.phi, otf.gradI(:,:,:,1), ...
		imat_part, jmat_part, phimat_part, 'linear', 0);
	pimg_jac_y				= -interpn(otf.irel, otf.jrel, otf.phi, otf.gradI(:,:,:,2), ...
		imat_part, jmat_part, phimat_part, 'linear', 0);
	
	%% 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]);
	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]);
	pimg_row						= reshape(pimg_row  , [n_support   n_part]);
	pimg_col						= reshape(pimg_col  , [n_support   n_part]);
	pimg_jac_x						= reshape(pimg_jac_x, [n_support 1 n_part]);
	pimg_jac_y						= reshape(pimg_jac_y, [n_support 1 n_part]);
	pimg_jac						= cat(2, pimg_jac_x, pimg_jac_y);
		
	%% 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