function [gl_pos, gl_dist, ij_dist] = triangulate(camera, ij_pos)
	% [gl_pos, gl_dist, ij_dist] = triangulate(cae
	%
	% given pairs of points yA and yB in image coordinates 
	% on cameras A and B, triangulate the position of these particles
	% by minimising the least-squares distance in object space
	% between the particle and the back-projected rays from either camera
	%
	% inputs:
	%	camA			calibration struct for camera A
	%	camB			calibration struct for camera B
	%	yA				2xN array of position in A, undistorted image coords
	%	yB				2xN array of position in B, undistorted image coords
	%
	% output:
	%	gl_pos		3xN array of position in object space that is the best fit
	%					to the true particle position, given its two projections
	%	gl_dist		2xN array of distance in object space between particle
	%					and either ray
	%	ij_dist		2xN array of distance in image space between particle
	%					and its projection
	
	
	%% calculate ray vectors
	% these are unit vectors which correspond to the direction of rays
	% back-projected from the location of each point on each camera
	b_mat				= zeros(3, n_points, n_cameras);
	xc_mat			= zeros(3, 1, 1, n_cameras);
	for c = 1 : n_cameras
		% unit vector
		G						= camera(c).calib.G;
		R						= camera(c).calib.R;
		GR						= G*R;
		b						= GR \ [ij_dewarped(:,:,c); ones(1,n_points)];
		b_norm				= sqrt(sum(b.^2, 1));
		b_mat(:,:,c)		= bsxfun(@rdivide, b, b_norm);
		xc_mat(:,1,1,c)	= camera(c).calib.xc;
	end
	xc_mat			= repmat(xc_mat, [1 1 n_points 1]);
	
	%% find projection tensor
	% Pb_{ij} = d_{ij} - b_i b_j
	b_outer			= bsxfun(@times,	reshape(b_mat, [3 1 n_points n_cameras]), ...
												reshape(b_mat, [1 3 n_points n_cameras]));
	Pb_mat			= bsxfun(@minus, eye(3), b_outer);
	Pb_sum			= sum(Pb_mat, 4);
	
	%% project origin
	Pb_x0_mat		= mmx('mult', Pb_mat, xc_mat);
	Pb_x0_sum		= sum(Pb_x0_mat, 4);
	%% solve for 3D position
	gl_pos			= mmx('backslash', Pb_sum, Pb_x0_sum);
	gl_pos			= squeeze(gl_pos);
end
