function [bb_lim, A, b] = bounding_box(camera)
	% bbmax = bounding_box(calib)
	% determine the corners of the smallest axis aligned bounding box
	% which fit encapsulate the extent of the measurement volume
	%
	% bbmax is a 2 x 3 matrix containing the x, y and z limits of the
	% smallest bounding box which encapsulate the entire measurement volume
	
	%% 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);
	for c = 1 : n_cameras
		M_ary(:,:,c)= frustum(camera(c))';
	end
	
	%% solve 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);
	opts				= optimoptions('linprog', 'Display', 'off');
	for i = 1 : 3
		e_i			= [0 0 0]';
		e_i(i)		= 1;
		bb_lim(1,i)	= dot(e_i, linprog(+e_i, A, b, [], [], [], [], [], opts));
		bb_lim(2,i)	= dot(e_i, linprog(-e_i, A, b, [], [], [], [], [], opts));
	end
	
	%% test
	%{
	N					= 1E5;
	bb_size			= bb_lim(2,:) - bb_lim(1,:);
	gl_pos			= bsxfun(@plus, bb_lim(1,:)', diag(bb_size) * rand(3,N));
	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