triangulate cuts polygons into ears, not fans

changelog.d/triangulate-ears.md

@@ -0,0 +1,8 @@
+### Fixed
+
+- `triangulate` cut a polygon into a fan from its first corner, which covers a
+  concave polygon twice in part and leaves part uncovered (a U-shaped polygon of
+  area 6 gave triangles of area 9). It now cuts ears in the polygon's plane, as
+  VTK's `vtkTriangleFilter`: `n - 2` triangles that cover the polygon once, each
+  wound as it is. Found running VTK's Python tests on Tack (the concave polygons of
+  `vtkDiscreteFlyingEdges2D` and `vtkImageToPolyDataFilter`).

packages/tack-vis/src/tack/dataset/algorithms/_refine.py

@@ -150,16 +150,109 @@ def _fan_counts(loop_offsets, counts):
         counts[c] = max(loop_offsets[c + 1] - loop_offsets[c] - 2, 0)
 
 
+@tack.func
+def _plane_uv(p, drop):
+    """``p``'s coordinates in the plane that drops axis ``drop``, the other two in
+    cyclic order (so the plane's normal is the dropped axis, sign kept)."""
+    u = p[1]
+    v = p[2]
+    if drop == 1:
+        u = p[2]
+        v = p[0]
+    elif drop == 2:
+        u = p[0]
+        v = p[1]
+    return u, v
+
+
+@tack.func
+def _turn(x, loop_points, first, a, b, c, drop):
+    """Twice the signed area of the triangle of loop corners ``a``, ``b``, ``c`` in the
+    plane (``_plane_uv``)."""
+    au, av = _plane_uv(x[loop_points[first + a]], drop)
+    bu, bv = _plane_uv(x[loop_points[first + b]], drop)
+    cu, cv = _plane_uv(x[loop_points[first + c]], drop)
+    return (bu - au) * (cv - av) - (bv - av) * (cu - au)
+
+
 @tack.kernel
-def _fans(loop_offsets, loop_points, starts, connectivity, sources):
+def _ears(loop_offsets, loop_points, x, starts, after, before, connectivity, sources):
+    # Each polygon by ear cutting in its plane -- the plane of its Newell normal's
+    # largest component -- so a concave one is covered once (a fan from its first
+    # corner covers part of it twice and part not at all). An ear is a corner
+    # turning as the polygon does whose triangle holds no other corner; cutting it
+    # leaves a polygon of one corner fewer. ``after`` and ``before`` link the
+    # corners left, within the polygon's own range. Where no ear is left (a
+    # degenerate remainder) the rest is a fan, so every polygon of n corners gives
+    # n - 2 triangles, each wound as the polygon.
     for c in range(len(loop_offsets) - 1):
         first = loop_offsets[c]
+        n = loop_offsets[c + 1] - first
         at = starts[c]
-        for k in range(loop_offsets[c + 1] - first - 2):
-            connectivity[3 * (at + k)] = loop_points[first]
-            connectivity[3 * (at + k) + 1] = loop_points[first + k + 1]
-            connectivity[3 * (at + k) + 2] = loop_points[first + k + 2]
-            sources[at + k] = c
+        if n >= 3:
+            nx = 0.0
+            ny = 0.0
+            nz = 0.0
+            for k in range(n):
+                p = x[loop_points[first + k]]
+                q = x[loop_points[first + (k + 1) % n]]
+                nx += (p[1] - q[1]) * (p[2] + q[2])
+                ny += (p[2] - q[2]) * (p[0] + q[0])
+                nz += (p[0] - q[0]) * (p[1] + q[1])
+            drop = 2
+            sense = nz
+            if abs(nx) >= abs(ny) and abs(nx) >= abs(nz):
+                drop = 0
+                sense = nx
+            elif abs(ny) >= abs(nz):
+                drop = 1
+                sense = ny
+            sign = 1.0
+            if sense < 0.0:
+                sign = -1.0
+            for k in range(n):
+                after[first + k] = (k + 1) % n
+                before[first + k] = (k + n - 1) % n
+            left = n
+            i = 0
+            misses = 0
+            t = 0
+            while left > 3 and misses < left:
+                a = before[first + i]
+                b = after[first + i]
+                ear = sign * _turn(x, loop_points, first, a, i, b, drop) > 0.0
+                if ear:
+                    j = after[first + b]
+                    while ear and j != a:
+                        inside = (sign * _turn(x, loop_points, first, a, i, j, drop) > 0.0
+                                  and sign * _turn(x, loop_points, first, i, b, j, drop) > 0.0
+                                  and sign * _turn(x, loop_points, first, b, a, j, drop) > 0.0)
+                        if inside:
+                            ear = False
+                        j = after[first + j]
+                if ear:
+                    connectivity[3 * (at + t)] = loop_points[first + a]
+                    connectivity[3 * (at + t) + 1] = loop_points[first + i]
+                    connectivity[3 * (at + t) + 2] = loop_points[first + b]
+                    sources[at + t] = c
+                    t += 1
+                    after[first + a] = b
+                    before[first + b] = a
+                    left -= 1
+                    i = a
+                    misses = 0
+                else:
+                    i = b
+                    misses += 1
+            # The last triangle, or a fan of a degenerate remainder from corner i.
+            j = after[first + i]
+            while t < n - 2:
+                connectivity[3 * (at + t)] = loop_points[first + i]
+                connectivity[3 * (at + t) + 1] = loop_points[first + j]
+                connectivity[3 * (at + t) + 2] = loop_points[first + after[first + j]]
+                sources[at + t] = c
+                t += 1
+                j = after[first + j]
 
 
 @tack.kernel

packages/tack-vis/src/tack/dataset/algorithms/triangulate.py

@@ -1,23 +1,24 @@
 # SPDX-FileCopyrightText: Copyright (c) Kitware Inc.
 # SPDX-License-Identifier: BSD-3-Clause
 
-"""``triangulate``: every 2D cell as triangles: a quad or pixel in 2, a polygon in a fan from its first point; other cells are left out. The points are the same; cell fields go to each cell's triangles."""
+"""``triangulate``: every 2D cell as triangles: a quad or pixel in 2, a polygon by ear cutting; other cells are left out. The points are the same; cell fields go to each cell's triangles."""
 
 import tack
 from tack.algorithms.scan import exclusive_scan
-from tack.arrays import size_of
+from tack.arrays import materialize, size_of
 from tack.dataset import ids, shapes
-from tack.dataset.algorithms._refine import _fan_counts, _fans, _pieces_topology, _split
+from tack.dataset.algorithms._refine import _ears, _fan_counts, _pieces_topology, _split
 from tack.runtime.dispatch import on_inputs_device
 
 __all__ = ['triangulate']
 
 
 @on_inputs_device
 def triangulate(data, fields=None):
-    """Every 2D cell as triangles: a quad or pixel in 2, a polygon in a fan from its
-    first point; other cells are left out. The points are the same; cell fields
-    go to each cell's triangles."""
+    """Every 2D cell as triangles: a quad or pixel in 2, a polygon of ``n`` corners in
+    ``n - 2`` by ear cutting in its plane (as VTK's vtkTriangleFilter), so a concave
+    one is covered once, each triangle wound as the polygon; other cells are left
+    out. The points are the same; cell fields go to each cell's triangles."""
     from tack.dataset.ghosts import visible
 
     data = visible(data)
@@ -37,5 +38,9 @@ def triangulate(data, fields=None):
     connectivity = tack.empty((3 * count,), odt)
     sources = tack.empty((count,), data.id_dtype)
     if count:
-        _fans(loop_offsets, loop_points, starts, connectivity, sources)
+        # The corners left of each polygon as it is cut, linked within its range.
+        after = tack.empty((size_of(loop_points),), tack.i32)
+        before = tack.empty((size_of(loop_points),), tack.i32)
+        _ears(loop_offsets, loop_points, materialize(data.geometry.values), starts, after,
+              before, connectivity, sources)
     return _pieces_topology(data, connectivity, sources, count, shapes.TRIANGLE, 3, fields)

packages/tack-vis/tests/test_refinement_and_sources.py

@@ -153,8 +153,45 @@ def test_triangulate_keeps_the_area(backend, cell_type, per_cell):
                                                                    rel=1e-5)
 
 
+#: A U, concave: a fan from its first corner covers 9 of its area 6.
+_U = [(0, 0), (4, 0), (4, 2), (3, 2), (3, 1), (1, 1), (1, 2), (0, 2)]
+
+
+def _polygons(loops):
+    """Polygons of the given corner loops, ``(x, y, z)`` each, as one dataset."""
+    points = np.array([p for loop in loops for p in loop], np.float32)
+    offsets = np.cumsum([0] + [len(loop) for loop in loops]).astype(np.int64)
+    topology = td.PolygonalTopology(tack.array(offsets), tack.array(np.arange(
+        len(points), dtype=np.int64)), num_points=len(points))
+    return td.DataSet(topology, points)
+
+
+@pytest.mark.parametrize("loop", [
+    [(x, y, 0.0) for x, y in _U],                     # in xy
+    [(x, y, 0.0) for x, y in _U[::-1]],               # wound the other way
+    [(0.0, x, y) for x, y in _U],                     # in yz
+    [(x, y, 0.25 * x) for x, y in _U],                # in a tilted plane
+    [(0, 0, 0), (1, 0, 0), (2, 0, 0), (2, 2, 0), (0, 2, 0)],    # three corners in line
+])
+def test_triangulate_cuts_concave_polygons_into_ears(backend, loop):
+    """VTK's vtkTriangleFilter's answer: n - 2 triangles covering the polygon once,
+    each wound as it is, where a fan from the first corner covered a concave one
+    twice in part (the U's 9 for its 6)."""
+    data = _polygons([loop, [(5, 0, 0), (6, 0, 0), (6, 1, 0), (5, 1, 0)]])
+    out = td.triangulate(data)
+    assert out.num_cells == (len(loop) - 2) + 2
+    area = td.cell_geometry(out).measure.values.numpy()
+    exact = td.cell_geometry(data).measure.values.numpy()
+    np.testing.assert_allclose([area[:len(loop) - 2].sum(), area[len(loop) - 2:].sum()],
+                               exact, rtol=1e-6)
+    normals = td.surface_normals(out).fields["cell_normals"].values.numpy()
+    polygon_normal = td.surface_normals(data).fields["cell_normals"].values.numpy()[0]
+    wound = normals[:len(loop) - 2][area[:len(loop) - 2] > 1e-9]
+    np.testing.assert_allclose(wound, np.tile(polygon_normal, (len(wound), 1)), atol=1e-6)
+
+
 @needs_vtk
-def test_triangulate_polygons_as_fans(backend):
+def test_triangulate_polygons_as_ears(backend):
     data = td.as_polygons(_plane())
     out = td.triangulate(data)
     loops, _ = data.topology.loops()