Question [Algorithm Question] Mesh-plane intersection to edge loop polygons

This isn't Unity-specific, but rather a general computational geometry question. I've found some resources that do about 90% of what I need, but not 100%

Given a mesh (triangles + vertex positions), I need to get the 2D polygons where it intersects the XY plane. There could be more than one (for example, consider a torus), and they might not be closed (consider a single triangle crossing the plane). However, at least for the ones that are closed, I need the complete silhouette/cross-section of intersection as (possibly convex) counter-clockwise 2D polygons. Eg a sphere which crosses the plane needs to generate a single ellipsoid, no matter how awkward the geometry of the sphere. It needs to handle arbitrary meshes, including ones with holes/self-intersections, etc.

Finding the *points* of intersection consists of finding edges, then lerping along them; it's fairly trivial. Getting a convex hull of those can be done and then you have a "good enough" cross section for a lot of use cases. I need to actually trace them and get the polylines, though.

The next obvious algorithm (alluded to elsewhere) is to trace connecting edges. That is, if an edge crosses the plane, you look for the next edge that connects back to the first one. However, because a single triangle can have a fan of points it connects to, this will result in many discontinuities. Just put the Unity sphere at the origin and try that; it will result in 4 lines. This is good enough for mesh slicing, since you can just connect these edges back up; such discontinuites won't matter, but it's not sufficient to get an actual 2D silhouette. The disconnected paths must be joined.

My intuition says this is an NP-hard problem, since it may reduce to the longest path problem on a graph (equivalent to finding a Hamiltonian cycle)*. Those are well-studied for DAGs, but an edge graph is not that. I can construct an adjacency matrix from edges which cross the plane and maybe with enough "elbow grease performance" (bit packing, etc) I can get the brute force solution fast enough for reasonably sized meshes. (* Or maybe not, since finding simple/elementary cycles is much more tractable, eg Johnson's algorithm)

I've also found a C++ project https://github.com/intents-software/mesh-plane-intersection . That seems like it gets the short edge loops and tries to connect them, which is also suggested in this paper, although the description in the paper is quite vague:

Spoiler

After the intersection boundary was found by the al-
gorithm from Section 3.1, we can fix the mesh. We
propose a novel method for repairing the mesh by us-
ing the detected intersection boundary. For each in-
tersection segment found during the neighbor tracing,
we store the pair of triangles which formed the in-
tersection. This information is used during the repair
step to fix the connectivity and topology of the mesh.
For each triangle participating in the intersection,
a new valid polygon is created by cutting the incon-
sistent parts off (see Figure 9). To create the polygon,
we first need to identify all the polylines that inter-
sect it. We insert the polylines into an auxiliary data
structure which helps to correctly trace the polygon
as shown in Figure 10. We start with an empty list
and insert the vertices of the triangle. We then insert
the first point of each intersecting polyline into the
list at the appropriate position corresponding to the
line segment where the point is located in the orig-
inal triangle. Finally, we connect the end points of
the polylines with the vertices that come after them in
the polygon boundary. Then we trace the boundary of
the polygon: starting with any polyline, we trace the
boundary following the pointers we set up in the aux-
iliary structure until we get back to the starting point.
The points we visited form the boundary of the new
polygon. If any polyline is left unvisited, the origi-
nal triangle needs to be cut into more pieces. We re-
peat this process until all the polygons are found (Fig-
ure 9). The order of the processing does not affect the
outcome of the algorithm.

I'm going to try to port it, but it's fairly undocumented and my C++ is rusty. My fear is that it makes assumptions about the mesh layouts that may not be true when dealing with arbitrary non-manifold meshes. (EDIT: at the very least, it fails for meshes which have vertices lying directly on the plane... for example, the default unity sphere at the origin)

Everything else is closed source; eg you can do it in Maya.

Has anyone successfully done this in Unity for convex, possibly disconnected meshes?

Spoiler

Just so this post isn't completely without any Unity context, here's some code that extracts all triangles and vertex positions from a GameObject hierarchy, which is the first step in this sort of algorithm:

Code (CSharp):

1. /// <summary>
2. /// Flattens an object hierarchy and gets all the post-transformation vertices and indices for MeshRenderers
3. /// in the hierarchy.
4. /// </summary>
5. private static (float3[] vertices, int3[] triangles) extractVertexSoup(GameObject prefab)
6. {
7.     List<Mesh> uniqueMeshes = new();
8.     Dictionary<int, int> iidToUniqueIndex = new();
9.     List<(int meshIndex, float4x4 transform)> indexAndTransforms = new();
10.     foreach(GameObject go in traverse(prefab))
11.     {
12.         if(go.TryGetComponent(out MeshFilter mf))
13.         {
14.             Mesh mesh = mf.sharedMesh;
15.             if(mesh)
16.             {
17.                 int iid = mesh.GetInstanceID();
18.                 if(!iidToUniqueIndex.TryGetValue(iid, out int meshIndex))
19.                 {
20.                     meshIndex = uniqueMeshes.Count;
21.                     iidToUniqueIndex.Add(iid, meshIndex);
22.                     uniqueMeshes.Add(mesh);
23.                 }
24.                 indexAndTransforms.Add((meshIndex, go.transform.localToWorldMatrix));
25.             }
26.         }
27.     }
28.  
29.     if(uniqueMeshes.Count == 0)
30.         return (null, null);
31.  
32.     (float3[] vertices, int3[] triangles)[] untransformedMeshData = new (float3[], int3[])[uniqueMeshes.Count];
33.     using(Mesh.MeshDataArray mda = Mesh.AcquireReadOnlyMeshData(uniqueMeshes))
34.     {
35.         for(int iMesh = 0; iMesh < untransformedMeshData.Length; ++iMesh)
36.         {
37.             Mesh.MeshData md = mda[iMesh];
38.             VertexAttributeAccessor<float3> verticesIn = new(md, VertexAttribute.Position);
39.             int nVertices = verticesIn.length;
40.             float3[] vertices = new float3[nVertices];
41.             for(int iVertex = 0; iVertex < nVertices; ++iVertex)
42.                 vertices[iVertex] = verticesIn[iVertex];
43.  
44.             int3[] triangles;
45.             if(md.indexFormat == IndexFormat.UInt16)
46.             {
47.                 NativeArray<ushort> indices = md.GetIndexData<ushort>();
48.                 int nIndices = indices.Length, nTriangles = nIndices / 3;
49.                 if(nIndices % 3 != 0)
50.                     throw new Exception($"Mesh {uniqueMeshes[iMesh].name} has {nIndices} indices, which is not a multiple of 3. Only triangular meshes are supported");
51.                 triangles = new int3[nTriangles];
52.                 for(int iTriangle = 0; iTriangle < nTriangles; ++iTriangle)
53.                     triangles[iTriangle] = new(indices[iTriangle * 3], indices[iTriangle * 3 + 1], indices[iTriangle * 3 + 2]);
54.             }
55.             else
56.             {
57.                 NativeArray<int> indices = md.GetIndexData<int>();
58.                 int nIndices = indices.Length, nTriangles = nIndices / 3;
59.                 if(nIndices % 3 != 0)
60.                     throw new Exception($"Mesh {uniqueMeshes[iMesh].name} has {nIndices} indices, which is not a multiple of 3. Only triangular meshes are supported");
61.                 triangles = new int3[nTriangles];
62.                 for(int iTriangle = 0; iTriangle < nTriangles; ++iTriangle)
63.                     triangles[iTriangle] = new(indices[iTriangle * 3], indices[iTriangle * 3 + 1], indices[iTriangle * 3 + 2]);
64.             }
65.  
66.             untransformedMeshData[iMesh] = (vertices, triangles);
67.         }
68.     }
69.  
70.     int nInputs = indexAndTransforms.Count;
71.     int totalVertices = 0, totalTriangles = 0;
72.     for(int iInput = 0; iInput < nInputs; ++iInput)
73.     {
74.         (float3[] vertices, int3[] triangles) = untransformedMeshData[indexAndTransforms[iInput].meshIndex];
75.         totalVertices += vertices.Length;
76.         totalTriangles += triangles.Length;
77.     }
78.  
79.     float3[] verticesOut = new float3[totalVertices];
80.     int3[] trianglesOut = new int3[totalTriangles];
81.     totalVertices = 0;
82.     totalTriangles = 0;
83.     for(int iInput = 0; iInput < nInputs; ++iInput)
84.     {
85.         (int meshIndex, float4x4 transform) = indexAndTransforms[iInput];
86.         (float3[] vertices, int3[] triangles) = untransformedMeshData[meshIndex];
87.         int3 indexOffset = new(totalVertices, totalVertices, totalVertices);
88.         for(int iVertex = 0; iVertex < vertices.Length; ++iVertex)
89.             verticesOut[totalVertices + iVertex] = math.transform(transform, vertices[iVertex]);
90.         for(int iTriangle = 0; iTriangle < triangles.Length; ++iTriangle)
91.             trianglesOut[totalTriangles + iTriangle] = triangles[iTriangle] + indexOffset;
92.         totalVertices += vertices.Length;
93.         totalTriangles += triangles.Length;
94.     }
95.  
96.     if(verticesOut.Length < 3 || trianglesOut.Length < 1)
97.         return (null, null);
98.  
99.     return (verticesOut, trianglesOut);
100.  
101.     static IEnumerable<GameObject> traverse(GameObject go)
102.     {
103.         yield return go;
104.         Transform t = go.transform;
105.         int nChildren = t.childCount;
106.         for(int i = 0 ; i < nChildren; ++i)
107.             foreach(GameObject c in traverse(t.GetChild(i).gameObject))
108.                 yield return c;
109.     }
110. }
111.  
112. /// <summary>
113. /// Merges vertices that are identical (or within a small distance), quantizes them all to an integer grid,
114. /// and modifies the input triangles array to point to the new vertices. This is needed because meshes often
115. /// have 2 vertices at the same position with different normals or UVs. It also helps to have everything on
116. /// an integer grid to avoid any floating point errors or instability.
117. /// </summary>
118. private static int3[] quantizeAndMergeCoincident(float3[] verticesIn, int3[] triangles)
119. {
120.     List<int3> verticesOut = new(verticesIn.Length);
121.     Dictionary<int3, int> remap = new(verticesIn.Length);
122.     int[] oldToNew = new int[verticesIn.Length];
123.     for(int iOld = 0; iOld < verticesIn.Length; ++iOld)
124.     {
125.         float3 p = verticesIn[iOld];
126.         int3 q = (int3) math.round(p * QUANTIZE_FACTOR);
127.         if(!remap.TryGetValue(q, out int iNew))
128.         {
129.             iNew = verticesOut.Count;
130.             verticesOut.Add(q);
131.             remap.Add(q, iNew);
132.         }
133.         oldToNew[iOld] = iNew;
134.     }
135.  
136.     for(int iTriangle = 0; iTriangle < triangles.Length; ++iTriangle)
137.     {
138.         static int lookupOrClampInvalid(int i, int[] r) => i >= 0 && i < r.Length ? r[i] : 0;
139.         int3 t = triangles[iTriangle];
140.         t.x = lookupOrClampInvalid(t.x, oldToNew);
141.         t.y = lookupOrClampInvalid(t.y, oldToNew);
142.         t.z = lookupOrClampInvalid(t.z, oldToNew);
143.         triangles[iTriangle] = t;
144.     }
145.  
146.     return verticesOut.ToArray();
147. }
148.  

 

So here's a more practical example of a tricky case. Here's a stool (positioned slightly inwards):



And here's some edge loops from my attempt at porting the above C++ library:



Which isn't a loop at all.

I plan to eventually use Clipper to combine polygons, however those polygons must be closed if I am to combine them. I know this problem may be intractable for all possible meshes, but I'd like to cover as many possible cases as I can.

Note that I am already quantizing all vertices to an integer grid, and combining coincident vertices.

 

So, you want the two black paths, did I understand that correctly?


Are all lines connected properly? Do you have a vertex where two edges touch?
E.g. is this line made up of 4 segments?


If yes, you could remove all interior line segments to get what you want. You'd "just" have to find out if a line segment is an interior segment.

 

Eventually, yes, that would be what I want for this case (a polygon with a hole). But a more realistic decomposition, which I could then pass into clipper to combine them for the final polygon would be several closed loops that could be unioned. So in the above example, I'd need to combine the yellow and blue lines and the pink line, like so, while leaving all the other paths. Then they can be unioned for a full polygon, triangulated, etc, etc.



But of course I need to do this programmatically, and with arbitrarily complex meshes.

Not yet. In fact, I think the chair mesh is just a union of cubes, so there aren't vertices on both sides at the point the lines touch.

An option would be to identify vertices that fall directly on another edge, and then insert new vertices at the points they meet, thereby splitting those edges. That might help if doing the "longest cycle" approach. If I did that, then (along with the vertex quantization to int grid thing), I think I *could* make that a guarantee. The question is if that would provide a stepping stone to a general solution.

That's kinda putting the cart before the horse, though. I don't have a closed polygon yet, I'm trying to construct one. So there's no such thing as an "interior segment", just a lot of lines.

 

This is basically constructive solid geometry (CSG) what you are asking for. That's not trivial, indeed.

One way to approach this would be by converting the geometry into a BSP tree, do your boolean operations and then convert the BSP tree back to polygons. There are libraries that can do this like CGAL (the most unintuitive piece of software I have ever used which I do NOT recommend using). Writing something like this is very challenging.

Do you need polygons as output or would it be enough to see the result on the screen? If the latter is enough, there is an easy algorithm for screen-space CSG which uses the depth buffer. Signed distance fields might also be an option.

 

That's a possibility, although a rather heavyweight one. I feel like this *should* be possible with just the tri-mesh. But maybe that's the best way to handle all those weird cases like self-intersection, degenerate geometry, etc -- just convert it into a format that doesn't allow those things.

Sadly, this is for generating geometry, not just visuals. My goal is to come up with a "good enough" set of polygons that I can hand them to Clipper, which will do the polygon boolean operations, and then from there into Triangle.NET to get the triangulation.

Thanks for the help, though. Your idea of splitting edges where another vertex meets them is likely going to be a requirement for any graph-based approach.