Add option to deduplicate edges in edge triangulation

src/PartSegCore_compiled_backend/triangulate.pyx

@@ -13,8 +13,6 @@ from libcpp cimport bool
 from libcpp.unordered_set cimport unordered_set
 from libcpp.utility cimport pair
 from libcpp.vector cimport vector
-from libcpp.unordered_map cimport unordered_map
-from libcpp.algorithm cimport sort
 
 
 cdef extern from "triangulation/point.hpp" namespace "partsegcore::point":
@@ -96,6 +94,7 @@ cdef extern from "triangulation/triangulate.hpp" namespace "partsegcore::triangu
     pair[vector[Triangle], vector[Point]] triangulate_polygon_face(const vector[Point]& polygon) except + nogil
     pair[vector[Triangle], vector[Point]] triangulate_polygon_face(const vector[vector[Point]]& polygon_list) except + nogil
     PathTriangulation triangulate_path_edge(const vector[Point]& path, bool closed, float limit, bool bevel) except + nogil
+    vector[vector[Point]] split_polygon_on_repeated_edges(const vector[Point]& polygon) except + nogil
 
 
 ctypedef fused float_types:
@@ -429,14 +428,15 @@ def triangulate_path_edge_numpy(cnp.ndarray[cnp.float32_t, ndim=2] path, bool cl
     )
 
 
-def triangulate_polygon_with_edge_numpy_li(polygon_li: list[np.ndarray]) -> tuple[tuple[np.ndarray, np.ndarray], tuple[np.ndarray, np.ndarray, np.ndarray]]:
+def triangulate_polygon_with_edge_numpy_li(polygon_li: list[np.ndarray], split_edges: bool=False) -> tuple[tuple[np.ndarray, np.ndarray], tuple[np.ndarray, np.ndarray, np.ndarray]]:
     """ Triangulate polygon"""
     cdef vector[Point] polygon_vector
-    cdef vector[vector[Point]] polygon_vector_list
+    cdef vector[vector[Point]] polygon_vector_list, edge_split_list
     cdef Point p1, p2
     cdef pair[vector[Triangle], vector[Point]] triangulation_result
     cdef vector[PathTriangulation] edge_result
     cdef cnp.ndarray[cnp.uint32_t, ndim=2] triangles, edge_triangles
+    cdef size_t triangle_count = 0
 
     cdef cnp.ndarray[cnp.float32_t, ndim=2] points, edge_offsets, edges_centers, polygon
     cdef size_t i, j, len_path, edge_triangle_count, edge_center_count, edge_triangle_index, edge_center_index
@@ -457,8 +457,14 @@ def triangulate_polygon_with_edge_numpy_li(polygon_li: list[np.ndarray]) -> tupl
         if polygon_vector.size() > 1 and polygon_vector.front() == polygon_vector.back():
             polygon_vector.pop_back()
         polygon_vector_list.push_back(polygon_vector)
-        with cython.nogil:
-            edge_result.push_back(triangulate_path_edge(polygon_vector, True, 3.0, False))
+        if split_edges:
+            with cython.nogil:
+                edge_split_list = split_polygon_on_repeated_edges(polygon_vector)
+                for edge_li in edge_split_list:
+                    edge_result.push_back(triangulate_path_edge(edge_li, True, 3.0, False))
+        else:
+            with cython.nogil:
+                edge_result.push_back(triangulate_path_edge(polygon_vector, True, 3.0, False))
 
     with cython.nogil:
         triangulation_result = triangulate_polygon_face(polygon_vector_list)
@@ -470,7 +476,7 @@ def triangulate_polygon_with_edge_numpy_li(polygon_li: list[np.ndarray]) -> tupl
         triangles[i, 1] = triangulation_result.first[i].y
         triangles[i, 2] = triangulation_result.first[i].z
 
-    points =np.empty((triangulation_result.second.size(), 2), dtype=np.float32)
+    points = np.empty((triangulation_result.second.size(), 2), dtype=np.float32)
     for i in range(triangulation_result.second.size()):
         points[i, 0] = triangulation_result.second[i].x
         points[i, 1] = triangulation_result.second[i].y
@@ -489,10 +495,11 @@ def triangulate_polygon_with_edge_numpy_li(polygon_li: list[np.ndarray]) -> tupl
     edge_center_index = 0
     for i in range(edge_result.size()):
         for j in range(edge_result[i].triangles.size()):
-            edge_triangles[edge_triangle_index, 0] = edge_result[i].triangles[j].x
-            edge_triangles[edge_triangle_index, 1] = edge_result[i].triangles[j].y
-            edge_triangles[edge_triangle_index, 2] = edge_result[i].triangles[j].z
+            edge_triangles[edge_triangle_index, 0] = edge_result[i].triangles[j].x + triangle_count
+            edge_triangles[edge_triangle_index, 1] = edge_result[i].triangles[j].y + triangle_count
+            edge_triangles[edge_triangle_index, 2] = edge_result[i].triangles[j].z + triangle_count
             edge_triangle_index += 1
+        triangle_count += edge_result[i].centers.size()
 
         for j in range(edge_result[i].centers.size()):
             edges_centers[edge_center_index, 0] = edge_result[i].centers[j].x
@@ -511,3 +518,20 @@ def triangulate_polygon_with_edge_numpy_li(polygon_li: list[np.ndarray]) -> tupl
         edge_triangles,
     )
     )
+
+
+def split_polygon_on_repeated_edges_py(polygon: Sequence[Sequence[float]]) -> list[list[tuple[float, float]]]:
+    """ Split polygon on repeated edges"""
+    cdef vector[Point] polygon_vector
+    cdef vector[vector[Point]] result
+    cdef Point p1, p2
+
+    polygon_vector.reserve(len(polygon))
+    for point in polygon:
+        polygon_vector.push_back(Point(point[0], point[1]))
+
+    result = split_polygon_on_repeated_edges(polygon_vector)
+    return [
+        [(point.x, point.y) for point in polygon_vector]
+        for polygon_vector in result
+    ]

src/PartSegCore_compiled_backend/triangulation/triangulate.hpp

@@ -2,9 +2,8 @@
 #define PARTSEGCORE_TRIANGULATE_H
 
 #include <algorithm>
-#include <cmath>
 #include <map>
-#include <memory>
+#include <memory>  // memory header is required on linux, and not on macos
 #include <set>
 #include <sstream>
 #include <unordered_map>
@@ -1288,6 +1287,87 @@ inline PathTriangulation triangulate_path_edge(
   return result;
 }
 
+/**
+ * Represents an edge in a graph structure used for polygon processing.
+ * Each edge contains a reference to its opposite point and a flag to track
+ * if it has been visited during graph traversal.
+ */
+struct GraphEdge {
+  point::Point opposite_point;
+  bool visited;
+  explicit GraphEdge(point::Point p) : opposite_point(p), visited(false) {}
+};
+
+/**
+ * Represents a node in a graph structure used for polygon processing.
+ * Each node contains its edges, a sub-index for traversal tracking,
+ * and a visited flag for graph traversal.
+ */
+struct GraphNode {
+  std::vector<GraphEdge> edges;
+  std::size_t sub_index;
+  bool visited;
+
+  GraphNode() : sub_index(0), visited(false) {}
+};
+
+/**
+ * Splits a polygon into sub-polygons by identifying and removing edges that
+ * appear more than once in the polygon's edge list.
+ *
+ * This function processes the given polygon and separates it wherever an
+ * edge is repeated. It generates a collection of sub-polygons such that each
+ * resulting sub-polygon contains unique edges. This operation can help to
+ * resolve ambiguities in complex or self-intersecting polygons.
+ *
+ * @param polygon The input polygon represented as a list of edges.
+ *
+ * @return A vector of sub-polygons, where each sub-polygon is free of repeated
+ * edges.
+ */
+inline std::vector<std::vector<point::Point>> split_polygon_on_repeated_edges(
+    const std::vector<point::Point> &polygon) {
+  auto edges_dedup = calc_dedup_edges({polygon});
+  std::vector<std::vector<point::Point>> result;
+  point::Segment segment;
+
+  std::unordered_set edges_set(edges_dedup.begin(), edges_dedup.end());
+  std::unordered_map<point::Point, GraphNode> edges_map;
+  for (std::size_t i = 0; i < polygon.size() - 1; i++) {
+    segment = {polygon[i], polygon[(i + 1)]};
+    if (edges_set.count(segment) > 0) {
+      edges_map[polygon[i]].edges.emplace_back(polygon[i + 1]);
+    }
+  }
+  segment = {polygon.back(), polygon.front()};
+  if (edges_set.count(segment) > 0) {
+    edges_map[polygon.back()].edges.emplace_back(polygon.front());
+  }
+  for (auto &edge : edges_map) {
+    if (edge.second.visited) {
+      continue;
+    }
+    edge.second.visited = true;
+    std::vector<point::Point> new_polygon;
+    new_polygon.push_back(edge.first);
+    auto *current_edge = &edge.second;
+    while (current_edge->sub_index < current_edge->edges.size()) {
+      auto *prev = current_edge;
+      auto next_point =
+          current_edge->edges[current_edge->sub_index].opposite_point;
+      current_edge = &edges_map.at(next_point);
+      prev->sub_index++;
+      current_edge->visited = true;
+      new_polygon.push_back(next_point);
+    }
+    while (new_polygon.front() == new_polygon.back()) {
+      new_polygon.pop_back();
+    }
+    result.push_back(new_polygon);
+  }
+  return result;
+}
+
 }  // namespace partsegcore::triangulation
 
 #endif  // PARTSEGCORE_TRIANGULATE_H

src/tests/test_triangulate.py

@@ -10,6 +10,7 @@
     on_segment,
     orientation,
     segment_left_to_right_comparator,
+    split_polygon_on_repeated_edges_py,
     triangle_convex_polygon,
     triangulate_monotone_polygon_py,
     triangulate_path_edge_numpy,
@@ -661,3 +662,32 @@ def test_triangulate_polygon_with_edge_numpy_li(polygon, expected):
     triangles_ = _renumerate_triangles(polygon, points, triangles)
     assert triangles_ == expected
     assert centers.shape == offsets.shape
+
+
+def test_split_polygon_on_repeated_edges_py_no_split():
+    res = split_polygon_on_repeated_edges_py([[0, 0], [0, 1], [1, 1], [1, 0]])
+    assert len(res) == 1
+    assert len(res[0]) == 4
+    idx = res[0].index((0, 0))
+    assert res[0][idx:] + res[0][:idx] == [(0, 0), (0, 1), (1, 1), (1, 0)]
+
+
+def test_split_polygon_on_repeated_edges_py_square_in_square():
+    res = split_polygon_on_repeated_edges_py(
+        [[0, 0], [0, 5], [1, 5], [1, 1], [9, 1], [9, 9], [1, 9], [1, 5], [0, 5], [0, 10], [10, 10], [10, 0]]
+    )
+    assert len(res) == 2
+    assert len(res[0]) == 5
+    assert len(res[1]) == 5
+    # idx = res[0].index((0, 0))
+    # assert res[0][idx:] + res[0][:idx] == [(0, 0), (0, 1), (1, 1), (1, 0)]
+
+
+@pytest.mark.parametrize(('split_edges', 'triangles'), [(True, 20), (False, 24)])
+def test_splitting_edges(split_edges, triangles):
+    polygon = np.array(
+        [[0, 0], [0, 5], [1, 5], [1, 1], [9, 1], [9, 9], [1, 9], [1, 5], [0, 5], [0, 10], [10, 10], [10, 0]],
+        dtype=np.float32,
+    )
+    triangles_ = triangulate_polygon_with_edge_numpy_li([polygon], split_edges=split_edges)[1][2]
+    assert len(triangles_) == triangles