[CURA-9178] Port external Voronoi fixes

include/SkeletalTrapezoidation.h

@@ -138,6 +138,9 @@ class SkeletalTrapezoidation
         {}
     };
 
+    bool detectMissingVoronoiVertex(const vd_t& voronoi_diagram, std::vector<Point>& points, const std::vector<SkeletalTrapezoidation::Segment>& segments);
+    bool isVoronoiDiagramPlanarAngle(const vd_t& voronoi_diagram);
+
     /*!
      * Compute the skeletal trapezoidation decomposition of the input shape.
      * 

include/utils/VoronoiUtils.h

@@ -32,6 +32,9 @@ class VoronoiUtils
 
     static Point p(const vd_t::vertex_type* node);
 
+    static bool isFinite(const vd_t::vertex_type* const v);
+    static bool hasFiniteEndpoints(const vd_t::edge_type* const edge);
+
     static bool isSourcePoint(Point p, const vd_t::cell_type& cell, const std::vector<Point>& points, const std::vector<Segment>& segments, coord_t snap_dist = 10);
 
     static coord_t getDistance(Point p, const vd_t::cell_type& cell, const std::vector<Point>& points, const std::vector<Segment>& segments);

src/SkeletalTrapezoidation.cpp

@@ -376,6 +376,123 @@ SkeletalTrapezoidation::SkeletalTrapezoidation(const Polygons& polys,
     constructFromPolygons(polys);
 }
 
+bool SkeletalTrapezoidation::detectMissingVoronoiVertex(const vd_t& voronoi_diagram, std::vector<Point>& points, const std::vector<SkeletalTrapezoidation::Segment>& segments)
+{
+    // Defensive programming: Handle any points at infinity, for supposedly 'finite' edges, as missing as well.
+    for (const auto& edge : voronoi_diagram.edges())
+    {
+        if (edge.is_finite() && ! VoronoiUtils::hasFiniteEndpoints(&edge))
+        {
+            return true;
+        }
+    }
+
+    // Now for the actual mising vertices:
+    for (vd_t::cell_type cell : voronoi_diagram.cells())
+    {
+        if (!cell.incident_edge())
+        {
+            continue; // There is no spoon
+        }
+
+        if (cell.contains_segment())
+        {
+            const SkeletalTrapezoidation::Segment& source_segment = VoronoiUtils::getSourceSegment(cell, points, segments);
+            const Point from = source_segment.from();
+            const Point to = source_segment.to();
+
+            // Find starting edge
+            // Find end edge
+            bool seen_possible_start = false;
+            bool after_start = false;
+            bool ending_edge_is_set_before_start = false;
+            VoronoiUtils::vd_t::edge_type* starting_vd_edge = nullptr;
+            VoronoiUtils::vd_t::edge_type* ending_vd_edge = nullptr;
+            VoronoiUtils::vd_t::edge_type* edge = cell.incident_edge();
+            do
+            {
+                if (edge->is_infinite() || ! VoronoiUtils::hasFiniteEndpoints(edge))
+                {
+                    continue;
+                }
+
+                const Point& v0 = VoronoiUtils::p(edge->vertex0());
+                const Point& v1 = VoronoiUtils::p(edge->vertex1());
+
+                assert(!(v0 == to && v1 == from));
+                if (v0 == to && !after_start)
+                { // Use the last edge which starts in source_segment.to
+                    starting_vd_edge = edge;
+                    seen_possible_start = true;
+                }
+                else if (seen_possible_start)
+                {
+                    after_start = true;
+                }
+
+                if (v1 == from && (!ending_vd_edge || ending_edge_is_set_before_start))
+                {
+                    ending_edge_is_set_before_start = !after_start;
+                    ending_vd_edge = edge;
+                }
+            } while (edge = edge->next(), edge != cell.incident_edge());
+
+            if (!starting_vd_edge || !ending_vd_edge || starting_vd_edge == ending_vd_edge)
+            {
+                return true;
+            }
+        }
+    }
+
+    return false;
+}
+
+bool SkeletalTrapezoidation::isVoronoiDiagramPlanarAngle(const vd_t& voronoi_diagram)
+{
+    for (const vd_t::vertex_type& vertex : voronoi_diagram.vertices())
+    {
+        std::vector<const vd_t::edge_type*> edges;
+        const vd_t::edge_type* edge = vertex.incident_edge();
+
+        do
+        {
+            // NOTE: Currently, it's not known if these degenaracies can also affect parabolic segments. They're not processed at the moment.
+            if (edge->is_finite() && edge->is_linear() && VoronoiUtils::hasFiniteEndpoints(edge))
+            {
+                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 vd_t::edge_type* prev_edge = edge_it == edges.begin() ? edges.back() : *std::prev(edge_it);
+                const vd_t::edge_type* curr_edge = *edge_it;
+                const vd_t::edge_type* next_edge = std::next(edge_it) == edges.end() ? edges.front() : *std::next(edge_it);
+
+                const bool isCCW =
+                    LinearAlg2D::isInsideCorner
+                    (
+                        VoronoiUtils::p(prev_edge->vertex0()),
+                        VoronoiUtils::p(prev_edge->vertex1()),
+                        VoronoiUtils::p(curr_edge->vertex1()),
+                        VoronoiUtils::p(next_edge->vertex1())
+                    );
+                if (! isCCW)
+                {
+                    return false;
+                }
+            }
+        }
+    }
+
+    return true;
+}
+
 void SkeletalTrapezoidation::constructFromPolygons(const Polygons& polys)
 {
     vd_edge_to_he_edge.clear();
@@ -393,10 +510,72 @@ void SkeletalTrapezoidation::constructFromPolygons(const Polygons& polys)
         }
     }
 
-    vd_t vonoroi_diagram;
-    construct_voronoi(segments.begin(), segments.end(), &vonoroi_diagram);
+    vd_t voronoi_diagram;
+    construct_voronoi(segments.begin(), segments.end(), &voronoi_diagram);
 
-    for (vd_t::cell_type cell : vonoroi_diagram.cells())
+    // Try to detect cases when some Voronoi vertex is missing.
+    // When any Voronoi vertex is missing, rotate input polygon and try again.
+    constexpr double fix_angle = 0.5;  // Choose an off-kilter (non-Pi/non-Tau based) angle to counteract rotational symmetry.
+    const bool has_missing_voronoi_vertex = detectMissingVoronoiVertex(voronoi_diagram, points, segments);
+    const bool is_voronoi_planar = isVoronoiDiagramPlanarAngle(voronoi_diagram);
+    const bool is_voronoi_misconstructed = has_missing_voronoi_vertex || ! is_voronoi_planar;
+    std::unordered_map<Point, Point> vertex_mapping;  // NOTE: Should maybe add a functor to specify specialized hash as 3rd template parameter.
+    Polygons polys_copy; // ATTN!: Don't move inside the scope below, 'segments' and 'vertex_mapping' above both can reference points out of it!
+    if (is_voronoi_misconstructed)
+    {
+        if (has_missing_voronoi_vertex)
+        {
+            spdlog::debug("Detected missing Voronoi vertex, input polygons will be rotated back and forth.");
+        }
+        if (! is_voronoi_planar)
+        {
+            spdlog::debug("Detected non-planar Voronoi diagram, input polygons will be rotated back and forth.");
+        }
+
+        polys_copy = polys;
+        const auto rot_matrix = LinearAlg2D::rotateAround(Point(0, 0), fix_angle);
+        for (auto& poly : polys_copy)
+        {
+            for (auto& pt : poly)
+            {
+                pt = rot_matrix.apply(pt);
+            }
+        }
+
+        assert(polys_copy.size() == polys.size());
+        for (size_t poly_idx = 0; poly_idx < polys.size(); ++poly_idx)
+        {
+            assert(polys_copy[poly_idx].size() == polys[poly_idx].size());
+            for (size_t point_idx = 0; point_idx < polys[poly_idx].size(); ++point_idx)
+            {
+                vertex_mapping.insert({ polys[poly_idx][point_idx], polys_copy[poly_idx][point_idx] });
+            }
+        }
+
+        segments.clear();
+        for (size_t poly_idx = 0; poly_idx < polys_copy.size(); poly_idx++)
+        {
+            for (size_t point_idx = 0; point_idx < polys_copy[poly_idx].size(); point_idx++)
+            {
+                segments.emplace_back(&polys_copy, poly_idx, point_idx);
+            }
+        }
+
+        voronoi_diagram.clear();
+        construct_voronoi(segments.begin(), segments.end(), &voronoi_diagram);
+        assert(!detectMissingVoronoiVertex(voronoi_diagram, points, segments));
+        assert(isVoronoiDiagramPlanarAngle(voronoi_diagram));
+        if (detectMissingVoronoiVertex(voronoi_diagram, points, segments))
+        {
+            spdlog::error("Detected missing Voronoi vertex even after the rotation of input.");
+        }
+        if (! isVoronoiDiagramPlanarAngle(voronoi_diagram))
+        {
+            spdlog::error("Detected non-planar Voronoi diagram even after the rotation of input.");
+        }
+    }
+
+    for (vd_t::cell_type cell : voronoi_diagram.cells())
     {
         if (! cell.incident_edge())
         { // There is no spoon
@@ -449,6 +628,23 @@ void SkeletalTrapezoidation::constructFromPolygons(const Polygons& polys)
         prev_edge->to->data.distance_to_boundary = 0;
     }
 
+    if (is_voronoi_misconstructed)
+    {
+        const auto inv_rot_matrix = LinearAlg2D::rotateAround(Point(0, 0), -fix_angle);
+        for (node_t& node : graph.nodes)
+        {
+            // If a mapping exists between a rotated point and an original point, use this mapping. Otherwise, rotate a point in the opposite direction.
+            if (auto node_it = vertex_mapping.find(node.p); node_it != vertex_mapping.end())
+            {
+                node.p = node_it->second;
+            }
+            else
+            {
+                node.p = inv_rot_matrix.apply(node.p);
+            }
+        }
+    }
+
     separatePointyQuadEndNodes();
 
     graph.collapseSmallEdges();
@@ -1701,7 +1897,8 @@ void SkeletalTrapezoidation::generateJunctions(ptr_vector_t<BeadingPropagation>&
         for (junction_idx = (std::max(size_t(1), beading->toolpath_locations.size()) - 1) / 2; junction_idx < num_junctions; junction_idx--)
         {
             coord_t bead_R = beading->toolpath_locations[junction_idx];
-            if (bead_R <= start_R)
+            // Adding a small epsilon (+1) to resolve an edge-case caused by rounding errors. (Would result in missing middle line.)
+            if (bead_R <= start_R + 1)
             { // Junction coinciding with start node is used in this function call
                 break;
             }

src/WallToolPaths.cpp

@@ -66,6 +66,7 @@ const std::vector<VariableWidthLines>& WallToolPaths::generate()
     PolygonUtils::fixSelfIntersections(epsilon_offset, prepared_outline);
     prepared_outline.removeDegenerateVerts();
     prepared_outline.removeSmallAreas(small_area_length * small_area_length, false);
+    prepared_outline = prepared_outline.unionPolygons();
 
     if (prepared_outline.area() <= 0)
     {

src/utils/VoronoiUtils.cpp

@@ -1,6 +1,7 @@
 // Copyright (c) 2022 Ultimaker B.V.
 // CuraEngine is released under the terms of the AGPLv3 or higher
 
+#include <cmath>
 #include <optional>
 #include <stack>
 
@@ -19,6 +20,16 @@ Point VoronoiUtils::p(const vd_t::vertex_type* node)
     return Point(x + 0.5 - (x < 0), y + 0.5 - (y < 0)); // Round to nearest integer coordinates.
 }
 
+bool VoronoiUtils::isFinite(const vd_t::vertex_type* const v)
+{
+    return v != nullptr && std::isfinite(v->x()) && std::isfinite(v->y());
+}
+
+bool VoronoiUtils::hasFiniteEndpoints(const vd_t::edge_type* const edge)
+{
+    return isFinite(edge->vertex0()) && isFinite(edge->vertex1());
+}
+
 bool VoronoiUtils::isSourcePoint(Point p, const vd_t::cell_type& cell, const std::vector<Point>& points, const std::vector<Segment>& segments, coord_t snap_dist)
 {
     if (cell.contains_point())