bounding_polyhedron.m

function [A_hull,b_hull] = bounding_polyhedron(tri, vert)
	% [A,b] = bounding_polyhedron(tri,vert)
	% 
	% determine the system of linear constraints
	%			Ax <= b
	% that describe the inside of the convex polyhedron described by
	% tri and vert. That is, if x lies within the polyhedron, Ax<=b
	%
	% vert is a 3 x n_vert array of vertex coordinates
	% tri	is a n_face x 3 array of indices, describing the triangular
	% faces of the polyhedron
	
	n_vert		= size(vert, 2);
	n_face		= size(tri, 1);
	
	%% determine surface normals
	% vertices of each triangular face
	A				= vert(:, tri(:,1));
	B				= vert(:, tri(:,2));
	C				= vert(:, tri(:,3));
	D				= mean(vert, 2);
	DA				= bsxfun(@minus, A, D);
	% surface normal
	n				= cross(B-A, C-A, 1);
	n				= bsxfun(@rdivide, n, sqrt(sum(n.^2,1)));
	% sign of surface normal is outwards positive
	% D is inside polyhedron so (A-D).n >= 0
	S				= sign(sum(DA.*n, 1));
	n				= bsxfun(@times, S, n);
	
	%% equation of plane: x.n + d = 0
	d				= -sum(n.*A, 1);
	
	%% determine linear constraints
	% each equation of plane has been chosen so that
	% n.x + d = 0
	% and if inside
	% n.x + d > 0
	% since n points outwards
	A_hull		= n';
	b_hull		= -d';
end