Update Arachne algorithm from PrusaSlicer 2.9 (#7819)

See individual commits for change details

src/libslic3r/Arachne/BeadingStrategy/BeadingStrategy.cpp

@@ -1,10 +1,8 @@
 //Copyright (c) 2022 Ultimaker B.V.
 //CuraEngine is released under the terms of the AGPLv3 or higher.
 
-#include <cassert>
-
 #include "BeadingStrategy.hpp"
-#include "Point.hpp"
+#include "libslic3r/Point.hpp"
 
 namespace Slic3r::Arachne
 {

src/libslic3r/Arachne/BeadingStrategy/BeadingStrategy.hpp

@@ -4,9 +4,13 @@
 #ifndef BEADING_STRATEGY_H
 #define BEADING_STRATEGY_H
 
+#include <math.h>
 #include <memory>
+#include <string>
+#include <vector>
+#include <cmath>
 
-#include "../../libslic3r.h"
+#include "libslic3r/libslic3r.h"
 
 namespace Slic3r::Arachne
 {

src/libslic3r/Arachne/BeadingStrategy/BeadingStrategyFactory.cpp

@@ -3,51 +3,55 @@
 
 #include "BeadingStrategyFactory.hpp"
 
+#include <boost/log/trivial.hpp>
+#include <memory>
+#include <utility>
+
 #include "LimitedBeadingStrategy.hpp"
 #include "WideningBeadingStrategy.hpp"
 #include "DistributedBeadingStrategy.hpp"
 #include "RedistributeBeadingStrategy.hpp"
 #include "OuterWallInsetBeadingStrategy.hpp"
+#include "libslic3r/Arachne/BeadingStrategy/BeadingStrategy.hpp"
 
-#include <limits>
-#include <boost/log/trivial.hpp>
+namespace Slic3r::Arachne {
 
-namespace Slic3r::Arachne
+BeadingStrategyPtr BeadingStrategyFactory::makeStrategy(const coord_t preferred_bead_width_outer,
+                                                        const coord_t preferred_bead_width_inner,
+                                                        const coord_t preferred_transition_length,
+                                                        const float   transitioning_angle,
+                                                        const bool    print_thin_walls,
+                                                        const coord_t min_bead_width,
+                                                        const coord_t min_feature_size,
+                                                        const double  wall_split_middle_threshold,
+                                                        const double  wall_add_middle_threshold,
+                                                        const coord_t max_bead_count,
+                                                        const coord_t outer_wall_offset,
+                                                        const int     inward_distributed_center_wall_count,
+                                                        const double  minimum_variable_line_ratio)
 {
+    // Handle a special case when there is just one external perimeter.
+    // Because big differences in bead width for inner and other perimeters cause issues with current beading strategies.
+    const coord_t      optimal_width = max_bead_count <= 2 ? preferred_bead_width_outer : preferred_bead_width_inner;
+    BeadingStrategyPtr ret = std::make_unique<DistributedBeadingStrategy>(optimal_width, preferred_transition_length, transitioning_angle,
+                                                                          wall_split_middle_threshold, wall_add_middle_threshold,
+                                                                          inward_distributed_center_wall_count);
 
-BeadingStrategyPtr BeadingStrategyFactory::makeStrategy(
-    const coord_t preferred_bead_width_outer,
-    const coord_t preferred_bead_width_inner,
-    const coord_t preferred_transition_length,
-    const float transitioning_angle,
-    const bool print_thin_walls,
-    const coord_t min_bead_width,
-    const coord_t min_feature_size,
-    const double wall_split_middle_threshold,
-    const double wall_add_middle_threshold,
-    const coord_t max_bead_count,
-    const coord_t outer_wall_offset,
-    const int inward_distributed_center_wall_count,
-    const double minimum_variable_line_ratio
-)
-{
-    BeadingStrategyPtr ret = std::make_unique<DistributedBeadingStrategy>(preferred_bead_width_inner, preferred_transition_length, transitioning_angle, wall_split_middle_threshold, wall_add_middle_threshold, inward_distributed_center_wall_count);
-    BOOST_LOG_TRIVIAL(debug) << "Applying the Redistribute meta-strategy with outer-wall width = " << preferred_bead_width_outer << ", inner-wall width = " << preferred_bead_width_inner << ".";
+    BOOST_LOG_TRIVIAL(trace) << "Applying the Redistribute meta-strategy with outer-wall width = " << preferred_bead_width_outer << ", inner-wall width = " << preferred_bead_width_inner << ".";
     ret = std::make_unique<RedistributeBeadingStrategy>(preferred_bead_width_outer, minimum_variable_line_ratio, std::move(ret));
 
     if (print_thin_walls) {
-        BOOST_LOG_TRIVIAL(debug) << "Applying the Widening Beading meta-strategy with minimum input width " << min_feature_size << " and minimum output width " << min_bead_width << ".";
+        BOOST_LOG_TRIVIAL(trace) << "Applying the Widening Beading meta-strategy with minimum input width " << min_feature_size << " and minimum output width " << min_bead_width << ".";
         ret = std::make_unique<WideningBeadingStrategy>(std::move(ret), min_feature_size, min_bead_width);
     }
     // Orca: we allow negative outer_wall_offset here
-    if (outer_wall_offset != 0) 
-    {
-        BOOST_LOG_TRIVIAL(debug) << "Applying the OuterWallOffset meta-strategy with offset = " << outer_wall_offset << ".";
+    if (outer_wall_offset != 0) {
+        BOOST_LOG_TRIVIAL(trace) << "Applying the OuterWallOffset meta-strategy with offset = " << outer_wall_offset << ".";
         ret = std::make_unique<OuterWallInsetBeadingStrategy>(outer_wall_offset, std::move(ret));
     }
 
-    //Apply the LimitedBeadingStrategy last, since that adds a 0-width marker wall which other beading strategies shouldn't touch.
-    BOOST_LOG_TRIVIAL(debug) << "Applying the Limited Beading meta-strategy with maximum bead count = " << max_bead_count << ".";
+    // Apply the LimitedBeadingStrategy last, since that adds a 0-width marker wall which other beading strategies shouldn't touch.
+    BOOST_LOG_TRIVIAL(trace) << "Applying the Limited Beading meta-strategy with maximum bead count = " << max_bead_count << ".";
     ret = std::make_unique<LimitedBeadingStrategy>(max_bead_count, std::move(ret));
     return ret;
 }

src/libslic3r/Arachne/BeadingStrategy/BeadingStrategyFactory.hpp

@@ -4,8 +4,12 @@
 #ifndef BEADING_STRATEGY_FACTORY_H
 #define BEADING_STRATEGY_FACTORY_H
 
+#include <math.h>
+#include <cmath>
+
 #include "BeadingStrategy.hpp"
 #include "../../Point.hpp"
+#include "libslic3r/libslic3r.h"
 
 namespace Slic3r::Arachne
 {

src/libslic3r/Arachne/BeadingStrategy/DistributedBeadingStrategy.cpp

@@ -1,7 +1,12 @@
 // Copyright (c) 2022 Ultimaker B.V.
 // CuraEngine is released under the terms of the AGPLv3 or higher.
 #include <numeric>
+#include <algorithm>
+#include <vector>
+#include <cassert>
+
 #include "DistributedBeadingStrategy.hpp"
+#include "libslic3r/Arachne/BeadingStrategy/BeadingStrategy.hpp"
 
 namespace Slic3r::Arachne
 {

src/libslic3r/Arachne/BeadingStrategy/DistributedBeadingStrategy.hpp

@@ -5,6 +5,7 @@
 #define DISTRIBUTED_BEADING_STRATEGY_H
 
 #include "BeadingStrategy.hpp"
+#include "libslic3r/libslic3r.h"
 
 namespace Slic3r::Arachne
 {

src/libslic3r/Arachne/BeadingStrategy/LimitedBeadingStrategy.cpp

@@ -1,11 +1,14 @@
 //Copyright (c) 2022 Ultimaker B.V.
 //CuraEngine is released under the terms of the AGPLv3 or higher.
 
-#include <cassert>
 #include <boost/log/trivial.hpp>
+#include <cassert>
+#include <utility>
+#include <cstddef>
 
 #include "LimitedBeadingStrategy.hpp"
-#include "Point.hpp"
+#include "libslic3r/Point.hpp"
+#include "libslic3r/Arachne/BeadingStrategy/BeadingStrategy.hpp"
 
 namespace Slic3r::Arachne
 {

src/libslic3r/Arachne/BeadingStrategy/LimitedBeadingStrategy.hpp

@@ -4,7 +4,10 @@
 #ifndef LIMITED_BEADING_STRATEGY_H
 #define LIMITED_BEADING_STRATEGY_H
 
+#include <string>
+
 #include "BeadingStrategy.hpp"
+#include "libslic3r/libslic3r.h"
 
 namespace Slic3r::Arachne
 {

src/libslic3r/Arachne/BeadingStrategy/OuterWallInsetBeadingStrategy.cpp

@@ -4,6 +4,9 @@
 #include "OuterWallInsetBeadingStrategy.hpp"
 
 #include <algorithm>
+#include <utility>
+
+#include "libslic3r/Arachne/BeadingStrategy/BeadingStrategy.hpp"
 
 namespace Slic3r::Arachne
 {

src/libslic3r/Arachne/BeadingStrategy/OuterWallInsetBeadingStrategy.hpp

@@ -4,7 +4,10 @@
 #ifndef OUTER_WALL_INSET_BEADING_STRATEGY_H
 #define OUTER_WALL_INSET_BEADING_STRATEGY_H
 
+#include <string>
+
 #include "BeadingStrategy.hpp"
+#include "libslic3r/libslic3r.h"
 
 namespace Slic3r::Arachne
 {

src/libslic3r/Arachne/BeadingStrategy/RedistributeBeadingStrategy.cpp

@@ -5,6 +5,9 @@
 
 #include <algorithm>
 #include <numeric>
+#include <utility>
+
+#include "libslic3r/Arachne/BeadingStrategy/BeadingStrategy.hpp"
 
 namespace Slic3r::Arachne
 {

src/libslic3r/Arachne/BeadingStrategy/RedistributeBeadingStrategy.hpp

@@ -4,7 +4,10 @@
 #ifndef REDISTRIBUTE_DISTRIBUTED_BEADING_STRATEGY_H
 #define REDISTRIBUTE_DISTRIBUTED_BEADING_STRATEGY_H
 
+#include <string>
+
 #include "BeadingStrategy.hpp"
+#include "libslic3r/libslic3r.h"
 
 namespace Slic3r::Arachne
 {

src/libslic3r/Arachne/BeadingStrategy/WideningBeadingStrategy.cpp

@@ -3,6 +3,11 @@
 
 #include "WideningBeadingStrategy.hpp"
 
+#include <algorithm>
+#include <utility>
+
+#include "libslic3r/Arachne/BeadingStrategy/BeadingStrategy.hpp"
+
 namespace Slic3r::Arachne
 {
 
@@ -24,17 +29,16 @@ WideningBeadingStrategy::Beading WideningBeadingStrategy::compute(coord_t thickn
     if (thickness < optimal_width) {
         Beading ret;
         ret.total_thickness = thickness;
-        if (thickness >= min_input_width)
-        {
+        if (thickness >= min_input_width) {
             ret.bead_widths.emplace_back(std::max(thickness, min_output_width));
             ret.toolpath_locations.emplace_back(thickness / 2);
-        } else {
+            ret.left_over = 0;
+        } else
             ret.left_over = thickness;
-        }
+
         return ret;
-    } else {
+    } else
         return parent->compute(thickness, bead_count);
-    }
 }
 
 coord_t WideningBeadingStrategy::getOptimalThickness(coord_t bead_count) const

src/libslic3r/Arachne/BeadingStrategy/WideningBeadingStrategy.hpp

@@ -4,7 +4,11 @@
 #ifndef WIDENING_BEADING_STRATEGY_H
 #define WIDENING_BEADING_STRATEGY_H
 
+#include <string>
+#include <vector>
+
 #include "BeadingStrategy.hpp"
+#include "libslic3r/libslic3r.h"
 
 namespace Slic3r::Arachne
 {

src/libslic3r/Arachne/PerimeterOrder.cpp

@@ -0,0 +1,280 @@
+#include <stack>
+#include <algorithm>
+#include <cmath>
+
+#include "PerimeterOrder.hpp"
+#include "libslic3r/Arachne/utils/ExtrusionJunction.hpp"
+#include "libslic3r/Point.hpp"
+
+namespace Slic3r::Arachne::PerimeterOrder {
+
+using namespace Arachne;
+
+static size_t get_extrusion_lines_count(const Perimeters &perimeters) {
+    size_t extrusion_lines_count = 0;
+    for (const Perimeter &perimeter : perimeters)
+        extrusion_lines_count += perimeter.size();
+
+    return extrusion_lines_count;
+}
+
+static PerimeterExtrusions get_sorted_perimeter_extrusions_by_area(const Perimeters &perimeters) {
+    PerimeterExtrusions sorted_perimeter_extrusions;
+    sorted_perimeter_extrusions.reserve(get_extrusion_lines_count(perimeters));
+
+    for (const Perimeter &perimeter : perimeters) {
+        for (const ExtrusionLine &extrusion_line : perimeter) {
+            if (extrusion_line.empty())
+                continue; // This shouldn't ever happen.
+
+            const BoundingBox bbox    = get_extents(extrusion_line);
+            // Be aware that Arachne produces contours with clockwise orientation and holes with counterclockwise orientation.
+            const double      area    = std::abs(extrusion_line.area());
+            const Polygon     polygon = extrusion_line.is_closed ? to_polygon(extrusion_line) : Polygon{};
+
+            sorted_perimeter_extrusions.emplace_back(extrusion_line, area, polygon, bbox);
+        }
+    }
+
+    // Open extrusions have an area equal to zero, so sorting based on the area ensures that open extrusions will always be before closed ones.
+    std::sort(sorted_perimeter_extrusions.begin(), sorted_perimeter_extrusions.end(),
+              [](const PerimeterExtrusion &l, const PerimeterExtrusion &r) { return l.area < r.area; });
+
+    return sorted_perimeter_extrusions;
+}
+
+// Functions fill adjacent_perimeter_extrusions field for every PerimeterExtrusion by pointers to PerimeterExtrusions that contain or are inside this PerimeterExtrusion.
+static void construct_perimeter_extrusions_adjacency_graph(PerimeterExtrusions &sorted_perimeter_extrusions) {
+    // Construct a graph (defined using adjacent_perimeter_extrusions field) where two PerimeterExtrusion are adjacent when one is inside the other.
+    std::vector<bool> root_candidates(sorted_perimeter_extrusions.size(), false);
+    for (PerimeterExtrusion &perimeter_extrusion : sorted_perimeter_extrusions) {
+        const size_t perimeter_extrusion_idx = &perimeter_extrusion - sorted_perimeter_extrusions.data();
+
+        if (!perimeter_extrusion.is_closed()) {
+            root_candidates[perimeter_extrusion_idx] = true;
+            continue;
+        }
+
+        for (PerimeterExtrusion &root_candidate : sorted_perimeter_extrusions) {
+            const size_t root_candidate_idx = &root_candidate - sorted_perimeter_extrusions.data();
+
+            if (!root_candidates[root_candidate_idx])
+                continue;
+
+            if (perimeter_extrusion.bbox.contains(root_candidate.bbox) && perimeter_extrusion.polygon.contains(root_candidate.extrusion.junctions.front().p)) {
+                perimeter_extrusion.adjacent_perimeter_extrusions.emplace_back(&root_candidate);
+                root_candidate.adjacent_perimeter_extrusions.emplace_back(&perimeter_extrusion);
+                root_candidates[root_candidate_idx] = false;
+            }
+        }
+
+        root_candidates[perimeter_extrusion_idx] = true;
+    }
+}
+
+// Perform the depth-first search to assign the nearest external perimeter for every PerimeterExtrusion.
+// When some PerimeterExtrusion is achievable from more than one external perimeter, then we choose the
+// one that comes from a contour.
+static void assign_nearest_external_perimeter(PerimeterExtrusions &sorted_perimeter_extrusions) {
+    std::stack<PerimeterExtrusion *> stack;
+    for (PerimeterExtrusion &perimeter_extrusion : sorted_perimeter_extrusions) {
+        if (perimeter_extrusion.is_external_perimeter()) {
+            perimeter_extrusion.depth                      = 0;
+            perimeter_extrusion.nearest_external_perimeter = &perimeter_extrusion;
+            stack.push(&perimeter_extrusion);
+        }
+    }
+
+    while (!stack.empty()) {
+        PerimeterExtrusion *current_extrusion = stack.top();
+        stack.pop();
+
+        for (PerimeterExtrusion *adjacent_extrusion : current_extrusion->adjacent_perimeter_extrusions) {
+            const size_t adjacent_extrusion_depth = current_extrusion->depth + 1;
+            // Update depth when the new depth is smaller or when we can achieve the same depth from a contour.
+            // This will ensure that the internal perimeter will be extruded before the outer external perimeter
+            // when there are two external perimeters and one internal.
+            if (adjacent_extrusion_depth < adjacent_extrusion->depth) {
+                adjacent_extrusion->nearest_external_perimeter = current_extrusion->nearest_external_perimeter;
+                adjacent_extrusion->depth                      = adjacent_extrusion_depth;
+                stack.push(adjacent_extrusion);
+            } else if (adjacent_extrusion_depth == adjacent_extrusion->depth && !adjacent_extrusion->nearest_external_perimeter->is_contour() && current_extrusion->is_contour()) {
+                adjacent_extrusion->nearest_external_perimeter = current_extrusion->nearest_external_perimeter;
+                stack.push(adjacent_extrusion);
+            }
+        }
+    }
+}
+
+inline Point get_end_position(const ExtrusionLine &extrusion) {
+    if (extrusion.is_closed)
+        return extrusion.junctions[0].p; // We ended where we started.
+    else
+        return extrusion.junctions.back().p; // Pick the other end from where we started.
+}
+
+// Returns ordered extrusions.
+static std::vector<const PerimeterExtrusion *> ordered_perimeter_extrusions_to_minimize_distances(Point current_position, std::vector<const PerimeterExtrusion *> extrusions) {
+    // Ensure that open extrusions will be placed before the closed one.
+    std::sort(extrusions.begin(), extrusions.end(),
+              [](const PerimeterExtrusion *l, const PerimeterExtrusion *r) -> bool { return l->is_closed() < r->is_closed(); });
+
+    std::vector<const PerimeterExtrusion *> ordered_extrusions;
+    std::vector<bool> already_selected(extrusions.size(), false);
+    while (ordered_extrusions.size() < extrusions.size()) {
+        double nearest_distance_sqr  = std::numeric_limits<double>::max();
+        size_t nearest_extrusion_idx = 0;
+        bool   is_nearest_closed     = false;
+
+        for (size_t extrusion_idx = 0; extrusion_idx < extrusions.size(); ++extrusion_idx) {
+            if (already_selected[extrusion_idx])
+                continue;
+
+            const ExtrusionLine &extrusion_line           = extrusions[extrusion_idx]->extrusion;
+            const Point         &extrusion_start_position = extrusion_line.junctions.front().p;
+            const double         distance_sqr             = (current_position - extrusion_start_position).cast<double>().squaredNorm();
+            if (distance_sqr < nearest_distance_sqr) {
+                if (extrusion_line.is_closed || (!extrusion_line.is_closed && nearest_distance_sqr == std::numeric_limits<double>::max()) || (!extrusion_line.is_closed && !is_nearest_closed)) {
+                    nearest_extrusion_idx = extrusion_idx;
+                    nearest_distance_sqr  = distance_sqr;
+                    is_nearest_closed     = extrusion_line.is_closed;
+                }
+            }
+        }
+
+        already_selected[nearest_extrusion_idx]     = true;
+        const PerimeterExtrusion *nearest_extrusion = extrusions[nearest_extrusion_idx];
+        current_position                            = get_end_position(nearest_extrusion->extrusion);
+        ordered_extrusions.emplace_back(nearest_extrusion);
+    }
+
+    return ordered_extrusions;
+}
+
+struct GroupedPerimeterExtrusions
+{
+    GroupedPerimeterExtrusions() = delete;
+    explicit GroupedPerimeterExtrusions(const PerimeterExtrusion *external_perimeter_extrusion)
+        : external_perimeter_extrusion(external_perimeter_extrusion) {}
+
+    std::vector<const PerimeterExtrusion *> extrusions;
+    const PerimeterExtrusion               *external_perimeter_extrusion  = nullptr;
+};
+
+// Returns vector of indexes that represent the order of grouped extrusions in grouped_extrusions.
+static std::vector<size_t> order_of_grouped_perimeter_extrusions_to_minimize_distances(Point current_position, std::vector<GroupedPerimeterExtrusions> grouped_extrusions) {
+    // Ensure that holes will be placed before contour and open extrusions before the closed one.
+    std::sort(grouped_extrusions.begin(), grouped_extrusions.end(), [](const GroupedPerimeterExtrusions &l, const GroupedPerimeterExtrusions &r) -> bool {
+        return (l.external_perimeter_extrusion->is_contour() <  r.external_perimeter_extrusion->is_contour()) ||
+               (l.external_perimeter_extrusion->is_contour() == r.external_perimeter_extrusion->is_contour()  && l.external_perimeter_extrusion->is_closed() < r.external_perimeter_extrusion->is_closed());
+    });
+
+    const size_t holes_cnt = std::count_if(grouped_extrusions.begin(), grouped_extrusions.end(), [](const GroupedPerimeterExtrusions &grouped_extrusions) {
+        return !grouped_extrusions.external_perimeter_extrusion->is_contour();
+    });
+
+    std::vector<size_t> grouped_extrusions_order;
+    std::vector<bool>   already_selected(grouped_extrusions.size(), false);
+    while (grouped_extrusions_order.size() < grouped_extrusions.size()) {
+        double nearest_distance_sqr           = std::numeric_limits<double>::max();
+        size_t nearest_grouped_extrusions_idx = 0;
+        bool   is_nearest_closed              = false;
+
+        // First we order all holes and then we start ordering contours.
+        const size_t grouped_extrusion_end = grouped_extrusions_order.size() < holes_cnt ? holes_cnt: grouped_extrusions.size();
+        for (size_t grouped_extrusion_idx = 0; grouped_extrusion_idx < grouped_extrusion_end; ++grouped_extrusion_idx) {
+            if (already_selected[grouped_extrusion_idx])
+                continue;
+
+            const ExtrusionLine &external_perimeter_extrusion_line = grouped_extrusions[grouped_extrusion_idx].external_perimeter_extrusion->extrusion;
+            const Point         &extrusion_start_position          = external_perimeter_extrusion_line.junctions.front().p;
+            const double         distance_sqr                      = (current_position - extrusion_start_position).cast<double>().squaredNorm();
+            if (distance_sqr < nearest_distance_sqr) {
+                if (external_perimeter_extrusion_line.is_closed || (!external_perimeter_extrusion_line.is_closed && nearest_distance_sqr == std::numeric_limits<double>::max()) || (!external_perimeter_extrusion_line.is_closed && !is_nearest_closed)) {
+                    nearest_grouped_extrusions_idx = grouped_extrusion_idx;
+                    nearest_distance_sqr           = distance_sqr;
+                    is_nearest_closed              = external_perimeter_extrusion_line.is_closed;
+                }
+            }
+        }
+
+        grouped_extrusions_order.emplace_back(nearest_grouped_extrusions_idx);
+        already_selected[nearest_grouped_extrusions_idx]             = true;
+        const GroupedPerimeterExtrusions &nearest_grouped_extrusions = grouped_extrusions[nearest_grouped_extrusions_idx];
+        const ExtrusionLine              &last_extrusion_line        = nearest_grouped_extrusions.extrusions.back()->extrusion;
+        current_position                                             = get_end_position(last_extrusion_line);
+    }
+
+    return grouped_extrusions_order;
+}
+
+static PerimeterExtrusions extract_ordered_perimeter_extrusions(const PerimeterExtrusions &sorted_perimeter_extrusions, const bool external_perimeters_first) {
+    // Extrusions are ordered inside each group.
+    std::vector<GroupedPerimeterExtrusions> grouped_extrusions;
+
+    std::stack<const PerimeterExtrusion *> stack;
+    std::vector<bool>                      visited(sorted_perimeter_extrusions.size(), false);
+    for (const PerimeterExtrusion &perimeter_extrusion : sorted_perimeter_extrusions) {
+        if (!perimeter_extrusion.is_external_perimeter())
+            continue;
+
+        stack.push(&perimeter_extrusion);
+        visited.assign(sorted_perimeter_extrusions.size(), false);
+
+        grouped_extrusions.emplace_back(&perimeter_extrusion);
+        while (!stack.empty()) {
+            const PerimeterExtrusion *current_extrusion     = stack.top();
+            const size_t              current_extrusion_idx = current_extrusion - sorted_perimeter_extrusions.data();
+
+            stack.pop();
+            visited[current_extrusion_idx] = true;
+
+            if (current_extrusion->nearest_external_perimeter == &perimeter_extrusion) {
+                grouped_extrusions.back().extrusions.emplace_back(current_extrusion);
+            }
+
+            std::vector<const PerimeterExtrusion *> available_candidates;
+            for (const PerimeterExtrusion *adjacent_extrusion : current_extrusion->adjacent_perimeter_extrusions) {
+                const size_t adjacent_extrusion_idx = adjacent_extrusion - sorted_perimeter_extrusions.data();
+                if (!visited[adjacent_extrusion_idx] && !adjacent_extrusion->is_external_perimeter() && adjacent_extrusion->nearest_external_perimeter == &perimeter_extrusion) {
+                    available_candidates.emplace_back(adjacent_extrusion);
+                }
+            }
+
+            if (available_candidates.size() == 1) {
+                stack.push(available_candidates.front());
+            } else if (available_candidates.size() > 1) {
+                // When there is more than one available candidate, then order candidates to minimize distances between
+                // candidates and also to minimize the distance from the current_position.
+                std::vector<const PerimeterExtrusion *> adjacent_extrusions = ordered_perimeter_extrusions_to_minimize_distances(Point::Zero(), available_candidates);
+                for (auto extrusion_it = adjacent_extrusions.rbegin(); extrusion_it != adjacent_extrusions.rend(); ++extrusion_it) {
+                    stack.push(*extrusion_it);
+                }
+            }
+        }
+
+        if (!external_perimeters_first)
+            std::reverse(grouped_extrusions.back().extrusions.begin(), grouped_extrusions.back().extrusions.end());
+    }
+
+    const std::vector<size_t> grouped_extrusion_order = order_of_grouped_perimeter_extrusions_to_minimize_distances(Point::Zero(), grouped_extrusions);
+
+    PerimeterExtrusions ordered_extrusions;
+    for (size_t order_idx : grouped_extrusion_order) {
+        for (const PerimeterExtrusion *perimeter_extrusion : grouped_extrusions[order_idx].extrusions)
+            ordered_extrusions.emplace_back(*perimeter_extrusion);
+    }
+
+    return ordered_extrusions;
+}
+
+// FIXME: From the point of better patch planning, it should be better to do ordering when we have generated all extrusions (for now, when G-Code is exported).
+// FIXME: It would be better to extract the adjacency graph of extrusions from the SkeletalTrapezoidation graph.
+PerimeterExtrusions ordered_perimeter_extrusions(const Perimeters &perimeters, const bool external_perimeters_first) {
+    PerimeterExtrusions sorted_perimeter_extrusions = get_sorted_perimeter_extrusions_by_area(perimeters);
+    construct_perimeter_extrusions_adjacency_graph(sorted_perimeter_extrusions);
+    assign_nearest_external_perimeter(sorted_perimeter_extrusions);
+    return extract_ordered_perimeter_extrusions(sorted_perimeter_extrusions, external_perimeters_first);
+}
+
+} // namespace Slic3r::Arachne::PerimeterOrder

src/libslic3r/Arachne/PerimeterOrder.hpp

@@ -0,0 +1,51 @@
+#ifndef slic3r_GCode_PerimeterOrder_hpp_
+#define slic3r_GCode_PerimeterOrder_hpp_
+
+#include <stddef.h>
+#include <limits>
+#include <vector>
+#include <cstddef>
+
+#include "libslic3r/Arachne/utils/ExtrusionLine.hpp"
+#include "libslic3r/BoundingBox.hpp"
+#include "libslic3r/Polygon.hpp"
+
+namespace Slic3r::Arachne::PerimeterOrder {
+
+// Data structure stores ExtrusionLine (closed and open) together with additional data.
+struct PerimeterExtrusion
+{
+    explicit PerimeterExtrusion(const Arachne::ExtrusionLine &extrusion, const double area, const Polygon &polygon, const BoundingBox &bbox)
+        : extrusion(extrusion), area(area), polygon(polygon), bbox(bbox) {}
+
+    Arachne::ExtrusionLine             extrusion;
+    // Absolute value of the area of the polygon. The value is always non-negative, even for holes.
+    double                             area = 0;
+
+    // Polygon is non-empty only for closed extrusions.
+    Polygon                            polygon;
+    BoundingBox                        bbox;
+
+    std::vector<PerimeterExtrusion *>  adjacent_perimeter_extrusions;
+
+    // How far is this perimeter from the nearest external perimeter. Contour is always preferred over holes.
+    size_t                             depth                      = std::numeric_limits<size_t>::max();
+    PerimeterExtrusion                *nearest_external_perimeter = nullptr;
+
+    // Returns if ExtrusionLine is a contour or a hole.
+    bool is_contour() const { return extrusion.is_contour(); }
+
+    // Returns if ExtrusionLine is closed or opened.
+    bool is_closed() const { return extrusion.is_closed; }
+
+    // Returns if ExtrusionLine is an external or an internal perimeter.
+    bool is_external_perimeter() const { return extrusion.is_external_perimeter(); }
+};
+
+using PerimeterExtrusions = std::vector<PerimeterExtrusion>;
+
+PerimeterExtrusions ordered_perimeter_extrusions(const Perimeters &perimeters, bool external_perimeters_first);
+
+} // namespace Slic3r::Arachne::PerimeterOrder
+
+#endif // slic3r_GCode_Travels_hpp_

src/libslic3r/Arachne/SkeletalTrapezoidation.cpp

@@ -3,21 +3,27 @@
 
 #include "SkeletalTrapezoidation.hpp"
 
-#include <stack>
-#include <functional>
-#include <sstream>
-#include <queue>
-#include <functional>
 #include <boost/log/trivial.hpp>
+#include <boost/polygon/polygon.hpp>
+#include <queue>
+#include <algorithm>
+#include <cmath>
+#include <cstdint>
+#include <limits>
+#include <utility>
+#include <cassert>
+#include <cstdlib>
 
-#include "utils/linearAlg2D.hpp"
-#include "Utils.hpp"
-#include "SVG.hpp"
-#include "Geometry/VoronoiVisualUtils.hpp"
-#include "Geometry/VoronoiUtilsCgal.hpp"
-#include "../EdgeGrid.hpp"
+#include "libslic3r/Geometry/VoronoiUtils.hpp"
+#include "ankerl/unordered_dense.h"
+#include "libslic3r/Arachne/SkeletalTrapezoidationEdge.hpp"
+#include "libslic3r/Arachne/SkeletalTrapezoidationJoint.hpp"
+#include "libslic3r/Arachne/utils/ExtrusionJunction.hpp"
+#include "libslic3r/Arachne/utils/ExtrusionLine.hpp"
 
-#include "Geometry/VoronoiUtils.hpp"
+#ifndef NDEBUG
+    #include "libslic3r/EdgeGrid.hpp"
+#endif
 
 #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).
 
@@ -104,7 +110,7 @@ SkeletalTrapezoidation::node_t &SkeletalTrapezoidation::makeNode(const VD::verte
     }
 }
 
-void SkeletalTrapezoidation::transferEdge(Point from, Point to, const VD::edge_type &vd_edge, edge_t *&prev_edge, Point &start_source_point, Point &end_source_point, const std::vector<Segment> &segments) {
+void SkeletalTrapezoidation::transferEdge(const Point &from, const Point &to, const VD::edge_type &vd_edge, edge_t *&prev_edge, const Point &start_source_point, const Point &end_source_point, const std::vector<Segment> &segments) {
     auto he_edge_it = vd_edge_to_he_edge.find(vd_edge.twin());
     if (he_edge_it != vd_edge_to_he_edge.end())
     { // Twin segment(s) have already been made
@@ -152,8 +158,7 @@ void SkeletalTrapezoidation::transferEdge(Point from, Point to, const VD::edge_t
             assert(twin->prev->twin); // Back rib
             assert(twin->prev->twin->prev); // Prev segment along parabola
             
-            constexpr bool is_not_next_to_start_or_end = false; // Only ribs at the end of a cell should be skipped
-            graph.makeRib(prev_edge, start_source_point, end_source_point, is_not_next_to_start_or_end);
+            graph.makeRib(prev_edge, start_source_point, end_source_point);
         }
         assert(prev_edge);
     }
@@ -203,10 +208,8 @@ void SkeletalTrapezoidation::transferEdge(Point from, Point to, const VD::edge_t
             p0 = p1;
             v0 = v1;
             
-            if (p1_idx < discretized.size() - 1)
-            { // Rib for last segment gets introduced outside this function!
-                constexpr bool is_not_next_to_start_or_end = false; // Only ribs at the end of a cell should be skipped
-                graph.makeRib(prev_edge, start_source_point, end_source_point, is_not_next_to_start_or_end);
+            if (p1_idx < discretized.size() - 1) { // Rib for last segment gets introduced outside this function!
+                graph.makeRib(prev_edge, start_source_point, end_source_point);
             }
         }
         assert(prev_edge);
@@ -326,50 +329,6 @@ Points SkeletalTrapezoidation::discretize(const VD::edge_type& vd_edge, const st
     }
 }
 
-bool SkeletalTrapezoidation::computePointCellRange(const VD::cell_type &cell, Point &start_source_point, Point &end_source_point, const VD::edge_type *&starting_vd_edge, const VD::edge_type *&ending_vd_edge, const std::vector<Segment> &segments) {
-    if (cell.incident_edge()->is_infinite())
-        return false; //Infinite edges only occur outside of the polygon. Don't copy any part of this cell.
-
-    // Check if any point of the cell is inside or outside polygon
-    // Copy whole cell into graph or not at all
-
-    // If the cell.incident_edge()->vertex0() is far away so much that it doesn't even fit into Vec2i64, then there is no way that it will be inside the input polygon.
-    if (const VD::vertex_type &vert = *cell.incident_edge()->vertex0();
-        vert.x() >= double(std::numeric_limits<int64_t>::max()) || vert.x() <= double(std::numeric_limits<int64_t>::lowest()) ||
-        vert.y() >= double(std::numeric_limits<int64_t>::max()) || vert.y() <= double(std::numeric_limits<int64_t>::lowest()))
-        return false; // Don't copy any part of this cell
-
-    const Point              source_point       = Geometry::VoronoiUtils::get_source_point(cell, segments.begin(), segments.end());
-    const PolygonsPointIndex source_point_index = Geometry::VoronoiUtils::get_source_point_index(cell, segments.begin(), segments.end());
-    Vec2i64                  some_point         = Geometry::VoronoiUtils::to_point(cell.incident_edge()->vertex0());
-    if (some_point == source_point.cast<int64_t>())
-        some_point = Geometry::VoronoiUtils::to_point(cell.incident_edge()->vertex1());
-
-    //Test if the some_point is even inside the polygon.
-    //The edge leading out of a polygon must have an endpoint that's not in the corner following the contour of the polygon at that vertex.
-    //So if it's inside the corner formed by the polygon vertex, it's all fine.
-    //But if it's outside of the corner, it must be a vertex of the Voronoi diagram that goes outside of the polygon towards infinity.
-    if (!LinearAlg2D::isInsideCorner(source_point_index.prev().p(), source_point_index.p(), source_point_index.next().p(), some_point))
-        return false; // Don't copy any part of this cell
-
-    const VD::edge_type* vd_edge = cell.incident_edge();
-    do {
-        assert(vd_edge->is_finite());
-        if (Vec2i64 p1 = Geometry::VoronoiUtils::to_point(vd_edge->vertex1()); p1 == source_point.cast<int64_t>()) {
-            start_source_point = source_point;
-            end_source_point = source_point;
-            starting_vd_edge = vd_edge->next();
-            ending_vd_edge = vd_edge;
-        } else {
-            assert((Geometry::VoronoiUtils::to_point(vd_edge->vertex0()) == source_point.cast<int64_t>() || !vd_edge->is_secondary()) && "point cells must end in the point! They cannot cross the point with an edge, because collinear edges are not allowed in the input.");
-        }
-    }
-    while (vd_edge = vd_edge->next(), vd_edge != cell.incident_edge());
-    assert(starting_vd_edge && ending_vd_edge);
-    assert(starting_vd_edge != ending_vd_edge);
-    return true;
-}
-
 SkeletalTrapezoidation::SkeletalTrapezoidation(const Polygons& polys, const BeadingStrategy& beading_strategy,
                                                double transitioning_angle, coord_t discretization_step_size,
                                                coord_t transition_filter_dist, coord_t allowed_filter_deviation,
@@ -437,15 +396,20 @@ void SkeletalTrapezoidation::constructFromPolygons(const Polygons& polys)
         // Compute and store result in above variables
 
         if (cell.contains_point()) {
-            const bool keep_going = computePointCellRange(cell, start_source_point, end_source_point, starting_voronoi_edge, ending_voronoi_edge, segments);
-            if (!keep_going)
+            Geometry::PointCellRange<Point> cell_range = Geometry::VoronoiUtils::compute_point_cell_range(cell, segments.cbegin(), segments.cend());
+            start_source_point    = cell_range.source_point;
+            end_source_point      = cell_range.source_point;
+            starting_voronoi_edge = cell_range.edge_begin;
+            ending_voronoi_edge   = cell_range.edge_end;
+
+            if (!cell_range.is_valid())
                 continue;
         } else {
             assert(cell.contains_segment());
             Geometry::SegmentCellRange<Point> cell_range = Geometry::VoronoiUtils::compute_segment_cell_range(cell, segments.cbegin(), segments.cend());
             assert(cell_range.is_valid());
-            start_source_point    = cell_range.segment_start_point;
-            end_source_point      = cell_range.segment_end_point;
+            start_source_point    = cell_range.source_segment_start_point;
+            end_source_point      = cell_range.source_segment_end_point;
             starting_voronoi_edge = cell_range.edge_begin;
             ending_voronoi_edge   = cell_range.edge_end;
         }
@@ -462,8 +426,7 @@ void SkeletalTrapezoidation::constructFromPolygons(const Polygons& polys)
         node_t *starting_node                    = vd_node_to_he_node[starting_voronoi_edge->vertex0()];
         starting_node->data.distance_to_boundary = 0;
 
-        constexpr bool is_next_to_start_or_end = true;
-        graph.makeRib(prev_edge, start_source_point, end_source_point, is_next_to_start_or_end);
+        graph.makeRib(prev_edge, start_source_point, end_source_point);
         for (const VD::edge_type* vd_edge = starting_voronoi_edge->next(); vd_edge != ending_voronoi_edge; vd_edge = vd_edge->next()) {
             assert(vd_edge->is_finite());
             assert(Geometry::VoronoiUtils::is_in_range<coord_t>(*vd_edge));
@@ -471,7 +434,7 @@ void SkeletalTrapezoidation::constructFromPolygons(const Polygons& polys)
             Point v1 = Geometry::VoronoiUtils::to_point(vd_edge->vertex0()).cast<coord_t>();
             Point v2 = Geometry::VoronoiUtils::to_point(vd_edge->vertex1()).cast<coord_t>();
             transferEdge(v1, v2, *vd_edge, prev_edge, start_source_point, end_source_point, segments);
-            graph.makeRib(prev_edge, start_source_point, end_source_point, vd_edge->next() == ending_voronoi_edge);
+            graph.makeRib(prev_edge, start_source_point, end_source_point);
         }
 
         transferEdge(Geometry::VoronoiUtils::to_point(ending_voronoi_edge->vertex0()).cast<coord_t>(), end_source_point, *ending_voronoi_edge, prev_edge, start_source_point, end_source_point, segments);
@@ -493,6 +456,8 @@ void SkeletalTrapezoidation::constructFromPolygons(const Polygons& polys)
             edge.from->incident_edge = &edge;
 }
 
+using NodeSet = SkeletalTrapezoidation::NodeSet;
+
 void SkeletalTrapezoidation::separatePointyQuadEndNodes()
 {
     NodeSet visited_nodes;
@@ -1968,6 +1933,8 @@ void SkeletalTrapezoidation::addToolpathSegment(const ExtrusionJunction& from, c
 
 void SkeletalTrapezoidation::connectJunctions(ptr_vector_t<LineJunctions>& edge_junctions)
 {
+    using EdgeSet = ankerl::unordered_dense::set<edge_t*>;
+
     EdgeSet unprocessed_quad_starts(graph.edges.size() * 5 / 2);
     for (edge_t& edge : graph.edges)
     {

src/libslic3r/Arachne/SkeletalTrapezoidation.hpp

@@ -5,11 +5,11 @@
 #define SKELETAL_TRAPEZOIDATION_H
 
 #include <boost/polygon/voronoi.hpp>
-
+#include <ankerl/unordered_dense.h>
 #include <memory> // smart pointers
 #include <utility> // pair
-
-#include <ankerl/unordered_dense.h>
+#include <list>
+#include <vector>
 
 #include "utils/HalfEdgeGraph.hpp"
 #include "utils/PolygonsSegmentIndex.hpp"
@@ -20,6 +20,10 @@
 #include "libslic3r/Arachne/BeadingStrategy/BeadingStrategy.hpp"
 #include "SkeletalTrapezoidationGraph.hpp"
 #include "../Geometry/Voronoi.hpp"
+#include "libslic3r/Line.hpp"
+#include "libslic3r/Point.hpp"
+#include "libslic3r/Polygon.hpp"
+#include "libslic3r/libslic3r.h"
 
 //#define ARACHNE_DEBUG
 //#define ARACHNE_DEBUG_VORONOI
@@ -82,10 +86,7 @@ class SkeletalTrapezoidation
 
 public:
     using Segment  = PolygonsSegmentIndex;
-    using NodeSet = ankerl::unordered_dense::set<node_t *>;
-    using EdgeSet = ankerl::unordered_dense::set<edge_t *>;
-    using EdgeMap = ankerl::unordered_dense::map<const VD::edge_type *, edge_t *>;
-    using NodeMap = ankerl::unordered_dense::map<const VD::vertex_type *, node_t *>;
+    using NodeSet  = ankerl::unordered_dense::set<node_t*>;
 
     /*!
      * Construct a new trapezoidation problem to solve.
@@ -169,8 +170,8 @@ protected:
      * mapping each voronoi VD edge to the corresponding halfedge HE edge
      * In case the result segment is discretized, we map the VD edge to the *last* HE edge
      */
-    EdgeMap vd_edge_to_he_edge;
-    NodeMap vd_node_to_he_node;
+    ankerl::unordered_dense::map<const VD::edge_type *, edge_t *> vd_edge_to_he_edge;
+    ankerl::unordered_dense::map<const VD::vertex_type *, node_t *> vd_node_to_he_node;
     node_t &makeNode(const VD::vertex_type &vd_node, Point p); //!< Get the node which the VD node maps to, or create a new mapping if there wasn't any yet.
 
     /*!
@@ -182,7 +183,7 @@ protected:
      * Transfer an edge from the VD to the HE and perform discretization of parabolic edges (and vertex-vertex edges)
      * \p prev_edge serves as input and output. May be null as input.
      */
-    void transferEdge(Point from, Point to, const VD::edge_type &vd_edge, edge_t *&prev_edge, Point &start_source_point, Point &end_source_point, const std::vector<Segment> &segments);
+    void transferEdge(const Point &from, const Point &to, const VD::edge_type &vd_edge, edge_t *&prev_edge, const Point &start_source_point, const Point &end_source_point, const std::vector<Segment> &segments);
 
     /*!
      * Discretize a Voronoi edge that represents the medial axis of a vertex-
@@ -211,32 +212,6 @@ protected:
      */
     Points discretize(const VD::edge_type& segment, const std::vector<Segment>& segments);
 
-    /*!
-     * Compute the range of line segments that surround a cell of the skeletal
-     * graph that belongs to a point on the medial axis.
-     *
-     * This should only be used on cells that belong to a corner in the skeletal
-     * graph, e.g. triangular cells, not trapezoid cells.
-     *
-     * The resulting line segments is just the first and the last segment. They
-     * are linked to the neighboring segments, so you can iterate over the
-     * segments until you reach the last segment.
-     * \param cell The cell to compute the range of line segments for.
-     * \param[out] start_source_point The start point of the source segment of
-     * this cell.
-     * \param[out] end_source_point The end point of the source segment of this
-     * cell.
-     * \param[out] starting_vd_edge The edge of the Voronoi diagram where the
-     * loop around the cell starts.
-     * \param[out] ending_vd_edge The edge of the Voronoi diagram where the loop
-     * around the cell ends.
-     * \param points All vertices of the input Polygons.
-     * \param segments All edges of the input Polygons.
-     * /return Whether the cell is inside of the polygon. If it's outside of the
-     * polygon we should skip processing it altogether.
-     */
-    static bool computePointCellRange(const VD::cell_type &cell, Point &start_source_point, Point &end_source_point, const VD::edge_type *&starting_vd_edge, const VD::edge_type *&ending_vd_edge, const std::vector<Segment> &segments);
-
     /*!
      * For VD cells associated with an input polygon vertex, we need to separate the node at the end and start of the cell into two
      * That way we can reach both the quad_start and the quad_end from the [incident_edge] of the two new nodes

src/libslic3r/Arachne/SkeletalTrapezoidationGraph.cpp

@@ -2,14 +2,18 @@
 //CuraEngine is released under the terms of the AGPLv3 or higher.
 
 #include "SkeletalTrapezoidationGraph.hpp"
-#include "../Line.hpp"
+
 #include <ankerl/unordered_dense.h>
-
-
 #include <boost/log/trivial.hpp>
+#include <algorithm>
+#include <iostream>
+#include <cassert>
+#include <cinttypes>
 
-#include "utils/linearAlg2D.hpp"
 #include "../Line.hpp"
+#include "libslic3r/Arachne/SkeletalTrapezoidationEdge.hpp"
+#include "libslic3r/Arachne/SkeletalTrapezoidationJoint.hpp"
+#include "libslic3r/Point.hpp"
 
 namespace Slic3r::Arachne
 {
@@ -182,8 +186,8 @@ bool STHalfEdgeNode::isLocalMaximum(bool strict) const
 
 void SkeletalTrapezoidationGraph::collapseSmallEdges(coord_t snap_dist)
 {
-    ankerl::unordered_dense::map<edge_t*, std::list<edge_t>::iterator> edge_locator;
-    ankerl::unordered_dense::map<node_t*, std::list<node_t>::iterator> node_locator;
+    ankerl::unordered_dense::map<edge_t*, Edges::iterator> edge_locator;
+    ankerl::unordered_dense::map<node_t*, Nodes::iterator> node_locator;
     
     for (auto edge_it = edges.begin(); edge_it != edges.end(); ++edge_it)
     {
@@ -195,7 +199,7 @@ void SkeletalTrapezoidationGraph::collapseSmallEdges(coord_t snap_dist)
         node_locator.emplace(&*node_it, node_it);
     }
     
-    auto safelyRemoveEdge = [this, &edge_locator](edge_t* to_be_removed, std::list<edge_t>::iterator& current_edge_it, bool& edge_it_is_updated)
+    auto safelyRemoveEdge = [this, &edge_locator](edge_t* to_be_removed, Edges::iterator& current_edge_it, bool& edge_it_is_updated)
     {
         if (current_edge_it != edges.end()
             && to_be_removed == &*current_edge_it)
@@ -315,8 +319,7 @@ void SkeletalTrapezoidationGraph::collapseSmallEdges(coord_t snap_dist)
     }
 }
 
-void SkeletalTrapezoidationGraph::makeRib(edge_t*& prev_edge, Point start_source_point, Point end_source_point, bool is_next_to_start_or_end)
-{
+void SkeletalTrapezoidationGraph::makeRib(edge_t *&prev_edge, const Point &start_source_point, const Point &end_source_point) {
     Point p;
     Line(start_source_point, end_source_point).distance_to_infinite_squared(prev_edge->to->p, &p);
     coord_t dist = (prev_edge->to->p - p).cast<int64_t>().norm();

src/libslic3r/Arachne/SkeletalTrapezoidationGraph.hpp

@@ -5,10 +5,20 @@
 #define SKELETAL_TRAPEZOIDATION_GRAPH_H
 
 #include <optional>
+#include <utility>
 
 #include "utils/HalfEdgeGraph.hpp"
 #include "SkeletalTrapezoidationEdge.hpp"
 #include "SkeletalTrapezoidationJoint.hpp"
+#include "libslic3r/Arachne/utils/HalfEdge.hpp"
+#include "libslic3r/Arachne/utils/HalfEdgeNode.hpp"
+#include "libslic3r/libslic3r.h"
+
+namespace Slic3r
+{
+class Line;
+class Point;
+};
 
 namespace Slic3r::Arachne
 {
@@ -83,7 +93,7 @@ public:
      */
     void collapseSmallEdges(coord_t snap_dist = 5);
 
-    void makeRib(edge_t*& prev_edge, Point start_source_point, Point end_source_point, bool is_next_to_start_or_end);
+    void makeRib(edge_t*& prev_edge, const Point &start_source_point, const Point &end_source_point);
 
     /*!
      * Insert a node into the graph and connect it to the input polygon using ribs

src/libslic3r/Arachne/WallToolPaths.cpp

@@ -361,7 +361,7 @@ void removeSmallAreas(Polygons &thiss, const double min_area_size, const bool re
         }
     } else {
         // For each polygon, computes the signed area, move small outlines at the end of the vector and keep pointer on small holes
-        std::vector<Polygon> small_holes;
+        Polygons small_holes;
         for (auto it = thiss.begin(); it < new_end;) {
             if (double area = ClipperLib::Area(to_path(*it)); fabs(area) < min_area_size) {
                 if (area >= 0) {
@@ -782,98 +782,4 @@ bool WallToolPaths::removeEmptyToolPaths(std::vector<VariableWidthLines> &toolpa
     return toolpaths.empty();
 }
 
-/*!
-     * Get the order constraints of the insets when printing walls per region / hole.
-     * Each returned pair consists of adjacent wall lines where the left has an inset_idx one lower than the right.
-     *
-     * Odd walls should always go after their enclosing wall polygons.
-     *
-     * \param outer_to_inner Whether the wall polygons with a lower inset_idx should go before those with a higher one.
- */
-WallToolPaths::ExtrusionLineSet WallToolPaths::getRegionOrder(const std::vector<ExtrusionLine *> &input, const bool outer_to_inner)
-{
-    ExtrusionLineSet order_requirements;
-    // We build a grid where we map toolpath vertex locations to toolpaths,
-    // so that we can easily find which two toolpaths are next to each other,
-    // which is the requirement for there to be an order constraint.
-    //
-    // We use a PointGrid rather than a LineGrid to save on computation time.
-    // In very rare cases two insets might lie next to each other without having neighboring vertices, e.g.
-    //  \            .
-    //   |  /        .
-    //   | /         .
-    //   ||          .
-    //   | \         .
-    //   |  \        .
-    //  /            .
-    // However, because of how Arachne works this will likely never be the case for two consecutive insets.
-    // On the other hand one could imagine that two consecutive insets of a very large circle
-    // could be simplify()ed such that the remaining vertices of the two insets don't align.
-    // In those cases the order requirement is not captured,
-    // which means that the PathOrderOptimizer *might* result in a violation of the user set path order.
-    // This problem is expected to be not so severe and happen very sparsely.
-
-    coord_t max_line_w = 0u;
-    for (const ExtrusionLine *line : input) // compute max_line_w
-        for (const ExtrusionJunction &junction : *line)
-            max_line_w = std::max(max_line_w, junction.w);
-    if (max_line_w == 0u)
-        return order_requirements;
-
-    struct LineLoc
-    {
-        ExtrusionJunction    j;
-        const ExtrusionLine *line;
-    };
-    struct Locator
-    {
-        Point operator()(const LineLoc &elem) { return elem.j.p; }
-    };
-
-    // How much farther two verts may be apart due to corners.
-    // This distance must be smaller than 2, because otherwise
-    // we could create an order requirement between e.g.
-    // wall 2 of one region and wall 3 of another region,
-    // while another wall 3 of the first region would lie in between those two walls.
-    // However, higher values are better against the limitations of using a PointGrid rather than a LineGrid.
-    constexpr float diagonal_extension = 1.9f;
-    const auto      searching_radius   = coord_t(max_line_w * diagonal_extension);
-    using GridT                        = SparsePointGrid<LineLoc, Locator>;
-    GridT grid(searching_radius);
-
-    for (const ExtrusionLine *line : input)
-        for (const ExtrusionJunction &junction : *line) grid.insert(LineLoc{junction, line});
-    for (const std::pair<const SquareGrid::GridPoint, LineLoc> &pair : grid) {
-        const LineLoc       &lineloc_here = pair.second;
-        const ExtrusionLine *here         = lineloc_here.line;
-        Point                loc_here     = pair.second.j.p;
-        std::vector<LineLoc> nearby_verts = grid.getNearby(loc_here, searching_radius);
-        for (const LineLoc &lineloc_nearby : nearby_verts) {
-            const ExtrusionLine *nearby = lineloc_nearby.line;
-            if (nearby == here)
-                continue;
-            if (nearby->inset_idx == here->inset_idx)
-                continue;
-            if (nearby->inset_idx > here->inset_idx + 1)
-                continue; // not directly adjacent
-            if (here->inset_idx > nearby->inset_idx + 1)
-                continue; // not directly adjacent
-            if (!shorter_then(loc_here - lineloc_nearby.j.p, (lineloc_here.j.w + lineloc_nearby.j.w) / 2 * diagonal_extension))
-                continue; // points are too far away from each other
-            if (here->is_odd || nearby->is_odd) {
-                if (here->is_odd && !nearby->is_odd && nearby->inset_idx < here->inset_idx)
-                    order_requirements.emplace(std::make_pair(nearby, here));
-                if (nearby->is_odd && !here->is_odd && here->inset_idx < nearby->inset_idx)
-                    order_requirements.emplace(std::make_pair(here, nearby));
-            } else if ((nearby->inset_idx < here->inset_idx) == outer_to_inner) {
-                order_requirements.emplace(std::make_pair(nearby, here));
-            } else {
-                assert((nearby->inset_idx > here->inset_idx) == outer_to_inner);
-                order_requirements.emplace(std::make_pair(here, nearby));
-            }
-        }
-    }
-    return order_requirements;
-}
-
 } // namespace Slic3r::Arachne

src/libslic3r/Arachne/WallToolPaths.hpp

@@ -90,16 +90,6 @@ public:
 
     using ExtrusionLineSet = ankerl::unordered_dense::set<std::pair<const ExtrusionLine *, const ExtrusionLine *>, boost::hash<std::pair<const ExtrusionLine *, const ExtrusionLine *>>>;
 
-    /*!
-     * Get the order constraints of the insets when printing walls per region / hole.
-     * Each returned pair consists of adjacent wall lines where the left has an inset_idx one lower than the right.
-     *
-     * Odd walls should always go after their enclosing wall polygons.
-     *
-     * \param outer_to_inner Whether the wall polygons with a lower inset_idx should go before those with a higher one.
-     */
-    static ExtrusionLineSet getRegionOrder(const std::vector<ExtrusionLine *> &input, bool outer_to_inner);
-
 protected:
     /*!
      * Stitch the polylines together and form closed polygons.

src/libslic3r/Arachne/utils/ExtrusionJunction.cpp

@@ -1,18 +0,0 @@
-//Copyright (c) 2020 Ultimaker B.V.
-//CuraEngine is released under the terms of the AGPLv3 or higher.
-
-#include "ExtrusionJunction.hpp"
-
-namespace Slic3r::Arachne
-{
-
-bool ExtrusionJunction::operator ==(const ExtrusionJunction& other) const
-{
-    return p == other.p
-        && w == other.w
-        && perimeter_index == other.perimeter_index;
-}
-
-ExtrusionJunction::ExtrusionJunction(const Point p, const coord_t w, const coord_t perimeter_index) : p(p), w(w), perimeter_index(perimeter_index) {}
-
-}

src/libslic3r/Arachne/utils/ExtrusionJunction.hpp

@@ -37,9 +37,11 @@ struct ExtrusionJunction
 	 */
     size_t perimeter_index;
 
-    ExtrusionJunction(const Point p, const coord_t w, const coord_t perimeter_index);
+    ExtrusionJunction(const Point p, const coord_t w, const coord_t perimeter_index) : p(p), w(w), perimeter_index(perimeter_index) {}
 
-    bool operator==(const ExtrusionJunction& other) const;
+    bool operator==(const ExtrusionJunction &other) const {
+        return p == other.p && w == other.w && perimeter_index == other.perimeter_index;
+    }
 
     coord_t x() const { return p.x(); }
     coord_t y() const { return p.y(); }

src/libslic3r/Arachne/utils/ExtrusionLine.cpp

@@ -2,10 +2,21 @@
 //CuraEngine is released under the terms of the AGPLv3 or higher.
 
 #include <algorithm>
+#include <cmath>
+#include <cstdlib>
 
 #include "ExtrusionLine.hpp"
-#include "linearAlg2D.hpp"
 #include "../../VariableWidth.hpp"
+#include "libslic3r/Arachne/utils/ExtrusionJunction.hpp"
+#include "libslic3r/BoundingBox.hpp"
+#include "libslic3r/ExtrusionEntity.hpp"
+#include "libslic3r/Line.hpp"
+#include "libslic3r/Polygon.hpp"
+#include "libslic3r/Polyline.hpp"
+
+namespace Slic3r {
+class Flow;
+}  // namespace Slic3r
 
 namespace Slic3r::Arachne
 {
@@ -29,15 +40,6 @@ int64_t ExtrusionLine::getLength() const
     return len;
 }
 
-coord_t ExtrusionLine::getMinimalWidth() const
-{
-    return std::min_element(junctions.cbegin(), junctions.cend(),
-                            [](const ExtrusionJunction& l, const ExtrusionJunction& r)
-                            {
-                                return l.w < r.w;
-                            })->w;
-}
-
 void ExtrusionLine::simplify(const int64_t smallest_line_segment_squared, const int64_t allowed_error_distance_squared, const int64_t maximum_extrusion_area_deviation)
 {
     const size_t min_path_size = is_closed ? 3 : 2;
@@ -262,9 +264,10 @@ bool ExtrusionLine::is_contour() const
     return poly.is_clockwise();
 }
 
-double ExtrusionLine::area() const
-{
-    assert(this->is_closed);
+double ExtrusionLine::area() const {
+    if (!this->is_closed)
+        return 0.;
+
     double a = 0.;
     if (this->junctions.size() >= 3) {
         Vec2d p1 = this->junctions.back().p.cast<double>();
@@ -274,9 +277,25 @@ double ExtrusionLine::area() const
             p1 = p2;
         }
     }
+
     return 0.5 * a;
 }
 
+Points to_points(const ExtrusionLine &extrusion_line) {
+    Points points;
+    points.reserve(extrusion_line.junctions.size());
+    for (const ExtrusionJunction &junction : extrusion_line.junctions)
+        points.emplace_back(junction.p);
+    return points;
+}
+
+BoundingBox get_extents(const ExtrusionLine &extrusion_line) {
+    BoundingBox bbox;
+    for (const ExtrusionJunction &junction : extrusion_line.junctions)
+        bbox.merge(junction.p);
+    return bbox;
+}
+
 } // namespace Slic3r::Arachne
 
 namespace Slic3r {

src/libslic3r/Arachne/utils/ExtrusionLine.hpp

@@ -5,16 +5,28 @@
 #ifndef UTILS_EXTRUSION_LINE_H
 #define UTILS_EXTRUSION_LINE_H
 
+#include <clipper/clipper_z.hpp>
+#include <assert.h>
+#include <stddef.h>
+#include <stdint.h>
+#include <algorithm>
+#include <utility>
+#include <vector>
+#include <cassert>
+#include <cinttypes>
+#include <cstddef>
+
 #include "ExtrusionJunction.hpp"
 #include "../../Polyline.hpp"
 #include "../../Polygon.hpp"
 #include "../../BoundingBox.hpp"
 #include "../../ExtrusionEntity.hpp"
 #include "../../Flow.hpp"
-#include "../../../clipper/clipper_z.hpp"
+#include "libslic3r/Point.hpp"
 
 namespace Slic3r {
 class ThickPolyline;
+class Flow;
 }
 
 namespace Slic3r::Arachne
@@ -136,11 +148,6 @@ struct ExtrusionLine
         return ret;
     }
 
-    /*!
-     * Get the minimal width of this path
-     */
-    coord_t getMinimalWidth() const;
-
     /*!
      * Removes vertices of the ExtrusionLines to make sure that they are not too high
      * resolution.
@@ -192,6 +199,8 @@ struct ExtrusionLine
     bool is_contour() const;
 
     double area() const;
+
+    bool is_external_perimeter() const { return this->inset_idx == 0; }
 };
 
 template<class PathType>
@@ -218,6 +227,7 @@ static inline Slic3r::ThickPolyline to_thick_polyline(const PathType &path)
 static inline Polygon to_polygon(const ExtrusionLine &line)
 {
     Polygon out;
+    assert(line.is_closed);
     assert(line.junctions.size() >= 3);
     assert(line.junctions.front().p == line.junctions.back().p);
     out.points.reserve(line.junctions.size() - 1);
@@ -226,24 +236,11 @@ static inline Polygon to_polygon(const ExtrusionLine &line)
     return out;
 }
 
-static Points to_points(const ExtrusionLine &extrusion_line)
-{
-    Points points;
-    points.reserve(extrusion_line.junctions.size());
-    for (const ExtrusionJunction &junction : extrusion_line.junctions)
-        points.emplace_back(junction.p);
-    return points;
-}
+Points to_points(const ExtrusionLine &extrusion_line);
+
+BoundingBox get_extents(const ExtrusionLine &extrusion_line);
 
 #if 0
-static BoundingBox get_extents(const ExtrusionLine &extrusion_line)
-{
-    BoundingBox bbox;
-    for (const ExtrusionJunction &junction : extrusion_line.junctions)
-        bbox.merge(junction.p);
-    return bbox;
-}
-
 static BoundingBox get_extents(const std::vector<ExtrusionLine> &extrusion_lines)
 {
     BoundingBox bbox;
@@ -274,6 +271,8 @@ static std::vector<Points> to_points(const std::vector<const ExtrusionLine *> &e
 #endif
 
 using VariableWidthLines = std::vector<ExtrusionLine>; //<! The ExtrusionLines generated by libArachne
+using Perimeter = VariableWidthLines;
+using Perimeters = std::vector<Perimeter>;
 
 } // namespace Slic3r::Arachne
 

src/libslic3r/Arachne/utils/HalfEdgeGraph.hpp

@@ -21,8 +21,10 @@ class HalfEdgeGraph
 public:
     using edge_t = derived_edge_t;
     using node_t = derived_node_t;
-    std::list<edge_t> edges;
-    std::list<node_t> nodes;
+    using Edges = std::list<edge_t>;
+    using Nodes = std::list<node_t>;
+    Edges edges;
+    Nodes nodes;
 };
 
 } // namespace Slic3r::Arachne

src/libslic3r/Arachne/utils/PolygonsPointIndex.hpp

@@ -156,8 +156,6 @@ struct PathsPointIndexLocator
     }
 };
 
-using PolygonsPointIndexLocator = PathsPointIndexLocator<Polygons>;
-
 }//namespace Slic3r::Arachne
 
 namespace std

src/libslic3r/Arachne/utils/PolylineStitcher.cpp

@@ -2,7 +2,16 @@
 //CuraEngine is released under the terms of the AGPLv3 or higher.
 
 #include "PolylineStitcher.hpp"
+
 #include "ExtrusionLine.hpp"
+#include "libslic3r/Arachne/utils/PolygonsPointIndex.hpp"
+#include "libslic3r/Polygon.hpp"
+
+namespace Slic3r {
+namespace Arachne {
+struct ExtrusionJunction;
+}  // namespace Arachne
+}  // namespace Slic3r
 
 namespace Slic3r::Arachne {
 

src/libslic3r/Arachne/utils/PolylineStitcher.hpp

@@ -4,11 +4,20 @@
 #ifndef UTILS_POLYLINE_STITCHER_H
 #define UTILS_POLYLINE_STITCHER_H
 
+#include <stddef.h>
+#include <stdint.h>
+#include <cassert>
+#include <functional>
+#include <limits>
+#include <vector>
+#include <cinttypes>
+#include <cstddef>
+
 #include "SparsePointGrid.hpp"
 #include "PolygonsPointIndex.hpp"
 #include "../../Polygon.hpp"
-#include <unordered_set>
-#include <cassert>
+#include "libslic3r/Point.hpp"
+#include "libslic3r/libslic3r.h"
 
 namespace Slic3r::Arachne
 {

src/libslic3r/Arachne/utils/SparseGrid.hpp

@@ -6,7 +6,6 @@
 #define UTILS_SPARSE_GRID_H
 
 #include <cassert>
-#include <unordered_map>
 #include <vector>
 #include <functional>
 

src/libslic3r/Arachne/utils/SparseLineGrid.hpp

@@ -6,7 +6,6 @@
 #define UTILS_SPARSE_LINE_GRID_H
 
 #include <cassert>
-#include <unordered_map>
 #include <vector>
 #include <functional>
 

src/libslic3r/Arachne/utils/SparsePointGrid.hpp

@@ -6,7 +6,6 @@
 #define UTILS_SPARSE_POINT_GRID_H
 
 #include <cassert>
-#include <unordered_map>
 #include <vector>
 
 #include "SparseGrid.hpp"
@@ -40,16 +39,6 @@ public:
      */
     void insert(const Elem &elem);
 
-    /*!
-     * Get just any element that's within a certain radius of a point.
-     *
-     * Rather than giving a vector of nearby elements, this function just gives
-     * a single element, any element, in no particular order.
-     * \param query_pt The point to query for an object nearby.
-     * \param radius The radius of what is considered "nearby".
-     */
-    const ElemT *getAnyNearby(const Point &query_pt, coord_t radius);
-
 protected:
     using GridPoint = typename SparseGrid<ElemT>::GridPoint;
 
@@ -69,22 +58,6 @@ void SparsePointGrid<ElemT, Locator>::insert(const Elem &elem)
     SparseGrid<ElemT>::m_grid.emplace(grid_loc, elem);
 }
 
-template<class ElemT, class Locator>
-const ElemT *SparsePointGrid<ElemT, Locator>::getAnyNearby(const Point &query_pt, coord_t radius)
-{
-    const ElemT                              *ret          = nullptr;
-    const std::function<bool(const ElemT &)> &process_func = [&ret, query_pt, radius, this](const ElemT &maybe_nearby) {
-        if (shorter_then(m_locator(maybe_nearby) - query_pt, radius)) {
-            ret = &maybe_nearby;
-            return false;
-        }
-        return true;
-    };
-    SparseGrid<ElemT>::processNearby(query_pt, radius, process_func);
-
-    return ret;
-}
-
 } // namespace Slic3r::Arachne
 
 #endif // UTILS_SPARSE_POINT_GRID_H

src/libslic3r/Arachne/utils/SquareGrid.cpp

@@ -2,7 +2,10 @@
 //CuraEngine is released under the terms of the AGPLv3 or higher.
 
 #include "SquareGrid.hpp"
-#include "../../Point.hpp"
+
+#include <cassert>
+
+#include "libslic3r/Point.hpp"
 
 using namespace Slic3r::Arachne;
 

src/libslic3r/Arachne/utils/SquareGrid.hpp

@@ -4,12 +4,15 @@
 #ifndef UTILS_SQUARE_GRID_H
 #define UTILS_SQUARE_GRID_H
 
-#include "../../Point.hpp"
-
+#include <stdint.h>
 #include <cassert>
-#include <unordered_map>
 #include <vector>
 #include <functional>
+#include <utility>
+#include <cinttypes>
+
+#include "../../Point.hpp"
+#include "libslic3r/libslic3r.h"
 
 namespace Slic3r::Arachne {
 

src/libslic3r/Arachne/utils/linearAlg2D.hpp

@@ -9,63 +9,6 @@
 namespace Slic3r::Arachne::LinearAlg2D
 {
 
-/*!
-     * Test whether a point is inside a corner.
-     * Whether point \p query_point is left of the corner abc.
-     * Whether the \p query_point is in the circle half left of ab and left of bc, rather than to the right.
-     *
-     * Test whether the \p query_point is inside of a polygon w.r.t a single corner.
- */
-inline static bool isInsideCorner(const Point &a, const Point &b, const Point &c, const Vec2i64 &query_point)
-{
-    //     Visualisation for the algorithm below:
-    //
-    //                 query
-    //                   |
-    //                   |
-    //                   |
-    //    perp-----------b
-    //                  / \       (note that the lines
-    //                 /   \      AB and AC are normalized
-    //                /     \     to 10000 units length)
-    //               a       c
-    //
-
-    auto normal = [](const Point &p0, coord_t len) -> Point {
-        int64_t _len = p0.norm();
-        if (_len < 1)
-            return {len, 0};
-        return (p0.cast<int64_t>() * int64_t(len) / _len).cast<coord_t>();
-    };
-
-    auto rotate_90_degree_ccw = [](const Vec2d &p) -> Vec2d {
-        return {-p.y(), p.x()};
-    };
-
-    constexpr coord_t normal_length = 10000; //Create a normal vector of reasonable length in order to reduce rounding error.
-    const Point ba = normal(a - b, normal_length);
-    const Point bc = normal(c - b, normal_length);
-    const Vec2d bq = query_point.cast<double>() - b.cast<double>();
-    const Vec2d perpendicular = rotate_90_degree_ccw(bq); //The query projects to this perpendicular to coordinate 0.
-
-    const double project_a_perpendicular = ba.cast<double>().dot(perpendicular); //Project vertex A on the perpendicular line.
-    const double project_c_perpendicular = bc.cast<double>().dot(perpendicular); //Project vertex C on the perpendicular line.
-    if ((project_a_perpendicular > 0.) != (project_c_perpendicular > 0.)) //Query is between A and C on the projection.
-    {
-        return project_a_perpendicular > 0.; //Due to the winding order of corner ABC, this means that the query is inside.
-    }
-    else //Beyond either A or C, but it could still be inside of the polygon.
-    {
-        const double project_a_parallel = ba.cast<double>().dot(bq); //Project not on the perpendicular, but on the original.
-        const double project_c_parallel = bc.cast<double>().dot(bq);
-
-        //Either:
-        // * A is to the right of B (project_a_perpendicular > 0) and C is below A (project_c_parallel < project_a_parallel), or
-        // * A is to the left of B (project_a_perpendicular < 0) and C is above A (project_c_parallel > project_a_parallel).
-        return (project_c_parallel < project_a_parallel) == (project_a_perpendicular > 0.);
-    }
-}
-
 /*!
      * Returns the determinant of the 2D matrix defined by the the vectors ab and ap as rows.
      * 

src/libslic3r/CMakeLists.txt

@@ -419,7 +419,6 @@ set(lisbslic3r_sources
     Arachne/BeadingStrategy/WideningBeadingStrategy.hpp
     Arachne/BeadingStrategy/WideningBeadingStrategy.cpp
     Arachne/utils/ExtrusionJunction.hpp
-    Arachne/utils/ExtrusionJunction.cpp
     Arachne/utils/ExtrusionLine.hpp
     Arachne/utils/ExtrusionLine.cpp
     Arachne/utils/HalfEdge.hpp
@@ -434,6 +433,8 @@ set(lisbslic3r_sources
     Arachne/utils/PolygonsSegmentIndex.hpp
     Arachne/utils/PolylineStitcher.hpp
     Arachne/utils/PolylineStitcher.cpp
+    Arachne/PerimeterOrder.hpp
+    Arachne/PerimeterOrder.cpp
     Arachne/SkeletalTrapezoidation.hpp
     Arachne/SkeletalTrapezoidation.cpp
     Arachne/SkeletalTrapezoidationEdge.hpp

src/libslic3r/Geometry.cpp

@@ -839,4 +839,50 @@ TransformationSVD::TransformationSVD(const Transform3d& trafo)
     return curMat;
 }
 
+bool is_point_inside_polygon_corner(const Point &a, const Point &b, const Point &c, const Point &query_point) {
+    // Cast all input points into int64_t to prevent overflows when points are close to max values of coord_t.
+    const Vec2i64 a_i64           = a.cast<int64_t>();
+    const Vec2i64 b_i64           = b.cast<int64_t>();
+    const Vec2i64 c_i64           = c.cast<int64_t>();
+    const Vec2i64 query_point_i64 = query_point.cast<int64_t>();
+
+    // Shift all points to have a base in vertex B.
+    // Then construct normalized vectors to ensure that we will work with vectors with endpoints on the unit circle.
+    const Vec2d ba = (a_i64 - b_i64).cast<double>().normalized();
+    const Vec2d bc = (c_i64 - b_i64).cast<double>().normalized();
+    const Vec2d bq = (query_point_i64 - b_i64).cast<double>().normalized();
+
+    // Points A and C has to be different.
+    assert(ba != bc);
+
+    // Construct a normal for the vector BQ that points to the left side of the vector BQ.
+    const Vec2d bq_left_normal = perp(bq);
+
+    const double proj_a_on_bq_normal = ba.dot(bq_left_normal); // Project point A on the normal of BQ.
+    const double proj_c_on_bq_normal = bc.dot(bq_left_normal); // Project point C on the normal of BQ.
+    if ((proj_a_on_bq_normal > 0. && proj_c_on_bq_normal <= 0.) || (proj_a_on_bq_normal <= 0. && proj_c_on_bq_normal > 0.)) {
+        // Q is between points A and C or lies on one of those vectors (BA or BC).
+
+        // Based on the CCW order of polygons (contours) and order of corner ABC,
+        // when this condition is met, the query point is inside the corner.
+        return proj_a_on_bq_normal > 0.;
+    } else {
+        // Q isn't between points A and C, but still it can be inside the corner.
+
+        const double proj_a_on_bq = ba.dot(bq); // Project point A on BQ.
+        const double proj_c_on_bq = bc.dot(bq); // Project point C on BQ.
+
+        // The value of proj_a_on_bq_normal is the same when we project the vector BA on the normal of BQ.
+        // So we can say that the Q is on the right side of the vector BA when proj_a_on_bq_normal > 0, and
+        // that the Q is on the left side of the vector BA proj_a_on_bq_normal < 0.
+        // Also, the Q is on the right side of the bisector of oriented angle ABC when proj_c_on_bq < proj_a_on_bq, and
+        // the Q is on the left side of the bisector of oriented angle ABC when proj_c_on_bq > proj_a_on_bq.
+
+        // So the Q is inside the corner when one of the following conditions is met:
+        //  * The Q is on the right side of the vector BA, and the Q is on the right side of the bisector of the oriented angle ABC.
+        //  * The Q is on the left side of the vector BA, and the Q is on the left side of the bisector of the oriented angle ABC.
+        return (proj_a_on_bq_normal > 0. && proj_c_on_bq < proj_a_on_bq) || (proj_a_on_bq_normal <= 0. && proj_c_on_bq >= proj_a_on_bq);
+    }
+}
+
 }} // namespace Slic3r::Geometry

src/libslic3r/Geometry.hpp

@@ -545,6 +545,23 @@ inline bool is_rotation_ninety_degrees(const Vec3d &rotation)
 }
 
 Transformation mat_around_a_point_rotate(const Transformation& innMat, const Vec3d &pt, const Vec3d &axis, float rotate_theta_radian);
+
+/**
+ * Checks if a given point is inside a corner of a polygon.
+ *
+ * The corner of a polygon is defined by three points A, B, C in counterclockwise order.
+ *
+ * Adapted from CuraEngine LinearAlg2D::isInsideCorner by Tim Kuipers @BagelOrb
+ * and @Ghostkeeper.
+ *
+ * @param a The first point of the corner.
+ * @param b The second point of the corner (the common vertex of the two edges forming the corner).
+ * @param c The third point of the corner.
+ * @param query_point The point to be checked if is inside the corner.
+ * @return True if the query point is inside the corner, false otherwise.
+ */
+bool is_point_inside_polygon_corner(const Point &a, const Point &b, const Point &c, const Point &query_point);
+
 } } // namespace Slicer::Geometry
 
 #endif

src/libslic3r/Geometry/MedialAxis.cpp

@@ -1,9 +1,16 @@
 #include <boost/log/trivial.hpp>
 #include "MedialAxis.hpp"
 
-#include "clipper.hpp"
+#include <boost/log/trivial.hpp>
+#include <boost/polygon/polygon.hpp>
+#include <cassert>
+#include <cmath>
+
 #include "VoronoiOffset.hpp"
-#include "../ClipperUtils.hpp"
+#include "libslic3r/ClipperUtils.hpp"
+#include "libslic3r/ExPolygon.hpp"
+#include "libslic3r/Point.hpp"
+#include "libslic3r/libslic3r.h"
 
 #ifdef SLIC3R_DEBUG
 namespace boost { namespace polygon {

src/libslic3r/Geometry/MedialAxis.hpp

@@ -1,8 +1,15 @@
 #ifndef slic3r_Geometry_MedialAxis_hpp_
 #define slic3r_Geometry_MedialAxis_hpp_
 
+#include <stddef.h>
+#include <utility>
+#include <vector>
+#include <cstddef>
+
 #include "Voronoi.hpp"
 #include "../ExPolygon.hpp"
+#include "libslic3r/Line.hpp"
+#include "libslic3r/Polyline.hpp"
 
 namespace Slic3r::Geometry {
 

src/libslic3r/Geometry/Voronoi.cpp

@@ -1,11 +1,13 @@
 #include "Voronoi.hpp"
 
+#include <boost/log/trivial.hpp>
+#include <cassert>
+
 #include "libslic3r/Arachne/utils/PolygonsSegmentIndex.hpp"
 #include "libslic3r/Geometry/VoronoiUtils.hpp"
 #include "libslic3r/Geometry/VoronoiUtilsCgal.hpp"
 #include "libslic3r/MultiMaterialSegmentation.hpp"
-
-#include <boost/log/trivial.hpp>
+#include "libslic3r/Line.hpp"
 
 namespace Slic3r::Geometry {
 
@@ -35,6 +37,8 @@ VoronoiDiagram::construct_voronoi(const SegmentIterator segment_begin, const Seg
                 BOOST_LOG_TRIVIAL(warning) << "Detected Voronoi edge intersecting input segment, input polygons will be rotated back and forth.";
             } else if (m_issue_type == IssueType::FINITE_EDGE_WITH_NON_FINITE_VERTEX) {
                 BOOST_LOG_TRIVIAL(warning) << "Detected finite Voronoi vertex with non finite vertex, input polygons will be rotated back and forth.";
+            } else if (m_issue_type == IssueType::PARABOLIC_VORONOI_EDGE_WITHOUT_FOCUS_POINT) {
+                BOOST_LOG_TRIVIAL(warning) << "Detected parabolic Voronoi edges without focus point, input polygons will be rotated back and forth.";
             } else {
                 BOOST_LOG_TRIVIAL(error) << "Detected unknown Voronoi diagram issue, input polygons will be rotated back and forth.";
             }
@@ -48,6 +52,8 @@ VoronoiDiagram::construct_voronoi(const SegmentIterator segment_begin, const Seg
                     BOOST_LOG_TRIVIAL(error) << "Detected Voronoi edge intersecting input segment even after the rotation of input.";
                 } else if (m_issue_type == IssueType::FINITE_EDGE_WITH_NON_FINITE_VERTEX) {
                     BOOST_LOG_TRIVIAL(error) << "Detected finite Voronoi vertex with non finite vertex even after the rotation of input.";
+                } else if (m_issue_type == IssueType::PARABOLIC_VORONOI_EDGE_WITHOUT_FOCUS_POINT) {
+                    BOOST_LOG_TRIVIAL(error) << "Detected parabolic Voronoi edges without focus point even after the rotation of input.";
                 } else {
                     BOOST_LOG_TRIVIAL(error) << "Detected unknown Voronoi diagram issue even after the rotation of input.";
                 }
@@ -159,8 +165,8 @@ typename boost::polygon::enable_if<
     VoronoiDiagram::IssueType>::type
 VoronoiDiagram::detect_known_issues(const VoronoiDiagram &voronoi_diagram, SegmentIterator segment_begin, SegmentIterator segment_end)
 {
-    if (has_finite_edge_with_non_finite_vertex(voronoi_diagram)) {
-        return IssueType::FINITE_EDGE_WITH_NON_FINITE_VERTEX;
+    if (const IssueType edge_issue_type = detect_known_voronoi_edge_issues(voronoi_diagram); edge_issue_type != IssueType::NO_ISSUE_DETECTED) {
+        return edge_issue_type;
     } else if (const IssueType cell_issue_type = detect_known_voronoi_cell_issues(voronoi_diagram, segment_begin, segment_end); cell_issue_type != IssueType::NO_ISSUE_DETECTED) {
         return cell_issue_type;
     } else if (!VoronoiUtilsCgal::is_voronoi_diagram_planar_angle(voronoi_diagram, segment_begin, segment_end)) {
@@ -218,16 +224,20 @@ VoronoiDiagram::detect_known_voronoi_cell_issues(const VoronoiDiagram &voronoi_d
     return IssueType::NO_ISSUE_DETECTED;
 }
 
-bool VoronoiDiagram::has_finite_edge_with_non_finite_vertex(const VoronoiDiagram &voronoi_diagram)
+VoronoiDiagram::IssueType VoronoiDiagram::detect_known_voronoi_edge_issues(const VoronoiDiagram &voronoi_diagram)
 {
     for (const voronoi_diagram_type::edge_type &edge : voronoi_diagram.edges()) {
         if (edge.is_finite()) {
             assert(edge.vertex0() != nullptr && edge.vertex1() != nullptr);
             if (edge.vertex0() == nullptr || edge.vertex1() == nullptr || !VoronoiUtils::is_finite(*edge.vertex0()) || !VoronoiUtils::is_finite(*edge.vertex1()))
-                return true;
+                return IssueType::FINITE_EDGE_WITH_NON_FINITE_VERTEX;
+
+            if (edge.is_curved() && !edge.cell()->contains_point() && !edge.twin()->cell()->contains_point())
+                return IssueType::PARABOLIC_VORONOI_EDGE_WITHOUT_FOCUS_POINT;
         }
     }
-    return false;
+
+    return IssueType::NO_ISSUE_DETECTED;
 }
 
 template<typename SegmentIterator>

src/libslic3r/Geometry/Voronoi.hpp

@@ -1,8 +1,15 @@
 #ifndef slic3r_Geometry_Voronoi_hpp_
 #define slic3r_Geometry_Voronoi_hpp_
 
+#include <boost/polygon/polygon.hpp>
+#include <cstddef>
+#include <iterator>
+#include <vector>
+
 #include "../Line.hpp"
 #include "../Polyline.hpp"
+#include "libslic3r/Point.hpp"
+#include "libslic3r/libslic3r.h"
 
 #ifdef _MSC_VER
 // Suppress warning C4146 in OpenVDB: unary minus operator applied to unsigned type, result still unsigned
@@ -10,6 +17,12 @@
 #pragma warning(disable : 4146)
 #endif // _MSC_VER
 #include "boost/polygon/voronoi.hpp"
+
+namespace boost {
+namespace polygon {
+template <typename Segment> struct segment_traits;
+}  // namespace polygon
+}  // namespace boost
 #ifdef _MSC_VER
 #pragma warning(pop)
 #endif // _MSC_VER
@@ -44,7 +57,8 @@ public:
         MISSING_VORONOI_VERTEX,
         NON_PLANAR_VORONOI_DIAGRAM,
         VORONOI_EDGE_INTERSECTING_INPUT_SEGMENT,
-        UNKNOWN                                  // Repairs are disabled in the constructor.
+        PARABOLIC_VORONOI_EDGE_WITHOUT_FOCUS_POINT,
+        UNKNOWN                                     // Repairs are disabled in the constructor.
     };
 
     enum class State {
@@ -158,7 +172,10 @@ private:
         IssueType>::type
     detect_known_voronoi_cell_issues(const VoronoiDiagram &voronoi_diagram, SegmentIterator segment_begin, SegmentIterator segment_end);
 
-    static bool has_finite_edge_with_non_finite_vertex(const VoronoiDiagram &voronoi_diagram);
+    // Detect issues related to Voronoi edges, or that can be detected by iterating over Voronoi edges.
+    // The first type of issue that can be detected is a finite Voronoi edge with a non-finite vertex.
+    // The second type of issue that can be detected is a parabolic Voronoi edge without a focus point (produced by two segments).
+    static IssueType detect_known_voronoi_edge_issues(const VoronoiDiagram &voronoi_diagram);
 
     voronoi_diagram_type  m_voronoi_diagram;
     vertex_container_type m_vertices;

src/libslic3r/Geometry/VoronoiOffset.cpp

@@ -1,15 +1,21 @@
 // Polygon offsetting using Voronoi diagram prodiced by boost::polygon.
 
-#include "Geometry.hpp"
-#include "VoronoiOffset.hpp"
-#include "libslic3r.h"
-
 #include <cmath>
+#include <algorithm>
+#include <array>
+#include <limits>
+#include <tuple>
+#include <utility>
+#include <cassert>
+#include <cstdlib>
+
+#include "libslic3r/Geometry.hpp"
+#include "VoronoiOffset.hpp"
+#include "libslic3r/libslic3r.h"
+#include "libslic3r/Geometry/Voronoi.hpp"
 
 // #define VORONOI_DEBUG_OUT
 
-#include <boost/polygon/detail/voronoi_ctypes.hpp>
-
 #ifdef VORONOI_DEBUG_OUT
 #include <libslic3r/VoronoiVisualUtils.hpp>
 #endif

src/libslic3r/Geometry/VoronoiOffset.hpp

@@ -3,9 +3,15 @@
 #ifndef slic3r_VoronoiOffset_hpp_
 #define slic3r_VoronoiOffset_hpp_
 
-#include "../libslic3r.h"
+#include <boost/polygon/polygon.hpp>
+#include <cmath>
+#include <vector>
 
+#include "libslic3r/libslic3r.h"
 #include "Voronoi.hpp"
+#include "libslic3r/Line.hpp"
+#include "libslic3r/Point.hpp"
+#include "libslic3r/Polygon.hpp"
 
 namespace Slic3r {
 

src/libslic3r/Geometry/VoronoiUtils.cpp

@@ -1,10 +1,17 @@
 #include <boost/log/trivial.hpp>
-
-#include <Arachne/utils/PolygonsSegmentIndex.hpp>
-#include <MultiMaterialSegmentation.hpp>
+#include <libslic3r/Arachne/utils/PolygonsSegmentIndex.hpp>
+#include <libslic3r/MultiMaterialSegmentation.hpp>
+#include <libslic3r/Geometry.hpp>
+#include <cmath>
+#include <cstdint>
+#include <utility>
+#include <vector>
+#include <cassert>
+#include <cstdlib>
 
 #include "VoronoiUtils.hpp"
-#include "libslic3r.h"
+#include "libslic3r/Exception.hpp"
+#include "libslic3r/Line.hpp"
 
 namespace Slic3r::Geometry {
 
@@ -28,6 +35,7 @@ template SegmentCellRange<Point> VoronoiUtils::compute_segment_cell_range(const
 template SegmentCellRange<Point> VoronoiUtils::compute_segment_cell_range(const VoronoiDiagram::cell_type &, VD::SegmentIt, VD::SegmentIt);
 template SegmentCellRange<Point> VoronoiUtils::compute_segment_cell_range(const VoronoiDiagram::cell_type &, ColoredLinesConstIt, ColoredLinesConstIt);
 template SegmentCellRange<Point> VoronoiUtils::compute_segment_cell_range(const VoronoiDiagram::cell_type &, PolygonsSegmentIndexConstIt, PolygonsSegmentIndexConstIt);
+template PointCellRange<Point> VoronoiUtils::compute_point_cell_range(const VoronoiDiagram::cell_type &, PolygonsSegmentIndexConstIt, PolygonsSegmentIndexConstIt);
 template Points VoronoiUtils::discretize_parabola(const Point &, const Arachne::PolygonsSegmentIndex &, const Point &, const Point &, coord_t, float);
 template Arachne::PolygonsPointIndex VoronoiUtils::get_source_point_index(const VoronoiDiagram::cell_type &, PolygonsSegmentIndexConstIt, PolygonsSegmentIndexConstIt);
 
@@ -124,7 +132,7 @@ VoronoiUtils::discretize_parabola(const Point &source_point, const Segment &sour
     Point pxx;
     Line(a, b).distance_to_infinite_squared(source_point, &pxx);
     const Point   ppxx = pxx - source_point;
-    const coord_t d    = ppxx.norm();
+    const coord_t d    = ppxx.cast<int64_t>().norm();
 
     const Vec2d  rot           = perp(ppxx).cast<double>().normalized();
     const double rot_cos_theta = rot.x();
@@ -137,8 +145,8 @@ VoronoiUtils::discretize_parabola(const Point &source_point, const Segment &sour
     }
 
     const double marking_bound = atan(transitioning_angle * 0.5);
-    int64_t      msx           = -marking_bound * d; // projected marking_start
-    int64_t      mex           = marking_bound * d;  // projected marking_end
+    int64_t      msx           = -marking_bound * int64_t(d); // projected marking_start
+    int64_t      mex           = marking_bound * int64_t(d);  // projected marking_end
 
     const coord_t marking_start_end_h = msx * msx / (2 * d) + d / 2;
     Point         marking_start       = Point(coord_t(msx), marking_start_end_h).rotated(rot_cos_theta, rot_sin_theta) + pxx;
@@ -152,7 +160,7 @@ VoronoiUtils::discretize_parabola(const Point &source_point, const Segment &sour
     bool add_marking_start = msx * int64_t(dir) > int64_t(sx - px) * int64_t(dir) && msx * int64_t(dir) < int64_t(ex - px) * int64_t(dir);
     bool add_marking_end   = mex * int64_t(dir) > int64_t(sx - px) * int64_t(dir) && mex * int64_t(dir) < int64_t(ex - px) * int64_t(dir);
 
-    const Point apex     = Point(coord_t(0), coord_t(d / 2)).rotated(rot_cos_theta, rot_sin_theta) + pxx;
+    const Point apex     = Point(0, d / 2).rotated(rot_cos_theta, rot_sin_theta) + pxx;
     bool        add_apex = int64_t(sx - px) * int64_t(dir) < 0 && int64_t(ex - px) * int64_t(dir) > 0;
 
     assert(!add_marking_start || !add_marking_end || add_apex);
@@ -245,6 +253,62 @@ VoronoiUtils::compute_segment_cell_range(const VD::cell_type &cell, const Segmen
     return cell_range;
 }
 
+template<typename SegmentIterator>
+typename boost::polygon::enable_if<
+    typename boost::polygon::gtl_if<typename boost::polygon::is_segment_concept<
+        typename boost::polygon::geometry_concept<typename std::iterator_traits<SegmentIterator>::value_type>::type>::type>::type,
+    Geometry::PointCellRange<
+        typename boost::polygon::segment_point_type<typename std::iterator_traits<SegmentIterator>::value_type>::type>>::type
+VoronoiUtils::compute_point_cell_range(const VD::cell_type &cell, const SegmentIterator segment_begin, const SegmentIterator segment_end)
+{
+    using Segment        = typename std::iterator_traits<SegmentIterator>::value_type;
+    using Point          = typename boost::polygon::segment_point_type<Segment>::type;
+    using PointCellRange = PointCellRange<Point>;
+    using CoordType      = typename Point::coord_type;
+
+    const Point source_point = Geometry::VoronoiUtils::get_source_point(cell, segment_begin, segment_end);
+
+    // We want to ignore (by returning PointCellRange without assigned edge_begin and edge_end) cells outside the input polygon.
+    PointCellRange cell_range(source_point);
+
+    const VD::edge_type *edge = cell.incident_edge();
+    if (edge->is_infinite() || !is_in_range<CoordType>(*edge)) {
+        // Ignore infinite edges, because they only occur outside the polygon.
+        // Also ignore edges with endpoints that don't fit into CoordType, because such edges are definitely outside the polygon.
+        return cell_range;
+    }
+
+    const Arachne::PolygonsPointIndex source_point_idx = Geometry::VoronoiUtils::get_source_point_index(cell, segment_begin, segment_end);
+    const Point                       edge_v0          = Geometry::VoronoiUtils::to_point(edge->vertex0()).template cast<CoordType>();
+    const Point                       edge_v1          = Geometry::VoronoiUtils::to_point(edge->vertex1()).template cast<CoordType>();
+    const Point                       edge_query_point = (edge_v0 == source_point) ? edge_v1 : edge_v0;
+
+    // Check if the edge has another endpoint inside the corner of the polygon.
+    if (!Geometry::is_point_inside_polygon_corner(source_point_idx.prev().p(), source_point_idx.p(), source_point_idx.next().p(), edge_query_point)) {
+        // If the endpoint isn't inside the corner of the polygon, it means that
+        // the whole cell isn't inside the polygons, and we will ignore such cells.
+        return cell_range;
+    }
+
+    const Vec2i64 source_point_i64 = source_point.template cast<int64_t>();
+    edge = cell.incident_edge();
+    do {
+        assert(edge->is_finite());
+
+        if (Vec2i64 v1 = Geometry::VoronoiUtils::to_point(edge->vertex1()); v1 == source_point_i64) {
+            cell_range.edge_begin = edge->next();
+            cell_range.edge_end   = edge;
+        } else {
+            // FIXME @hejllukas: With Arachne, we don't support polygons with collinear edges,
+            //                   because with collinear edges we have to handle secondary edges.
+            //                   Such edges goes through the endpoints of the input segments.
+            assert((Geometry::VoronoiUtils::to_point(edge->vertex0()) == source_point_i64 || edge->is_primary()) && "Point cells must end in the point! They cannot cross the point with an edge, because collinear edges are not allowed in the input.");
+        }
+    } while (edge = edge->next(), edge != cell.incident_edge());
+
+    return cell_range;
+}
+
 Vec2i64 VoronoiUtils::to_point(const VD::vertex_type *vertex)
 {
     assert(vertex != nullptr);

src/libslic3r/Geometry/VoronoiUtils.hpp

@@ -1,8 +1,15 @@
 #ifndef slic3r_VoronoiUtils_hpp_
 #define slic3r_VoronoiUtils_hpp_
 
+#include <boost/polygon/polygon.hpp>
+#include <iterator>
+#include <limits>
+
 #include "libslic3r/Geometry/Voronoi.hpp"
 #include "libslic3r/Arachne/utils/PolygonsSegmentIndex.hpp"
+#include "libslic3r/Arachne/utils/PolygonsPointIndex.hpp"
+#include "libslic3r/Point.hpp"
+#include "libslic3r/libslic3r.h"
 
 using VD = Slic3r::Geometry::VoronoiDiagram;
 
@@ -11,15 +18,28 @@ namespace Slic3r::Geometry {
 // Represent trapezoid Voronoi cell around segment.
 template<typename PT> struct SegmentCellRange
 {
-    const PT             segment_start_point;  // The start point of the source segment of this cell.
-    const PT             segment_end_point;    // The end point of the source segment of this cell.
+    const PT             source_segment_start_point; // The start point of the source segment of this cell.
+    const PT             source_segment_end_point;   // The end point of the source segment of this cell.
+    const VD::edge_type *edge_begin = nullptr;       // The edge of the Voronoi diagram where the loop around the cell starts.
+    const VD::edge_type *edge_end   = nullptr;       // The edge of the Voronoi diagram where the loop around the cell ends.
+
+    SegmentCellRange() = delete;
+    explicit SegmentCellRange(const PT &source_segment_start_point, const PT &source_segment_end_point)
+        : source_segment_start_point(source_segment_start_point), source_segment_end_point(source_segment_end_point)
+    {}
+
+    bool is_valid() const { return edge_begin && edge_end && edge_begin != edge_end; }
+};
+
+// Represent trapezoid Voronoi cell around point.
+template<typename PT> struct PointCellRange
+{
+    const PT             source_point;  // The source point of this cell.
     const VD::edge_type *edge_begin = nullptr; // The edge of the Voronoi diagram where the loop around the cell starts.
     const VD::edge_type *edge_end   = nullptr; // The edge of the Voronoi diagram where the loop around the cell ends.
 
-    SegmentCellRange() = delete;
-    explicit SegmentCellRange(const PT &segment_start_point, const PT &segment_end_point)
-        : segment_start_point(segment_start_point), segment_end_point(segment_end_point)
-    {}
+    PointCellRange() = delete;
+    explicit PointCellRange(const PT &source_point) : source_point(source_point) {}
 
     bool is_valid() const { return edge_begin && edge_end && edge_begin != edge_end; }
 };
@@ -80,7 +100,7 @@ public:
      * are linked to the neighboring segments, so you can iterate over the
      * segments until you reach the last segment.
      *
-     * Adapted from CuraEngine VoronoiUtils::computePointCellRange by Tim Kuipers @BagelOrb,
+     * Adapted from CuraEngine VoronoiUtils::computeSegmentCellRange by Tim Kuipers @BagelOrb,
      * Jaime van Kessel @nallath, Remco Burema @rburema and @Ghostkeeper.
      *
      * @param cell The cell to compute the range of line segments for.
@@ -96,6 +116,33 @@ public:
             typename boost::polygon::segment_point_type<typename std::iterator_traits<SegmentIterator>::value_type>::type>>::type
     compute_segment_cell_range(const VD::cell_type &cell, SegmentIterator segment_begin, SegmentIterator segment_end);
 
+    /**
+     * Compute the range of line segments that surround a cell of the skeletal
+     * graph that belongs to a point on the medial axis.
+     *
+     * This should only be used on cells that belong to a corner in the skeletal
+     * graph, e.g. triangular cells, not trapezoid cells.
+     *
+     * The resulting line segments is just the first and the last segment. They
+     * are linked to the neighboring segments, so you can iterate over the
+     * segments until you reach the last segment.
+     *
+     * Adapted from CuraEngine VoronoiUtils::computePointCellRange by Tim Kuipers @BagelOrb
+     * Jaime van Kessel @nallath, Remco Burema @rburema and @Ghostkeeper.
+     *
+     * @param cell The cell to compute the range of line segments for.
+     * @param segment_begin Begin iterator for all edges of the input Polygons.
+     * @param segment_end End iterator for all edges of the input Polygons.
+     * @return Range of line segments that surround the cell.
+     */
+    template<typename SegmentIterator>
+    static typename boost::polygon::enable_if<
+        typename boost::polygon::gtl_if<typename boost::polygon::is_segment_concept<
+            typename boost::polygon::geometry_concept<typename std::iterator_traits<SegmentIterator>::value_type>::type>::type>::type,
+        Geometry::PointCellRange<
+            typename boost::polygon::segment_point_type<typename std::iterator_traits<SegmentIterator>::value_type>::type>>::type
+    compute_point_cell_range(const VD::cell_type &cell, SegmentIterator segment_begin, SegmentIterator segment_end);
+
     template<typename T> static bool is_in_range(double value)
     {
         return double(std::numeric_limits<T>::lowest()) <= value && value <= double(std::numeric_limits<T>::max());

src/libslic3r/Geometry/VoronoiUtilsCgal.cpp

@@ -3,13 +3,21 @@
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 #include <CGAL/Arr_segment_traits_2.h>
 #include <CGAL/Surface_sweep_2_algorithms.h>
+#include <boost/variant/get.hpp>
+#include <vector>
+#include <cassert>
 
 #include "libslic3r/Geometry/Voronoi.hpp"
 #include "libslic3r/Geometry/VoronoiUtils.hpp"
 #include "libslic3r/Arachne/utils/PolygonsSegmentIndex.hpp"
 #include "libslic3r/MultiMaterialSegmentation.hpp"
-
 #include "VoronoiUtilsCgal.hpp"
+#include "libslic3r/Line.hpp"
+#include "libslic3r/Point.hpp"
+
+namespace CGAL {
+class MP_Float;
+}  // namespace CGAL
 
 using VD = Slic3r::Geometry::VoronoiDiagram;
 

src/libslic3r/Geometry/VoronoiUtilsCgal.hpp

@@ -1,6 +1,9 @@
 #ifndef slic3r_VoronoiUtilsCgal_hpp_
 #define slic3r_VoronoiUtilsCgal_hpp_
 
+#include <boost/polygon/polygon.hpp>
+#include <iterator>
+
 #include "Voronoi.hpp"
 #include "../Arachne/utils/PolygonsSegmentIndex.hpp"
 

src/libslic3r/Geometry/VoronoiVisualUtils.hpp

@@ -1,3 +1,6 @@
+#ifndef slic3r_VoronoiVisualUtils_hpp_
+#define slic3r_VoronoiVisualUtils_hpp_
+
 #include <stack>
 
 #include <libslic3r/Geometry.hpp>
@@ -451,3 +454,5 @@ static inline void dump_voronoi_to_svg(
 }
 
 } // namespace Slic3r
+
+#endif // slic3r_VoronoiVisualUtils_hpp_