Fix of #8446: Non-planar Voronoi diagram.master_250

This is the follow-up to 63c66f4f189b6f3d3a3b95ddd158c66ee5166420. Detection of non-planar (degenerated) Voronoi diagrams was rewritten to check if all neighboring edges of the Voronoi vertex are CCW ordered.

6 files changed, 147 insertions, 55 deletions

diff --git a/src/libslic3r/Arachne/SkeletalTrapezoidation.cpp b/src/libslic3r/Arachne/SkeletalTrapezoidation.cpp
index d98cf7ba9..34d6d058d 100644
--- a/src/libslic3r/Arachne/SkeletalTrapezoidation.cpp
+++ b/src/libslic3r/Arachne/SkeletalTrapezoidation.cpp
@@ -17,12 +17,9 @@
 #include "Utils.hpp"
 #include "SVG.hpp"
 #include "Geometry/VoronoiVisualUtils.hpp"
+#include "Geometry/VoronoiUtilsCgal.hpp"
 #include "../EdgeGrid.hpp"
 
-#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
-#include <CGAL/Arr_segment_traits_2.h>
-#include <CGAL/Surface_sweep_2_algorithms.h>
-
 #define SKELETAL_TRAPEZOIDATION_BEAD_SEARCH_MAX 1000 //A limit to how long it'll keep searching for adjacent beads. Increasing will re-use beadings more often (saving performance), but search longer for beading (costing performance).
 
 //#define ARACHNE_DEBUG
@@ -391,40 +388,6 @@ SkeletalTrapezoidation::SkeletalTrapezoidation(const Polygons& polys, const Bead
     constructFromPolygons(polys);
 }
 
-using CGAL_Point   = CGAL::Exact_predicates_exact_constructions_kernel::Point_2;
-using CGAL_Segment = CGAL::Arr_segment_traits_2<CGAL::Exact_predicates_exact_constructions_kernel>::Curve_2;
-
-inline static CGAL_Point to_cgal_point(const Voronoi::VD::vertex_type &pt)
-{
-    return {pt.x(), pt.y()};
-}
-
-bool is_voronoi_diagram_planar(const Geometry::VoronoiDiagram &voronoi_diagram) {
-    assert(std::all_of(voronoi_diagram.edges().cbegin(), voronoi_diagram.edges().cend(),
-                       [](const Geometry::VoronoiDiagram::edge_type &edge) { return edge.color() == 0; }));
-
-    std::vector<CGAL_Segment> segments;
-    segments.reserve(voronoi_diagram.num_edges());
-
-    for (const Geometry::VoronoiDiagram::edge_type &edge : voronoi_diagram.edges()) {
-        if (edge.color() != 0)
-            continue;
-
-        if (edge.is_finite() && edge.is_linear()) {
-            segments.emplace_back(to_cgal_point(*edge.vertex0()), to_cgal_point(*edge.vertex1()));
-            edge.color(1);
-            assert(edge.twin() != nullptr);
-            edge.twin()->color(1);
-        }
-    }
-
-    for (const Geometry::VoronoiDiagram::edge_type &edge : voronoi_diagram.edges())
-        edge.color(0);
-
-    std::vector<CGAL_Point> intersections_pt;
-    CGAL::compute_intersection_points(segments.begin(), segments.end(), std::back_inserter(intersections_pt));
-    return intersections_pt.empty();
-}
 
 static bool detect_missing_voronoi_vertex(const Geometry::VoronoiDiagram &voronoi_diagram, const std::vector<SkeletalTrapezoidation::Segment> &segments) {
     for (VoronoiUtils::vd_t::cell_type cell : voronoi_diagram.cells()) {
@@ -506,7 +469,7 @@ inline static std::unordered_map<Point, Point, PointHash> try_to_fix_degenerated
     voronoi_diagram.clear();
     construct_voronoi(segments.begin(), segments.end(), &voronoi_diagram);
 
-    assert(is_voronoi_diagram_planar(voronoi_diagram));
+    assert(Geometry::VoronoiUtilsCgal::is_voronoi_diagram_planar_intersection(voronoi_diagram));
 
     return vertex_mapping;
 }
@@ -562,26 +525,38 @@ void SkeletalTrapezoidation::constructFromPolygons(const Polygons& polys)
 #endif
 
 #ifdef ARACHNE_DEBUG
-    assert(is_voronoi_diagram_planar(voronoi_diagram));
+    assert(is_voronoi_diagram_planar_intersection(voronoi_diagram));
 #endif
 
-    // Try to detect cases when some Voronoi vertex is missing.
-    // When any Voronoi vertex is missing, rotate input polygon and try again.
-    const bool                                  has_missing_voronoi_vertex = detect_missing_voronoi_vertex(voronoi_diagram, segments);
-    const double                                fix_angle                  = PI / 6;
+    // Try to detect cases when some Voronoi vertex is missing and when
+    // the Voronoi diagram is not planar.
+    // When any Voronoi vertex is missing, or the Voronoi diagram is not
+    // planar, rotate the input polygon and try again.
+    const bool   has_missing_voronoi_vertex = detect_missing_voronoi_vertex(voronoi_diagram, segments);
+    // Detection of non-planar Voronoi diagram detects at least GH issues #8474, #8514 and #8446.
+    const bool   is_voronoi_diagram_planar  = Geometry::VoronoiUtilsCgal::is_voronoi_diagram_planar_angle(voronoi_diagram);
+    const double fix_angle                  = PI / 6;
+
     std::unordered_map<Point, Point, PointHash> vertex_mapping;
     // polys_copy is referenced through items stored in the std::vector segments.
     Polygons                                    polys_copy = polys;
-    if (has_missing_voronoi_vertex) {
-        BOOST_LOG_TRIVIAL(warning) << "Detected missing Voronoi vertex, input polygons will be rotated back and forth.";
+    if (has_missing_voronoi_vertex || !is_voronoi_diagram_planar) {
+        if (has_missing_voronoi_vertex)
+            BOOST_LOG_TRIVIAL(warning) << "Detected missing Voronoi vertex, input polygons will be rotated back and forth.";
+        else if (!is_voronoi_diagram_planar)
+            BOOST_LOG_TRIVIAL(warning) << "Detected non-planar Voronoi diagram, input polygons will be rotated back and forth.";
+
         vertex_mapping = try_to_fix_degenerated_voronoi_diagram_by_rotation(voronoi_diagram, polys, polys_copy, segments, fix_angle);
 
         assert(!detect_missing_voronoi_vertex(voronoi_diagram, segments));
+        assert(Geometry::VoronoiUtilsCgal::is_voronoi_diagram_planar_angle(voronoi_diagram));
         if (detect_missing_voronoi_vertex(voronoi_diagram, segments))
             BOOST_LOG_TRIVIAL(error) << "Detected missing Voronoi vertex even after the rotation of input.";
+        else if (!Geometry::VoronoiUtilsCgal::is_voronoi_diagram_planar_angle(voronoi_diagram))
+            BOOST_LOG_TRIVIAL(error) << "Detected non-planar Voronoi diagram even after the rotation of input.";
     }
 
-    bool degenerated_voronoi_diagram = has_missing_voronoi_vertex;
+    bool degenerated_voronoi_diagram = has_missing_voronoi_vertex || !is_voronoi_diagram_planar;
 
 process_voronoi_diagram:
     assert(this->graph.edges.empty() && this->graph.nodes.empty() && this->vd_edge_to_he_edge.empty() && this->vd_node_to_he_node.empty());
@@ -644,8 +619,6 @@ process_voronoi_diagram:
     // When this degenerated Voronoi diagram is processed, the resulting half-edge structure contains some edges that don't have
     // a twin edge. Based on this, we created a fast mechanism that detects those causes and tries to recompute the Voronoi
     // diagram on slightly rotated input polygons that usually make the Voronoi generator generate a non-degenerated Voronoi diagram.
-    // FIXME Lukas H.: Replace has_missing_twin_edge with detection if the Voronoi diagram is planar using a custom algorithm based on the
-    //                 sweeping line algorithm. At least it is required to fix GH #8446.
     if (!degenerated_voronoi_diagram && has_missing_twin_edge(this->graph)) {
         BOOST_LOG_TRIVIAL(warning) << "Detected degenerated Voronoi diagram, input polygons will be rotated back and forth.";
         degenerated_voronoi_diagram = true;
@@ -655,7 +628,7 @@ process_voronoi_diagram:
         if (detect_missing_voronoi_vertex(voronoi_diagram, segments))
             BOOST_LOG_TRIVIAL(error) << "Detected missing Voronoi vertex after the rotation of input.";
 
-        assert(is_voronoi_diagram_planar(voronoi_diagram));
+        assert(Geometry::VoronoiUtilsCgal::is_voronoi_diagram_planar_intersection(voronoi_diagram));
 
         this->graph.edges.clear();
         this->graph.nodes.clear();
diff --git a/src/libslic3r/Arachne/SkeletalTrapezoidation.hpp b/src/libslic3r/Arachne/SkeletalTrapezoidation.hpp
index 321ace9d7..83065cf87 100644
--- a/src/libslic3r/Arachne/SkeletalTrapezoidation.hpp
+++ b/src/libslic3r/Arachne/SkeletalTrapezoidation.hpp
@@ -592,7 +592,5 @@ protected:
     void generateLocalMaximaSingleBeads();
 };
 
-bool is_voronoi_diagram_planar(const Geometry::VoronoiDiagram &voronoi_diagram);
-
 } // namespace Slic3r::Arachne
 #endif // VORONOI_QUADRILATERALIZATION_H
diff --git a/src/libslic3r/CMakeLists.txt b/src/libslic3r/CMakeLists.txt
index 498c01d87..41ad68db5 100644
--- a/src/libslic3r/CMakeLists.txt
+++ b/src/libslic3r/CMakeLists.txt
@@ -357,7 +357,7 @@ find_package(CGAL REQUIRED)
 cmake_policy(POP)
 
 add_library(libslic3r_cgal STATIC MeshBoolean.cpp MeshBoolean.hpp  TryCatchSignal.hpp
-    TryCatchSignal.cpp)
+    TryCatchSignal.cpp Geometry/VoronoiUtilsCgal.hpp Geometry/VoronoiUtilsCgal.cpp)
 target_include_directories(libslic3r_cgal PRIVATE ${CMAKE_CURRENT_BINARY_DIR})
 
 # Reset compile options of libslic3r_cgal. Despite it being linked privately, CGAL options
diff --git a/src/libslic3r/Geometry/VoronoiUtilsCgal.cpp b/src/libslic3r/Geometry/VoronoiUtilsCgal.cpp
new file mode 100644
index 000000000..caaf1ee9c
--- /dev/null
+++ b/src/libslic3r/Geometry/VoronoiUtilsCgal.cpp
@@ -0,0 +1,100 @@
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Arr_segment_traits_2.h>
+#include <CGAL/Surface_sweep_2_algorithms.h>
+
+#include "libslic3r/Geometry/Voronoi.hpp"
+
+#include "VoronoiUtilsCgal.hpp"
+
+using VD = Slic3r::Geometry::VoronoiDiagram;
+
+namespace Slic3r::Geometry {
+
+using CGAL_Point   = CGAL::Exact_predicates_exact_constructions_kernel::Point_2;
+using CGAL_Segment = CGAL::Arr_segment_traits_2<CGAL::Exact_predicates_exact_constructions_kernel>::Curve_2;
+
+inline static CGAL_Point to_cgal_point(const VD::vertex_type &pt) { return {pt.x(), pt.y()}; }
+
+// FIXME Lukas H.: Also includes parabolic segments.
+bool VoronoiUtilsCgal::is_voronoi_diagram_planar_intersection(const VD &voronoi_diagram)
+{
+    assert(std::all_of(voronoi_diagram.edges().cbegin(), voronoi_diagram.edges().cend(),
+                       [](const VD::edge_type &edge) { return edge.color() == 0; }));
+
+    std::vector<CGAL_Segment> segments;
+    segments.reserve(voronoi_diagram.num_edges());
+
+    for (const VD::edge_type &edge : voronoi_diagram.edges()) {
+        if (edge.color() != 0)
+            continue;
+
+        if (edge.is_finite() && edge.is_linear()) {
+            segments.emplace_back(to_cgal_point(*edge.vertex0()), to_cgal_point(*edge.vertex1()));
+            edge.color(1);
+            assert(edge.twin() != nullptr);
+            edge.twin()->color(1);
+        }
+    }
+
+    for (const VD::edge_type &edge : voronoi_diagram.edges())
+        edge.color(0);
+
+    std::vector<CGAL_Point> intersections_pt;
+    CGAL::compute_intersection_points(segments.begin(), segments.end(), std::back_inserter(intersections_pt));
+    return intersections_pt.empty();
+}
+
+static bool check_if_three_vectors_are_ccw(const CGAL_Point &common_pt, const CGAL_Point &pt_1, const CGAL_Point &pt_2, const CGAL_Point &test_pt) {
+    CGAL::Orientation orientation = CGAL::orientation(common_pt, pt_1, pt_2);
+    if (orientation == CGAL::Orientation::COLLINEAR) {
+        // The first two edges are collinear, so the third edge must be on the right side on the first of them.
+        return CGAL::orientation(common_pt, pt_1, test_pt) == CGAL::Orientation::RIGHT_TURN;
+    } else if (orientation == CGAL::Orientation::LEFT_TURN) {
+        // CCW oriented angle between vectors (common_pt, pt1) and (common_pt, pt2) is bellow PI.
+        // So we need to check if test_pt isn't between them.
+        CGAL::Orientation orientation1 = CGAL::orientation(common_pt, pt_1, test_pt);
+        CGAL::Orientation orientation2 = CGAL::orientation(common_pt, pt_2, test_pt);
+        return (orientation1 != CGAL::Orientation::LEFT_TURN || orientation2 != CGAL::Orientation::RIGHT_TURN);
+    } else {
+        assert(orientation == CGAL::Orientation::RIGHT_TURN);
+        // CCW oriented angle between vectors (common_pt, pt1) and (common_pt, pt2) is upper PI.
+        // So we need to check if test_pt is between them.
+        CGAL::Orientation orientation1 = CGAL::orientation(common_pt, pt_1, test_pt);
+        CGAL::Orientation orientation2 = CGAL::orientation(common_pt, pt_2, test_pt);
+        return (orientation1 == CGAL::Orientation::RIGHT_TURN || orientation2 == CGAL::Orientation::LEFT_TURN);
+    }
+}
+
+bool VoronoiUtilsCgal::is_voronoi_diagram_planar_angle(const VoronoiDiagram &voronoi_diagram)
+{
+    for (const VD::vertex_type &vertex : voronoi_diagram.vertices()) {
+        std::vector<const VD::edge_type *> edges;
+        const VD::edge_type               *edge = vertex.incident_edge();
+
+        do {
+            // FIXME Lukas H.: Also process parabolic segments.
+            if (edge->is_finite() && edge->is_linear())
+                edges.emplace_back(edge);
+
+            edge = edge->rot_next();
+        } while (edge != vertex.incident_edge());
+
+        // Checking for CCW make sense for three and more edges.
+        if (edges.size() > 2) {
+            for (auto edge_it = edges.begin() ; edge_it != edges.end(); ++edge_it) {
+                const Geometry::VoronoiDiagram::edge_type *prev_edge = edge_it == edges.begin() ? edges.back() : *std::prev(edge_it);
+                const Geometry::VoronoiDiagram::edge_type *curr_edge = *edge_it;
+                const Geometry::VoronoiDiagram::edge_type *next_edge = std::next(edge_it) == edges.end() ? edges.front() : *std::next(edge_it);
+
+                if (!check_if_three_vectors_are_ccw(to_cgal_point(*prev_edge->vertex0()), to_cgal_point(*prev_edge->vertex1()),
+                                                    to_cgal_point(*curr_edge->vertex1()), to_cgal_point(*next_edge->vertex1())))
+                    return false;
+            }
+        }
+    }
+
+    return true;
+}
+
+
+} // namespace Slic3r::Geometry
\ No newline at end of file
diff --git a/src/libslic3r/Geometry/VoronoiUtilsCgal.hpp b/src/libslic3r/Geometry/VoronoiUtilsCgal.hpp
new file mode 100644
index 000000000..897891bd9
--- /dev/null
+++ b/src/libslic3r/Geometry/VoronoiUtilsCgal.hpp
@@ -0,0 +1,21 @@
+#ifndef slic3r_VoronoiUtilsCgal_hpp_
+#define slic3r_VoronoiUtilsCgal_hpp_
+
+#include "Voronoi.hpp"
+
+namespace Slic3r::Geometry {
+class VoronoiDiagram;
+
+class VoronoiUtilsCgal
+{
+public:
+    // Check if the Voronoi diagram is planar using CGAL sweeping edge algorithm for enumerating all intersections between lines.
+    static bool is_voronoi_diagram_planar_intersection(const VoronoiDiagram &voronoi_diagram);
+
+    // Check if the Voronoi diagram is planar using verification that all neighboring edges are ordered CCW for each vertex.
+    static bool is_voronoi_diagram_planar_angle(const VoronoiDiagram &voronoi_diagram);
+
+};
+} // namespace Slic3r::Geometry
+
+#endif // slic3r_VoronoiUtilsCgal_hpp_
diff --git a/tests/libslic3r/test_voronoi.cpp b/tests/libslic3r/test_voronoi.cpp
index 5d8dce146..36323aa0a 100644
--- a/tests/libslic3r/test_voronoi.cpp
+++ b/tests/libslic3r/test_voronoi.cpp
@@ -5,7 +5,7 @@
 #include <libslic3r/Polyline.hpp>
 #include <libslic3r/EdgeGrid.hpp>
 #include <libslic3r/Geometry.hpp>
-#include <libslic3r/Arachne/SkeletalTrapezoidation.hpp>
+#include "libslic3r/Geometry/VoronoiUtilsCgal.hpp"
 
 #include <libslic3r/Geometry/VoronoiOffset.hpp>
 #include <libslic3r/Geometry/VoronoiVisualUtils.hpp>
@@ -2188,5 +2188,5 @@ TEST_CASE("Non-planar voronoi diagram", "[VoronoiNonPlanar]")
     dump_voronoi_to_svg(debug_out_path("voronoi-non-planar-out.svg").c_str(), vd, Points(), lines);
 #endif
 
-//    REQUIRE(Arachne::is_voronoi_diagram_planar(vd));
+//    REQUIRE(Geometry::VoronoiUtilsCgal::is_voronoi_diagram_planar_intersection(vd));
 }