diff --git a/changelog.d/triangulate-ears.md b/changelog.d/triangulate-ears.md new file mode 100644 index 000000000..2d0e182ab --- /dev/null +++ b/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`). diff --git a/packages/tack-vis/src/tack/dataset/algorithms/_refine.py b/packages/tack-vis/src/tack/dataset/algorithms/_refine.py index b5f2969a0..c6c4628c4 100644 --- a/packages/tack-vis/src/tack/dataset/algorithms/_refine.py +++ b/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 diff --git a/packages/tack-vis/src/tack/dataset/algorithms/triangulate.py b/packages/tack-vis/src/tack/dataset/algorithms/triangulate.py index c0cf2adce..76ed9fdd3 100644 --- a/packages/tack-vis/src/tack/dataset/algorithms/triangulate.py +++ b/packages/tack-vis/src/tack/dataset/algorithms/triangulate.py @@ -1,13 +1,13 @@ # 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'] @@ -15,9 +15,10 @@ @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) diff --git a/packages/tack-vis/tests/test_refinement_and_sources.py b/packages/tack-vis/tests/test_refinement_and_sources.py index 7ae87f376..3f4051b32 100644 --- a/packages/tack-vis/tests/test_refinement_and_sources.py +++ b/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()