Reworked constrained Delanay triangulation of polygons / expolygons

using CGAL CDT implementation:
Removed all the sets / maps, replaced with vectors and CDT vertex
intrusive indices.
Reworked the outside / inside classification using just the CDT
"constrained edge" attributes and a single queue.

Cherry pick commit 1648ae853d6c69a1118efbc694dadeb9965154ee

src/libslic3r/Triangulation.cpp

@@ -12,71 +12,99 @@ inline void insert_points(Points& points, const Polygon &polygon) {
     points.insert(points.end(), polygon.points.begin(), polygon.points.end());
 }
 
-inline void insert_edge(Triangulation::HalfEdges &edges, uint32_t &offset, const Polygon &polygon) {
+inline void insert_edges(Triangulation::HalfEdges &edges, uint32_t &offset, const Polygon &polygon) {
     const Points &pts = polygon.points;
     for (uint32_t i = 1; i < pts.size(); ++i) {
         uint32_t i2 = i + offset;
-        edges.insert({i2 - 1, i2});
+        edges.push_back({i2 - 1, i2});
     }
     uint32_t size = static_cast<uint32_t>(pts.size());
     // add connection from first to last point
-    edges.insert({offset + size - 1, offset});
+    edges.push_back({offset + size - 1, offset});
     offset += size;
 }
 
 } // namespace Private
 
-Triangulation::Indices Triangulation::triangulate(const Points &   points,
+Triangulation::Indices Triangulation::triangulate(const Points    &points,
-                                                  const HalfEdges &half_edges,
+                                                  const HalfEdges &constrained_half_edges)
-                                                  bool allow_opposite_edge)
 {
     // IMPROVE use int point insted of float !!!
 
     // use cgal triangulation
     using K    = CGAL::Exact_predicates_inexact_constructions_kernel;
-    using Itag = CGAL::Exact_predicates_tag;
+    using Vb   = CGAL::Triangulation_vertex_base_with_info_2<uint32_t, K>;
+    using Fb   = CGAL::Constrained_triangulation_face_base_2<K>;
+    using Tds  = CGAL::Triangulation_data_structure_2<Vb, Fb>;
     using CDT =
-        CGAL::Constrained_Delaunay_triangulation_2<K, CGAL::Default, Itag>;
+        CGAL::Constrained_Delaunay_triangulation_2<K, Tds, CGAL::Exact_predicates_tag>;
-    using Point = CDT::Point;
 
     // construct a constrained triangulation
     CDT cdt;
-    std::map<CDT::Vertex_handle, uint32_t> map;             // for indices
+    {
         std::vector<CDT::Vertex_handle> vertices_handle; // for constriants
         vertices_handle.reserve(points.size());
         for (const auto &p : points) {
-        Point cdt_p(p.x(), p.y());
-        auto handle = cdt.insert(cdt_p);
-        vertices_handle.push_back(handle);
-        // point index
             uint32_t pi = &p - &points.front();
-        map[handle]  = pi;
+            auto  handle = cdt.insert({ p.x(), p.y() });
+            handle->info() = pi;
+            vertices_handle.push_back(handle);
         }
-
+        // Constrain the triangulation.
-    // triangle can not contain forbiden edge
+        for (const HalfEdge &edge : constrained_half_edges)
-    for (const HalfEdge &edge : half_edges) {
+            cdt.insert_constraint(vertices_handle[edge.first], vertices_handle[edge.second]);
-        const CDT::Vertex_handle &vh1 = vertices_handle[edge.first];
-        const CDT::Vertex_handle &vh2 = vertices_handle[edge.second];
-        cdt.insert_constraint(vh1, vh2);
     }
 
     auto faces = cdt.finite_face_handles();
+
+    // Unmark constrained edges of outside faces.
+    size_t num_faces = 0;
+    for (CDT::Face_handle fh : faces) {
+        for (int i = 0; i < 3; ++ i) {
+            if (fh->is_constrained(i)) {
+                auto key = std::make_pair(fh->vertex((i + 2) % 3)->info(), fh->vertex((i + 1) % 3)->info());
+                if (auto it = std::lower_bound(constrained_half_edges.begin(), constrained_half_edges.end(), key); it != constrained_half_edges.end() && *it == key) {
+                    // This face contains a constrained edge and it is outside.
+                    for (int j = 0; j < 3; ++ j)
+                        fh->set_constraint(j, false);
+                    goto end;
+                }
+            }
+        }
+        ++ num_faces;
+    end:;
+    }
+
+    // Propagate inside the constrained regions.
+    std::vector<CDT::Face_handle> queue;
+    queue.reserve(num_faces);
+    auto inside = [](CDT::Face_handle &fh) { return fh->neighbor(0) != fh && (fh->is_constrained(0) || fh->is_constrained(1) || fh->is_constrained(2)); };
+    for (CDT::Face_handle seed : faces)
+        if (inside(seed)) {
+            // Seed fill to neighbor faces.
+            queue.emplace_back(seed);
+            while (! queue.empty()) {
+                CDT::Face_handle fh = queue.back();
+                queue.pop_back();
+                for (int i = 0; i < 3; ++ i)
+                    if (! fh->is_constrained(i)) {
+                        // Propagate along this edge.
+                        fh->set_constraint(i, true);
+                        CDT::Face_handle nh = fh->neighbor(i);
+                        bool was_inside = inside(nh);
+                        // Mark the other side of this edge.
+                        nh->set_constraint(nh->index(fh), true);
+                        if (! was_inside)
+                            queue.push_back(nh);
+                    }
+            }
+        }
+
     std::vector<Vec3i> indices;
-    indices.reserve(faces.size());
+    indices.reserve(num_faces);
-    for (CDT::Face_handle face : faces) {
+    for (CDT::Face_handle fh : faces)
-        // point indices
+        if (inside(fh))
-        std::array<uint32_t, 3> pi;
+            indices.emplace_back(fh->vertex(0)->info(), fh->vertex(1)->info(), fh->vertex(2)->info());
-        for (size_t i = 0; i < 3; ++i) pi[i] = map[face->vertex(i)];
-
-        // Do not use triangles with opposit edges
-        if (!allow_opposite_edge) {
-            if (half_edges.find(std::make_pair(pi[1], pi[0])) != half_edges.end()) continue;
-            if (half_edges.find(std::make_pair(pi[2], pi[1])) != half_edges.end()) continue;
-            if (half_edges.find(std::make_pair(pi[0], pi[2])) != half_edges.end()) continue;
-        }
-
-        indices.emplace_back(pi[0], pi[1], pi[2]);
-    }
     return indices;
 }
 
@@ -84,11 +112,11 @@ Triangulation::Indices Triangulation::triangulate(const Polygon &polygon)
 {
     const Points &pts = polygon.points;
     HalfEdges edges;
+    edges.reserve(pts.size());
     uint32_t offset = 0;
-    Private::insert_edge(edges, offset, polygon);
+    Private::insert_edges(edges, offset, polygon);
-    Triangulation::Indices indices = triangulate(pts, edges);
+    std::sort(edges.begin(), edges.end());
-    remove_outer(indices, edges);
+    return triangulate(pts, edges);
-    return indices;
 }
 
 Triangulation::Indices Triangulation::triangulate(const Polygons &polygons)
@@ -98,16 +126,16 @@ Triangulation::Indices Triangulation::triangulate(const Polygons &polygons)
     points.reserve(count);
 
     HalfEdges edges;
+    edges.reserve(count);
     uint32_t  offset = 0;
 
     for (const Polygon &polygon : polygons) {
         Private::insert_points(points, polygon);
-        Private::insert_edge(edges, offset, polygon);
+        Private::insert_edges(edges, offset, polygon);
     }
 
-    Triangulation::Indices indices = triangulate(points, edges);
+    std::sort(edges.begin(), edges.end());
-    remove_outer(indices, edges);
+    return triangulate(points, edges);
-    return indices;
 }
 
 Triangulation::Indices Triangulation::triangulate(const ExPolygons &expolygons)
@@ -117,105 +145,18 @@ Triangulation::Indices Triangulation::triangulate(const ExPolygons &expolygons)
     points.reserve(count);
 
     HalfEdges edges;
+    edges.reserve(count);
     uint32_t offset = 0;
 
     for (const ExPolygon &expolygon : expolygons) { 
         Private::insert_points(points, expolygon.contour);
-        Private::insert_edge(edges, offset, expolygon.contour);
+        Private::insert_edges(edges, offset, expolygon.contour);
         for (const Polygon &hole : expolygon.holes) {
             Private::insert_points(points, hole);
-            Private::insert_edge(edges, offset, hole);        
+            Private::insert_edges(edges, offset, hole);        
         }
     }        
 
-    Triangulation::Indices indices = triangulate(points, edges);
+    std::sort(edges.begin(), edges.end());
-    remove_outer(indices, edges);
+    return triangulate(points, edges);
-    return indices;
-}
-
-void Triangulation::remove_outer(Indices &indices, const HalfEdges &half_edges)
-{
-    // convert half edge to triangle index
-    std::map<HalfEdge, uint32_t> edge2triangle;
-    // triangles with all edges out of half_edge, candidate to remove
-    std::vector<uint32_t> triangles_to_check;
-    triangles_to_check.reserve(indices.size() / 3);
-
-    // triangle half edge
-    auto create_half_edge = [](const Vec3i& triangle, size_t index) {
-        size_t next_index = (index == 2) ? 0 : (index + 1);
-        return HalfEdge(triangle[index], triangle[next_index]);
-    };
-
-    // neighbor triangle half edge
-    auto create_opposit_half_edge = [](const Vec3i &triangle, size_t index) {
-        size_t prev_index = (index == 0) ? 2 : (index - 1);
-        return HalfEdge(triangle[index], triangle[prev_index]);
-    };
-
-    for (const auto &t : indices) {
-        uint32_t index     = &t - &indices.front();
-        bool     is_border = false;
-        for (size_t j = 0; j < 3; ++j) {
-            HalfEdge he = create_half_edge(t, j);
-            if (half_edges.find(he) != half_edges.end())
-                is_border = true;
-            else {
-                // half edge already exists
-                assert(edge2triangle.find(he) == edge2triangle.end());
-                edge2triangle[he] = index;
-            }
-        }
-        if (!is_border) triangles_to_check.push_back(index);
-    }
-
-    std::set<uint32_t>   remove;
-    std::queue<uint32_t> insert;
-    for (uint32_t index : triangles_to_check) {
-        auto it = remove.find(index);
-        if (it != remove.end()) continue; // already removed
-
-        bool is_edge = false;
-        const Vec3i &t = indices[index];
-        // check if exist opposit triangle
-        for (size_t j = 0; j < 3; ++j) {
-            HalfEdge he = create_opposit_half_edge(t, j);
-            auto it = edge2triangle.find(he);
-            if (it == edge2triangle.end()) {
-                is_edge = true;
-                break;
-            }
-            assert(it->second != index);
-        }
-        if (!is_edge) continue; // correct
-        // when there is no opposit edge, remove triangle and all neighbors recursive
-
-        insert.push(index);
-        do {
-            uint32_t triangle_index = insert.front();
-            insert.pop();
-            if (remove.find(triangle_index) != remove.end()) continue; // already removed
-            remove.insert(triangle_index);
-
-            for (size_t j = 0; j < 3; ++j) {
-                const Vec3i &t  = indices[triangle_index];
-                HalfEdge he = create_opposit_half_edge(t, j);
-                auto     it = edge2triangle.find(he);
-                if (it == edge2triangle.end()) continue; // edge triangle without neighbor
-                // process neighbor
-                uint32_t neighbor_index = it->second;
-                assert(neighbor_index != triangle_index);
-                insert.push(neighbor_index);
-            }
-        } while (!insert.empty());
-    }
-
-    // remove indices
-    std::vector<uint32_t> rem(remove.begin(), remove.end());
-    std::sort(rem.begin(), rem.end());
-    uint32_t offset = 0;
-    for (uint32_t i : rem) {
-        indices.erase(indices.begin() + (i - offset));
-        ++offset;
-    }
 }

src/libslic3r/Triangulation.hpp

@@ -16,7 +16,7 @@ public:
 
     // define oriented connection of 2 vertices(defined by its index)
     using HalfEdge  = std::pair<uint32_t, uint32_t>;
-    using HalfEdges = std::set<HalfEdge>;
+    using HalfEdges = std::vector<HalfEdge>;
     using Indices   = std::vector<Vec3i>;
 
     /// <summary>
@@ -26,12 +26,10 @@ public:
     /// </summary>
     /// <param name="points">Points to connect</param>
     /// <param name="edges">Constraint for edges, pair is from point(first) to
-    /// point(second)</param> 
+    /// point(second), sorted lexicographically</param> 
-    /// <param name="allow_opposite_edge">Flag for filtration result indices by opposit half edge</param>
     /// <returns>Triangles</returns>
-    static Indices triangulate(const Points &   points,
+    static Indices triangulate(const Points &points,
-                               const HalfEdges &half_edges,
+                               const HalfEdges &half_edges);
-                               bool allow_opposite_edge = false);
     static Indices triangulate(const Polygon &polygon);
     static Indices triangulate(const Polygons &polygons);
     static Indices triangulate(const ExPolygons &expolygons);

tests/libslic3r/test_triangulation.cpp

@@ -46,7 +46,7 @@ TEST_CASE("Triangulate rectangle with restriction on edge", "[Triangulation]")
     //                    0            1            2            3 
     Points points = {Point(1, 1), Point(2, 1), Point(2, 2), Point(1, 2)};
     Triangulation::HalfEdges edges1 = {{1, 3}};
-    std::vector<Vec3i> indices1 = Triangulation::triangulate(points, edges1, true);
+    std::vector<Vec3i> indices1 = Triangulation::triangulate(points, edges1);
 
     auto check = [](int i1, int i2, Vec3i t) -> bool {
         return true;
@@ -59,7 +59,7 @@ TEST_CASE("Triangulate rectangle with restriction on edge", "[Triangulation]")
     CHECK(check(i1, i2, indices1[1]));
 
     Triangulation::HalfEdges edges2 = {{0, 2}};
-    std::vector<Vec3i> indices2 = Triangulation::triangulate(points, edges2, true);
+    std::vector<Vec3i> indices2 = Triangulation::triangulate(points, edges2);
     REQUIRE(indices2.size() == 2);
     i1 = edges2.begin()->first;
     i2 = edges2.begin()->second;