function [bb_lim_global, A, b] = bounding_box3(camera, A_additional, b_additional)
	% 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
	%
	% optionally, add constraints A_additional*x <= b_additional
	%	A_additional	M x 3	
	%	b_additional	M x 1
	%	
	
	if nargin < 3
		A_additional = [];
		b_additional = [];
	end
	
	%% 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
	
	%% identify all possible sets of three cameras
	cam_perms		= nchoosek(1:n_cameras, 3);
	n_perms			= size(cam_perms, 1);
	opts			= optimoptions('linprog', 'Display', 'off');
	b_feasible		= true(1,n_perms);
	bb_lim_global	= [inf(1,3); -inf(1,3)];
	
	for iperm = 1 : n_perms
		%% define linear programming problem for this subset of cameras
		ic					= cam_perms(iperm, :);
		% to obtain volume bounds, we find the x that 
		% minimises c'x subject to Ax < b
		A					= reshape(-M_ary(1:3,:,ic), [3 12])';
		b					= reshape( M_ary(4  ,:,ic), [12 1]);
		% add additional constraints
		A					= [A; A_additional];
		b					= [b; b_additional];
		
		%% test feasibility of problem
		[x,fval,flag,output,lambda] = linprog(ones(3,1),A,b,[],[],[],[],opts);
		if flag <= 0
			%fprintf('%d %d %d not suitable\n', ic);
			b_feasible(iperm) = false;
			continue;
		end
		
		%% identify local set of bounds
		bb_lim				= zeros(2,3);
		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
	
		fprintf('%d %d %d is bounded [%+6.1f - %+6.1f, %+6.1f-%+6.1f, %+6.1f-%+6.1f]\n', ic, bb_lim);
		
		%% expand bounds
		bb_lim_global	= [min(bb_lim_global(1,:), bb_lim(1,:));
						   max(bb_lim_global(2,:), bb_lim(2,:))];
	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