function [sA, sB, n0] = back_project_plane(camA, camB, x_O, yA, yB)
	% s = back_project_plane(camA, camB, x0, yA, yB)
	% 
	% for a pair of cameras A and B, find the epipolar planes that points
	% in the image space of either camera lie on
	%
	% The epipolar planes are parameterised by:
	%	position of camera A
	%  position of camera B
	%  position s along a line L, satisfying x = x0 + s*n0
	%
	% sA is the set of s coordinates found for the (2 x NA) array of points yA
	% in the image space of camera A
	% sB is the set of s coordinates found for the (2 x NB) array of points yB
	% in the image space of camera B
	% 
	% The line L is defined by x0 and is normal to the epipolar plane at
	% x = x0.
	
	NA					= size(yA, 2);
	NB					= size(yB, 2);
	
	%% get parameters
	x_A				= camA.xc;
	x_B				= camB.xc;
	
	r_OA				= x_A - x_O;
	r_OB				= x_B - x_O;
	r_AB				= x_B - x_A;
	e_OA				= r_OA / norm(r_OA);
	e_OB				= r_OB / norm(r_OB);
	e_AB				= r_AB / norm(r_AB);
	n0					= cross(e_OA, e_OB);
    n0                  = n0 / norm(n0);
	
	%% exit if trivial
	if NA < 1 || NB < 1
		sA				= [];
		sB				= [];
		return;
	end

	%% find s coordinate for camera A
	% the epipolar plane is defined by three conditions:
	% xA.n + d = 0
	% xB.n + d = 0
	% t.n      = 0
	% t is the vector along ray back projected from image coordinate yA
	%
	% the unit vector n can be found by taking the cross product of unit
	% vectors constructed from t and rAB = xB-xA
	% the plane coefficient d then follows from the 
	
	% calculate unit vector of ray, e_t
	GRA_inv			= inv(camA.G * camA.R);
	GRB_inv			= inv(camB.G * camB.R);
	
	yA_tilde		= [yA; ones(1,NA)];
	yB_tilde		= [yB; ones(1,NB)];
	%tA_mat			= mmx('mult', GRA_inv, yA_tilde);
	%tB_mat			= mmx('mult', GRB_inv, yB_tilde);
	tA_mat			= pagemtimes(GRA_inv, yA_tilde);
	tB_mat			= pagemtimes(GRB_inv, yB_tilde);
	
	
	tA_norm			= sqrt(sum(tA_mat.^2,1));
	tB_norm			= sqrt(sum(tB_mat.^2,1));
	etA_mat			= bsxfun(@rdivide, tA_mat, tA_norm);
	etB_mat			= bsxfun(@rdivide, tB_mat, tB_norm);
	% calculate unit vector of plane, n using cross product
	nA_mat			= [etA_mat(2,:)*e_AB(3) - etA_mat(3,:)*e_AB(2);
							etA_mat(3,:)*e_AB(1) - etA_mat(1,:)*e_AB(3);
							etA_mat(1,:)*e_AB(2) - etA_mat(2,:)*e_AB(1)];
	nB_mat			= [etB_mat(2,:)*e_AB(3) - etB_mat(3,:)*e_AB(2);
							etB_mat(3,:)*e_AB(1) - etB_mat(1,:)*e_AB(3);
							etB_mat(1,:)*e_AB(2) - etB_mat(2,:)*e_AB(1)];
	% substituting in to get s, we find:
	% s = (x_A - x_O).n / n0.n = -r_OA.n / n0.n
	% calculate these products
% 	nA_dot_AO		= mmx('mult', -r_OA', nA_mat);
% 	nB_dot_AO		= mmx('mult', -r_OA', nB_mat);
% 	nA_dot_n0		= mmx('mult',    n0', nA_mat);
% 	nB_dot_n0		= mmx('mult',    n0', nB_mat);
	nA_dot_AO		= pagemtimes( -r_OA', nA_mat);
	nB_dot_AO		= pagemtimes( -r_OA', nB_mat);
	nA_dot_n0		= pagemtimes(    n0', nA_mat);
	nB_dot_n0		= pagemtimes(    n0', nB_mat);
	
	sA					= bsxfun(@rdivide, nA_dot_AO, nA_dot_n0);
	sB					= bsxfun(@rdivide, nB_dot_AO, nB_dot_n0);
end