Polygon triangulation supports collinear vertices

CHANGELOG.md

@@ -10,6 +10,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Added
 - Users can manually specify the background medium for a structure to be used for geometry gradient calculations by supplying `Structure.background_permittivity`. This is useful when there are overlapping structures or structures embedded in other mediums.
 
+### Changed
+- Triangulation of `PolySlab` polygons now supports polygons with collinear vertices.
+
 ### Fixed
 - Minor gradient direction and normalization fixes for polyslab, field monitors, and diffraction monitors in autograd.
 

tests/test_components/test_geometry.py

@@ -957,3 +957,9 @@ def test_update_from_bounds():
     for geom2d in geometries:
         with pytest.raises(NotImplementedError):
             geom_update = geom2d._update_from_bounds(bounds=new_bounds, axis=axis)
+
+
+def test_triangulation_with_collinear_vertices():
+    xr = np.linspace(0, 1, 6)
+    a = np.array([[x, -0.5] for x in xr] + [[x, 0.5] for x in xr[::-1]])
+    assert len(td.components.geometry.triangulation.triangulate(a)) == 10

tidy3d/components/geometry/triangulation.py

@@ -4,6 +4,7 @@
 import numpy as np
 import shapely
 
+from ...exceptions import Tidy3dError
 from ..types import ArrayFloat1D, ArrayFloat2D
 
 
@@ -17,8 +18,8 @@ class Vertex:
         Vertex coordinate.
     index : int
         Vertex index in the original polygon.
-    is_convex : bool = False
-        Flag indicating whether this is a convex vertex in the polygon.
+    convexity : float = 0.0
+        Value representing the convexity (> 0) or concavity (< 0) of the vertex in the polygon.
     is_ear : bool = False
         Flag indicating whether this is an ear of the polygon.
     """
@@ -27,12 +28,12 @@ class Vertex:
 
     index: int
 
-    is_convex: bool = False
+    convexity: float = 0.0
 
     is_ear: bool = False
 
 
-def update_convexity(vertices: List[Vertex], i: int) -> None:
+def update_convexity(vertices: List[Vertex], i: int) -> int:
     """Update the convexity of a vertex in a polygon.
 
     Parameters
@@ -41,15 +42,21 @@ def update_convexity(vertices: List[Vertex], i: int) -> None:
         Vertices of the polygon.
     i : int
         Index of the vertex to be updated.
+
+    Returns
+    -------
+    int
+        -1 if vertex was collinear, +1 if it became collinear, 0 otherwise.
     """
+    result = -1 if vertices[i].convexity == 0.0 else 0
     j = (i + 1) % len(vertices)
-    vertices[i].is_convex = (
-        np.cross(
-            vertices[i].coordinate - vertices[i - 1].coordinate,
-            vertices[j].coordinate - vertices[i].coordinate,
-        )
-        > 0
+    vertices[i].convexity = np.cross(
+        vertices[i].coordinate - vertices[i - 1].coordinate,
+        vertices[j].coordinate - vertices[i].coordinate,
     )
+    if vertices[i].convexity == 0.0:
+        result += 1
+    return result
 
 
 def is_inside(
@@ -87,10 +94,10 @@ def update_ear_flag(vertices: List[Vertex], i: int) -> None:
     h = (i - 1) % len(vertices)
     j = (i + 1) % len(vertices)
     triangle = (vertices[h].coordinate, vertices[i].coordinate, vertices[j].coordinate)
-    vertices[i].is_ear = vertices[i].is_convex and not any(
+    vertices[i].is_ear = vertices[i].convexity > 0 and not any(
         is_inside(v.coordinate, triangle)
         for k, v in enumerate(vertices)
-        if not (v.is_convex or k == h or k == i or k == j)
+        if not (v.convexity > 0 or k == h or k == i or k == j)
     )
 
 
@@ -114,31 +121,45 @@ def triangulate(vertices: ArrayFloat2D) -> List[Tuple[int, int, int]]:
     if not is_ccw:
         vertices.reverse()
 
+    collinears = 0
     for i in range(len(vertices)):
         update_convexity(vertices, i)
+        if vertices[i].convexity == 0.0:
+            collinears += 1
 
     for i in range(len(vertices)):
         update_ear_flag(vertices, i)
 
     triangles = []
 
+    ear_found = True
     while len(vertices) > 3:
+        if not ear_found:
+            raise Tidy3dError(
+                "Impossible to triangulate polygon. Verify that the polygon is valid."
+            )
+        ear_found = False
         i = 0
         while i < len(vertices):
             if vertices[i].is_ear:
-                j = (i + 1) % len(vertices)
-                triangles.append((vertices[i - 1].index, vertices[i].index, vertices[j].index))
-                vertices.pop(i)
-                if len(vertices) == 3:
-                    break
+                removed = vertices.pop(i)
                 h = (i - 1) % len(vertices)
                 j = i % len(vertices)
-                if not vertices[h].is_convex:
-                    update_convexity(vertices, h)
-                if not vertices[j].is_convex:
-                    update_convexity(vertices, j)
-                update_ear_flag(vertices, h)
-                update_ear_flag(vertices, j)
+                collinears += update_convexity(vertices, h)
+                collinears += update_convexity(vertices, j)
+                if collinears == len(vertices):
+                    # Undo removal because only collinear vertices remain
+                    vertices.insert(i, removed)
+                    collinears += update_convexity(vertices, (i - 1) % len(vertices))
+                    collinears += update_convexity(vertices, (i + 1) % len(vertices))
+                    i += 1
+                else:
+                    ear_found = True
+                    triangles.append((vertices[h].index, removed.index, vertices[j].index))
+                    update_ear_flag(vertices, h)
+                    update_ear_flag(vertices, j)
+                    if len(vertices) == 3:
+                        break
             else:
                 i += 1