function [bb_lim, R, A, b] = min_bounding_box(camera)
	% [bb_lim, R] = min_bounding_box(calib)
	% determine the corners of the smallest arbitrarily oriented bounding box
	% which fit encapsulate the extent of the measurement volume
	%
	% bb_lim is a 2 x 3 matrix containing the x, y and z limits of the
	% smallest bounding box which encapsulate the entire measurement volume
	% in a rotated coordinate system
	%
	% R is a rotation matrix describing the rotation of object space
	% coordinates into a coordinate system in which all points fit within
	% the limits 
	
	%% get view frustums
	% M = M_ary(:,p,c) are the coefficients for the equation of plane p of the
	% viewing frustum for camera c. 
	% i.e.
	%		M(1)x + M(2)y + M(3)z + M(4) >= 0
	% when a point lies within the measurement volume
	n_cameras		= length(camera);
	M_ary				= zeros(4,4,n_cameras);
	U_ary				= zeros(3,4,n_cameras);
	for c = 1 : n_cameras
		[M,U]			= frustum(camera(c));
		M_ary(:,:,c)= M';
		U_ary(:,:,c)= U;
	end
	
	%% get orientation of view frustum
	R_ary				= zeros(3,3,n_cameras);
	normalise		= @(x) x / norm(x);
	for c = 1 : n_cameras
		e1				= normalise(U_ary(:,1,c));
		e2				= normalise(U_ary(:,2,c) - U_ary(:,1,c));
		e3				= normalise(cross(e1,e2));		
		R_ary(:,:,c)= [e1, e2, e3];
	end
	
	%% setup linear programming problem
	% to obtain volume bounds, we find the x that 
	% minimises c'x subject to Ax < b
	A					= reshape(-M_ary(1:3,:,:), [3 4*n_cameras])';
	b					= reshape( M_ary(4  ,:,:), [4*n_cameras 1]);
	bb_lim			= zeros(2,3,n_cameras);
	bb_vol			= zeros(1,n_cameras);
	opts				= optimoptions('linprog', 'Display', 'off');
	
	%% try orientations corresponding to orientation of cameras
	for c = 1 : n_cameras
		for i = 1 : 3
			e_i			= R_ary(:,i,c);
			bb_lim(1,i,c)	= dot(e_i, linprog(+e_i, A, b, [], [], [], [], [], opts));
			bb_lim(2,i,c)	= dot(e_i, linprog(-e_i, A, b, [], [], [], [], [], opts));
		end
		bb_size		= bb_lim(2,:,c) - bb_lim(1,:,c);
		bb_vol(c)	= abs(prod(bb_size));
	end
	
	%% identify smallest volume bounding box
	[~,c_best]		= min(bb_vol);
	R					= R_ary(:,:,c_best);
	bb_lim			= bb_lim(:,:,c_best);
	
	%% test
	%{
	N					= 1E5;
	bb_size			= bb_lim(2,:) - bb_lim(1,:);
	bb_pos			= bsxfun(@plus, bb_lim(1,:)', diag(bb_size) * rand(3,N));
	gl_pos			= R*bb_pos;
	bb_dist			= bsxfun(@minus, A*gl_pos, b);
	b_inside			= all(bb_dist <= 0, 1);
	figure(1);
	clf;
	plot3(gl_pos(1, b_inside), gl_pos(2, b_inside), gl_pos(3, b_inside), 'g.');
	hold on;
	plot3(gl_pos(1,~b_inside), gl_pos(2,~b_inside), gl_pos(3,~b_inside), 'r.');
	axis equal tight xy;
	%}
end