How to calculate the area of intersection between a 3D volume and a 2D plane

Question

Hello if anyone can offer insight on how to solve my problem that would be great!

I am looking to calculate the area of intersection between a 3D volume and a 2D plane. 3D volume: defined by 6 points (will always be a 3D wedge resembling a thick triangle) 2D plane: defined by 3 points

any suggestions? are there any open source software out there that can perform this task?

thanks very much Chris

Answer 1

There is not going to be a magic bullet answer to this question; at some point you will just have to suck it up and consider all the cases. I once had to compute the intersection of a triangle and a circle... it was awful. Since your 3D shape has certain symmetries (like the fact that it is always a triangular prism) helps narrow down the possibilities tremendously.

1. Find which side of the plane each point lies.
2. If all points are on one side, done.
3. If one point is separated from the other 5 by the plane, the intersection is a triangle. Presumably, you have connectivity information to tell you which faces are incident on the vertex, and you can compute the 3 lines due to plane-plane intersections, and then the intersection triangle. Best do all this in a parameterization of the intersecting plane.
4. Two points are separated, either both are on the same triangular face of the prism, or they are two corresponding points on opposite triangular faces. In either case, the intersection is a quadrilateral.
5. For three points, they could all be the same triangular face or only two on it. You have two cases here, for a resulting triangle or hexagon.

Those are really the only high level cases, there are obviously symmetry reductions you need to apply (in case 5, the one point case is equivalent to the 2 point case on the other triangular facet).

My one suggestion to you is to pick a robust representation for your shapes, and to use robust predicates to perform the geometric tests. For the plane, it is best represented as a point on the plane, and a unit normal vector. The prism is represented by defining a orthonormal basis (triad) with one axis aligned with the extrusion direction of the prism. Let one vertex be at the origin, the other two on the triangular face be represented in the uv-coordinates of the triad, and then you only need to store its height, and global offset of its base vertex. Essentially, I'm thinking

struct plane{
    double p[3]; // point on plane
    double n[3]; // unit normal vector to plane
};
struct TriPrism{
    double basis[9]; // 3x3 orthogonal matrix of local coordinate frame (det = +1)
    // Stored columnwise, first two vectors are in the plane of the triangular face)
    // Last vector is parallel to extrusion direction, call the set [u, v, w]
    double base[3]; // global coordinates of base vertex (where the basis vectors are "rooted")
    double buv[2]; // the uv-coordinates of the second point on the triangular face
    double cuv[2]; // third point on triangular face
    // Assume that the "bottom" triangular face is formed by vertices (a,b,c)
    // with base being a, and (b-a) cross (c-a) directed along vector w (instead of against)
    double h; // height of prism ("top" triangular face is h*w offset from the "bottom" face)
};

This representation is robust to moving the prism around and should maintain high relative precision for all but the most pathological cases.

As for robust predicates, I highly recommend Shewchuk's robust predicates assuming you use ordinary floating point numbers. You would mainly use orient3d.

In terms of representing the final output polygon, choose a suitable parameterization of the plane. The best choice, in my opinion, is to first choose an orthonormal basis set in the plane, determined by Gram-Schmidt on the plane normal vector (see geom_maketriad3d here). Then, let the origin of that 2D parametric plane be the projection of the prism base point into the plane. This ensures your parameterization is actually rooted near where the resulting polygon will be, to ensure effective use of the bits of the floating point numbers involved. Do all the remaining computations using this parameterization if at all possible.

Generally, geometric computation is fraught with peril due to these silly numerical considerations (a slight error with a vertex near the plane in case 5 above could drastically change the result from a 4-gon to a 6-gon). I suggest you handle a few cases at a time, and try to visualize the result. I have a rather simple viewer program here, with a kind of postscript-like input language capability. You can dump the facets of the prism and draw the plane and also draw the resulting polygons and visually check them to see if they're right.

Addendum: I forgot that you originally wanted the area of the intersection polygon. It is trivial, but in case you don't know, once you have the collection of lines defining the edges of the polygon in the cut plane, you have to convert it to a vertex representation. Since you constructed the intersections from facetplane-and-plane intersections, you should have been able to keep track of the ordering of the lines, and hence the ordering of the edges going around the polygon. You only need to compute the intersection of successive pairs of these lines to get the vertices. Once you have the vertices, it's a simple matter to get the area (see e.g. geom_polygon_area2d here). If you worked entirely in uv-coordinates of the cut-plane, then you can feed them directly to such a function to get the area.

I should add that there is an obvious dumb approach, which is to pick a suitable large region in the cut-plane, and randomly sample points and checking to see if they're in the prism (which is cheap, since it's convex). Then you can calculate the area as a ratio of points inside vs. total points, times the area of the sampling region. If you really don't care about accuracy, then this will probably be faster, but otherwise, the analytic method should not be much slower, despite its combinatorial complexity.