I faced a problem of finding the material side of a polygon when I was working on slicing algorithms for 3D printing. This is not something new and all the slicing software need to know which side of the polygon is empty and which side is solid. An easy solution is to get the normal vector of triangles from the surface mesh (.stl). This is part of standard workflow in 3D printing. A STL file already has the normal direction of each triangles and the opposite of normal direction is the material side which needs to be filled. But I need something more general. Not every STL file is valid or sometimes I don’t even have a STL file to start with. So my problem statement can be rephrased to:
Given any set of valid polygons (i.e., closed loop polygon with no intersections and no non-manifold geometries.), how to determine which polygon is interior and exterior? Even further, how to know which polygons are inside and outside of another polygon?
In short, I need to output a data structure which unveils the relationships between each polygon and intuitively tells us which polygon is the outmost and innermost. So, simply getting the normal directions does bot help here. Not to mention I do not even have normals to start with. Based on my research, I put together a algorithm.
- Slice a surface mesh data with a plane to get undirected graph $G = (V, E)$. This optional if you already have $G$ to start with.
- Use union-find data structure to find disjoint sets $V^d = V_1^d, V_2^d, …, V_n^d$ where $V_i^d$ is a disjoint set of vertices.
- Each $V_i^d$ represents a closed polygon. If they are not already sorted, we can sort them in clockwise order and rewrite them to a graph which represents a closed polygon $P^d = P_1^d, P_2^d, …, P_n^d$.
-
Initialize an empty $n$ by $n$ matrix $M$ which reveals the relationships between each $P_i^d$. If $P_i^d$ encapsulates $P_j^d$, then $M[i][j] = 1$. If $P_i^d$ is inside of $P_j^d$, then $M[i][j] = -1$. All other cases will be $M[i][j] = 0$, meaning $P_i^d$ and $P_j^d$ have no relation. Pseudo code
1 2 3 4 5 6 7 8 9 10 11 12 13 14
Initialize M[n][n] for i = 0 to n - 1 for j = 0 to n - 1 if (i == j) M[i][j] = -1 continue if (M[i][j] == -1) continue if (polygon i is inside polygon j) M[i][j] = 1 M[j][i] = -1 else M[j][i] = 0 return M[n][n] - Build a $n$ by 3 matrix $A$ where $A[i][0]$ is the index, $A[i][i]$ is an integer which indicates whether $P_i^d$ is interior (-1) or exterior (1), and $A[i][2]$ is the number of polygons encapsulated. $A[i][1]$ can be calculated by multiplying all non-zero terms in $M[i][.]$. $A[i][2]$ can be calculated by adding all non-negative terms in $M[i][.]$.
- Sort $A$ by the number of polygons encapsulated, i.e. $A[i][2]$, in a descending order to get a new sorted matrix $A^s$.
Now the problem is solved. Really the key here is to build the matrix $M$. The $i$ row of $M$ represents the polygon $P_i^d$. Multiplying all the non-zero terms at row $i$ gives either 1/exterior or -1/interior. Adding all the non-negative terms at row $i$ gives us the number of polygons are inside the polygon $P_i^d$. Hence, we know the topological order of all the polygons. The first polygon in the sorted matrix $A^s$ is the outermost polygon and the last polygon is the innermost polygon.
This is a pretty straight forward algorithm. The only thing I abstract away is the part which tests whether if a point is inside polygon or not. The method is to draw a ray starting from the point and detect how many times the ray intersects with the polygon. If there is only one intersection, the point is inside the polygon. If there are two intersections, the point is outside of the polygon. Its implementation is not as easy as it sounds and there are also some edge cases to deal with. But it is not the scope of this article. The BOOST library already has it implemented. I am happy to walk away using BOOST library. Finally, an example is given in this repo.