Determine intersection between plane and polyhedron

Given a plane in 3D space and a polyhedron with also extends out in 3d space - assuming that the plane and polyhedron always intersect - how do I figure out the portion of the plane that is enclosed by the polyhedron?

[Torsten]

"portion of the plane" means "area of the plane within the polyhedron" ?

Answer 1

[John D'Errico]

Um, not too hard, as long as the polyhedron is convex, so perhaps a convex hull. If it is not convex, then it is more difficult. As always, difficulty is measured by your skill at coding, and your understanding of the relevant mathematics. For example...

xyz = rand(100,3);
H = convhulln(xyz);

Now I'll choose some simple plane. A plane is defined by a point in the plane, and the normal vector to the plane. Lets see. I'll pick the point as:

Pip = [0.2 0.3 0.4];
Pnorml = [1 -2 2];
Pnorml = Pnorml/norm(Pnorml); % Normalize to unit length

So pretty arbitrary. I do know the plane will intersect the polyhedron, but I just picked some arbitrary numbers.

First, transform the problem so the point in plane is the origin, then rotate the polyhedron to have the plane as two of the axes. We will undo that rotation later on.

xyzhat = xyz - Pip;
Pbasis = null(Pnorml);

uvw = xyzhat*[Pbasis,Pnorml.'];

Each row of the array uvw has the first two elements representing the projection into the plane of each vertex of the set, The third column represents the distance from the plane. A positive number indicates the point lies "above" the plane. Negative means the point lies below.

Next, we get the list of all edges in the polyhedron. They will each appear twice, so drop the dups.

edgelist = [H(:,[1 2]);H(:,[2 3]);H(:,[1 3])];
edgelist = sort(edgelist,2);
edgelist = unique(edgelist,'rows');

Next, we care only about those edges that cross the plane. Throw away those that do not.

k = uvw(edgelist(:,1),3).*uvw(edgelist(:,2),3)<=0;
edgelist = edgelist(k,:);

Find the intersection point in the plane of those edges we just found.

t = uvw(edgelist(:,1),3)./(uvw(edgelist(:,1),3) - uvw(edgelist(:,2),3));
uvinplane = uvw(edgelist(:,1),[1 2]).*(1-t) + uvw(edgelist(:,2),[1 2]).*t;

Next, compute the convex hull of those points. The planar intersection of a plane with a convex domain will also be convex.

planarhulledges = convhull(uvinplane(:,1),uvinplane(:,2));

Finally, transform these points back into the original domain.

xyzhull = uvinplane*Pbasis.' + Pip;
plot3(xyz(:,1),xyz(:,2),xyz(:,3),'b.')
hold on
plot3(xyzhull(:,1),xyzhull(:,2),xyzhull(:,3),'ro')
plot3(xyzhull(planarhulledges,1),xyzhull(planarhulledges,2),xyzhull(planarhulledges,3),'r-')
grid on
box on

I could do a better job in the graphics, plotting the facets of the original polyhedron as semi-transparent. I might plot the intersection plane as a filled patch. But you should get the idea, and I've done enough for you here.

And again, if the original polyhedron is not convex, then things get more complicated, but it is still doable.

Answer 2

[Matt J]

Use intersectionHull() from this FEX download,

https://www.mathworks.com/matlabcentral/fileexchange/30892-analyze-n-dimensional-convex-polyhedra?s_tid=srchtitle