From 59e14cb7520070a082117f0aa78edc639a1e7dab Mon Sep 17 00:00:00 2001 From: Vojtech Bubnik Date: Wed, 2 Mar 2022 14:47:54 +0100 Subject: [PATCH] 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 | 221 +++++++++---------------- src/libslic3r/Triangulation.hpp | 10 +- tests/libslic3r/test_triangulation.cpp | 4 +- 3 files changed, 87 insertions(+), 148 deletions(-) diff --git a/src/libslic3r/Triangulation.cpp b/src/libslic3r/Triangulation.cpp index 4ecdd43ae..8c075ea2e 100644 --- a/src/libslic3r/Triangulation.cpp +++ b/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(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, - const HalfEdges &half_edges, - bool allow_opposite_edge) +Triangulation::Indices Triangulation::triangulate(const Points &points, + const HalfEdges &constrained_half_edges) { // 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; + using Fb = CGAL::Constrained_triangulation_face_base_2; + using Tds = CGAL::Triangulation_data_structure_2; using CDT = - CGAL::Constrained_Delaunay_triangulation_2; - using Point = CDT::Point; + CGAL::Constrained_Delaunay_triangulation_2; // construct a constrained triangulation - CDT cdt; - std::map map; // for indices - std::vector 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; - } - - // triangle can not contain forbiden edge - for (const HalfEdge &edge : half_edges) { - const CDT::Vertex_handle &vh1 = vertices_handle[edge.first]; - const CDT::Vertex_handle &vh2 = vertices_handle[edge.second]; - cdt.insert_constraint(vh1, vh2); + CDT cdt; + { + std::vector vertices_handle; // for constriants + vertices_handle.reserve(points.size()); + for (const auto &p : points) { + uint32_t pi = &p - &points.front(); + auto handle = cdt.insert({ p.x(), p.y() }); + handle->info() = pi; + vertices_handle.push_back(handle); + } + // Constrain the triangulation. + for (const HalfEdge &edge : constrained_half_edges) + cdt.insert_constraint(vertices_handle[edge.first], vertices_handle[edge.second]); } auto faces = cdt.finite_face_handles(); - std::vector indices; - indices.reserve(faces.size()); - for (CDT::Face_handle face : faces) { - // point indices - std::array pi; - 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; + // 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 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); + } + } } - indices.emplace_back(pi[0], pi[1], pi[2]); - } + std::vector indices; + indices.reserve(num_faces); + for (CDT::Face_handle fh : faces) + if (inside(fh)) + indices.emplace_back(fh->vertex(0)->info(), fh->vertex(1)->info(), fh->vertex(2)->info()); 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); - Triangulation::Indices indices = triangulate(pts, edges); - remove_outer(indices, edges); - return indices; + Private::insert_edges(edges, offset, polygon); + std::sort(edges.begin(), edges.end()); + return triangulate(pts, edges); } 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); - remove_outer(indices, edges); - return indices; + std::sort(edges.begin(), edges.end()); + return triangulate(points, edges); } 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); - remove_outer(indices, edges); - return indices; + std::sort(edges.begin(), edges.end()); + return triangulate(points, edges); } - -void Triangulation::remove_outer(Indices &indices, const HalfEdges &half_edges) -{ - // convert half edge to triangle index - std::map edge2triangle; - // triangles with all edges out of half_edge, candidate to remove - std::vector 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 remove; - std::queue 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 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; - } -} \ No newline at end of file diff --git a/src/libslic3r/Triangulation.hpp b/src/libslic3r/Triangulation.hpp index ac847a8e4..1adf1b93c 100644 --- a/src/libslic3r/Triangulation.hpp +++ b/src/libslic3r/Triangulation.hpp @@ -16,7 +16,7 @@ public: // define oriented connection of 2 vertices(defined by its index) using HalfEdge = std::pair; - using HalfEdges = std::set; + using HalfEdges = std::vector; using Indices = std::vector; /// @@ -26,12 +26,10 @@ public: /// /// Points to connect /// Constraint for edges, pair is from point(first) to - /// point(second) - /// Flag for filtration result indices by opposit half edge + /// point(second), sorted lexicographically /// Triangles - static Indices triangulate(const Points & points, - const HalfEdges &half_edges, - bool allow_opposite_edge = false); + static Indices triangulate(const Points &points, + const HalfEdges &half_edges); static Indices triangulate(const Polygon &polygon); static Indices triangulate(const Polygons &polygons); static Indices triangulate(const ExPolygons &expolygons); diff --git a/tests/libslic3r/test_triangulation.cpp b/tests/libslic3r/test_triangulation.cpp index bc0f371d7..ed2bbef83 100644 --- a/tests/libslic3r/test_triangulation.cpp +++ b/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 indices1 = Triangulation::triangulate(points, edges1, true); + std::vector 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 indices2 = Triangulation::triangulate(points, edges2, true); + std::vector indices2 = Triangulation::triangulate(points, edges2); REQUIRE(indices2.size() == 2); i1 = edges2.begin()->first; i2 = edges2.begin()->second;