function camera = pinhole_fit_zhang(camera, setup)
	%% iterate over views to build set of constraints on intrinsics
	n_views			= length(camera.view);
	V_ary			= cell(n_views, 1);
	H_ary			= cell(n_views, 1);
	for iv = 1 : n_views
		uv_mark			= camera.view(iv).uv_markers;
		ij_mark			= camera.view(iv).ij_markers;
		n_mark			= size(ij_mark,2);

		%% compute the homography
		%H				= dlt33_zissermann([ij_mark; ones(1,n_mark)], [uv_mark; ones(1,n_mark)]);
		H				= dlt33([ij_mark; ones(1,n_mark)], [uv_mark; ones(1,n_mark)]);
		H_ary{iv}		= H;
		%% compute the matrix V
		% equation (7) of Zhang 1999
		v				= cell(2,2);
		for i = 1 : 2
			for j = i : 2
				v{i,j}	= row_constraint(H(:,i), H(:,j));
				%% check
% 				b = randn(6,1);
% 				B = [b(1) b(2) b(4); b(2) b(3) b(5); b(4) b(5) b(6)];
% 				disp(v{i,j}*b - H(:,i).'*B*H(:,j))
			end
			
		end
		% equation (8)
		V_ary{iv}		= [v{1,2}; (v{1,1}-v{2,2})];
	end

	%% formulate a set of constraints upon B (equation 9)
	V					= cell2mat(V_ary);
	
	%% formulate additional constraints for zero skewness and principal point
	if ~setup.b_estimate_skew
		%% add additional skewness constraint, b(2) = 0
		V				= [V; 0 1 0 0 0 0];
		
		%% additional constraint if no shift of principal point
		if ~setup.b_tilt_shift
			% no tilt shift adapter
			% assume the principal point is at the centre of the image
			u0				= mean(camera.i_axis);
			v0				= mean(camera.j_axis);
			V				= [V; u0 0 0 1 0 0; 0 0 v0 0 1 0];
		end
	end
	
	%% formulate additional constraints based on orthogonal translation of plate
	% the plate is constrained to translate orthogonally (in Z) to the world
	% system. this poses an additional set of constraints upon the
	% intrinsic matrix
	for i = 1 : n_views
		for j = (i+1) : n_views
			Hi			= H_ary{i};
			Hj			= H_ary{j};
			c			= Hi(:,3)/norm(Hi(:,2)) - Hj(:,3)/norm(Hj(:,2));
			V			= [V;	row_constraint(Hi(:,1),c);
								row_constraint(Hi(:,2),c);
								row_constraint(Hj(:,1),c);
								row_constraint(Hj(:,2),c)];
		end
	end
		
	%% solve equation 9 in least squares sense
	[~,~,b]				= svds(V, 1, 'smallest');
	% obtain the matrix B from equation 6
	B					= [b(1) b(2) b(4); b(2) b(3) b(5); b(4) b(5) b(6)];

	%% extract the intrinsic matrix 
	if ~setup.b_estimate_skew && ~setup.b_tilt_shift
		alpha			= sqrt(1/abs(B(1,1)));
		beta			= sqrt(1/abs(B(2,2)));
		c				= 0;
	else
		% appendix B of Zhang 
		v0				= (B(1,2)*B(1,3) - B(1,1)*B(2,3))/(B(1,1)*B(2,2)-B(1,2)^2);
		lambda			= B(3,3) - (B(1,3)^2 + v0*(B(1,2)*B(1,3) - B(1,1)*B(2,3)))/B(1,1);
		% hack
		%alpha			= sqrt(lambda/B(1,1));
		%beta			= sqrt(lambda*B(1,1)/(B(1,1)*B(2,2)-B(1,2)^2));
		alpha			= sqrt(abs(lambda/B(1,1)));
		beta			= sqrt(abs(lambda*B(1,1)/(B(1,1)*B(2,2)-B(1,2)^2)));
		c				= -B(1,2)*alpha^2*beta/lambda;
		u0				= c*v0/alpha - B(1,3)*alpha^2/lambda;
	end
	G					= [alpha c u0; 0 beta v0; 0 0 1];
	
	%% estimate camera pose
	t_ary				= zeros([3 n_views]);
	R_ary				= zeros([3 3 n_views]);
	for v = 1 : n_views
		H				= H_ary{v};
		r1				= G \ H(:,1);
		r2				= G \ H(:,2);
		lambda1			= 1 / norm(r1);
		lambda2			= 1 / norm(r2);
		r1				= r1 * lambda1; % using this normalisation ensures unit vector
		r2				= r2 * lambda2;	
		r3				= cross(r1, r2);% ensures right handed system
		t				= sqrt(lambda1*lambda2) * G \ H(:,3);
		R				= [r1 r2 r3];
		R_ary(:,:,v)	= R;
		t_ary(:,v)		= t;
	end
	
	%% estimate the camera orientation
	% we take advantage of the fact that the plate orientation does not change
	% across views to estimate the orientation R
	C_r3				= zeros(12,3);
	C_r2				= zeros(7,3);
	for v = 1 : n_views
		C_r3(v  ,:)		= G\H_ary{v}(:,1); % r3 is orthogonla to r1
		C_r3(v+6,:)		= G\H_ary{v}(:,2); % r3 is orthogonal to r2
		C_r2(v  ,:)		= C_r3(v,:)/norm(C_r3(v,:)); % r2 is orthogonal to r1
	end
	[~,~,r3]			= svds(C_r3, 1, 'smallest');
	C_r2(n_views+1,:)	= r3.';
	[~,~,r2]			= svds(C_r2, 1, 'smallest');
	r1					= cross(r2, r3);
	R					= [r1 r2 r3];
	
	%% estimate camera pose with MATLAB
	cparams				= cameraParameters('IntrinsicMatrix', G');
	[gl_mark,ij_mark]	= plate_to_world(camera, setup);
	[R, xc]				= estimateWorldCameraPose(ij_mark', gl_mark', cparams);
	%R					= Rt';
	xc					= xc(:);
	P					= G*R*[eye(3) -xc(:)];
	
	%% store in structure
	calib				= struct();
	calib.P				= P;
	calib.G				= G;
	calib.xc			= xc;
	calib.R				= R;
	calib.ic			= [u0 v0];
	calib.kr			= [0 0];
	camera.calib		= calib;
	
	%% test
	ij_proj				= pinhole_model(gl_mark, calib, true);
	ij_res				= ij_mark - ij_proj;
	disp(rms(ij_res,2));
end

function V = row_constraint(a,c)
	% write a row constraint V*b = 0 
	% representing a.' * B * c = 0
	% with the symmetric matrix B represented as
	% B		= [b1 b2 b4;
	%          b2 b3 b5;
	%		   b4 b5 b6]
	V	= [	a(1)*c(1);
			a(1)*c(2)+a(2)*c(1);
			a(2)*c(2);
			a(1)*c(3)+a(3)*c(1);
			a(2)*c(3)+a(3)*c(2);
			a(3)*c(3)].';
end