function H = dlt33_zissermann(xk, yk)
	% compute homography between two homogeneous vectors
	% xk ~ H yk
	% Input:
	%	xk		3xN
	%	yk		3xN
	
	%% number of similarity relations we have
	N			= size(xk, 2);
	assert(size(yk,2) == N);
	assert(size(xk,1) == 3);
	assert(size(yk,1) == 3);
	
	%% set up the reduced measurement matrix A
	% A is a 2n x 9 matrix where individual 2x9 blocks Ai satisfy
	%	Ai h = 0
	% where h is a 12 x 1 vector containing the entries of the homography
	A			= cell(N, 1);
	for k = 1 : N
		x		= xk(:, k);
		y		= yk(:, k);
		% Equation (4.3) 
		A{k}	= [zeros(1,3), -x(3)*y', x(2)*y'; x(3)*y', zeros(1,3), -x(1)*y'];
	end
	A			= cell2mat(A);
	
	%% solution
	% least squares solution
	[~,~,h_vec]	= svds(A, 1, 'smallest');
	% ordering of h is along rows first (Equation 4.2)
	H			= reshape(h_vec, [3 3])'; 
end