Computational geometry with Matlab and MPT

For computational geometry problems there is a rather splendid tool called the Multi Parametric Toolbox for Matlab Here are a few examples of using this toolbox for a variety of purposes using 3D convex polyhedra represented as polytope arrays. In the following examples, P is always a domain bounded 3D polytope array

	V=volume(P);
	% return a vector of the number of sides in each polytope
	[H,K]=double(P);
	[numsides, Cncols] = cellfun(@size, H);
	% return a vector of polytope radii
	[xyz,Rad] = chebyball(P);
	% find the extrema of each polytope in P
	for ip=1:length(P)
	E(ip)=extreme(P(ip));
	end
	%Slice through a 3-dimensional volume. Here, I'm slicing through 3D polytope array P to return the number of polytopes in each slice (N) and the mean diameter of polytopes in each slice (M). Slices are made at arbitrary values in the space [0,1]
	M=[]; N=[];
	for ii=[0:.01:.2,.25:.05:.8]
	[Pcut, sliced] = slice(P,2,ii);
	[xc,rc] = chebyball(Pcut); M=[M;mean(rc.*2)]; N=[N;length(rc)];
	end %ii
	% here I find the indices where of polytopes whose centres are in the space bounded by x=min and y=max
	index=[];
	for ip=1:length(P)
	E=extreme(P(ip));
	if(any((r(1,ip)+E(:,1))<min) any((r(2,ip)+e(:,2))="" ||=""><min) any((r(3,ip)+e(:,3))="" ||=""><min)) continue="" elseif(any((r(1,ip)+e(:,1))="">max) || any((r(2,ip)+E(:,2))&gt;max) || any((r(3,ip)+E(:,3))&gt;max))
	continue
	else
	index=[index;ip];
	end
	% Here I reduce each 3D polytope in array P by factor fct to give new polytope array Q
	Q = [];
	for ip=1:length(P)
	E=extreme(P(ip));
	ey=E(:,2); ex=E(:,1); ez=E(:,3);
	ex = fct*ex+(1-fct)*mean(ex); % expand x by factor fct
	ey = fct*ey+(1-fct)*mean(ey); % expand y
	ez = fct*ez+(1-fct)*mean(ez); % expand y
	Er=[ex,ey,ez];
	Q = [Q,hull(Er)];
	end
	% Finally, here's a piece of code which will a) convert a 3D polytope array into a 3D matrix, then b) display slices of that volume (left side, back, and the centre line in both the x and y orientation) % first, slice P at small increments and build up matrix im
	im=zeros(501,601,length([0:.01:1]));
	counter=1;
	for i=[0:.01:1]
	Pcut=slice(Q,2,i);
	plot(Pcut)
	axis([0 1 0 1])
	axis off
	print -dtiffnocompression tmp.tiff
	I=double(rgb2gray(imread('tmp.tiff')));
	Ib=(rescale(I,0,1)~=1);
	Ib=Ib(200:700,400:1000);
	im(:,:,counter)=Ib;
	counter=counter+1;
	end
	close all
	delete tmp.tiff
	% resize by 50% (helps speed things up)
	im_50percent=imageresize(im,.5,.5);
	im=im_50percent;
	% create a grid
	[x,y,z]=meshgrid([1:size(im,2)],[1:size(im,1)],[1:size(im,3)]);
	xmin = min(x(:));
	ymin = min(y(:));
	zmin = min(z(:));
	xmax = max(x(:));
	ymax = max(y(:));
	zmax = max(z(:));
	% get the first slice
	hslice = surf(linspace(xmin,xmax,100),...
	linspace(ymin,ymax,100),...
	zeros(100));
	rotate(hslice,[-1,0,0],-90)
	xd = get(hslice,'XData');
	yd = get(hslice,'YData');
	zd = rescale(get(hslice,'ZData'),zmin,zmax);
	delete(hslice)
	% then build more slices on top
	h = slice(x,y,z,im,xd,yd,zd);
	set(h,'FaceColor','flat',...
	'EdgeColor','none',...
	'DiffuseStrength',5)
	colormap gray
	hold on
	hslice = surf(linspace(xmin,xmax,100),...
	linspace(ymin,ymax,100),...
	zeros(100));
	rotate(hslice,[0,-1,0],-90)
	xd = get(hslice,'XData');
	yd = get(hslice,'YData');
	zd = rescale(get(hslice,'ZData'),zmin,zmax);
	delete(hslice)
	h = slice(x,y,z,im,xd,yd,zd);
	set(h,'FaceColor','flat',...
	'EdgeColor','none',...
	'DiffuseStrength',5)
	hx = slice(x,y,z,im,xmax,[],[]);
	set(hx,'FaceColor','flat','EdgeColor','none')
	hy = slice(x,y,z,im,[],ymax,[]);
	set(hy,'FaceColor','flat','EdgeColor','none')
	% hz = slice(x,y,z,im,round(xmax/2),[],[]);
	% set(hz,'FaceColor','flat','EdgeColor','none')
	axis tight
	box on
	camproj perspective
	lightangle(-45,45)
	set(gcf,'Renderer','zbuffer')
	warning on all

view raw mpt_examples.m hosted with ❤ by GitHub