An algorithm to find the area of intersection between a convex polygon and a 3D polyhedron?

Suppose I have a set of points in 3D which are all co-planar, and which describe the vertices of a convex polygon. I know the coordinates of all of these vertices, I know the unit normal to the plane of the polygon, and I know the order (in a right-hand sense with respect to the normal) in which the vertices go around the polygon.

I also have a polyhedron, which is in general non-convex, and the surface of which is made of triangles. It is a triangle mesh like the one shown here. It describes a completely enclosed volume (no holes in the surface). I know the coordinates of all the vertices of all the triangles that make up the polyhedron. For a given vertex, I know the list of which triangles it is a vertex of.

One other thing I know, which may or may not be useful, is that the polygon lies "in the vicinity of" a particular pair of vertices of the polyhedron, say the ones with indices $i$ and $j$. So if there were to some kind of a search of the surface of the polyhedron, that would be the place to start.

My task is to find the area of the intersection between the polygon and the polyhedron. I believe this is a well-defined quantity. I need an algorithm to calculate this, hopefully efficiently but at this point I'll take anything.

My thoughts:

I can classify each polygon vertex as either inside or outside the polyhedron. I've written a ray-casting algorithm that computes this fairly accurately and efficiently. That means that each edge of the polygon can be classified as either "both vertices inside", "both vertices outside", or "one in, one out." I can calculate the point of intersection between a given polygon edge and the polyhedron by checking a bunch of nearby triangles.

But that's not enough. For instance, just because both vertices of an edge are inside doesn't mean the whole edge is inside. A "finger" of the polyhedron might intersect an edge without touching either of its endpoints. Furthermore, the polyhedron might intersect the polygon in such a way that it doesn't touch any of the edges at all. Clearly, one must also calculate the intersections between the triangles' edges and the polygon as well.

It gets a bit too complicated for me at that point, and I'm wondering if there is some established algorithm for this?

Perhaps something like: Project the polyhedron onto the plane containing the polygon and then... I don't know.

• If you already have the ray-casting algorithm, why not cast a ray from a vertex along its edge and find whether the first intersection is before the distance to the opposite vertex? Then you can find the intersections on each edge. – Dan Hulme Jun 2 '17 at 14:34
• That doesn't help you where a finger of the mesh pokes through the polygon, though. Perhaps you're better off starting by classifying the mesh's vertices by which side of the plane of the polygon they're on. – Dan Hulme Jun 2 '17 at 14:47
• @DanHulme Once I've determined which mesh vertices are on which side of the plane, what would I do with that information? – John Barber Jun 2 '17 at 15:11
• Use it to find mesh edges that cross the plane and intersect them with the plane to find all the crossing points. You can use the crossing points to make a new polygon which is the section of the plane and the polyhedron. (You need to use the topology of the polyhedron to know what order to connect the edges together, because this new polygon is non-convex.) Now just intersect the two polygons in the usual way. I'm making this up as I go along but it sounds plausible. – Dan Hulme Jun 2 '17 at 17:06
• @DanHulme So it seems to me that there are two kinds of crossing points: Where an edge of the polygon crosses a triangle, and where an edge of a triangle crosses a polygon. Both kinds can be found by some kind of exhaustive search. I guess the next step would be to throw away vertices of the original polygon that were known to be outside the polyhedron. How do I use the topology of the polyhedron to order the resulting bunch of points? – John Barber Jun 2 '17 at 17:31

There are different numerical approximations you could use:

A simple solution is to use brute-force Monte Carlo integration. Distribute $N$ random points on the polygon and calculate the number of points inside the 3D polyhedron $N_i$ using ray tracing. The area inside the polyhedron is $A_i=A\frac{N_i}{N}$, where $A$ is the area of the polygon.

To improve convergence of the above algorithm, you can distribute the points on grid on the polygon and jitter the sample positions within each grid cell. This is called stratified sampling.

Further improvement would be to perform the MC integration in 1D instead of 2D. Pick a sample line $L$ that's on the polygon plane and a sample direction $D$ that's on the plane perpendicular to $L$. Select sample line segment $L_s$ that bounds the polygon within $D$-directed lines (min & max of the polygon projection on $L$ for tight fit). For each sample on $L_s$ first calculate the line segment $S$ that intersects the polygon in direction $D$. Next calculate the length of $S$ that's inside the polyhedron $S_i$. You can calculate $S_i$ using ray tracing, by starting from one end of $S$ and marching along $S$ to intersection points of the segment with the polyhedron, summing up the segments of $S$ inside the polyhedron. More specifically, you need to first check if the start point on $S$ is within polyhedron and set the "inside" state accordingly. Then for each intersection flip the state and add the length from previous intersection to $S_i$ if the state changes from inside to outside. The approximation of the polygon area inside the 3D polyhedron is: $$A_i=\frac{|L_s|}{N}\sum_{i=1}^N S_i$$

If you have GL (or equivalent) available, the easiest way is probably to set up your projection matrix so that the plane of the polygon is the near clipping plane, draw the polygon into the stencil buffer, and then draw the polyhedron such that inside faces output a 1 fragment. (You could do this by flipping the winding of the polyhedron and turning on back-face culling, or you could check the normal in the fragment shader.) Then the sum of the white pixels in your framebuffer is the intersection area.

• I'm afraid I don't what any of what you just said means. I have almost no background in computer graphics. I'm learning this as I go over the last year or so. I'm writing this code from scratch in C++ for a computational physics application. For reasons I won't go into, this has to run in an environment with almost nothing available. Imagine that this has to be done on a non-network-connected machine in 1997. So straight up from scratch with no special libraries. – John Barber Jun 2 '17 at 15:06
• It comes down to the idea you suggested about projecting the mesh into the plane of the polygon. The near clipping plane part cuts off the parts of the mesh "in front of" the plane. Only drawing the inside faces means you'll only see the inside of the mesh where this cut has made a hole, giving you an image of the inside of the polyhedron. Finally, the stencil buffer intersects that image with your polygon. If you had a 3D library available it would be easy to write (because you don't need to do all the maths yourself) and fast (run on the GPU). – Dan Hulme Jun 2 '17 at 17:10
• What is "the near clipping plane part"? What do you mean when you say "inside faces"? Does that mean triangles which partially overlap the polygon in projection? What is the "stencil buffer"? As to using a 3D library... no can do. Some networks are VERY selective about what code can and can't be stored on it. I don't mind doing the math of implementing an algorithm. I've had a lot of practice at that. I just need to understand the geometry of the algorithm. – John Barber Jun 2 '17 at 17:27