function P = dlt33(xk, yk)
	% solve the set of similarity relations
	% xk ~ Pyk
	% 
	% for xk = [x1k x2k x3k]' is a 3x1 vector
	%     yk = [y1k y2k y3k]' is a 3x1 vector
	% 
	% to find the 3x3 matrix P
	% subject to |P| = 1, where |P| is the Froebenius norm
	
	% set up the problem as solving
	% xk' * H * P * yk = 0
	%
	% leveraging xk'*H*xk = 0
	% however, for xk on R^3, the dimension of H_m is 3
	% i.e. there are three basis-matrices H_m that can satisfy xk'*H*xk = 0
	% these are:
	% H1 = [0 0 0; 0 0 -1; 0 1 0]
	% H2 = [0 0 1; 0 0 0; -1 0 0]
	% H3 = [0 -1 0; 1 0 0; 0 0 0]
	%
	% so we get three sets of equations
	% bk' * p = [0 0 0]'
	%
	% 0 = P2_1*x3*y1 + P2_2*x3*y2 + P2_3*x3*y3 - P3_1*x2*y1 - P3_2*x2*y2 - P3_3*x2*y3
	% 0 = P3_1*x1*y1 - P1_2*x3*y2 - P1_3*x3*y3 - P1_1*x3*y1 + P3_2*x1*y2 + P3_3*x1*y3
	% 0 = P1_1*x2*y1 + P1_2*x2*y2 + P1_3*x2*y3 - P2_1*x1*y1 - P2_2*x1*y2 - P2_3*x1*y3
	
	% number of similarity relations we have
	N			= size(xk, 2);
	assert(size(yk,2) == N);
	
	% set up matrix B
	% with rows b_k
	% that satisfy b_k * p = 0
	B			= zeros(N*3, 9);
	for k = 1 : N
		x		= xk(:, k);
		y		= yk(:, k);
		B(3*(k-1)+1, :) = [  0; x(3)*y(1); -x(2)*y(1);  0; x(3)*y(2); -x(2)*y(2);  0; x(3)*y(3); -x(2)*y(3)];
		B(3*(k-1)+2, :) = [ -x(3)*y(1); 0; +x(1)*y(1); -x(3)*y(2); 0; +x(1)*y(2); -x(3)*y(3); 0; +x(1)*y(3)];
		B(3*(k-1)+3, :) = [ +x(2)*y(1); -x(1)*y(1); 0; +x(2)*y(2); -x(1)*y(2); 0; +x(2)*y(3); -x(1)*y(3); 0];
	end
	
	% now our solution p is the solution of
	% B*p = 0
	% i.e. p is in the null space of B
	% but due to errors, this is not exact, so we seek a solution in the
	% least squares sense
	% 
	% i.e. minimise |Bp| subject to |p| = 1
	% this can be solved with SVD
	%
	% p is the unit singular vector of B corresponding to the last column of
	% V, when the elements in S are ranked in decreasing size
	% which, helpfully, matlab does by default
	[U,S,V]	= svd(B);
	sigma		= diag(S);
	p			= V(:, end);
	
	P			= reshape(p, [3 3]);
end
	
		
	 