function [ij, jacobian] = pinhole_model(xyz, model, b_distort)
	% ij = pinhole_model(xyz, model, b_distort)
	%
	% use pinhole model specified in struct model
	% to map points in object space xyz
	% onto points on image plane
	%
	% xyz is 3 x N
	% ij  is 2 x N
	%
	% default is to apply radial disortion, but this can be turned off
	% by setting b_distort =false
	
	if nargin < 3
		b_distort = true;
	end
	
	%% model parameters
	P			= model.P;
	ic			= model.ic(1);
	jc			= model.ic(2);
	kr			= model.kr;
	
	%% pinhole projection onto distorted coordinates
	N			= size(xyz, 2); 
	lambda	=   P(3,1:3) * xyz + P(3,4);
	ij			= [(P(1,1:3) * xyz + P(1,4)) ./ lambda;
					(P(2,1:3) * xyz + P(2,4)) ./ lambda];
	%% memory management
	clear xyz;
	if nargout < 2
		clear lambda;
	end
				
	%% apply distortion		
	if b_distort
		% undistort
		rsq		= (ij(1,:) - ic).^2 + (ij(2,:) - jc).^2;
		rscale	= (1 + kr(1)*rsq + kr(2)*rsq.^2);
		yprime	= bsxfun(@minus, ij, [ic; jc]);
		ij(1,:)	= ic + yprime(1,:).*rscale;
		ij(2,:)	= jc + yprime(2,:).*rscale;
	end
	
	%% compute jacobian if requested
	if nargout > 1
		% Jacobian in undistorted coordinates
		Tij			= P(1:2, 1:3);
		Lambda		= P(3, 1:3);
		dyj_dxk		= bsxfun(@minus, Tij, bsxfun(@times, reshape(ij,[2 1 N]), Lambda));
		dyj_dxk		= bsxfun(@rdivide, dyj_dxk, shiftdim(lambda,-1));

		% correct Jacobian for distortion
		if b_distort
			yiyj_prime	= bsxfun(@times,	reshape(yprime, [2 1 N]), ...
													reshape(yprime, [1 2 N]));
			dypi_dyj		= bsxfun(@times, eye(2), shiftdim(rscale,-1)) ...
							+ bsxfun(@times, yiyj_prime, kr(1) + 2*kr(2)*shiftdim(rsq,-1));
			clear rsq;
			%jacobian		= mmx('mult', dypi_dyj, dyj_dxk);
			jacobian		= pagemtimes(dypi_dyj, dyj_dxk);
		else
			jacobian	= dyj_dxk;
		end
	end
end
