function [pimg_ary, pimg_jac, pimg_row, pimg_col] = otf_evaluate_diffraction_bicubic(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
	%  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

	%% 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);
	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(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
	                    
	%% define grid for integration of diffraction limited profile	
	% define the grid of neighbouring pixel edges
    irel_edge               = (0 : S(1)) - S(1)/2;
    jrel_edge               = (0 : S(2)) - S(2)/2;
	irel					= (0 : S(1)-1) - floor(S(1)/2);
	jrel					= (0 : S(2)-1) - floor(S(2)/2);						
	[imat,jmat]				= ndgrid(irel, jrel);
    [imat_edge,jmat_edge]   = ndgrid(irel_edge, jrel_edge);
	
	%% 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 
	im0_pos					= round(im_pos);
	igrid_part				= bsxfun(@plus, im0_pos(1,:), imat(:));
	jgrid_part				= bsxfun(@plus, im0_pos(2,:), jmat(:));
	% y^ = (y - Y) / sigma
	%    = (y - Y0 - (Y-Y0) / sigma
	dim_pos					= im_pos - im0_pos;
	irel_part				= bsxfun(@minus, imat_edge(:), dim_pos(1,:));
	jrel_part				= bsxfun(@minus, jmat_edge(:), dim_pos(2,:));
    
	r_part                  = sqrt(irel_part.^2 + jrel_part.^2);
    phi_part                = 2*pi*(Z_0 - Z)/wavelength*(D/(2*f))^2;
	phi_part				= repmat(phi_part, [size(irel_part,1) 1]);
    psi_part                = pi * (pitch/wavelength) * (D/f) * r_part;

	% evaluate particle image and its jacobian on finely spaced grid
	% result is of size n_fine x n_part
	I_part                  = interpn(otf.phi, otf.psi, otf.I, phi_part, psi_part, 'linear', 0);
	I_part_dr               = pi * (pitch/wavelength) * (D/f) * ...
        interpn(otf.phi, otf.phi, otf.gradI(:,:,2), phi_part, psi_part, 'linear', 0);
	I_part_d2r              = (pi * (pitch/wavelength) * (D/f))^2 * ...
        interpn(otf.phi, otf.phi, otf.d2Idphi, phi_part, psi_part, 'linear', 0);
	
    % jacobian with respect to image space coordinates 
	I_part_dx               = -I_part_dr .* irel_part ./ r_part;
	I_part_dy               = -I_part_dr .* jrel_part ./ r_part;
	I_part_dxdy             = +I_part_d2r.* irel_part .* jrel_part ./ r_part.^2;
    
    %% integrate image and its jacobian
    % use bicubic integration method
    % see integrate_bicubic.m
    % and maths_bicubic.m
	
	% reshape
	I_part					= reshape(I_part, [S+1, n_part]);
	I_part_dx				= reshape(I_part_dx, [S+1, n_part]);
	I_part_dy				= reshape(I_part_dy, [S+1, n_part]);
	I_part_dxdy				= reshape(I_part_dxdy, [S+1, n_part]);
	% get values at corners
    f00                     = I_part(1:end-1,1:end-1,:);
    f01                     = I_part(1:end-1,2:end  ,:);
    f10                     = I_part(2:end  ,1:end-1,:);
    f11                     = I_part(2:end  ,2:end  ,:);
	fx00                    = I_part_dx(1:end-1,1:end-1,:);
    fx01                    = I_part_dx(1:end-1,2:end  ,:);
    fx10                    = I_part_dx(2:end  ,1:end-1,:);
    fx11                    = I_part_dx(2:end  ,2:end  ,:);
	fy00                    = I_part_dy(1:end-1,1:end-1,:);
    fy01                    = I_part_dy(1:end-1,2:end  ,:);
    fy10                    = I_part_dy(2:end  ,1:end-1,:);
    fy11                    = I_part_dy(2:end  ,2:end  ,:);
    fxy00                   = I_part_dxdy(1:end-1,1:end-1,:);
    fxy01                   = I_part_dxdy(1:end-1,2:end  ,:);
    fxy10                   = I_part_dxdy(2:end  ,1:end-1,:);
    fxy11                   = I_part_dxdy(2:end  ,2:end  ,:); 
    % integrate
    pimg_ary				= delta_i*delta_j*(f00/4 + f01/4 + f10/4 + f11/4 + fx00/24 + fx01/24 - fx10/24 - fx11/24 + fxy00/144 - fxy01/144 - fxy10/144 + fxy11/144 + fy00/24 - fy01/24 + fy10/24 - fy11/24);
    % jacobian
    pimg_jac_x				= delta_j*(f10/2 - f01/2 - f00/2 + f11/2 - fy00/12 + fy01/12 + fy10/12 - fy11/12);
    pimg_jac_y				= delta_i*(f01/2 - f00/2 - f10/2 + f11/2 - fx00/12 + fx01/12 + fx10/12 - fx11/12);
	% flatten back
	pimg_ary				= reshape(pimg_ary, [n_support n_part]);
	pimg_jac_x				= reshape(pimg_jac_x, [n_support n_part]);
	pimg_jac_y				= reshape(pimg_jac_y, [n_support n_part]);
	
    %% 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