Initial port of Arachne from Cura.

src/libslic3r/Arachne/BeadingStrategy/BeadingStrategy.cpp

@@ -0,0 +1,62 @@
+//Copyright (c) 2020 Ultimaker B.V.
+//CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#include <cassert>
+
+#include "BeadingStrategy.hpp"
+#include "Point.hpp"
+
+namespace Slic3r::Arachne
+{
+
+BeadingStrategy::BeadingStrategy(coord_t optimal_width, coord_t default_transition_length, float transitioning_angle)
+    : optimal_width(optimal_width)
+    , default_transition_length(default_transition_length)
+    , transitioning_angle(transitioning_angle)
+{
+    name = "Unknown";
+}
+
+coord_t BeadingStrategy::getTransitioningLength(coord_t lower_bead_count) const
+{
+    if (lower_bead_count == 0)
+    {
+        return scaled<coord_t>(0.01);
+    }
+    return default_transition_length;
+}
+
+float BeadingStrategy::getTransitionAnchorPos(coord_t lower_bead_count) const
+{
+    coord_t lower_optimum = getOptimalThickness(lower_bead_count);
+    coord_t transition_point = getTransitionThickness(lower_bead_count);
+    coord_t upper_optimum = getOptimalThickness(lower_bead_count + 1);
+    return 1.0 - float(transition_point - lower_optimum) / float(upper_optimum - lower_optimum);
+}
+
+std::vector<coord_t> BeadingStrategy::getNonlinearThicknesses(coord_t lower_bead_count) const
+{
+    return std::vector<coord_t>();
+}
+
+std::string BeadingStrategy::toString() const
+{
+    return name;
+}
+
+coord_t BeadingStrategy::getDefaultTransitionLength() const
+{
+    return default_transition_length;
+}
+
+coord_t BeadingStrategy::getOptimalWidth() const
+{
+    return optimal_width;
+}
+
+double BeadingStrategy::getTransitioningAngle() const
+{
+    return transitioning_angle;
+}
+
+} // namespace Slic3r::Arachne

src/libslic3r/Arachne/BeadingStrategy/BeadingStrategy.hpp

@@ -0,0 +1,112 @@
+// Copyright (c) 2021 Ultimaker B.V.
+// CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#ifndef BEADING_STRATEGY_H
+#define BEADING_STRATEGY_H
+
+#include <memory>
+
+#include "../../libslic3r.h"
+
+namespace Slic3r::Arachne
+{
+
+template<typename T> constexpr T pi_div(const T div) { return static_cast<T>(M_PI) / div; }
+
+/*!
+ * Mostly virtual base class template.
+ * 
+ * Strategy for covering a given (constant) horizontal model thickness with a number of beads.
+ * 
+ * The beads may have different widths.
+ * 
+ * TODO: extend with printing order?
+ */
+class BeadingStrategy
+{
+public:
+    /*!
+     * The beading for a given horizontal model thickness.
+     */
+    struct Beading
+    {
+        coord_t total_thickness;
+        std::vector<coord_t> bead_widths; //! The line width of each bead from the outer inset inward
+        std::vector<coord_t> toolpath_locations; //! The distance of the toolpath location of each bead from the outline
+        coord_t left_over; //! The distance not covered by any bead; gap area.
+    };
+
+    BeadingStrategy(coord_t optimal_width, coord_t default_transition_length, float transitioning_angle = pi_div(3));
+
+    virtual ~BeadingStrategy() {}
+
+    /*!
+     * Retrieve the bead widths with which to cover a given thickness.
+     * 
+     * Requirement: Given a constant \p bead_count the output of each bead width must change gradually along with the \p thickness.
+     * 
+     * \note The \p bead_count might be different from the \ref BeadingStrategy::optimal_bead_count
+     */
+    virtual Beading compute(coord_t thickness, coord_t bead_count) const = 0;
+
+    /*!
+     * The ideal thickness for a given \param bead_count
+     */
+    virtual coord_t getOptimalThickness(coord_t bead_count) const = 0;
+
+    /*!
+     * The model thickness at which \ref BeadingStrategy::optimal_bead_count transitions from \p lower_bead_count to \p lower_bead_count + 1
+     */
+    virtual coord_t getTransitionThickness(coord_t lower_bead_count) const = 0;
+
+    /*!
+     * The number of beads should we ideally usefor a given model thickness
+     */
+    virtual coord_t getOptimalBeadCount(coord_t thickness) const = 0;
+
+    /*!
+     * The length of the transitioning region along the marked / significant regions of the skeleton.
+     * 
+     * Transitions are used to smooth out the jumps in integer bead count; the jumps turn into ramps with some incline defined by their length.
+     */
+    virtual coord_t getTransitioningLength(coord_t lower_bead_count) const;
+
+    /*!
+     * The fraction of the transition length to put between the lower end of the transition and the point where the unsmoothed bead count jumps.
+     * 
+     * Transitions are used to smooth out the jumps in integer bead count; the jumps turn into ramps which could be positioned relative to the jump location.
+     */
+    virtual float getTransitionAnchorPos(coord_t lower_bead_count) const;
+
+    /*!
+     * Get the locations in a bead count region where \ref BeadingStrategy::compute exhibits a bend in the widths.
+     * Ordered from lower thickness to higher.
+     * 
+     * This is used to insert extra support bones into the skeleton, so that the resulting beads in long trapezoids don't linearly change between the two ends.
+     */
+    virtual std::vector<coord_t> getNonlinearThicknesses(coord_t lower_bead_count) const;
+    
+    virtual std::string toString() const;
+
+    coord_t getOptimalWidth() const;
+    coord_t getDefaultTransitionLength() const;
+    double  getTransitioningAngle() const;
+
+protected:
+    std::string name;
+
+    coord_t optimal_width; //! Optimal bead width, nominal width off the walls in 'ideal' circumstances.
+
+    coord_t default_transition_length; //! The length of the region to smoothly transfer between bead counts
+
+    /*!
+     * The maximum angle between outline segments smaller than which we are going to add transitions
+     * Equals 180 - the "limit bisector angle" from the paper
+     */
+    double  transitioning_angle;
+};
+
+using BeadingStrategyPtr = std::unique_ptr<BeadingStrategy>;
+
+} // namespace Slic3r::Arachne
+#endif // BEADING_STRATEGY_H

src/libslic3r/Arachne/BeadingStrategy/BeadingStrategyFactory.cpp

@@ -0,0 +1,97 @@
+//Copyright (c) 2021 Ultimaker B.V.
+//CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#include "BeadingStrategyFactory.hpp"
+
+#include "LimitedBeadingStrategy.hpp"
+#include "CenterDeviationBeadingStrategy.hpp"
+#include "WideningBeadingStrategy.hpp"
+#include "DistributedBeadingStrategy.hpp"
+#include "RedistributeBeadingStrategy.hpp"
+#include "OuterWallInsetBeadingStrategy.hpp"
+
+#include <limits>
+#include <boost/log/trivial.hpp>
+
+namespace Slic3r::Arachne
+{
+
+coord_t getWeightedAverage(const coord_t preferred_bead_width_outer, const coord_t preferred_bead_width_inner, const coord_t max_bead_count)
+{
+    if(max_bead_count > preferred_bead_width_outer - preferred_bead_width_inner)
+    {
+        //The difference between outer and inner bead width would be spread out across so many lines that rounding errors would destroy the difference.
+        //Also catches the case of max_bead_count being "infinite" (max integer).
+        return (preferred_bead_width_outer + preferred_bead_width_inner) / 2;
+    }
+    if (max_bead_count > 2)
+    {
+        return ((preferred_bead_width_outer * 2) + preferred_bead_width_inner * (max_bead_count - 2)) / max_bead_count;
+    }
+    if (max_bead_count <= 0)
+    {
+        return preferred_bead_width_inner;
+    }
+    return preferred_bead_width_outer;
+}
+
+BeadingStrategyPtr BeadingStrategyFactory::makeStrategy
+(
+    const BeadingStrategyType type,
+    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_width
+)
+{
+    using std::make_unique;
+    using std::move;
+    const coord_t bar_preferred_wall_width = getWeightedAverage(preferred_bead_width_outer, preferred_bead_width_inner, max_bead_count);
+    BeadingStrategyPtr ret;
+    switch (type)
+    {
+        case BeadingStrategyType::Center:
+            ret = make_unique<CenterDeviationBeadingStrategy>(bar_preferred_wall_width, transitioning_angle, wall_split_middle_threshold, wall_add_middle_threshold);
+            break;
+        case BeadingStrategyType::Distributed:
+            ret = make_unique<DistributedBeadingStrategy>(bar_preferred_wall_width, preferred_transition_length, transitioning_angle, wall_split_middle_threshold, wall_add_middle_threshold, std::numeric_limits<int>::max());
+            break;
+        case BeadingStrategyType::InwardDistributed:
+            ret = make_unique<DistributedBeadingStrategy>(bar_preferred_wall_width, preferred_transition_length, transitioning_angle, wall_split_middle_threshold, wall_add_middle_threshold, inward_distributed_center_wall_count);
+            break;
+        default:
+            BOOST_LOG_TRIVIAL(error) << "Cannot make strategy!";
+            return nullptr;
+    }
+    
+    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 << ".";
+        ret = make_unique<WideningBeadingStrategy>(move(ret), min_feature_size, min_bead_width);
+    }
+    if (max_bead_count > 0)
+    {
+        BOOST_LOG_TRIVIAL(debug) << "Applying the Redistribute meta-strategy with outer-wall width = " << preferred_bead_width_outer << ", inner-wall width = " << preferred_bead_width_inner;
+        ret = make_unique<RedistributeBeadingStrategy>(preferred_bead_width_outer, preferred_bead_width_inner, minimum_variable_line_width, 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 << ".";
+        ret = make_unique<LimitedBeadingStrategy>(max_bead_count, move(ret));
+    }
+    
+    if (outer_wall_offset > 0)
+    {
+        BOOST_LOG_TRIVIAL(debug) << "Applying the OuterWallOffset meta-strategy with offset = " << outer_wall_offset << ".";
+        ret = make_unique<OuterWallInsetBeadingStrategy>(outer_wall_offset, move(ret));
+    }
+    return ret;
+}
+} // namespace Slic3r::Arachne

src/libslic3r/Arachne/BeadingStrategy/BeadingStrategyFactory.hpp

@@ -0,0 +1,36 @@
+// Copyright (c) 2021 Ultimaker B.V.
+// CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#ifndef BEADING_STRATEGY_FACTORY_H
+#define BEADING_STRATEGY_FACTORY_H
+
+#include "BeadingStrategy.hpp"
+#include "../../Point.hpp"
+
+namespace Slic3r::Arachne
+{
+
+class BeadingStrategyFactory
+{
+public:
+    static BeadingStrategyPtr makeStrategy
+    (
+        const BeadingStrategyType type,
+        const coord_t preferred_bead_width_outer = scaled<coord_t>(0.0005),
+        const coord_t preferred_bead_width_inner = scaled<coord_t>(0.0005),
+        const coord_t preferred_transition_length = scaled<coord_t>(0.0004),
+        const float transitioning_angle = M_PI / 4.0,
+        const bool print_thin_walls = false,
+        const coord_t min_bead_width = 0,
+        const coord_t min_feature_size = 0,
+        const double wall_split_middle_threshold = 0.5,
+        const double wall_add_middle_threshold = 0.5,
+        const coord_t max_bead_count = 0,
+        const coord_t outer_wall_offset = 0,
+        const int inward_distributed_center_wall_count = 2,
+        const double minimum_variable_line_width = 0.5
+    );
+};
+
+} // namespace Slic3r::Arachne
+#endif // BEADING_STRATEGY_FACTORY_H

src/libslic3r/Arachne/BeadingStrategy/CenterDeviationBeadingStrategy.cpp

@@ -0,0 +1,88 @@
+// Copyright (c) 2021 Ultimaker B.V.
+// CuraEngine is released under the terms of the AGPLv3 or higher.
+#include <algorithm>
+
+#include "CenterDeviationBeadingStrategy.hpp"
+
+namespace Slic3r::Arachne
+{
+CenterDeviationBeadingStrategy::CenterDeviationBeadingStrategy(const coord_t pref_bead_width,
+                                                               const double transitioning_angle,
+                                                               const double wall_split_middle_threshold,
+                                                               const double wall_add_middle_threshold)
+  : BeadingStrategy(pref_bead_width, pref_bead_width / 2, transitioning_angle),
+    minimum_line_width_split(pref_bead_width * wall_split_middle_threshold),
+    minimum_line_width_add(pref_bead_width * wall_add_middle_threshold)
+{
+    name = "CenterDeviationBeadingStrategy";
+}
+
+CenterDeviationBeadingStrategy::Beading CenterDeviationBeadingStrategy::compute(coord_t thickness, coord_t bead_count) const
+{
+    Beading ret;
+
+    ret.total_thickness = thickness;
+    if (bead_count > 0)
+    {
+        // Set the bead widths
+        ret.bead_widths = std::vector<coord_t>(static_cast<size_t>(bead_count), optimal_width);
+        const coord_t optimal_thickness = getOptimalThickness(bead_count);
+        const coord_t diff_thickness = thickness - optimal_thickness; //Amount of deviation. Either spread out over the middle 2 lines, or concentrated in the center line.
+        const size_t center_bead_idx = ret.bead_widths.size() / 2;
+        if (bead_count % 2 == 0) // Even lines
+        {
+            const coord_t inner_bead_widths = optimal_width + diff_thickness / 2;
+            if (inner_bead_widths < minimum_line_width_add)
+            {
+                return compute(thickness, bead_count - 1);
+            }
+            ret.bead_widths[center_bead_idx - 1] = inner_bead_widths;
+            ret.bead_widths[center_bead_idx] = inner_bead_widths;
+        }
+        else // Uneven lines
+        {
+            const coord_t inner_bead_widths = optimal_width + diff_thickness;
+            if (inner_bead_widths < minimum_line_width_split)
+            {
+                return compute(thickness, bead_count - 1);
+            }
+            ret.bead_widths[center_bead_idx] = inner_bead_widths;
+        }
+
+        // Set the center line location of the bead toolpaths.
+        ret.toolpath_locations.resize(ret.bead_widths.size());
+        ret.toolpath_locations.front() = ret.bead_widths.front() / 2;
+        for (size_t bead_idx = 1; bead_idx < ret.bead_widths.size(); ++bead_idx)
+        {
+            ret.toolpath_locations[bead_idx] =
+              ret.toolpath_locations[bead_idx - 1] + (ret.bead_widths[bead_idx] + ret.bead_widths[bead_idx - 1]) / 2;
+        }
+        ret.left_over = 0;
+    }
+    else
+    {
+        ret.left_over = thickness;
+    }
+
+    return ret;
+}
+
+coord_t CenterDeviationBeadingStrategy::getOptimalThickness(coord_t bead_count) const
+{
+    return bead_count * optimal_width;
+}
+
+coord_t CenterDeviationBeadingStrategy::getTransitionThickness(coord_t lower_bead_count) const
+{
+    return lower_bead_count * optimal_width + (lower_bead_count % 2 == 1 ? minimum_line_width_split : minimum_line_width_add);
+}
+
+coord_t CenterDeviationBeadingStrategy::getOptimalBeadCount(coord_t thickness) const
+{
+    const coord_t naive_count = thickness / optimal_width; // How many lines we can fit in for sure.
+    const coord_t remainder = thickness - naive_count * optimal_width; // Space left after fitting that many lines.
+    const coord_t minimum_line_width = naive_count % 2 == 1 ? minimum_line_width_split : minimum_line_width_add;
+    return naive_count + (remainder > minimum_line_width); // If there's enough space, fit an extra one.
+}
+
+} // namespace Slic3r::Arachne

src/libslic3r/Arachne/BeadingStrategy/CenterDeviationBeadingStrategy.hpp

@@ -0,0 +1,42 @@
+// Copyright (c) 2021 Ultimaker B.V.
+// CuraEngine is released under the terms of the AGPLv3 or higher.
+
+
+#ifndef CENTER_DEVIATION_BEADING_STRATEGY_H
+#define CENTER_DEVIATION_BEADING_STRATEGY_H
+
+#include "BeadingStrategy.hpp"
+
+namespace Slic3r::Arachne
+{
+
+/*!
+ * This beading strategy makes the deviation in the thickness of the part
+ * entirely compensated by the innermost wall.
+ *
+ * The outermost walls all use the ideal width, as far as possible.
+ */
+class CenterDeviationBeadingStrategy : public BeadingStrategy
+{
+  private:
+    // For uneven numbers of lines: Minimum line width for which the middle line will be split into two lines.
+    coord_t minimum_line_width_split;
+
+    // For even numbers of lines: Minimum line width for which a new middle line will be added between the two innermost lines.
+    coord_t minimum_line_width_add;
+
+  public:
+    CenterDeviationBeadingStrategy(coord_t pref_bead_width,
+                                   double transitioning_angle,
+                                   double wall_split_middle_threshold,
+                                   double wall_add_middle_threshold);
+
+    ~CenterDeviationBeadingStrategy() override{};
+    Beading compute(coord_t thickness, coord_t bead_count) const override;
+    coord_t getOptimalThickness(coord_t bead_count) const override;
+    coord_t getTransitionThickness(coord_t lower_bead_count) const override;
+    coord_t getOptimalBeadCount(coord_t thickness) const override;
+};
+
+} // namespace Slic3r::Arachne
+#endif // CENTER_DEVIATION_BEADING_STRATEGY_H

src/libslic3r/Arachne/BeadingStrategy/DistributedBeadingStrategy.cpp

@@ -0,0 +1,110 @@
+// Copyright (c) 2021 Ultimaker B.V.
+// CuraEngine is released under the terms of the AGPLv3 or higher.
+#include <numeric>
+#include "DistributedBeadingStrategy.hpp"
+
+namespace Slic3r::Arachne
+{
+
+DistributedBeadingStrategy::DistributedBeadingStrategy(const coord_t optimal_width,
+                                const coord_t default_transition_length,
+                                const double transitioning_angle,
+                                const double wall_split_middle_threshold,
+                                const double wall_add_middle_threshold,
+                                const int distribution_radius)
+    : BeadingStrategy(optimal_width, default_transition_length, transitioning_angle)
+    , wall_split_middle_threshold(wall_split_middle_threshold)
+    , wall_add_middle_threshold(wall_add_middle_threshold)
+{
+    if(distribution_radius >= 2)
+    {
+        one_over_distribution_radius_squared = 1.0f / (distribution_radius - 1) * 1.0f / (distribution_radius - 1);
+    }
+    else
+    {
+        one_over_distribution_radius_squared = 1.0f / 1 * 1.0f / 1;
+    }
+    name = "DistributedBeadingStrategy";
+}
+
+DistributedBeadingStrategy::Beading DistributedBeadingStrategy::compute(coord_t thickness, coord_t bead_count) const
+{
+    Beading ret;
+
+    ret.total_thickness = thickness;
+    if (bead_count > 2)
+    {
+        const coord_t to_be_divided = thickness - bead_count * optimal_width;
+        const float middle = static_cast<float>(bead_count - 1) / 2;
+
+        const auto getWeight = [middle, this](coord_t bead_idx)
+        {
+            const float dev_from_middle = bead_idx - middle;
+            return std::max(0.0f, 1.0f - one_over_distribution_radius_squared * dev_from_middle * dev_from_middle);
+        };
+
+        std::vector<float> weights;
+        weights.resize(bead_count);
+        for (coord_t bead_idx = 0; bead_idx < bead_count; bead_idx++)
+        {
+            weights[bead_idx] = getWeight(bead_idx);
+        }
+
+        const float total_weight = std::accumulate(weights.cbegin(), weights.cend(), 0.f);
+        for (coord_t bead_idx = 0; bead_idx < bead_count; bead_idx++)
+        {
+            const float weight_fraction = weights[bead_idx] / total_weight;
+            const coord_t splitup_left_over_weight = to_be_divided * weight_fraction;
+            const coord_t width = optimal_width + splitup_left_over_weight;
+            if (bead_idx == 0)
+            {
+                ret.toolpath_locations.emplace_back(width / 2);
+            }
+            else
+            {
+                ret.toolpath_locations.emplace_back(ret.toolpath_locations.back() + (ret.bead_widths.back() + width) / 2);
+            }
+            ret.bead_widths.emplace_back(width);
+        }
+        ret.left_over = 0;
+    }
+    else if (bead_count == 2)
+    {
+        const coord_t outer_width = thickness / 2;
+        ret.bead_widths.emplace_back(outer_width);
+        ret.bead_widths.emplace_back(outer_width);
+        ret.toolpath_locations.emplace_back(outer_width / 2);
+        ret.toolpath_locations.emplace_back(thickness - outer_width / 2);
+        ret.left_over = 0;
+    }
+    else if (bead_count == 1)
+    {
+        const coord_t outer_width = thickness;
+        ret.bead_widths.emplace_back(outer_width);
+        ret.toolpath_locations.emplace_back(outer_width / 2);
+        ret.left_over = 0;
+    }
+    else
+    {
+        ret.left_over = thickness;
+    }
+
+    return ret;
+}
+
+coord_t DistributedBeadingStrategy::getOptimalThickness(coord_t bead_count) const
+{
+    return bead_count * optimal_width;
+}
+
+coord_t DistributedBeadingStrategy::getTransitionThickness(coord_t lower_bead_count) const
+{
+    return lower_bead_count * optimal_width + optimal_width * (lower_bead_count % 2 == 1 ? wall_split_middle_threshold : wall_add_middle_threshold);
+}
+
+coord_t DistributedBeadingStrategy::getOptimalBeadCount(coord_t thickness) const
+{
+    return (thickness + optimal_width / 2) / optimal_width;
+}
+
+} // namespace Slic3r::Arachne

src/libslic3r/Arachne/BeadingStrategy/DistributedBeadingStrategy.hpp

@@ -0,0 +1,48 @@
+// Copyright (c) 2021 Ultimaker B.V.
+// CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#ifndef DISTRIBUTED_BEADING_STRATEGY_H
+#define DISTRIBUTED_BEADING_STRATEGY_H
+
+#include "BeadingStrategy.hpp"
+
+namespace Slic3r::Arachne
+{
+
+/*!
+ * This beading strategy chooses a wall count that would make the line width
+ * deviate the least from the optimal line width, and then distributes the lines
+ * evenly among the thickness available.
+ */
+class DistributedBeadingStrategy : public BeadingStrategy
+{
+protected:
+     // For uneven numbers of lines: Minimum factor of the optimal width for which the middle line will be split into two lines.
+    double wall_split_middle_threshold;
+    
+    // For even numbers of lines: Minimum factor of the optimal width for which a new middle line will be added between the two innermost lines.
+    double wall_add_middle_threshold;
+
+    float one_over_distribution_radius_squared; // (1 / distribution_radius)^2
+
+public:
+    /*!
+    * \param distribution_radius the radius (in number of beads) over which to distribute the discrepancy between the feature size and the optimal thickness
+    */
+    DistributedBeadingStrategy( const coord_t optimal_width,
+                                const coord_t default_transition_length,
+                                const double transitioning_angle,
+                                const double wall_split_middle_threshold,
+                                const double wall_add_middle_threshold,
+                                const int distribution_radius);
+
+    virtual ~DistributedBeadingStrategy() override {}
+
+    Beading compute(coord_t thickness, coord_t bead_count) const override;
+    coord_t getOptimalThickness(coord_t bead_count) const override;
+    coord_t getTransitionThickness(coord_t lower_bead_count) const override;
+    coord_t getOptimalBeadCount(coord_t thickness) const override;
+};
+
+} // namespace Slic3r::Arachne
+#endif // DISTRIBUTED_BEADING_STRATEGY_H

src/libslic3r/Arachne/BeadingStrategy/LimitedBeadingStrategy.cpp

@@ -0,0 +1,136 @@
+//Copyright (c) 2020 Ultimaker B.V.
+//CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#include <cassert>
+#include <boost/log/trivial.hpp>
+
+#include "LimitedBeadingStrategy.hpp"
+#include "Point.hpp"
+
+namespace Slic3r::Arachne
+{
+
+std::string LimitedBeadingStrategy::toString() const
+{
+    return std::string("LimitedBeadingStrategy+") + parent->toString();
+}
+
+coord_t LimitedBeadingStrategy::getTransitioningLength(coord_t lower_bead_count) const
+{
+    return parent->getTransitioningLength(lower_bead_count);
+}
+
+float LimitedBeadingStrategy::getTransitionAnchorPos(coord_t lower_bead_count) const
+{
+    return parent->getTransitionAnchorPos(lower_bead_count);
+}
+
+LimitedBeadingStrategy::LimitedBeadingStrategy(const coord_t max_bead_count, BeadingStrategyPtr parent)
+    : BeadingStrategy(parent->getOptimalWidth(), /*default_transition_length=*/-1, parent->getTransitioningAngle())
+    , max_bead_count(max_bead_count)
+    , parent(std::move(parent))
+{
+    if (max_bead_count % 2 == 1)
+    {
+        BOOST_LOG_TRIVIAL(warning) << "LimitedBeadingStrategy with odd bead count is odd indeed!";
+    }
+}
+
+LimitedBeadingStrategy::Beading LimitedBeadingStrategy::compute(coord_t thickness, coord_t bead_count) const
+{
+    if (bead_count <= max_bead_count)
+    {
+        Beading ret = parent->compute(thickness, bead_count);
+        bead_count = ret.toolpath_locations.size();
+
+        if (bead_count % 2 == 0 && bead_count == max_bead_count)
+        {
+            const coord_t innermost_toolpath_location = ret.toolpath_locations[max_bead_count / 2 - 1];
+            const coord_t innermost_toolpath_width = ret.bead_widths[max_bead_count / 2 - 1];
+            ret.toolpath_locations.insert(ret.toolpath_locations.begin() + max_bead_count / 2, innermost_toolpath_location + innermost_toolpath_width / 2);
+            ret.bead_widths.insert(ret.bead_widths.begin() + max_bead_count / 2, 0);
+        }
+        return ret;
+    }
+    assert(bead_count == max_bead_count + 1);
+    if(bead_count != max_bead_count + 1)
+    {
+        BOOST_LOG_TRIVIAL(warning) << "Too many beads! " << bead_count << " != " << max_bead_count + 1;
+    }
+
+    coord_t optimal_thickness = parent->getOptimalThickness(max_bead_count);
+    Beading ret = parent->compute(optimal_thickness, max_bead_count);
+    bead_count = ret.toolpath_locations.size();
+    ret.left_over += thickness - ret.total_thickness;
+    ret.total_thickness = thickness;
+    
+    // Enforce symmetry
+    if (bead_count % 2 == 1)
+    {
+        ret.toolpath_locations[bead_count / 2] = thickness / 2;
+        ret.bead_widths[bead_count / 2] = thickness - optimal_thickness;
+    }
+    for (coord_t bead_idx = 0; bead_idx < (bead_count + 1) / 2; bead_idx++)
+    {
+        ret.toolpath_locations[bead_count - 1 - bead_idx] = thickness - ret.toolpath_locations[bead_idx];
+    }
+
+    //Create a "fake" inner wall with 0 width to indicate the edge of the walled area.
+    //This wall can then be used by other structures to e.g. fill the infill area adjacent to the variable-width walls.
+    coord_t innermost_toolpath_location = ret.toolpath_locations[max_bead_count / 2 - 1];
+    coord_t innermost_toolpath_width = ret.bead_widths[max_bead_count / 2 - 1];
+    ret.toolpath_locations.insert(ret.toolpath_locations.begin() + max_bead_count / 2, innermost_toolpath_location + innermost_toolpath_width / 2);
+    ret.bead_widths.insert(ret.bead_widths.begin() + max_bead_count / 2, 0);
+
+    //Symmetry on both sides. Symmetry is guaranteed since this code is stopped early if the bead_count <= max_bead_count, and never reaches this point then.
+    const size_t opposite_bead = bead_count - (max_bead_count / 2 - 1);
+    innermost_toolpath_location = ret.toolpath_locations[opposite_bead];
+    innermost_toolpath_width = ret.bead_widths[opposite_bead];
+    ret.toolpath_locations.insert(ret.toolpath_locations.begin() + opposite_bead, innermost_toolpath_location - innermost_toolpath_width / 2);
+    ret.bead_widths.insert(ret.bead_widths.begin() + opposite_bead, 0);
+
+    return ret;
+}
+
+coord_t LimitedBeadingStrategy::getOptimalThickness(coord_t bead_count) const
+{
+    if (bead_count <= max_bead_count)
+    {
+        return parent->getOptimalThickness(bead_count);
+    }
+    assert(false);
+    return scaled<coord_t>(1000.); // 1 meter (Cura was returning 10 meter)
+}
+
+coord_t LimitedBeadingStrategy::getTransitionThickness(coord_t lower_bead_count) const
+{
+    if (lower_bead_count < max_bead_count)
+    {
+        return parent->getTransitionThickness(lower_bead_count);
+    }
+    if (lower_bead_count == max_bead_count)
+    {
+        return parent->getOptimalThickness(lower_bead_count + 1) - scaled<coord_t>(0.01);
+    }
+    assert(false);
+    return scaled<coord_t>(900.); // 0.9 meter;
+}
+
+coord_t LimitedBeadingStrategy::getOptimalBeadCount(coord_t thickness) const
+{
+    coord_t parent_bead_count = parent->getOptimalBeadCount(thickness);
+    if (parent_bead_count <= max_bead_count)
+    {
+        return parent->getOptimalBeadCount(thickness);
+    }
+    else if (parent_bead_count == max_bead_count + 1)
+    {
+        if (thickness < parent->getOptimalThickness(max_bead_count + 1) - scaled<coord_t>(0.01))
+            return max_bead_count;
+        else 
+            return max_bead_count + 1;
+    }
+    else return max_bead_count + 1;
+}
+
+} // namespace Slic3r::Arachne

src/libslic3r/Arachne/BeadingStrategy/LimitedBeadingStrategy.hpp

@@ -0,0 +1,49 @@
+//Copyright (c) 2020 Ultimaker B.V.
+//CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#ifndef LIMITED_BEADING_STRATEGY_H
+#define LIMITED_BEADING_STRATEGY_H
+
+#include "BeadingStrategy.hpp"
+
+namespace Slic3r::Arachne
+{
+
+/*!
+ * This is a meta-strategy that can be applied on top of any other beading
+ * strategy, which limits the thickness of the walls to the thickness that the
+ * lines can reasonably print.
+ *
+ * The width of the wall is limited to the maximum number of contours times the
+ * maximum width of each of these contours.
+ *
+ * If the width of the wall gets limited, this strategy outputs one additional
+ * bead with 0 width. This bead is used to denote the limits of the walled area.
+ * Other structures can then use this border to align their structures to, such
+ * as to create correctly overlapping infill or skin, or to align the infill
+ * pattern to any extra infill walls.
+ */
+class LimitedBeadingStrategy : public BeadingStrategy
+{
+public:
+    LimitedBeadingStrategy(const coord_t max_bead_count, BeadingStrategyPtr parent);
+    
+    virtual ~LimitedBeadingStrategy() override = default;
+    
+    Beading compute(coord_t thickness, coord_t bead_count) const override;
+    coord_t getOptimalThickness(coord_t bead_count) const override;
+    coord_t getTransitionThickness(coord_t lower_bead_count) const override;
+    coord_t getOptimalBeadCount(coord_t thickness) const override;
+    virtual std::string toString() const override;
+    
+    coord_t getTransitioningLength(coord_t lower_bead_count) const override;
+    
+    float getTransitionAnchorPos(coord_t lower_bead_count) const override;
+
+protected:
+    const coord_t max_bead_count;
+    const BeadingStrategyPtr parent;
+};
+
+} // namespace Slic3r::Arachne
+#endif // LIMITED_DISTRIBUTED_BEADING_STRATEGY_H

src/libslic3r/Arachne/BeadingStrategy/OuterWallInsetBeadingStrategy.cpp

@@ -0,0 +1,62 @@
+//Copyright (c) 2020 Ultimaker B.V.
+//CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#include "OuterWallInsetBeadingStrategy.hpp"
+
+#include <algorithm>
+
+namespace Slic3r::Arachne
+{
+OuterWallInsetBeadingStrategy::OuterWallInsetBeadingStrategy(coord_t outer_wall_offset, BeadingStrategyPtr parent) :
+            BeadingStrategy(parent->getOptimalWidth(), parent->getDefaultTransitionLength(), parent->getTransitioningAngle()),
+            parent(std::move(parent)),
+            outer_wall_offset(outer_wall_offset)
+{
+    name = "OuterWallOfsetBeadingStrategy";
+}
+
+
+coord_t OuterWallInsetBeadingStrategy::getOptimalThickness(coord_t bead_count) const
+{
+    return parent->getOptimalThickness(bead_count);
+}
+
+coord_t OuterWallInsetBeadingStrategy::getTransitionThickness(coord_t lower_bead_count) const
+{
+    return parent->getTransitionThickness(lower_bead_count);
+}
+
+coord_t OuterWallInsetBeadingStrategy::getOptimalBeadCount(coord_t thickness) const
+{
+    return parent->getOptimalBeadCount(thickness);
+}
+
+coord_t OuterWallInsetBeadingStrategy::getTransitioningLength(coord_t lower_bead_count) const
+{
+    return parent->getTransitioningLength(lower_bead_count);
+}
+
+std::string OuterWallInsetBeadingStrategy::toString() const
+{
+    return std::string("OuterWallOfsetBeadingStrategy+") + parent->toString();
+}
+
+BeadingStrategy::Beading OuterWallInsetBeadingStrategy::compute(coord_t thickness, coord_t bead_count) const
+{
+    Beading ret = parent->compute(thickness, bead_count);
+
+    // Actual count and thickness as represented by extant walls. Don't count any potential zero-width 'signaling' walls.
+    bead_count = std::count_if(ret.bead_widths.begin(), ret.bead_widths.end(), [](const coord_t width) { return width > 0; });
+
+    // No need to apply any inset if there is just a single wall.
+    if (bead_count < 2)
+    {
+        return ret;
+    }
+
+    // Actually move the outer wall inside. Ensure that the outer wall never goes beyond the middle line.
+    ret.toolpath_locations[0] = std::min(ret.toolpath_locations[0] + outer_wall_offset, thickness / 2);
+    return ret;
+}
+
+} // namespace Slic3r::Arachne

src/libslic3r/Arachne/BeadingStrategy/OuterWallInsetBeadingStrategy.hpp

@@ -0,0 +1,35 @@
+//Copyright (c) 2020 Ultimaker B.V.
+//CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#ifndef OUTER_WALL_INSET_BEADING_STRATEGY_H
+#define OUTER_WALL_INSET_BEADING_STRATEGY_H
+
+#include "BeadingStrategy.hpp"
+
+namespace Slic3r::Arachne
+{
+    /*
+     * This is a meta strategy that allows for the outer wall to be inset towards the inside of the model. 
+     */
+    class OuterWallInsetBeadingStrategy : public BeadingStrategy
+    {
+    public:
+        OuterWallInsetBeadingStrategy(coord_t outer_wall_offset, BeadingStrategyPtr parent);
+         
+        virtual ~OuterWallInsetBeadingStrategy() = default;
+
+        Beading compute(coord_t thickness, coord_t bead_count) const override;
+        
+        coord_t getOptimalThickness(coord_t bead_count) const override;
+        coord_t getTransitionThickness(coord_t lower_bead_count) const override;
+        coord_t getOptimalBeadCount(coord_t thickness) const override;
+        coord_t getTransitioningLength(coord_t lower_bead_count) const override;
+        
+        std::string toString() const override;
+        
+    private:
+        BeadingStrategyPtr parent;
+        coord_t outer_wall_offset;
+    };
+} // namespace Slic3r::Arachne
+#endif // OUTER_WALL_INSET_BEADING_STRATEGY_H

src/libslic3r/Arachne/BeadingStrategy/RedistributeBeadingStrategy.cpp

@@ -0,0 +1,180 @@
+//Copyright (c) 2021 Ultimaker B.V.
+//CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#include "RedistributeBeadingStrategy.hpp"
+
+#include <algorithm>
+#include <numeric>
+
+namespace Slic3r::Arachne
+{
+
+RedistributeBeadingStrategy::RedistributeBeadingStrategy(   const coord_t optimal_width_outer,
+                                                            const coord_t optimal_width_inner,
+                                                            const double minimum_variable_line_width,
+                                                            BeadingStrategyPtr parent) :
+        BeadingStrategy(parent->getOptimalWidth(), parent->getDefaultTransitionLength(), parent->getTransitioningAngle()),
+        parent(std::move(parent)),
+        optimal_width_outer(optimal_width_outer),
+        optimal_width_inner(optimal_width_inner),
+        minimum_variable_line_width(minimum_variable_line_width)
+{
+    name = "RedistributeBeadingStrategy";
+}
+
+coord_t RedistributeBeadingStrategy::getOptimalThickness(coord_t bead_count) const
+{
+    const coord_t inner_bead_count = bead_count > 2 ? bead_count - 2 : 0;
+    const coord_t outer_bead_count = bead_count - inner_bead_count;
+
+    return parent->getOptimalThickness(inner_bead_count) + optimal_width_outer * outer_bead_count;
+}
+
+coord_t RedistributeBeadingStrategy::getTransitionThickness(coord_t lower_bead_count) const
+{
+    return parent->getTransitionThickness(lower_bead_count);
+}
+
+coord_t RedistributeBeadingStrategy::getOptimalBeadCount(coord_t thickness) const
+{
+    return parent->getOptimalBeadCount(thickness);
+}
+
+coord_t RedistributeBeadingStrategy::getTransitioningLength(coord_t lower_bead_count) const
+{
+    return parent->getTransitioningLength(lower_bead_count);
+}
+
+float RedistributeBeadingStrategy::getTransitionAnchorPos(coord_t lower_bead_count) const
+{
+    return parent->getTransitionAnchorPos(lower_bead_count);
+}
+
+std::string RedistributeBeadingStrategy::toString() const
+{
+    return std::string("RedistributeBeadingStrategy+") + parent->toString();
+}
+
+BeadingStrategy::Beading RedistributeBeadingStrategy::compute(coord_t thickness, coord_t bead_count) const
+{
+    Beading ret;
+    if (bead_count > 2)
+    {
+        const coord_t inner_transition_width = optimal_width_inner * minimum_variable_line_width;
+        const coord_t outer_bead_width =
+            getOptimalOuterBeadWidth(thickness, optimal_width_outer, inner_transition_width);
+
+        // Outer wall is locked in size and position for wall regions of 3 and higher which have at least a
+        // thickness equal to two times the optimal outer width and the minimal inner wall width.
+        const coord_t virtual_thickness = thickness - outer_bead_width * 2;
+        const coord_t virtual_bead_count = bead_count - 2;
+
+        // Calculate the beads and widths of the inner walls only
+        ret = parent->compute(virtual_thickness, virtual_bead_count);
+
+        // Insert the outer beads
+        ret.bead_widths.insert(ret.bead_widths.begin(), outer_bead_width);
+        ret.bead_widths.emplace_back(outer_bead_width);
+    }
+    else
+    {
+        ret = parent->compute(thickness, bead_count);
+    }
+
+    // Filter out beads that violate the minimum inner wall widths and recompute if necessary
+    const coord_t outer_transition_width = optimal_width_inner * minimum_variable_line_width;
+    const bool removed_inner_beads = validateInnerBeadWidths(ret, outer_transition_width);
+    if (removed_inner_beads)
+    {
+        ret = compute(thickness, bead_count - 1);
+    }
+
+    // Ensure that the positions of the beads are distributed over the thickness
+    resetToolPathLocations(ret, thickness);
+
+    return ret;
+}
+
+coord_t RedistributeBeadingStrategy::getOptimalOuterBeadWidth(const coord_t thickness, const coord_t optimal_width_outer, const coord_t minimum_width_inner)
+{
+    const coord_t total_outer_optimal_width = optimal_width_outer * 2;
+    coord_t outer_bead_width = thickness / 2;
+    if (total_outer_optimal_width < thickness)
+    {
+        if (total_outer_optimal_width + minimum_width_inner > thickness)
+        {
+            outer_bead_width -= minimum_width_inner / 2;
+        }
+        else
+        {
+            outer_bead_width = optimal_width_outer;
+        }
+    }
+    return outer_bead_width;
+}
+
+void RedistributeBeadingStrategy::resetToolPathLocations(BeadingStrategy::Beading& beading, const coord_t thickness)
+{
+    const size_t bead_count = beading.bead_widths.size();
+    beading.toolpath_locations.resize(bead_count);
+
+    if (bead_count < 1)
+    {
+        beading.toolpath_locations.resize(0);
+        beading.total_thickness = thickness;
+        beading.left_over = thickness;
+        return;
+    }
+
+    // Update the first half of the toolpath-locations with the updated bead-widths (starting from 0, up to half):
+    coord_t last_coord = 0;
+    coord_t last_width = 0;
+    for (size_t i_location = 0; i_location < bead_count / 2; ++i_location)
+    {
+        beading.toolpath_locations[i_location] = last_coord + (last_width + beading.bead_widths[i_location]) / 2;
+        last_coord = beading.toolpath_locations[i_location];
+        last_width = beading.bead_widths[i_location];
+    }
+
+    // Handle the position of any middle wall (note that the width will already have been set correctly):
+    if (bead_count % 2 == 1)
+    {
+        beading.toolpath_locations[bead_count / 2] = thickness / 2;
+    }
+
+    // Update the last half of the toolpath-locations with the updated bead-widths (starting from thickness, down to half):
+    last_coord = thickness;
+    last_width = 0;
+    for (size_t i_location = bead_count - 1; i_location >= bead_count - (bead_count / 2); --i_location)
+    {
+        beading.toolpath_locations[i_location] = last_coord - (last_width + beading.bead_widths[i_location]) / 2;
+        last_coord = beading.toolpath_locations[i_location];
+        last_width = beading.bead_widths[i_location];
+    }
+
+    // Ensure correct total and left over thickness
+    beading.total_thickness = thickness;
+    beading.left_over = thickness - std::accumulate(beading.bead_widths.cbegin(), beading.bead_widths.cend(), static_cast<coord_t>(0));
+}
+
+bool RedistributeBeadingStrategy::validateInnerBeadWidths(BeadingStrategy::Beading& beading, const coord_t minimum_width_inner)
+{
+    // Filter out bead_widths that violate the transition width and recalculate if needed
+    const size_t unfiltered_beads = beading.bead_widths.size();
+    if(unfiltered_beads <= 2) //Outer walls are exempt. If there are 2 walls the range below will be empty. If there is 1 or 0 walls it would be invalid.
+    {
+        return false;
+    }
+    auto inner_begin = std::next(beading.bead_widths.begin());
+    auto inner_end = std::prev(beading.bead_widths.end());
+    beading.bead_widths.erase(
+        std::remove_if(inner_begin, inner_end,
+        [&minimum_width_inner](const coord_t width)
+        {
+            return width < minimum_width_inner;
+        }),
+        inner_end);
+    return unfiltered_beads != beading.bead_widths.size();
+    }
+
+} // namespace Slic3r::Arachne

src/libslic3r/Arachne/BeadingStrategy/RedistributeBeadingStrategy.hpp

@@ -0,0 +1,98 @@
+//Copyright (c) 2020 Ultimaker B.V.
+//CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#ifndef REDISTRIBUTE_DISTRIBUTED_BEADING_STRATEGY_H
+#define REDISTRIBUTE_DISTRIBUTED_BEADING_STRATEGY_H
+
+#include "BeadingStrategy.hpp"
+
+namespace Slic3r::Arachne
+{
+    /*!
+     * A meta-beading-strategy that takes outer and inner wall widths into account.
+     *
+     * The outer wall will try to keep a constant width by only applying the beading strategy on the inner walls. This
+     * ensures that this outer wall doesn't react to changes happening to inner walls. It will limit print artifacts on
+     * the surface of the print. Although this strategy technically deviates from the original philosophy of the paper.
+     * It will generally results in better prints because of a smoother motion and less variation in extrusion width in
+     * the outer walls.
+     *
+     * If the thickness of the model is less then two times the optimal outer wall width and once the minimum inner wall
+     * width it will keep the minimum inner wall at a minimum constant and vary the outer wall widths symmetrical. Until
+     * The thickness of the model is that of at least twice the optimal outer wall width it will then use two
+     * symmetrical outer walls only. Until it transitions into a single outer wall. These last scenario's are always
+     * symmetrical in nature, disregarding the user specified strategy.
+     */
+    class RedistributeBeadingStrategy : public BeadingStrategy
+    {
+    public:
+        /*!
+         * /param optimal_width_outer         Outer wall width, guaranteed to be the actual (save rounding errors) at a
+         *                                    bead count if the parent strategies' optimum bead width is a weighted
+         *                                    average of the outer and inner walls at that bead count.
+         * /param optimal_width_outer         Inner wall width, guaranteed to be the actual (save rounding errors) at a
+         *                                    bead count if the parent strategies' optimum bead width is a weighted
+         *                                    average of the outer and inner walls at that bead count.
+         * /param minimum_variable_line_width Minimum factor that the variable line might deviate from the optimal width.
+         */
+        RedistributeBeadingStrategy(const coord_t optimal_width_outer,
+                                    const coord_t optimal_width_inner,
+                                    const double minimum_variable_line_width,
+                                    BeadingStrategyPtr parent);
+
+        virtual ~RedistributeBeadingStrategy() override = default;
+
+        Beading compute(coord_t thickness, coord_t bead_count) const override;
+
+        coord_t getOptimalThickness(coord_t bead_count) const override;
+        coord_t getTransitionThickness(coord_t lower_bead_count) const override;
+        coord_t getOptimalBeadCount(coord_t thickness) const override;
+        coord_t getTransitioningLength(coord_t lower_bead_count) const override;
+        float getTransitionAnchorPos(coord_t lower_bead_count) const override;
+
+        std::string toString() const override;
+
+    protected:
+        /*!
+         * Determine the outer bead width.
+         *
+         * According to the following logic:
+         *  - If the thickness of the model is more then twice the optimal outer bead width and the minimum inner bead
+         *  width it will return the optimal outer bead width.
+         *  - If the thickness is less then twice the optimal outer bead width and the minimum inner bead width, but
+         *  more them twice the optimal outer bead with it will return the optimal bead width minus half the inner bead
+         *  width.
+         *  - If the thickness is less then twice the optimal  outer bead width it will return half the thickness as
+         *  outer bead width
+         *
+         * \param thickness Thickness of the total beads.
+         * \param optimal_width_outer User specified optimal outer bead width.
+         * \param minimum_width_inner Inner bead width times the minimum variable line width.
+         * \return The outer bead width.
+         */
+        static coord_t getOptimalOuterBeadWidth(coord_t thickness, coord_t optimal_width_outer, coord_t minimum_width_inner);
+
+        /*!
+         * Moves the beads towards the outer edges of thickness and ensures that the outer walls are locked in location
+         * \param beading The beading instance.
+         * \param thickness The thickness of the bead.
+         */
+        static void resetToolPathLocations(Beading& beading, coord_t thickness);
+
+        /*!
+         * Filters and validates the beads, to ensure that all inner beads are at least the minimum bead width.
+         *
+         * \param beading The beading instance.
+         * \param minimum_width_inner Inner bead width times the minimum variable line width.
+         * \return true if beads are removed.
+         */
+        static bool validateInnerBeadWidths(Beading& beading, coord_t minimum_width_inner);
+
+        BeadingStrategyPtr parent;
+        coord_t optimal_width_outer;
+        coord_t optimal_width_inner;
+        double minimum_variable_line_width;
+    };
+
+} // namespace Slic3r::Arachne
+#endif // INWARD_DISTRIBUTED_BEADING_STRATEGY_H

src/libslic3r/Arachne/BeadingStrategy/WideningBeadingStrategy.cpp

@@ -0,0 +1,89 @@
+//Copyright (c) 2020 Ultimaker B.V.
+//CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#include "WideningBeadingStrategy.hpp"
+
+namespace Slic3r::Arachne
+{
+
+WideningBeadingStrategy::WideningBeadingStrategy(BeadingStrategyPtr parent, const coord_t min_input_width, const coord_t min_output_width)
+    : BeadingStrategy(parent->getOptimalWidth(), /*default_transition_length=*/-1, parent->getTransitioningAngle())
+    , parent(std::move(parent))
+    , min_input_width(min_input_width)
+    , min_output_width(min_output_width)
+{
+}
+
+std::string WideningBeadingStrategy::toString() const
+{
+    return std::string("Widening+") + parent->toString();
+}
+
+WideningBeadingStrategy::Beading WideningBeadingStrategy::compute(coord_t thickness, coord_t bead_count) const
+{
+    if (thickness < optimal_width)
+    {
+        Beading ret;
+        ret.total_thickness = thickness;
+        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 = thickness;
+        }
+        return ret;
+    }
+    else
+    {
+        return parent->compute(thickness, bead_count);
+    }
+}
+
+coord_t WideningBeadingStrategy::getOptimalThickness(coord_t bead_count) const
+{
+    return parent->getOptimalThickness(bead_count);
+}
+
+coord_t WideningBeadingStrategy::getTransitionThickness(coord_t lower_bead_count) const
+{
+    if (lower_bead_count == 0)
+    {
+        return min_input_width;
+    }
+    else
+    {
+        return parent->getTransitionThickness(lower_bead_count);
+    }
+}
+
+coord_t WideningBeadingStrategy::getOptimalBeadCount(coord_t thickness) const
+{
+    if (thickness < min_input_width) return 0;
+    coord_t ret = parent->getOptimalBeadCount(thickness);
+    if (thickness >= min_input_width && ret < 1) return 1;
+    return ret;
+}
+
+coord_t WideningBeadingStrategy::getTransitioningLength(coord_t lower_bead_count) const 
+{
+    return parent->getTransitioningLength(lower_bead_count);
+}
+
+float WideningBeadingStrategy::getTransitionAnchorPos(coord_t lower_bead_count) const
+{
+    return parent->getTransitionAnchorPos(lower_bead_count);
+}
+
+std::vector<coord_t> WideningBeadingStrategy::getNonlinearThicknesses(coord_t lower_bead_count) const
+{
+    std::vector<coord_t> ret;
+    ret.emplace_back(min_output_width);
+    std::vector<coord_t> pret = parent->getNonlinearThicknesses(lower_bead_count);
+    ret.insert(ret.end(), pret.begin(), pret.end());
+    return ret;
+}
+
+} // namespace Slic3r::Arachne

src/libslic3r/Arachne/BeadingStrategy/WideningBeadingStrategy.hpp

@@ -0,0 +1,46 @@
+//Copyright (c) 2020 Ultimaker B.V.
+//CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#ifndef WIDENING_BEADING_STRATEGY_H
+#define WIDENING_BEADING_STRATEGY_H
+
+#include "BeadingStrategy.hpp"
+
+namespace Slic3r::Arachne
+{
+
+/*!
+ * This is a meta-strategy that can be applied on any other beading strategy. If
+ * the part is thinner than a single line, this strategy adjusts the part so
+ * that it becomes the minimum thickness of one line.
+ *
+ * This way, tiny pieces that are smaller than a single line will still be
+ * printed.
+ */
+class WideningBeadingStrategy : public BeadingStrategy
+{
+public:
+    /*!
+     * Takes responsibility for deleting \param parent
+     */
+    WideningBeadingStrategy(BeadingStrategyPtr parent, const coord_t min_input_width, const coord_t min_output_width);
+
+    virtual ~WideningBeadingStrategy() override = default;
+
+    virtual Beading compute(coord_t thickness, coord_t bead_count) const override;
+    virtual coord_t getOptimalThickness(coord_t bead_count) const override;
+    virtual coord_t getTransitionThickness(coord_t lower_bead_count) const override;
+    virtual coord_t getOptimalBeadCount(coord_t thickness) const override;
+    virtual coord_t getTransitioningLength(coord_t lower_bead_count) const override;
+    virtual float getTransitionAnchorPos(coord_t lower_bead_count) const override;
+    virtual std::vector<coord_t> getNonlinearThicknesses(coord_t lower_bead_count) const override;
+    virtual std::string toString() const override;
+
+protected:
+    BeadingStrategyPtr parent;
+    const coord_t min_input_width;
+    const coord_t min_output_width;
+};
+
+} // namespace Slic3r::Arachne
+#endif // WIDENING_BEADING_STRATEGY_H

src/libslic3r/Arachne/SkeletalTrapezoidation.hpp

@@ -0,0 +1,597 @@
+//Copyright (c) 2020 Ultimaker B.V.
+//CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#ifndef SKELETAL_TRAPEZOIDATION_H
+#define SKELETAL_TRAPEZOIDATION_H
+
+#include <boost/polygon/voronoi.hpp>
+
+#include <memory> // smart pointers
+#include <unordered_map>
+#include <utility> // pair
+
+#include "utils/HalfEdgeGraph.hpp"
+#include "utils/PolygonsSegmentIndex.hpp"
+#include "utils/ExtrusionJunction.hpp"
+#include "utils/ExtrusionLine.hpp"
+#include "SkeletalTrapezoidationEdge.hpp"
+#include "SkeletalTrapezoidationJoint.hpp"
+#include "libslic3r/Arachne/BeadingStrategy/BeadingStrategy.hpp"
+#include "SkeletalTrapezoidationGraph.hpp"
+
+namespace Slic3r::Arachne
+{
+
+/*!
+ * Main class of the dynamic beading strategies.
+ * 
+ * The input polygon region is decomposed into trapezoids and represented as a half-edge data-structure.
+ * 
+ * We determine which edges are 'central' accordinding to the transitioning_angle of the beading strategy,
+ * and determine the bead count for these central regions and apply them outward when generating toolpaths. [oversimplified]
+ * 
+ * The method can be visually explained as generating the 3D union of cones surface on the outline polygons,
+ * and changing the heights along central regions of that surface so that they are flat. 
+ * For more info, please consult the paper "A framework for adaptive width control of dense contour-parallel toolpaths in fused
+deposition modeling" by Kuipers et al. 
+ * This visual explanation aid explains the use of "upward", "lower" etc,
+ * i.e. the radial distance and/or the bead count are used as heights of this visualization, there is no coordinate called 'Z'.
+ *
+ * TODO: split this class into two:
+ * 1. Class for generating the decomposition and aux functions for performing updates
+ * 2. Class for editing the structure for our purposes.
+ */
+class SkeletalTrapezoidation
+{
+    using pos_t = double;
+    using vd_t = boost::polygon::voronoi_diagram<pos_t>;
+    using graph_t = SkeletalTrapezoidationGraph;
+    using edge_t = STHalfEdge;
+    using node_t = STHalfEdgeNode;
+    using Beading = BeadingStrategy::Beading;
+    using BeadingPropagation = SkeletalTrapezoidationJoint::BeadingPropagation;
+    using TransitionMiddle = SkeletalTrapezoidationEdge::TransitionMiddle;
+    using TransitionEnd = SkeletalTrapezoidationEdge::TransitionEnd;
+
+    template<typename T>
+    using ptr_vector_t = std::vector<std::shared_ptr<T>>;
+
+    double  transitioning_angle; //!< How pointy a region should be before we apply the method. Equals 180* - limit_bisector_angle
+    coord_t discretization_step_size; //!< approximate size of segments when parabolic VD edges get discretized (and vertex-vertex edges)
+    coord_t transition_filter_dist; //!< Filter transition mids (i.e. anchors) closer together than this
+    coord_t beading_propagation_transition_dist; //!< When there are different beadings propagated from below and from above, use this transitioning distance
+    static constexpr coord_t central_filter_dist = scaled<coord_t>(0.02); //!< Filter areas marked as 'central' smaller than this
+    static constexpr coord_t snap_dist = scaled<coord_t>(0.02); //!< Generic arithmatic inaccuracy. Only used to determine whether a transition really needs to insert an extra edge.
+
+    /*!
+     * The strategy to use to fill a certain shape with lines.
+     *
+     * Various BeadingStrategies are available that differ in which lines get to
+     * print at their optimal width, where the play is being compensated, and
+     * how the joints are handled where we transition to different numbers of
+     * lines.
+     */
+    const BeadingStrategy& beading_strategy;
+
+public:
+    using Segment = PolygonsSegmentIndex;
+
+    /*!
+     * Construct a new trapezoidation problem to solve.
+     * \param polys The shapes to fill with walls.
+     * \param beading_strategy The strategy to use to fill these shapes.
+     * \param transitioning_angle Where we transition to a different number of
+     * walls, how steep should this transition be? A lower angle means that the
+     * transition will be longer.
+     * \param discretization_step_size Since g-code can't represent smooth
+     * transitions in line width, the line width must change with discretized
+     * steps. This indicates how long the line segments between those steps will
+     * be.
+     * \param transition_filter_dist The minimum length of transitions.
+     * Transitions shorter than this will be considered for dissolution.
+     * \param beading_propagation_transition_dist When there are different
+     * beadings propagated from below and from above, use this transitioning
+     * distance.
+     */
+    SkeletalTrapezoidation(const Polygons& polys, 
+                           const BeadingStrategy& beading_strategy,
+                           double transitioning_angle
+    , coord_t discretization_step_size = scaled<coord_t>(0.0008)
+    , coord_t transition_filter_dist = scaled<coord_t>(0.001)
+    , coord_t beading_propagation_transition_dist = scaled<coord_t>(0.0004));
+
+    /*!
+     * A skeletal graph through the polygons that we need to fill with beads.
+     *
+     * The skeletal graph represents the medial axes through each part of the
+     * polygons, and the lines from these medial axes towards each vertex of the
+     * polygons. The graph can be used to see what the width is of a polygon in
+     * each place and where the width transitions.
+     */
+    graph_t graph;
+
+    /*!
+     * Generate the paths that the printer must extrude, to print the outlines
+     * in the input polygons.
+     * \param filter_outermost_central_edges Some edges are "central" but still
+     * touch the outside of the polygon. If enabled, don't treat these as
+     * "central" but as if it's a obtuse corner. As a result, sharp corners will
+     * no longer end in a single line but will just loop.
+     */
+    void generateToolpaths(VariableWidthPaths& generated_toolpaths, bool filter_outermost_central_edges = false);
+
+protected:
+    /*!
+     * Auxiliary for referencing one transition along an edge which may contain multiple transitions
+     */
+    struct TransitionMidRef
+    {
+        edge_t* edge;
+        std::list<TransitionMiddle>::iterator transition_it;
+        TransitionMidRef(edge_t* edge, std::list<TransitionMiddle>::iterator transition_it)
+            : edge(edge)
+            , transition_it(transition_it)
+        {}
+    };
+
+    /*!
+     * Compute the skeletal trapezoidation decomposition of the input shape.
+     * 
+     * Compute the Voronoi Diagram (VD) and transfer all inside edges into our half-edge (HE) datastructure.
+     * 
+     * The algorithm is currently a bit overcomplicated, because the discretization of parabolic edges is performed at the same time as all edges are being transfered,
+     * which means that there is no one-to-one mapping from VD edges to HE edges.
+     * Instead we map from a VD edge to the last HE edge.
+     * This could be cimplified by recording the edges which should be discretized and discretizing the mafterwards.
+     * 
+     * Another complication arises because the VD uses floating logic, which can result in zero-length segments after rounding to integers.
+     * We therefore collapse edges and their whole cells afterwards.
+     */
+    void constructFromPolygons(const Polygons& polys);
+
+    /*!
+     * 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
+     */
+    std::unordered_map<vd_t::edge_type*, edge_t*> vd_edge_to_he_edge;
+    std::unordered_map<vd_t::vertex_type*, node_t*> vd_node_to_he_node;
+    node_t& makeNode(vd_t::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.
+
+    /*!
+     * (Eventual) returned 'polylines per index' result (from generateToolpaths):
+     */
+    VariableWidthPaths* p_generated_toolpaths;
+
+    /*!
+     * 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, vd_t::edge_type& vd_edge, edge_t*& prev_edge, Point& start_source_point, Point& end_source_point, const std::vector<Segment>& segments);
+
+    /*!
+     * Discretize a Voronoi edge that represents the medial axis of a vertex-
+     * line region or vertex-vertex region into small segments that can be
+     * considered to have a straight medial axis and a linear line width
+     * transition.
+     *
+     * The medial axis between a point and a line is a parabola. The rest of the
+     * algorithm doesn't want to have to deal with parabola, so this discretises
+     * the parabola into straight line segments. This is necessary if there is a
+     * sharp inner corner (acts as a point) that comes close to a straight edge.
+     *
+     * The medial axis between a point and a point is a straight line segment.
+     * However the distance from the medial axis to either of those points draws
+     * a parabola as you go along the medial axis. That means that the resulting
+     * line width along the medial axis would not be linearly increasing or
+     * linearly decreasing, but needs to take the shape of a parabola. Instead,
+     * we'll break this edge up into tiny line segments that can approximate the
+     * parabola with tiny linear increases or decreases in line width.
+     * \param segment The variable-width Voronoi edge to discretize.
+     * \param points All vertices of the original Polygons to fill with beads.
+     * \param segments All line segments of the original Polygons to fill with
+     * beads.
+     * \return A number of coordinates along the edge where the edge is broken
+     * up into discrete pieces.
+     */
+    std::vector<Point> discretize(const vd_t::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.
+     */
+    bool computePointCellRange(vd_t::cell_type& cell, Point& start_source_point, Point& end_source_point, vd_t::edge_type*& starting_vd_edge, vd_t::edge_type*& ending_vd_edge, const std::vector<Segment>& segments);
+
+    /*!
+     * Compute the range of line segments that surround a cell of the skeletal
+     * graph that belongs to a line segment of the medial axis.
+     *
+     * This should only be used on cells that belong to a central line segment
+     * of the skeletal graph, e.g. trapezoid cells, not triangular 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.
+     */
+    void computeSegmentCellRange(vd_t::cell_type& cell, Point& start_source_point, Point& end_source_point, vd_t::edge_type*& starting_vd_edge, vd_t::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
+     * Otherwise if node.incident_edge = quad_start you couldnt reach quad_end.twin by normal iteration (i.e. it = it.twin.next)
+     */
+    void separatePointyQuadEndNodes();
+
+
+    // ^ init | v transitioning
+
+    void updateIsCentral(); // Update the "is_central" flag for each edge based on the transitioning_angle
+
+    /*!
+     * Filter out small central areas.
+     * 
+     * Only used to get rid of small edges which get marked as central because
+     * of rounding errors because the region is so small.
+     */
+    void filterCentral(coord_t max_length);
+
+    /*!
+     * Filter central areas connected to starting_edge recursively.
+     * \return Whether we should unmark this section marked as central, on the
+     * way back out of the recursion.
+     */
+    bool filterCentral(edge_t* starting_edge, coord_t traveled_dist, coord_t max_length);
+
+    /*!
+     * Unmark the outermost edges directly connected to the outline, as not
+     * being central.
+     * 
+     * Only used to emulate some related literature.
+     * 
+     * The paper shows that this function is bad for the stability of the framework.
+     */
+    void filterOuterCentral();
+
+    /*!
+     * Set bead count in central regions based on the optimal_bead_count of the
+     * beading strategy.
+     */
+    void updateBeadCount();
+
+    /*!
+     * Add central regions and set bead counts where there is an end of the
+     * central area and when traveling upward we get to another region with the
+     * same bead count.
+     */
+    void filterNoncentralRegions();
+
+    /*!
+     * Add central regions and set bead counts for a particular edge and all of
+     * its adjacent edges.
+     *
+     * Recursive subroutine for \ref filterNoncentralRegions().
+     * \return Whether to set the bead count on the way back
+     */
+    bool filterNoncentralRegions(edge_t* to_edge, coord_t bead_count, coord_t traveled_dist, coord_t max_dist);
+
+    /*!
+     * Generate middle points of all transitions on edges.
+     *
+     * The transition middle points are saved in the graph itself. They are also
+     * returned via the output parameter.
+     * \param[out] edge_transitions A list of transitions that were generated.
+     */
+    void generateTransitionMids(ptr_vector_t<std::list<TransitionMiddle>>& edge_transitions);
+
+    /*!
+     * Removes some transition middle points.
+     *
+     * Transitions can be removed if there are multiple intersecting transitions
+     * that are too close together. If transitions have opposite effects, both
+     * are removed.
+     */
+    void filterTransitionMids();
+
+    /*!
+     * Merge transitions that are too close together.
+     * \param edge_to_start Edge pointing to the node from which to start
+     * traveling in all directions except along \p edge_to_start .
+     * \param origin_transition The transition for which we are checking nearby
+     * transitions.
+     * \param traveled_dist The distance traveled before we came to
+     * \p edge_to_start.to .
+     * \param going_up Whether we are traveling in the upward direction as seen
+     * from the \p origin_transition. If this doesn't align with the direction
+     * according to the R diff on a consecutive edge we know there was a local
+     * optimum.
+     * \return Whether the origin transition should be dissolved.
+     */
+    std::list<TransitionMidRef> dissolveNearbyTransitions(edge_t* edge_to_start, TransitionMiddle& origin_transition, coord_t traveled_dist, coord_t max_dist, bool going_up);
+
+    /*!
+     * Spread a certain bead count over a region in the graph.
+     * \param edge_to_start One edge of the region to spread the bead count in.
+     * \param from_bead_count All edges with this bead count will be changed.
+     * \param to_bead_count The new bead count for those edges.
+     */
+    void dissolveBeadCountRegion(edge_t* edge_to_start, coord_t from_bead_count, coord_t to_bead_count);
+
+    /*!
+     * Change the bead count if the given edge is at the end of a central
+     * region.
+     *
+     * This is necessary to provide a transitioning bead count to the edges of a
+     * central region to transition more smoothly from a high bead count in the
+     * central region to a lower bead count at the edge.
+     * \param edge_to_start One edge from a zone that needs to be filtered.
+     * \param traveled_dist The distance along the edges we've traveled so far.
+     * \param max_distance Don't filter beyond this range.
+     * \param replacing_bead_count The new bead count for this region.
+     * \return ``true`` if the bead count of this edge was changed.
+     */
+    bool filterEndOfCentralTransition(edge_t* edge_to_start, coord_t traveled_dist, coord_t max_dist, coord_t replacing_bead_count);
+
+    /*!
+     * Generate the endpoints of all transitions for all edges in the graph.
+     * \param[out] edge_transition_ends The resulting transition endpoints.
+     */
+    void generateAllTransitionEnds(ptr_vector_t<std::list<TransitionEnd>>& edge_transition_ends);
+
+    /*!
+     * Also set the rest values at nodes in between the transition ends
+     */
+    void applyTransitions(ptr_vector_t<std::list<TransitionEnd>>& edge_transition_ends);
+
+    /*!
+     * Create extra edges along all edges, where it needs to transition from one
+     * bead count to another.
+     *
+     * For example, if an edge of the graph goes from a bead count of 6 to a
+     * bead count of 1, it needs to generate 5 places where the beads around
+     * this line transition to a lower bead count. These are the "ribs". They
+     * reach from the edge to the border of the polygon. Where the beads hit
+     * those ribs the beads know to make a transition.
+     */
+    void generateTransitioningRibs();
+
+    /*!
+     * Generate the endpoints of a specific transition midpoint.
+     * \param edge The edge to create transitions on.
+     * \param mid_R The radius of the transition middle point.
+     * \param transition_lower_bead_count The bead count at the lower end of the
+     * transition.
+     * \param[out] edge_transition_ends A list of endpoints to add the new
+     * endpoints to.
+     */
+    void generateTransitionEnds(edge_t& edge, coord_t mid_R, coord_t transition_lower_bead_count, ptr_vector_t<std::list<TransitionEnd>>& edge_transition_ends);
+
+    /*!
+     * Compute a single endpoint of a transition.
+     * \param edge The edge to generate the endpoint for.
+     * \param start_pos The position where the transition starts.
+     * \param end_pos The position where the transition ends on the other side.
+     * \param transition_half_length The distance to the transition middle
+     * point.
+     * \param start_rest The gap between the start of the transition and the
+     * starting endpoint, as ratio of the inner bead width at the high end of
+     * the transition.
+     * \param end_rest The gap between the end of the transition and the ending
+     * endpoint, as ratio of the inner bead width at the high end of the
+     * transition.
+     * \param transition_lower_bead_count The bead count at the lower end of the
+     * transition.
+     * \param[out] edge_transition_ends The list to put the resulting endpoints
+     * in.
+     * \return Whether the given edge is going downward (i.e. towards a thinner
+     * region of the polygon).
+     */
+    bool generateTransitionEnd(edge_t& edge, coord_t start_pos, coord_t end_pos, coord_t transition_half_length, double start_rest, double end_rest, coord_t transition_lower_bead_count, ptr_vector_t<std::list<TransitionEnd>>& edge_transition_ends);
+
+    /*!
+     * Determines whether an edge is going downwards or upwards in the graph.
+     *
+     * An edge is said to go "downwards" if it's going towards a narrower part
+     * of the polygon. The notion of "downwards" comes from the conical
+     * representation of the graph, where the polygon is filled with a cone of
+     * maximum radius.
+     *
+     * This function works by recursively checking adjacent edges until the edge
+     * is reached.
+     * \param outgoing The edge to check.
+     * \param traveled_dist The distance traversed so far.
+     * \param transition_half_length The radius of the transition width.
+     * \param lower_bead_count The bead count at the lower end of the edge.
+     * \return ``true`` if this edge is going down, or ``false`` if it's going
+     * up.
+     */
+    bool isGoingDown(edge_t* outgoing, coord_t traveled_dist, coord_t transition_half_length, coord_t lower_bead_count) const;
+
+    /*!
+     * Determines whether this edge marks the end of the central region.
+     * \param edge The edge to check.
+     * \return ``true`` if this edge goes from a central region to a non-central
+     * region, or ``false`` in every other case (central to central, non-central
+     * to non-central, non-central to central, or end-of-the-line).
+     */
+    bool isEndOfCentral(const edge_t& edge) const;
+
+    /*!
+     * Create extra ribs in the graph where the graph contains a parabolic arc
+     * or a straight between two inner corners.
+     *
+     * There might be transitions there as the beads go through a narrow
+     * bottleneck in the polygon.
+     */
+    void generateExtraRibs();
+
+    // ^ transitioning ^
+
+    /*!
+     * It's useful to know when the paths get back to the consumer, to (what part of) a polygon the paths 'belong'.
+     * A single polygon without a hole is one region, a polygon with (a) hole(s) has 2 regions.
+     */
+    void markRegions();
+
+    // v toolpath generation v
+
+    /*!
+     * \param[out] segments the generated segments
+     */
+    void generateSegments();
+
+    /*!
+     * From a quad (a group of linked edges in one cell of the Voronoi), find
+     * the edge that is furthest away from the border of the polygon.
+     * \param quad_start_edge The first edge of the quad.
+     * \return The edge of the quad that is furthest away from the border.
+     */
+    edge_t* getQuadMaxRedgeTo(edge_t* quad_start_edge);
+
+    /*!
+     * Propagate beading information from nodes that are closer to the edge
+     * (low radius R) to nodes that are farther from the edge (high R).
+     * 
+     * only propagate from nodes with beading info upward to nodes without beading info
+     * 
+     * Edges are sorted by their radius, so that we can do a depth-first walk
+     * without employing a recursive algorithm.
+     * 
+     * In upward propagated beadings we store the distance traveled, so that we can merge these beadings with the downward propagated beadings in \ref propagateBeadingsDownward(.)
+     * 
+     * \param upward_quad_mids all upward halfedges of the inner skeletal edges (not directly connected to the outline) sorted on their highest [distance_to_boundary]. Higher dist first.
+     */
+    void propagateBeadingsUpward(std::vector<edge_t*>& upward_quad_mids, ptr_vector_t<BeadingPropagation>& node_beadings);
+
+    /*!
+     * propagate beading info from higher R nodes to lower R nodes
+     * 
+     * merge with upward propagated beadings if they are encountered
+     * 
+     * don't transfer to nodes which lie on the outline polygon
+     * 
+     * edges are sorted so that we can do a depth-first walk without employing a recursive algorithm
+     * 
+     * \param upward_quad_mids all upward halfedges of the inner skeletal edges (not directly connected to the outline) sorted on their highest [distance_to_boundary]. Higher dist first.
+     */
+    void propagateBeadingsDownward(std::vector<edge_t*>& upward_quad_mids, ptr_vector_t<BeadingPropagation>& node_beadings);
+
+    /*!
+     * Subroutine of \ref propagateBeadingsDownward(std::vector<edge_t*>&, ptr_vector_t<BeadingPropagation>&)
+     */
+    void propagateBeadingsDownward(edge_t* edge_to_peak, ptr_vector_t<BeadingPropagation>& node_beadings);
+
+    /*!
+     * Find a beading in between two other beadings.
+     *
+     * This creates a new beading. With this we can find the coordinates of the
+     * endpoints of the actual line segments to draw.
+     *
+     * The parameters \p left and \p right are not actually always left or right
+     * but just arbitrary directions to visually indicate the difference.
+     * \param left One of the beadings to interpolate between.
+     * \param ratio_left_to_whole The position within the two beadings to sample
+     * an interpolation. Should be a ratio between 0 and 1.
+     * \param right One of the beadings to interpolate between.
+     * \param switching_radius The bead radius at which we switch from the left
+     * beading to the merged beading, if the beadings have a different number of
+     * beads.
+     * \return The beading at the interpolated location.
+     */
+    Beading interpolate(const Beading& left, double ratio_left_to_whole, const Beading& right, coord_t switching_radius) const;
+
+    /*!
+     * Subroutine of \ref interpolate(const Beading&, Ratio, const Beading&, coord_t)
+     *
+     * This creates a new Beading between two beadings, assuming that both have
+     * the same number of beads.
+     * \param left One of the beadings to interpolate between.
+     * \param ratio_left_to_whole The position within the two beadings to sample
+     * an interpolation. Should be a ratio between 0 and 1.
+     * \param right One of the beadings to interpolate between.
+     * \return The beading at the interpolated location.
+     */
+    Beading interpolate(const Beading& left, double ratio_left_to_whole, const Beading& right) const;
+
+    /*!
+     * Get the beading at a certain node of the skeletal graph, or create one if
+     * it doesn't have one yet.
+     *
+     * This is a lazy get.
+     * \param node The node to get the beading from.
+     * \param node_beadings A list of all beadings for nodes.
+     * \return The beading of that node.
+     */
+    std::shared_ptr<BeadingPropagation> getOrCreateBeading(node_t* node, ptr_vector_t<BeadingPropagation>& node_beadings);
+
+    /*!
+     * In case we cannot find the beading of a node, get a beading from the
+     * nearest node.
+     * \param node The node to attempt to get a beading from. The actual node
+     * that the returned beading is from may be a different, nearby node.
+     * \param max_dist The maximum distance to search for.
+     * \return A beading for the node, or ``nullptr`` if there is no node nearby
+     * with a beading.
+     */
+    std::shared_ptr<BeadingPropagation> getNearestBeading(node_t* node, coord_t max_dist);
+
+    /*!
+     * generate junctions for each bone
+     * \param edge_to_junctions junctions ordered high R to low R
+     */
+    void generateJunctions(ptr_vector_t<BeadingPropagation>& node_beadings, ptr_vector_t<LineJunctions>& edge_junctions);
+
+    /*!
+     * add a new toolpath segment, defined between two extrusion-juntions
+     */
+    void addToolpathSegment(const ExtrusionJunction& from, const ExtrusionJunction& to, bool is_odd, bool force_new_path);
+
+    /*!
+     * connect junctions in each quad
+     */
+    void connectJunctions(ptr_vector_t<LineJunctions>& edge_junctions);
+
+    /*!
+     * Genrate small segments for local maxima where the beading would only result in a single bead
+     */
+    void generateLocalMaximaSingleBeads();
+
+    /*!
+     * Extract region information from the junctions, for easier access to that info directly from the lines.
+     */
+    void liftRegionInfoToLines();
+};
+
+} // namespace Slic3r::Arachne
+#endif // VORONOI_QUADRILATERALIZATION_H

src/libslic3r/Arachne/SkeletalTrapezoidationEdge.hpp

@@ -0,0 +1,142 @@
+//Copyright (c) 2021 Ultimaker B.V.
+//CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#ifndef SKELETAL_TRAPEZOIDATION_EDGE_H
+#define SKELETAL_TRAPEZOIDATION_EDGE_H
+
+#include <memory> // smart pointers
+#include <list>
+#include <vector>
+
+#include "utils/ExtrusionJunction.hpp"
+
+namespace Slic3r::Arachne
+{
+
+class SkeletalTrapezoidationEdge
+{
+private:
+    enum class Central { UNKNOWN = -1, NO, YES };
+
+public:
+    /*!
+     * Representing the location along an edge where the anchor position of a transition should be placed.
+     */
+    struct TransitionMiddle
+    {
+        coord_t pos; // Position along edge as measure from edge.from.p
+        int lower_bead_count;
+        TransitionMiddle(coord_t pos, int lower_bead_count)
+            : pos(pos), lower_bead_count(lower_bead_count)
+        {}
+    };
+
+    /*!
+     * Represents the location along an edge where the lower or upper end of a transition should be placed.
+     */
+    struct TransitionEnd
+    {
+        coord_t pos; // Position along edge as measure from edge.from.p, where the edge is always the half edge oriented from lower to higher R
+        int lower_bead_count;
+        bool is_lower_end; // Whether this is the ed of the transition with lower bead count
+        TransitionEnd(coord_t pos, int lower_bead_count, bool is_lower_end)
+            : pos(pos), lower_bead_count(lower_bead_count), is_lower_end(is_lower_end)
+        {}
+    };
+
+    enum class EdgeType
+    {
+        NORMAL = 0, // from voronoi diagram
+        EXTRA_VD = 1, // introduced to voronoi diagram in order to make the gMAT
+        TRANSITION_END = 2 // introduced to voronoi diagram in order to make the gMAT
+    };
+    EdgeType type;
+
+    SkeletalTrapezoidationEdge()
+    : SkeletalTrapezoidationEdge(EdgeType::NORMAL)
+    {}
+    SkeletalTrapezoidationEdge(const EdgeType& type)
+    : type(type)
+    , is_central(Central::UNKNOWN)
+    , region(0)
+    {}
+
+    bool isCentral() const
+    {
+        assert(is_central != Central::UNKNOWN);
+        return is_central == Central::YES;
+    }
+    void setIsCentral(bool b)
+    {
+        is_central = b ? Central::YES : Central::NO;
+    }
+    bool centralIsSet() const
+    {
+        return is_central != Central::UNKNOWN;
+    }
+
+    size_t getRegion() const
+    {
+        assert(region != 0);
+        return region;
+    }
+    void setRegion(const size_t& r)
+    {
+        assert(region == 0);
+        region = r;
+    }
+    bool regionIsSet() const
+    {
+        return region > 0;
+    }
+
+    bool hasTransitions(bool ignore_empty = false) const
+    {
+        return transitions.use_count() > 0 && (ignore_empty || ! transitions.lock()->empty());
+    }
+    void setTransitions(std::shared_ptr<std::list<TransitionMiddle>> storage)
+    {
+        transitions = storage;
+    }
+    std::shared_ptr<std::list<TransitionMiddle>> getTransitions()
+    {
+        return transitions.lock();
+    }
+
+    bool hasTransitionEnds(bool ignore_empty = false) const
+    {
+        return transition_ends.use_count() > 0 && (ignore_empty || ! transition_ends.lock()->empty());
+    }
+    void setTransitionEnds(std::shared_ptr<std::list<TransitionEnd>> storage)
+    {
+        transition_ends = storage;
+    }
+    std::shared_ptr<std::list<TransitionEnd>> getTransitionEnds()
+    {
+        return transition_ends.lock();
+    }
+
+    bool hasExtrusionJunctions(bool ignore_empty = false) const
+    {
+        return extrusion_junctions.use_count() > 0 && (ignore_empty || ! extrusion_junctions.lock()->empty());
+    }
+    void setExtrusionJunctions(std::shared_ptr<LineJunctions> storage)
+    {
+        extrusion_junctions = storage;
+    }
+    std::shared_ptr<LineJunctions> getExtrusionJunctions()
+    {
+        return extrusion_junctions.lock();
+    }
+
+private:
+    Central is_central; //! whether the edge is significant; whether the source segments have a sharp angle; -1 is unknown
+    size_t region; //! what 'region' this edge is in ... if the originating polygon has no holes, there's one region -- useful for later algorithms that need to know where the paths came from 
+
+    std::weak_ptr<std::list<TransitionMiddle>> transitions;
+    std::weak_ptr<std::list<TransitionEnd>> transition_ends;
+    std::weak_ptr<LineJunctions> extrusion_junctions;
+};
+
+} // namespace Slic3r::Arachne
+#endif // SKELETAL_TRAPEZOIDATION_EDGE_H

src/libslic3r/Arachne/SkeletalTrapezoidationGraph.cpp

@@ -0,0 +1,496 @@
+//Copyright (c) 2020 Ultimaker B.V.
+//CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#include "SkeletalTrapezoidationGraph.hpp"
+#include <unordered_map>
+
+#include <boost/log/trivial.hpp>
+
+#include "utils/linearAlg2D.hpp"
+#include "../Line.hpp"
+
+namespace Slic3r::Arachne
+{
+
+STHalfEdge::STHalfEdge(SkeletalTrapezoidationEdge data) : HalfEdge(data) {}
+
+bool STHalfEdge::canGoUp(bool strict) const
+{
+    if (to->data.distance_to_boundary > from->data.distance_to_boundary)
+    {
+        return true;
+    }
+    if (to->data.distance_to_boundary < from->data.distance_to_boundary || strict)
+    {
+        return false;
+    }
+
+    // Edge is between equidistqant verts; recurse!
+    for (edge_t* outgoing = next; outgoing != twin; outgoing = outgoing->twin->next)
+    {
+        if (outgoing->canGoUp())
+        {
+            return true;
+        }
+        assert(outgoing->twin); if (!outgoing->twin) return false;
+        assert(outgoing->twin->next); if (!outgoing->twin->next) return true; // This point is on the boundary?! Should never occur
+    }
+    return false;
+}
+
+bool STHalfEdge::isUpward() const
+{
+    if (to->data.distance_to_boundary > from->data.distance_to_boundary)
+    {
+        return true;
+    }
+    if (to->data.distance_to_boundary < from->data.distance_to_boundary)
+    {
+        return false;
+    }
+
+    // Equidistant edge case:
+    std::optional<coord_t> forward_up_dist = this->distToGoUp();
+    std::optional<coord_t> backward_up_dist = twin->distToGoUp();
+    if (forward_up_dist && backward_up_dist)
+    {
+        return forward_up_dist < backward_up_dist;
+    }
+
+    if (forward_up_dist)
+    {
+        return true;
+    }
+
+    if (backward_up_dist)
+    {
+        return false;
+    }
+    return to->p < from->p; // Arbitrary ordering, which returns the opposite for the twin edge
+}
+
+std::optional<coord_t> STHalfEdge::distToGoUp() const
+{
+    if (to->data.distance_to_boundary > from->data.distance_to_boundary)
+    {
+        return 0;
+    }
+    if (to->data.distance_to_boundary < from->data.distance_to_boundary)
+    {
+        return std::optional<coord_t>();
+    }
+
+    // Edge is between equidistqant verts; recurse!
+    std::optional<coord_t> ret;
+    for (edge_t* outgoing = next; outgoing != twin; outgoing = outgoing->twin->next)
+    {
+        std::optional<coord_t> dist_to_up = outgoing->distToGoUp();
+        if (dist_to_up)
+        {
+            if (ret)
+            {
+                ret = std::min(*ret, *dist_to_up);
+            }
+            else
+            {
+                ret = dist_to_up;
+            }
+        }
+        assert(outgoing->twin); if (!outgoing->twin) return std::optional<coord_t>();
+        assert(outgoing->twin->next); if (!outgoing->twin->next) return 0; // This point is on the boundary?! Should never occur
+    }
+    if (ret)
+    {
+        ret = *ret + (to->p - from->p).cast<int64_t>().norm();
+    }
+    return ret;
+}
+
+STHalfEdge* STHalfEdge::getNextUnconnected()
+{
+    edge_t* result = static_cast<STHalfEdge*>(this);
+    while (result->next)
+    {
+        result = result->next;
+        if (result == this)
+        {
+            return nullptr;
+        }
+    }
+    return result->twin;
+}
+
+STHalfEdgeNode::STHalfEdgeNode(SkeletalTrapezoidationJoint data, Point p) : HalfEdgeNode(data, p) {}
+
+bool STHalfEdgeNode::isMultiIntersection()
+{
+    int odd_path_count = 0;
+    edge_t* outgoing = this->incident_edge;
+    do
+    {
+        if (outgoing->data.isCentral())
+        {
+            odd_path_count++;
+        }
+    } while (outgoing = outgoing->twin->next, outgoing != this->incident_edge);
+    return odd_path_count > 2;
+}
+
+bool STHalfEdgeNode::isCentral() const
+{
+    edge_t* edge = incident_edge;
+    do
+    {
+        if (edge->data.isCentral())
+        {
+            return true;
+        }
+        assert(edge->twin); if (!edge->twin) return false;
+    } while (edge = edge->twin->next, edge != incident_edge);
+    return false;
+}
+
+bool STHalfEdgeNode::isLocalMaximum(bool strict) const
+{
+    if (data.distance_to_boundary == 0)
+    {
+        return false;
+    }
+
+    edge_t* edge = incident_edge;
+    do
+    {
+        if (edge->canGoUp(strict))
+        {
+            return false;
+        }
+        assert(edge->twin); if (!edge->twin) return false;
+
+        if (!edge->twin->next)
+        { // This point is on the boundary
+            return false;
+        }
+    } while (edge = edge->twin->next, edge != incident_edge);
+    return true;
+}
+
+void SkeletalTrapezoidationGraph::fixNodeDuplication()
+{
+    for (auto node_it = nodes.begin(); node_it != nodes.end();)
+    {
+        node_t* replacing_node = nullptr;
+        for (edge_t* outgoing = node_it->incident_edge; outgoing != node_it->incident_edge; outgoing = outgoing->twin->next)
+        {
+            assert(outgoing);
+            if (outgoing->from != &*node_it)
+            {
+                replacing_node = outgoing->from;
+            }
+            if (outgoing->twin->to != &*node_it)
+            {
+                replacing_node = outgoing->twin->to;
+            }
+        }
+        if (replacing_node)
+        {
+            for (edge_t* outgoing = node_it->incident_edge; outgoing != node_it->incident_edge; outgoing = outgoing->twin->next)
+            {
+                outgoing->twin->to = replacing_node;
+                outgoing->from = replacing_node;
+            }
+            node_it = nodes.erase(node_it);
+        }
+        else
+        {
+            ++node_it;
+        }
+    }
+}
+    
+void SkeletalTrapezoidationGraph::collapseSmallEdges(coord_t snap_dist)
+{
+    std::unordered_map<edge_t*, std::list<edge_t>::iterator> edge_locator;
+    std::unordered_map<node_t*, std::list<node_t>::iterator> node_locator;
+    
+    for (auto edge_it = edges.begin(); edge_it != edges.end(); ++edge_it)
+    {
+        edge_locator.emplace(&*edge_it, edge_it);
+    }
+    
+    for (auto node_it = nodes.begin(); node_it != nodes.end(); ++node_it)
+    {
+        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)
+    {
+        if (current_edge_it != edges.end()
+            && to_be_removed == &*current_edge_it)
+        {
+            current_edge_it = edges.erase(current_edge_it);
+            edge_it_is_updated = true;
+        }
+        else
+        {
+            edges.erase(edge_locator[to_be_removed]);
+        }
+    };
+
+    auto should_collapse = [snap_dist](node_t* a, node_t* b) 
+    { 
+        return shorter_then(a->p - b->p, snap_dist);
+    };
+        
+    for (auto edge_it = edges.begin(); edge_it != edges.end();)
+    {
+        if (edge_it->prev)
+        {
+            edge_it++;
+            continue;
+        }
+        
+        edge_t* quad_start = &*edge_it;
+        edge_t* quad_end = quad_start; while (quad_end->next) quad_end = quad_end->next;
+        edge_t* quad_mid = (quad_start->next == quad_end)? nullptr : quad_start->next;
+
+        bool edge_it_is_updated = false;
+        if (quad_mid && should_collapse(quad_mid->from, quad_mid->to))
+        {
+            assert(quad_mid->twin);
+            if(!quad_mid->twin)
+            {
+                BOOST_LOG_TRIVIAL(warning) << "Encountered quad edge without a twin.";
+                continue; //Prevent accessing unallocated memory.
+            }
+            int count = 0;
+            for (edge_t* edge_from_3 = quad_end; edge_from_3 && edge_from_3 != quad_mid->twin; edge_from_3 = edge_from_3->twin->next)
+            {
+                edge_from_3->from = quad_mid->from;
+                edge_from_3->twin->to = quad_mid->from;
+                if (count > 50)
+                {
+                    std::cerr << edge_from_3->from->p << " - " << edge_from_3->to->p << '\n';
+                }
+                if (++count > 1000) 
+                {
+                    break;
+                }
+            }
+
+            // o-o > collapse top
+            // | |
+            // | |
+            // | |
+            // o o
+            if (quad_mid->from->incident_edge == quad_mid)
+            {
+                if (quad_mid->twin->next)
+                {
+                    quad_mid->from->incident_edge = quad_mid->twin->next;
+                }
+                else
+                {
+                    quad_mid->from->incident_edge = quad_mid->prev->twin;
+                }
+            }
+            
+            nodes.erase(node_locator[quad_mid->to]);
+
+            quad_mid->prev->next = quad_mid->next;
+            quad_mid->next->prev = quad_mid->prev;
+            quad_mid->twin->next->prev = quad_mid->twin->prev;
+            quad_mid->twin->prev->next = quad_mid->twin->next;
+
+            safelyRemoveEdge(quad_mid->twin, edge_it, edge_it_is_updated);
+            safelyRemoveEdge(quad_mid, edge_it, edge_it_is_updated);
+        }
+
+        //  o-o
+        //  | | > collapse sides
+        //  o o
+        if ( should_collapse(quad_start->from, quad_end->to) && should_collapse(quad_start->to, quad_end->from))
+        { // Collapse start and end edges and remove whole cell
+
+            quad_start->twin->to = quad_end->to;
+            quad_end->to->incident_edge = quad_end->twin;
+            if (quad_end->from->incident_edge == quad_end)
+            {
+                if (quad_end->twin->next)
+                {
+                    quad_end->from->incident_edge = quad_end->twin->next;
+                }
+                else
+                {
+                    quad_end->from->incident_edge = quad_end->prev->twin;
+                }
+            }
+            nodes.erase(node_locator[quad_start->from]);
+
+            quad_start->twin->twin = quad_end->twin;
+            quad_end->twin->twin = quad_start->twin;
+            safelyRemoveEdge(quad_start, edge_it, edge_it_is_updated);
+            safelyRemoveEdge(quad_end, edge_it, edge_it_is_updated);
+        }
+        // If only one side had collapsable length then the cell on the other side of that edge has to collapse
+        // if we would collapse that one edge then that would change the quad_start and/or quad_end of neighboring cells
+        // this is to do with the constraint that !prev == !twin.next
+
+        if (!edge_it_is_updated)
+        {
+            edge_it++;
+        }
+    }
+}
+
+void SkeletalTrapezoidationGraph::makeRib(edge_t*& prev_edge, Point start_source_point, Point end_source_point, bool is_next_to_start_or_end)
+{
+    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();
+    prev_edge->to->data.distance_to_boundary = dist;
+    assert(dist >= 0);
+
+    nodes.emplace_front(SkeletalTrapezoidationJoint(), p);
+    node_t* node = &nodes.front();
+    node->data.distance_to_boundary = 0;
+    
+    edges.emplace_front(SkeletalTrapezoidationEdge(SkeletalTrapezoidationEdge::EdgeType::EXTRA_VD));
+    edge_t* forth_edge = &edges.front();
+    edges.emplace_front(SkeletalTrapezoidationEdge(SkeletalTrapezoidationEdge::EdgeType::EXTRA_VD));
+    edge_t* back_edge = &edges.front();
+    
+    prev_edge->next = forth_edge;
+    forth_edge->prev = prev_edge;
+    forth_edge->from = prev_edge->to;
+    forth_edge->to = node;
+    forth_edge->twin = back_edge;
+    back_edge->twin = forth_edge;
+    back_edge->from = node;
+    back_edge->to = prev_edge->to;
+    node->incident_edge = back_edge;
+    
+    prev_edge = back_edge;
+}
+
+std::pair<SkeletalTrapezoidationGraph::edge_t*, SkeletalTrapezoidationGraph::edge_t*> SkeletalTrapezoidationGraph::insertRib(edge_t& edge, node_t* mid_node)
+{
+    edge_t* edge_before = edge.prev;
+    edge_t* edge_after = edge.next;
+    node_t* node_before = edge.from;
+    node_t* node_after = edge.to;
+    
+    Point p = mid_node->p;
+
+    const Line source_segment = getSource(edge);
+    Point      px;
+    source_segment.distance_to_squared(p, &px);
+    coord_t dist = (p - px).cast<int64_t>().norm();
+    assert(dist > 0);
+    mid_node->data.distance_to_boundary = dist;
+    mid_node->data.transition_ratio = 0; // Both transition end should have rest = 0, because at the ends a whole number of beads fits without rest
+
+    nodes.emplace_back(SkeletalTrapezoidationJoint(), px);
+    node_t* source_node = &nodes.back();
+    source_node->data.distance_to_boundary = 0;
+
+    edge_t* first = &edge;
+    edges.emplace_back(SkeletalTrapezoidationEdge());
+    edge_t* second = &edges.back();
+    edges.emplace_back(SkeletalTrapezoidationEdge(SkeletalTrapezoidationEdge::EdgeType::TRANSITION_END));
+    edge_t* outward_edge = &edges.back();
+    edges.emplace_back(SkeletalTrapezoidationEdge(SkeletalTrapezoidationEdge::EdgeType::TRANSITION_END));
+    edge_t* inward_edge = &edges.back();
+
+    if (edge_before)
+    {
+        edge_before->next = first;
+    }
+    first->next = outward_edge;
+    outward_edge->next = nullptr;
+    inward_edge->next = second;
+    second->next = edge_after;
+
+    if (edge_after) 
+    {
+        edge_after->prev = second;
+    }
+    second->prev = inward_edge;
+    inward_edge->prev = nullptr;
+    outward_edge->prev = first;
+    first->prev = edge_before;
+
+    first->to = mid_node;
+    outward_edge->to = source_node;
+    inward_edge->to = mid_node;
+    second->to = node_after;
+
+    first->from = node_before;
+    outward_edge->from = mid_node;
+    inward_edge->from = source_node;
+    second->from = mid_node;
+
+    node_before->incident_edge = first;
+    mid_node->incident_edge = outward_edge;
+    source_node->incident_edge = inward_edge;
+    if (edge_after) 
+    {
+        node_after->incident_edge = edge_after;
+    }
+
+    first->data.setIsCentral(true);
+    outward_edge->data.setIsCentral(false); // TODO verify this is always the case.
+    inward_edge->data.setIsCentral(false);
+    second->data.setIsCentral(true);
+
+    outward_edge->twin = inward_edge;
+    inward_edge->twin = outward_edge;
+
+    first->twin = nullptr; // we don't know these yet!
+    second->twin = nullptr;
+
+    assert(second->prev->from->data.distance_to_boundary == 0);
+
+    return std::make_pair(first, second);
+}
+
+SkeletalTrapezoidationGraph::edge_t* SkeletalTrapezoidationGraph::insertNode(edge_t* edge, Point mid, coord_t mide_node_bead_count)
+{
+    edge_t* last_edge_replacing_input = edge;
+
+    nodes.emplace_back(SkeletalTrapezoidationJoint(), mid);
+    node_t* mid_node = &nodes.back();
+
+    edge_t* twin = last_edge_replacing_input->twin;
+    last_edge_replacing_input->twin = nullptr;
+    twin->twin = nullptr;
+    std::pair<edge_t*, edge_t*> left_pair = insertRib(*last_edge_replacing_input, mid_node);
+    std::pair<edge_t*, edge_t*> right_pair = insertRib(*twin, mid_node);
+    edge_t* first_edge_replacing_input = left_pair.first;
+    last_edge_replacing_input = left_pair.second;
+    edge_t* first_edge_replacing_twin = right_pair.first;
+    edge_t* last_edge_replacing_twin = right_pair.second;
+
+    first_edge_replacing_input->twin = last_edge_replacing_twin;
+    last_edge_replacing_twin->twin = first_edge_replacing_input;
+    last_edge_replacing_input->twin = first_edge_replacing_twin;
+    first_edge_replacing_twin->twin = last_edge_replacing_input;
+
+    mid_node->data.bead_count = mide_node_bead_count;
+
+    return last_edge_replacing_input;
+}
+
+Line SkeletalTrapezoidationGraph::getSource(const edge_t &edge) const
+{
+    const edge_t *from_edge = &edge;
+    while (from_edge->prev)
+        from_edge = from_edge->prev;
+
+    const edge_t *to_edge = &edge;
+    while (to_edge->next)
+        to_edge = to_edge->next;
+
+    return Line(from_edge->from->p, to_edge->to->p);
+}
+
+}

src/libslic3r/Arachne/SkeletalTrapezoidationGraph.hpp

@@ -0,0 +1,106 @@
+//Copyright (c) 2020 Ultimaker B.V.
+//CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#ifndef SKELETAL_TRAPEZOIDATION_GRAPH_H
+#define SKELETAL_TRAPEZOIDATION_GRAPH_H
+
+#include <optional>
+
+#include "utils/HalfEdgeGraph.hpp"
+#include "SkeletalTrapezoidationEdge.hpp"
+#include "SkeletalTrapezoidationJoint.hpp"
+
+namespace Slic3r::Arachne
+{
+
+class STHalfEdgeNode;
+
+class STHalfEdge : public HalfEdge<SkeletalTrapezoidationJoint, SkeletalTrapezoidationEdge, STHalfEdgeNode, STHalfEdge>
+{
+    using edge_t = STHalfEdge;
+    using node_t = STHalfEdgeNode;
+public:
+    STHalfEdge(SkeletalTrapezoidationEdge data);
+
+    /*!
+        * Check (recursively) whether there is any upward edge from the distance_to_boundary of the from of the \param edge
+        *
+        * \param strict Whether equidistant edges can count as a local maximum
+        */
+    bool canGoUp(bool strict = false) const;
+
+    /*!
+        * Check whether the edge goes from a lower to a higher distance_to_boundary.
+        * Effectively deals with equidistant edges by looking beyond this edge.
+        */
+    bool isUpward() const;
+
+    /*!
+        * Calculate the traversed distance until we meet an upward edge.
+        * Useful for calling on edges between equidistant points.
+        *
+        * If we can go up then the distance includes the length of the \param edge
+        */
+    std::optional<coord_t> distToGoUp() const;
+
+    STHalfEdge* getNextUnconnected();
+};
+
+class STHalfEdgeNode : public HalfEdgeNode<SkeletalTrapezoidationJoint, SkeletalTrapezoidationEdge, STHalfEdgeNode, STHalfEdge>
+{
+    using edge_t = STHalfEdge;
+    using node_t = STHalfEdgeNode;
+public:
+    STHalfEdgeNode(SkeletalTrapezoidationJoint data, Point p);
+
+    bool isMultiIntersection();
+
+    bool isCentral() const;
+
+    /*!
+        * Check whether this node has a locally maximal distance_to_boundary
+        *
+        * \param strict Whether equidistant edges can count as a local maximum
+        */
+    bool isLocalMaximum(bool strict = false) const;
+};
+
+class SkeletalTrapezoidationGraph: public HalfEdgeGraph<SkeletalTrapezoidationJoint, SkeletalTrapezoidationEdge, STHalfEdgeNode, STHalfEdge>
+{
+    using edge_t = STHalfEdge;
+    using node_t = STHalfEdgeNode;
+public:
+    void fixNodeDuplication();
+
+    /*!
+     * If an edge is too small, collapse it and its twin and fix the surrounding edges to ensure a consistent graph.
+     * 
+     * Don't collapse support edges, unless we can collapse the whole quad.
+     * 
+     * o-,
+     * |  "-o
+     * |    | > Don't collapse this edge only.
+     * o    o
+     */
+    void collapseSmallEdges(coord_t snap_dist = 5000);
+
+    void makeRib(edge_t*& prev_edge, Point start_source_point, Point end_source_point, bool is_next_to_start_or_end);
+
+    /*!
+     * Insert a node into the graph and connect it to the input polygon using ribs
+     * 
+     * \return the last edge which replaced [edge], which points to the same [to] node
+     */
+    edge_t* insertNode(edge_t* edge, Point mid, coord_t mide_node_bead_count);
+
+    /*!
+     * Return the first and last edge of the edges replacing \p edge pointing to the same node
+     */
+    std::pair<edge_t*, edge_t*> insertRib(edge_t& edge, node_t* mid_node);
+
+protected:
+    Line getSource(const edge_t& edge) const;
+};
+
+}
+#endif 

src/libslic3r/Arachne/SkeletalTrapezoidationJoint.hpp

@@ -0,0 +1,60 @@
+//Copyright (c) 2020 Ultimaker B.V.
+//CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#ifndef SKELETAL_TRAPEZOIDATION_JOINT_H
+#define SKELETAL_TRAPEZOIDATION_JOINT_H
+
+#include <memory> // smart pointers
+
+#include "libslic3r/Arachne/BeadingStrategy/BeadingStrategy.hpp"
+
+namespace Slic3r::Arachne
+{
+
+class SkeletalTrapezoidationJoint
+{
+    using Beading = BeadingStrategy::Beading;
+public:
+    struct BeadingPropagation
+    {
+        Beading beading;
+        coord_t dist_to_bottom_source;
+        coord_t dist_from_top_source;
+        bool is_upward_propagated_only;
+        BeadingPropagation(const Beading& beading)
+            : beading(beading)
+            , dist_to_bottom_source(0)
+            , dist_from_top_source(0)
+            , is_upward_propagated_only(false)
+        {}
+    };
+
+    coord_t distance_to_boundary;
+    coord_t bead_count;
+    float transition_ratio; //! The distance near the skeleton to leave free because this joint is in the middle of a transition, as a fraction of the inner bead width of the bead at the higher transition.
+    SkeletalTrapezoidationJoint()
+    : distance_to_boundary(-1)
+    , bead_count(-1)
+    , transition_ratio(0)
+    {}
+
+    bool hasBeading() const
+    {
+        return beading.use_count() > 0;
+    }
+    void setBeading(std::shared_ptr<BeadingPropagation> storage)
+    {
+        beading = storage;
+    }
+    std::shared_ptr<BeadingPropagation> getBeading()
+    {
+        return beading.lock();
+    }
+
+private:
+
+    std::weak_ptr<BeadingPropagation> beading;
+};
+
+} // namespace Slic3r::Arachne
+#endif // SKELETAL_TRAPEZOIDATION_JOINT_H

src/libslic3r/Arachne/WallToolPaths.cpp

@@ -0,0 +1,864 @@
+// Copyright (c) 2020 Ultimaker B.V.
+// CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#include <algorithm> //For std::partition_copy and std::min_element.
+#include <unordered_set>
+
+#include "WallToolPaths.hpp"
+
+#include "SkeletalTrapezoidation.hpp"
+#include "../ClipperUtils.hpp"
+#include "Arachne/utils/linearAlg2D.hpp"
+#include "EdgeGrid.hpp"
+#include "utils/SparseLineGrid.hpp"
+#include "Geometry.hpp"
+
+namespace Slic3r::Arachne
+{
+
+WallToolPaths::WallToolPaths(const Polygons& outline, const coord_t nominal_bead_width, const size_t inset_count, const coord_t wall_0_inset,
+                             const PrintConfig &print_config)
+    : outline(outline)
+    , bead_width_0(nominal_bead_width)
+    , bead_width_x(nominal_bead_width)
+    , inset_count(inset_count)
+    , wall_0_inset(wall_0_inset)
+    , strategy_type(print_config.beading_strategy_type.value)
+    , print_thin_walls(Slic3r::Arachne::fill_outline_gaps)
+    , min_feature_size(scaled<coord_t>(print_config.min_feature_size.value))
+    , min_bead_width(scaled<coord_t>(print_config.min_bead_width.value))
+    , small_area_length(static_cast<double>(nominal_bead_width) / 2.)
+    , toolpaths_generated(false)
+    , print_config(print_config)
+{
+}
+
+WallToolPaths::WallToolPaths(const Polygons& outline, const coord_t bead_width_0, const coord_t bead_width_x,
+                             const size_t inset_count, const coord_t wall_0_inset, const PrintConfig &print_config)
+    : outline(outline)
+    , bead_width_0(bead_width_0)
+    , bead_width_x(bead_width_x)
+    , inset_count(inset_count)
+    , wall_0_inset(wall_0_inset)
+    , strategy_type(print_config.beading_strategy_type.value)
+    , print_thin_walls(Slic3r::Arachne::fill_outline_gaps)
+    , min_feature_size(scaled<coord_t>(print_config.min_feature_size.value))
+    , min_bead_width(scaled<coord_t>(print_config.min_bead_width.value))
+    , small_area_length(static_cast<double>(bead_width_0) / 2.)
+    , toolpaths_generated(false)
+    , print_config(print_config)
+{
+}
+
+void simplify(Polygon &thiss, const int64_t smallest_line_segment_squared, const int64_t allowed_error_distance_squared)
+{
+    if (thiss.size() < 3)
+    {
+        thiss.points.clear();
+        return;
+    }
+    if (thiss.size() == 3)
+    {
+        return;
+    }
+
+    Polygon new_path;
+    Point previous = thiss.points.back();
+    Point previous_previous = thiss.points.at(thiss.points.size() - 2);
+    Point current = thiss.points.at(0);
+
+    /* When removing a vertex, we check the height of the triangle of the area
+     being removed from the original polygon by the simplification. However,
+     when consecutively removing multiple vertices the height of the previously
+     removed vertices w.r.t. the shortcut path changes.
+     In order to not recompute the new height value of previously removed
+     vertices we compute the height of a representative triangle, which covers
+     the same amount of area as the area being cut off. We use the Shoelace
+     formula to accumulate the area under the removed segments. This works by
+     computing the area in a 'fan' where each of the blades of the fan go from
+     the origin to one of the segments. While removing vertices the area in
+     this fan accumulates. By subtracting the area of the blade connected to
+     the short-cutting segment we obtain the total area of the cutoff region.
+     From this area we compute the height of the representative triangle using
+     the standard formula for a triangle area: A = .5*b*h
+     */
+    int64_t accumulated_area_removed = int64_t(previous.x()) * int64_t(current.y()) - int64_t(previous.y()) * int64_t(current.x()); // Twice the Shoelace formula for area of polygon per line segment.
+
+    for (size_t point_idx = 0; point_idx < thiss.points.size(); point_idx++)
+    {
+        current = thiss.points.at(point_idx % thiss.points.size());
+
+        //Check if the accumulated area doesn't exceed the maximum.
+        Point next;
+        if (point_idx + 1 < thiss.points.size())
+        {
+            next = thiss.points.at(point_idx + 1);
+        }
+        else if (point_idx + 1 == thiss.points.size() && new_path.size() > 1)
+        { // don't spill over if the [next] vertex will then be equal to [previous]
+            next = new_path[0]; //Spill over to new polygon for checking removed area.
+        }
+        else
+        {
+            next = thiss.points.at((point_idx + 1) % thiss.points.size());
+        }
+        const int64_t removed_area_next = int64_t(current.x()) * int64_t(next.y()) - int64_t(current.y()) * int64_t(next.x()); // Twice the Shoelace formula for area of polygon per line segment.
+        const int64_t negative_area_closing = int64_t(next.x()) * int64_t(previous.y()) - int64_t(next.y()) * int64_t(previous.x()); // area between the origin and the short-cutting segment
+        accumulated_area_removed += removed_area_next;
+
+        const int64_t length2 = (current - previous).cast<int64_t>().squaredNorm();
+        if (length2 < scaled<int64_t>(25.))
+        {
+            // We're allowed to always delete segments of less than 5 micron.
+            continue;
+        }
+
+        const int64_t area_removed_so_far = accumulated_area_removed + negative_area_closing; // close the shortcut area polygon
+        const int64_t base_length_2 = (next - previous).cast<int64_t>().squaredNorm();
+
+        if (base_length_2 == 0) //Two line segments form a line back and forth with no area.
+        {
+            continue; //Remove the vertex.
+        }
+        //We want to check if the height of the triangle formed by previous, current and next vertices is less than allowed_error_distance_squared.
+        //1/2 L = A           [actual area is half of the computed shoelace value] // Shoelace formula is .5*(...) , but we simplify the computation and take out the .5
+        //A = 1/2 * b * h     [triangle area formula]
+        //L = b * h           [apply above two and take out the 1/2]
+        //h = L / b           [divide by b]
+        //h^2 = (L / b)^2     [square it]
+        //h^2 = L^2 / b^2     [factor the divisor]
+        const int64_t height_2 = double(area_removed_so_far) * double(area_removed_so_far) / double(base_length_2);
+        if ((height_2 <= Slic3r::sqr(scaled<coord_t>(0.005)) //Almost exactly colinear (barring rounding errors).
+             && Line::distance_to_infinite(current, previous, next) <= scaled<double>(0.005))) // make sure that height_2 is not small because of cancellation of positive and negative areas
+        {
+            continue;
+        }
+
+        if (length2 < smallest_line_segment_squared
+            && height_2 <= allowed_error_distance_squared) // removing the vertex doesn't introduce too much error.)
+        {
+            const int64_t next_length2 = (current - next).cast<int64_t>().squaredNorm();
+            if (next_length2 > smallest_line_segment_squared)
+            {
+                // Special case; The next line is long. If we were to remove this, it could happen that we get quite noticeable artifacts.
+                // We should instead move this point to a location where both edges are kept and then remove the previous point that we wanted to keep.
+                // By taking the intersection of these two lines, we get a point that preserves the direction (so it makes the corner a bit more pointy).
+                // We just need to be sure that the intersection point does not introduce an artifact itself.
+                Point intersection_point;
+                bool has_intersection = Line(previous_previous, previous).intersection_infinite(Line(current, next), &intersection_point);
+                if (!has_intersection
+                    || Line::distance_to_infinite_squared(intersection_point, previous, current) > double(allowed_error_distance_squared)
+                    || (intersection_point - previous).cast<int64_t>().squaredNorm() > smallest_line_segment_squared  // The intersection point is way too far from the 'previous'
+                    || (intersection_point - next).cast<int64_t>().squaredNorm() > smallest_line_segment_squared)     // and 'next' points, so it shouldn't replace 'current'
+                {
+                    // We can't find a better spot for it, but the size of the line is more than 5 micron.
+                    // So the only thing we can do here is leave it in...
+                }
+                else
+                {
+                    // New point seems like a valid one.
+                    current = intersection_point;
+                    // If there was a previous point added, remove it.
+                    if(!new_path.empty())
+                    {
+                        new_path.points.pop_back();
+                        previous = previous_previous;
+                    }
+                }
+            }
+            else
+            {
+                continue; //Remove the vertex.
+            }
+        }
+        //Don't remove the vertex.
+        accumulated_area_removed = removed_area_next; // so that in the next iteration it's the area between the origin, [previous] and [current]
+        previous_previous = previous;
+        previous = current; //Note that "previous" is only updated if we don't remove the vertex.
+        new_path.points.push_back(current);
+    }
+
+    thiss = new_path;
+}
+
+/*!
+     * Removes vertices of the polygons to make sure that they are not too high
+     * resolution.
+     *
+     * This removes points which are connected to line segments that are shorter
+     * than the `smallest_line_segment`, unless that would introduce a deviation
+     * in the contour of more than `allowed_error_distance`.
+     *
+     * Criteria:
+     * 1. Never remove a vertex if either of the connceted segments is larger than \p smallest_line_segment
+     * 2. Never remove a vertex if the distance between that vertex and the final resulting polygon would be higher than \p allowed_error_distance
+     * 3. The direction of segments longer than \p smallest_line_segment always
+     * remains unaltered (but their end points may change if it is connected to
+     * a small segment)
+     *
+     * Simplify uses a heuristic and doesn't neccesarily remove all removable
+     * vertices under the above criteria, but simplify may never violate these
+     * criteria. Unless the segments or the distance is smaller than the
+     * rounding error of 5 micron.
+     *
+     * Vertices which introduce an error of less than 5 microns are removed
+     * anyway, even if the segments are longer than the smallest line segment.
+     * This makes sure that (practically) colinear line segments are joined into
+     * a single line segment.
+     * \param smallest_line_segment Maximal length of removed line segments.
+     * \param allowed_error_distance If removing a vertex introduces a deviation
+     * from the original path that is more than this distance, the vertex may
+     * not be removed.
+ */
+void simplify(Polygons &thiss, const int64_t smallest_line_segment = scaled<coord_t>(0.01), const int64_t allowed_error_distance = scaled<coord_t>(0.005))
+{
+    const int64_t allowed_error_distance_squared = int64_t(allowed_error_distance) * int64_t(allowed_error_distance);
+    const int64_t smallest_line_segment_squared = int64_t(smallest_line_segment) * int64_t(smallest_line_segment);
+    for (size_t p = 0; p < thiss.size(); p++)
+    {
+        simplify(thiss[p], smallest_line_segment_squared, allowed_error_distance_squared);
+        if (thiss[p].size() < 3)
+        {
+            thiss.erase(thiss.begin() + p);
+            p--;
+        }
+    }
+}
+
+/*!
+ * Locator to extract a line segment out of a \ref PolygonsPointIndex
+ */
+struct PolygonsPointIndexSegmentLocator
+{
+    std::pair<Point, Point> operator()(const PolygonsPointIndex &val) const
+    {
+        const Polygon &poly           = (*val.polygons)[val.poly_idx];
+        const Point    start          = poly[val.point_idx];
+        unsigned int   next_point_idx = (val.point_idx + 1) % poly.size();
+        const Point    end            = poly[next_point_idx];
+        return std::pair<Point, Point>(start, end);
+    }
+};
+
+typedef SparseLineGrid<PolygonsPointIndex, PolygonsPointIndexSegmentLocator> LocToLineGrid;
+std::unique_ptr<LocToLineGrid>                                               createLocToLineGrid(const Polygons &polygons, int square_size)
+{
+    unsigned int n_points = 0;
+    for (const auto &poly : polygons)
+        n_points += poly.size();
+
+    auto ret = std::make_unique<LocToLineGrid>(square_size, n_points);
+
+    for (unsigned int poly_idx = 0; poly_idx < polygons.size(); poly_idx++)
+        for (unsigned int point_idx = 0; point_idx < polygons[poly_idx].size(); point_idx++)
+            ret->insert(PolygonsPointIndex(&polygons, poly_idx, point_idx));
+    return ret;
+}
+
+/* Note: Also tries to solve for near-self intersections, when epsilon >= 1
+ */
+void fixSelfIntersections(const coord_t epsilon, Polygons &thiss)
+{
+    if (epsilon < 1) {
+        ClipperLib::SimplifyPolygons(ClipperUtils::PolygonsProvider(thiss));
+        return;
+    }
+
+    const int64_t half_epsilon = (epsilon + 1) / 2;
+
+    // Points too close to line segments should be moved a little away from those line segments, but less than epsilon,
+    //   so at least half-epsilon distance between points can still be guaranteed.
+    constexpr coord_t grid_size  = scaled<coord_t>(2.);
+    auto              query_grid = createLocToLineGrid(thiss, grid_size);
+
+    const auto    move_dist         = std::max<int64_t>(2L, half_epsilon - 2);
+    const int64_t half_epsilon_sqrd = half_epsilon * half_epsilon;
+
+    const size_t n = thiss.size();
+    for (size_t poly_idx = 0; poly_idx < n; poly_idx++) {
+        const size_t pathlen = thiss[poly_idx].size();
+        for (size_t point_idx = 0; point_idx < pathlen; ++point_idx) {
+            Point &pt = thiss[poly_idx][point_idx];
+            for (const auto &line : query_grid->getNearby(pt, epsilon)) {
+                const size_t line_next_idx = (line.point_idx + 1) % thiss[line.poly_idx].size();
+                if (poly_idx == line.poly_idx && (point_idx == line.point_idx || point_idx == line_next_idx))
+                    continue;
+
+                const Line segment(thiss[line.poly_idx][line.point_idx], thiss[line.poly_idx][line_next_idx]);
+                Point      segment_closest_point;
+                segment.distance_to_squared(pt, &segment_closest_point);
+
+                if (half_epsilon_sqrd >= (pt - segment_closest_point).cast<int64_t>().squaredNorm()) {
+                    const Point  &other = thiss[poly_idx][(point_idx + 1) % pathlen];
+                    const Vec2i64 vec   = (LinearAlg2D::pointIsLeftOfLine(other, segment.a, segment.b) > 0 ? segment.b - segment.a : segment.a - segment.b).cast<int64_t>();
+                    assert(Slic3r::sqr(double(vec.x())) < double(std::numeric_limits<int64_t>::max()));
+                    assert(Slic3r::sqr(double(vec.y())) < double(std::numeric_limits<int64_t>::max()));
+                    const int64_t len   = vec.norm();
+                    pt.x() += (-vec.y() * move_dist) / len;
+                    pt.y() += (vec.x() * move_dist) / len;
+                }
+            }
+        }
+    }
+
+    ClipperLib::SimplifyPolygons(ClipperUtils::PolygonsProvider(thiss));
+}
+
+/*!
+     * Removes overlapping consecutive line segments which don't delimit a positive area.
+ */
+void removeDegenerateVerts(Polygons &thiss)
+{
+    for (unsigned int poly_idx = 0; poly_idx < thiss.size(); poly_idx++) {
+        Polygon &poly = thiss[poly_idx];
+        Polygon  result;
+
+        auto isDegenerate = [](const Point &last, const Point &now, const Point &next) {
+            Vec2i64 last_line = (now - last).cast<int64_t>();
+            Vec2i64 next_line = (next - now).cast<int64_t>();
+            return last_line.dot(next_line) == -1 * last_line.norm() * next_line.norm();
+        };
+        bool isChanged = false;
+        for (unsigned int idx = 0; idx < poly.size(); idx++) {
+            const Point &last = (result.size() == 0) ? poly.back() : result.back();
+            if (idx + 1 == poly.size() && result.size() == 0) { break; }
+            Point &next = (idx + 1 == poly.size()) ? result[0] : poly[idx + 1];
+            if (isDegenerate(last, poly[idx], next)) { // lines are in the opposite direction
+                // don't add vert to the result
+                isChanged = true;
+                while (result.size() > 1 && isDegenerate(result[result.size() - 2], result.back(), next)) { result.points.pop_back(); }
+            } else {
+                result.points.emplace_back(poly[idx]);
+            }
+        }
+
+        if (isChanged) {
+            if (result.size() > 2) {
+                poly = result;
+            } else {
+                thiss.erase(thiss.begin() + poly_idx);
+                poly_idx--; // effectively the next iteration has the same poly_idx (referring to a new poly which is not yet processed)
+            }
+        }
+    }
+}
+
+void removeSmallAreas(Polygons &thiss, const double min_area_size, const bool remove_holes)
+{
+    auto to_path = [](const Polygon &poly) -> ClipperLib::Path {
+        ClipperLib::Path out;
+        for (const Point &pt : poly.points)
+            out.emplace_back(ClipperLib::cInt(pt.x()), ClipperLib::cInt(pt.y()));
+        return out;
+    };
+
+    auto new_end = thiss.end();
+    if(remove_holes)
+    {
+        for(auto it = thiss.begin(); it < new_end; it++)
+        {
+            // All polygons smaller than target are removed by replacing them with a polygon from the back of the vector
+            if(fabs(ClipperLib::Area(to_path(*it))) < min_area_size)
+            {
+                new_end--;
+                *it = std::move(*new_end);
+                it--; // wind back the iterator such that the polygon just swaped in is checked next
+            }
+        }
+    }
+    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;
+        for(auto it = thiss.begin(); it < new_end; it++) {
+            double area = ClipperLib::Area(to_path(*it));
+            if (fabs(area) < min_area_size)
+            {
+                if(area >= 0)
+                {
+                    new_end--;
+                    if(it < new_end) {
+                        std::swap(*new_end, *it);
+                        it--;
+                    }
+                    else
+                    { // Don't self-swap the last Path
+                        break;
+                    }
+                }
+                else
+                {
+                    small_holes.push_back(*it);
+                }
+            }
+        }
+
+        // Removes small holes that have their first point inside one of the removed outlines
+        // Iterating in reverse ensures that unprocessed small holes won't be moved
+        const auto removed_outlines_start = new_end;
+        for(auto hole_it = small_holes.rbegin(); hole_it < small_holes.rend(); hole_it++)
+        {
+            for(auto outline_it = removed_outlines_start; outline_it < thiss.end() ; outline_it++)
+            {
+                if(Polygon(*outline_it).contains(*hole_it->begin())) {
+                    new_end--;
+                    *hole_it = std::move(*new_end);
+                    break;
+                }
+            }
+        }
+    }
+    thiss.resize(new_end-thiss.begin());
+}
+
+void removeColinearEdges(Polygon &poly, const double max_deviation_angle)
+{
+    // TODO: Can be made more efficient (for example, use pointer-types for process-/skip-indices, so we can swap them without copy).
+    size_t num_removed_in_iteration = 0;
+    do {
+        num_removed_in_iteration = 0;
+        std::vector<bool> process_indices(poly.points.size(), true);
+
+        bool go = true;
+        while (go) {
+            go = false;
+
+            const auto  &rpath   = poly;
+            const size_t pathlen = rpath.size();
+            if (pathlen <= 3)
+                return;
+
+            std::vector<bool> skip_indices(poly.points.size(), false);
+
+            Polygon new_path;
+            for (size_t point_idx = 0; point_idx < pathlen; ++point_idx) {
+                // Don't iterate directly over process-indices, but do it this way, because there are points _in_ process-indices that should nonetheless
+                // be skipped:
+                if (!process_indices[point_idx]) {
+                    new_path.points.push_back(rpath[point_idx]);
+                    continue;
+                }
+
+                // Should skip the last point for this iteration if the old first was removed (which can be seen from the fact that the new first was skipped):
+                if (point_idx == (pathlen - 1) && skip_indices[0]) {
+                    skip_indices[new_path.size()] = true;
+                    go                            = true;
+                    new_path.points.push_back(rpath[point_idx]);
+                    break;
+                }
+
+                const Point &prev = rpath[(point_idx - 1 + pathlen) % pathlen];
+                const Point &pt   = rpath[point_idx];
+                const Point &next = rpath[(point_idx + 1) % pathlen];
+
+                float angle = LinearAlg2D::getAngleLeft(prev, pt, next); // [0 : 2 * pi]
+                if (angle >= float(M_PI)) { angle -= float(M_PI); }                    // map [pi : 2 * pi] to [0 : pi]
+
+                // Check if the angle is within limits for the point to 'make sense', given the maximum deviation.
+                // If the angle indicates near-parallel segments ignore the point 'pt'
+                if (angle > max_deviation_angle && angle < M_PI - max_deviation_angle) {
+                    new_path.points.push_back(pt);
+                } else if (point_idx != (pathlen - 1)) {
+                    // Skip the next point, since the current one was removed:
+                    skip_indices[new_path.size()] = true;
+                    go                            = true;
+                    new_path.points.push_back(next);
+                    ++point_idx;
+                }
+            }
+            poly = new_path;
+            num_removed_in_iteration += pathlen - poly.points.size();
+
+            process_indices.clear();
+            process_indices.insert(process_indices.end(), skip_indices.begin(), skip_indices.end());
+        }
+    } while (num_removed_in_iteration > 0);
+}
+
+void removeColinearEdges(Polygons &thiss, const double max_deviation_angle = 0.0005)
+{
+    for (int p = 0; p < int(thiss.size()); p++) {
+        removeColinearEdges(thiss[p], max_deviation_angle);
+        if (thiss[p].size() < 3) {
+            thiss.erase(thiss.begin() + p);
+            p--;
+        }
+    }
+}
+
+const VariableWidthPaths& WallToolPaths::generate()
+{
+    const coord_t smallest_segment = Slic3r::Arachne::meshfix_maximum_resolution;
+    const coord_t allowed_distance = Slic3r::Arachne::meshfix_maximum_deviation;
+    const coord_t epsilon_offset = (allowed_distance / 2) - 1;
+    const double  transitioning_angle = Geometry::deg2rad(this->print_config.wall_transition_angle.value);
+    constexpr coord_t discretization_step_size = scaled<coord_t>(0.8);
+
+    // Simplify outline for boost::voronoi consumption. Absolutely no self intersections or near-self intersections allowed:
+    // TODO: Open question: Does this indeed fix all (or all-but-one-in-a-million) cases for manifold but otherwise possibly complex polygons?
+    Polygons prepared_outline = offset(offset(offset(outline, -epsilon_offset), epsilon_offset * 2), -epsilon_offset);
+    simplify(prepared_outline, smallest_segment, allowed_distance);
+    fixSelfIntersections(epsilon_offset, prepared_outline);
+    removeDegenerateVerts(prepared_outline);
+    removeColinearEdges(prepared_outline, 0.005);
+    // Removing collinear edges may introduce self intersections, so we need to fix them again
+    fixSelfIntersections(epsilon_offset, prepared_outline);
+    removeDegenerateVerts(prepared_outline);
+    removeSmallAreas(prepared_outline, small_area_length * small_area_length, false);
+
+    if (area(prepared_outline) > 0)
+    {
+        const coord_t wall_transition_length = scaled<coord_t>(this->print_config.wall_transition_length.value);
+        const double wall_split_middle_threshold = this->print_config.wall_split_middle_threshold.value / 100.;  // For an uneven nr. of lines: When to split the middle wall into two.
+        const double wall_add_middle_threshold = this->print_config.wall_add_middle_threshold.value / 100.;      // For an even nr. of lines: When to add a new middle in between the innermost two walls.
+        const int wall_distribution_count = this->print_config.wall_distribution_count.value;
+        const size_t max_bead_count = (inset_count < std::numeric_limits<coord_t>::max() / 2) ? 2 * inset_count : std::numeric_limits<coord_t>::max();
+        const auto beading_strat = BeadingStrategyFactory::makeStrategy
+            (
+                strategy_type,
+                bead_width_0,
+                bead_width_x,
+                wall_transition_length,
+                transitioning_angle,
+                print_thin_walls,
+                min_bead_width,
+                min_feature_size,
+                wall_split_middle_threshold,
+                wall_add_middle_threshold,
+                max_bead_count,
+                wall_0_inset,
+                wall_distribution_count
+            );
+        const coord_t transition_filter_dist = scaled<coord_t>(this->print_config.wall_transition_filter_distance.value);
+        SkeletalTrapezoidation wall_maker
+        (
+            prepared_outline,
+            *beading_strat,
+            beading_strat->getTransitioningAngle(),
+            discretization_step_size,
+            transition_filter_dist,
+            wall_transition_length
+        );
+        wall_maker.generateToolpaths(toolpaths);
+        computeInnerContour();
+    }
+    simplifyToolPaths(toolpaths);
+
+    removeEmptyToolPaths(toolpaths);
+    assert(std::is_sorted(toolpaths.cbegin(), toolpaths.cend(),
+                          [](const VariableWidthLines& l, const VariableWidthLines& r)
+                          {
+                              return l.front().inset_idx < r.front().inset_idx;
+                          }) && "WallToolPaths should be sorted from the outer 0th to inner_walls");
+    toolpaths_generated = true;
+    return toolpaths;
+}
+
+void WallToolPaths::simplifyToolPaths(VariableWidthPaths& toolpaths/*, const Settings& settings*/)
+{
+    for (size_t toolpaths_idx = 0; toolpaths_idx < toolpaths.size(); ++toolpaths_idx)
+    {
+        const int64_t maximum_resolution = Slic3r::Arachne::meshfix_maximum_resolution;
+        const int64_t maximum_deviation = Slic3r::Arachne::meshfix_maximum_deviation;
+        const int64_t maximum_extrusion_area_deviation = Slic3r::Arachne::meshfix_maximum_extrusion_area_deviation; // unit: μm²
+        for (auto& line : toolpaths[toolpaths_idx])
+        {
+            line.simplify(maximum_resolution * maximum_resolution, maximum_deviation * maximum_deviation, maximum_extrusion_area_deviation);
+        }
+    }
+}
+
+const VariableWidthPaths& WallToolPaths::getToolPaths()
+{
+    if (!toolpaths_generated)
+    {
+        return generate();
+    }
+    return toolpaths;
+}
+
+void WallToolPaths::computeInnerContour()
+{
+    //We'll remove all 0-width paths from the original toolpaths and store them separately as polygons.
+    VariableWidthPaths actual_toolpaths;
+    actual_toolpaths.reserve(toolpaths.size()); //A bit too much, but the correct order of magnitude.
+    VariableWidthPaths contour_paths;
+    contour_paths.reserve(toolpaths.size() / inset_count);
+    std::partition_copy(toolpaths.begin(), toolpaths.end(), std::back_inserter(actual_toolpaths), std::back_inserter(contour_paths),
+        [](const VariableWidthLines& path)
+        {
+            for(const ExtrusionLine& line : path)
+            {
+                for(const ExtrusionJunction& junction : line.junctions)
+                {
+                    return junction.w != 0; //On the first actual junction, decide: If it's got 0 width, this is a contour. Otherwise it is an actual toolpath.
+                }
+            }
+            return true; //No junctions with any vertices? Classify it as a toolpath then.
+        });
+    if (! actual_toolpaths.empty())
+    {
+        toolpaths = std::move(actual_toolpaths); //Filtered out the 0-width paths.
+    }
+    else
+    {
+        toolpaths.clear();
+    }
+
+    //Now convert the contour_paths to Polygons to denote the inner contour of the walled areas.
+    inner_contour.clear();
+
+    //We're going to have to stitch these paths since not all walls may be closed contours.
+    //Since these walls have 0 width they should theoretically be closed. But there may be rounding errors.
+    const coord_t minimum_line_width = bead_width_0 / 2;
+    stitchContours(contour_paths, minimum_line_width, inner_contour);
+
+    //The output walls from the skeletal trapezoidation have no known winding order, especially if they are joined together from polylines.
+    //They can be in any direction, clockwise or counter-clockwise, regardless of whether the shapes are positive or negative.
+    //To get a correct shape, we need to make the outside contour positive and any holes inside negative.
+    //This can be done by applying the even-odd rule to the shape. This rule is not sensitive to the winding order of the polygon.
+    //The even-odd rule would be incorrect if the polygon self-intersects, but that should never be generated by the skeletal trapezoidation.
+    inner_contour = union_(inner_contour);
+}
+
+const Polygons& WallToolPaths::getInnerContour()
+{
+    if (!toolpaths_generated && inset_count > 0)
+    {
+        generate();
+    }
+    else if(inset_count == 0)
+    {
+        return outline;
+    }
+    return inner_contour;
+}
+
+bool WallToolPaths::removeEmptyToolPaths(VariableWidthPaths& toolpaths)
+{
+    toolpaths.erase(std::remove_if(toolpaths.begin(), toolpaths.end(), [](const VariableWidthLines& lines)
+                                   {
+                                       return lines.empty();
+                                   }), toolpaths.end());
+    return toolpaths.empty();
+}
+
+void WallToolPaths::stitchContours(const VariableWidthPaths& input, const coord_t stitch_distance, Polygons& output)
+{
+    // Create a bucket grid to find endpoints that are close together.
+    struct ExtrusionLineStartLocator
+    {
+        const Point *operator()(const ExtrusionLine *line) { return &line->junctions.front().p; }
+    };
+    struct ExtrusionLineEndLocator
+    {
+        const Point *operator()(const ExtrusionLine *line) { return &line->junctions.back().p; }
+    };
+
+    // Only find endpoints closer than minimum_line_width, so we can't ever accidentally make crossing contours.
+    ClosestPointInRadiusLookup<const ExtrusionLine*, ExtrusionLineStartLocator> line_starts(coord_t(stitch_distance * std::sqrt(2.)));
+    ClosestPointInRadiusLookup<const ExtrusionLine*, ExtrusionLineEndLocator> line_ends(coord_t(stitch_distance * std::sqrt(2.)));
+
+    auto get_search_bbox = [](const Point &pt, const coord_t radius) -> BoundingBox {
+        const Point min_grid((pt - Point(radius, radius)) / radius);
+        const Point max_grid((pt + Point(radius, radius)) / radius);
+        return {min_grid * radius, (max_grid + Point(1, 1)) * radius - Point(1, 1)};
+    };
+
+    for (const VariableWidthLines &path : input) {
+        for (const ExtrusionLine &line : path) {
+            line_starts.insert(&line);
+            line_ends.insert(&line);
+        }
+    }
+    //Then go through all lines and construct chains of polylines if the endpoints are nearby.
+    std::unordered_set<const ExtrusionLine*> processed_lines; //Track which lines were already processed to not process them twice.
+    for(const VariableWidthLines& path : input)
+    {
+        for(const ExtrusionLine& line : path)
+        {
+            if(processed_lines.find(&line) != processed_lines.end()) //We already added this line before. It got added as a nearby line.
+            {
+                continue;
+            }
+            //We'll create a chain of polylines that get joined together. We can add polylines on both ends!
+            std::deque<const ExtrusionLine*> chain;
+            std::deque<bool> is_reversed; //Lines could need to be inserted in reverse. Must coincide with the `chain` deque.
+            const ExtrusionLine* nearest = &line; //At every iteration, add the polyline that joins together most closely.
+            bool nearest_reverse = false; //Whether the next line to insert must be inserted in reverse.
+            bool nearest_before = false; //Whether the next line to insert must be inserted in the front of the chain.
+            while(nearest)
+            {
+                if(processed_lines.find(nearest) != processed_lines.end())
+                {
+                    break; //Looping. This contour is already processed.
+                }
+                processed_lines.insert(nearest);
+                if(nearest_before)
+                {
+                    chain.push_front(nearest);
+                    is_reversed.push_front(nearest_reverse);
+                }
+                else
+                {
+                    chain.push_back(nearest);
+                    is_reversed.push_back(nearest_reverse);
+                }
+
+                //Find any nearby lines to attach. Look on both ends of our current chain and find both ends of polylines.
+                const Point chain_start = is_reversed.front() ? chain.front()->junctions.back().p : chain.front()->junctions.front().p;
+                const Point chain_end = is_reversed.back() ? chain.back()->junctions.front().p : chain.back()->junctions.back().p;
+
+                std::vector<std::pair<const ExtrusionLine *const *, double>> starts_near_start = line_starts.find_all(chain_start);
+                std::vector<std::pair<const ExtrusionLine *const *, double>> ends_near_start   = line_ends.find_all(chain_start);
+                std::vector<std::pair<const ExtrusionLine *const *, double>> starts_near_end   = line_starts.find_all(chain_end);
+                std::vector<std::pair<const ExtrusionLine *const *, double>> ends_near_end     = line_ends.find_all(chain_end);
+
+                nearest = nullptr;
+                int64_t nearest_dist2 = std::numeric_limits<int64_t>::max();
+                for (const auto &candidate_ptr : starts_near_start) {
+                    const ExtrusionLine* candidate = *candidate_ptr.first;
+                    if(processed_lines.find(candidate) != processed_lines.end())
+                        continue; //Already processed this line before. It's linked to something else.
+
+                    if (const int64_t dist2 = (candidate->junctions.front().p - chain_start).cast<int64_t>().squaredNorm(); dist2 < nearest_dist2) {
+                        nearest         = candidate;
+                        nearest_dist2   = dist2;
+                        nearest_reverse = true;
+                        nearest_before  = true;
+                    }
+                }
+                for (const auto &candidate_ptr : ends_near_start) {
+                    const ExtrusionLine* candidate = *candidate_ptr.first;
+                    if(processed_lines.find(candidate) != processed_lines.end())
+                        continue;
+
+                    if (const int64_t dist2 = (candidate->junctions.back().p - chain_start).cast<int64_t>().squaredNorm(); dist2 < nearest_dist2) {
+                        nearest         = candidate;
+                        nearest_dist2   = dist2;
+                        nearest_reverse = false;
+                        nearest_before  = true;
+                    }
+                }
+                for (const auto &candidate_ptr : starts_near_end) {
+                    const ExtrusionLine* candidate = *candidate_ptr.first;
+                    if(processed_lines.find(candidate) != processed_lines.end())
+                        continue; //Already processed this line before. It's linked to something else.
+
+                    if (const int64_t dist2 = (candidate->junctions.front().p - chain_start).cast<int64_t>().squaredNorm(); dist2 < nearest_dist2) {
+                        nearest         = candidate;
+                        nearest_dist2   = dist2;
+                        nearest_reverse = false;
+                        nearest_before  = false;
+                    }
+                }
+                for (const auto &candidate_ptr : ends_near_end) {
+                    const ExtrusionLine* candidate = *candidate_ptr.first;
+                    if (processed_lines.find(candidate) != processed_lines.end())
+                        continue;
+
+                    if (const int64_t dist2 = (candidate->junctions.back().p - chain_start).cast<int64_t>().squaredNorm(); dist2 < nearest_dist2) {
+                        nearest         = candidate;
+                        nearest_dist2   = dist2;
+                        nearest_reverse = true;
+                        nearest_before  = false;
+                    }
+                }
+            }
+
+            //Now serialize the entire chain into one polygon.
+            output.emplace_back();
+            for (size_t i = 0; i < chain.size(); ++i) {
+                if(!is_reversed[i])
+                    for (const ExtrusionJunction& junction : chain[i]->junctions)
+                        output.back().points.emplace_back(junction.p);
+                else
+                    for (auto junction = chain[i]->junctions.rbegin(); junction != chain[i]->junctions.rend(); ++junction)
+                        output.back().points.emplace_back(junction->p);
+            }
+        }
+    }
+}
+
+size_t getOuterRegionId(const Arachne::VariableWidthPaths& toolpaths, size_t& out_max_region_id)
+{
+    // Polygons show up here one by one, so there are always only a) the outer lines and b) the lines that are part of the holes.
+    // Therefore, the outer-regions' lines will always have the region-id that is larger then all of the other ones.
+
+    // First, build the bounding boxes:
+    std::map<size_t, BoundingBox> region_ids_to_bboxes; // Use a sorted map, ordered by region_id, so that we can find the largest region_id quickly.
+    for (const Arachne::VariableWidthLines &path : toolpaths) {
+        for (const Arachne::ExtrusionLine &line : path) {
+            BoundingBox &aabb =
+                region_ids_to_bboxes[line.region_id]; // Empty AABBs are default initialized when region_ids are encountered for the first time.
+            for (const auto &junction : line.junctions) aabb.merge(junction.p);
+        }
+    }
+
+    // Then, the largest of these will be the one that's needed for the outer region, the others' all belong to hole regions:
+    BoundingBox outer_bbox;
+    size_t      outer_region_id = 0; // Region-ID 0 is reserved for 'None'.
+    for (const auto &region_id_bbox_pair : region_ids_to_bboxes) {
+        if (region_id_bbox_pair.second.contains(outer_bbox)) {
+            outer_bbox      = region_id_bbox_pair.second;
+            outer_region_id = region_id_bbox_pair.first;
+        }
+    }
+
+    // Maximum Region-ID (using the ordering of the map)
+    out_max_region_id = region_ids_to_bboxes.empty() ? 0 : region_ids_to_bboxes.rbegin()->first;
+    return outer_region_id;
+}
+
+Arachne::BinJunctions variableWidthPathToBinJunctions(const Arachne::VariableWidthPaths& toolpaths, const bool pack_regions_by_inset, const bool center_last, std::set<size_t>* p_bins_with_index_zero_insets)
+{
+    // Find the largest inset-index:
+    size_t max_inset_index = 0;
+    for (const Arachne::VariableWidthLines &path : toolpaths)
+        max_inset_index = std::max(path.front().inset_idx, max_inset_index);
+
+    // Find which regions are associated with the outer-outer walls (which region is the one the rest is holes inside of):
+    size_t       max_region_id   = 0;
+    const size_t outer_region_id = getOuterRegionId(toolpaths, max_region_id);
+
+    //Since we're (optionally!) splitting off in the outer and inner regions, it may need twice as many bins as inset-indices.
+    //Add two extra bins for the center-paths, if they need to be stored separately. One bin for inner and one for outer walls.
+    const size_t max_bin = (pack_regions_by_inset ? (max_region_id * 2) + 2 : (max_inset_index + 1) * 2) + center_last * 2;
+    Arachne::BinJunctions insets(max_bin + 1);
+    for (const Arachne::VariableWidthLines &path : toolpaths) {
+        if (path.empty()) // Don't bother printing these.
+            continue;
+
+        const size_t inset_index = path.front().inset_idx;
+
+        // Convert list of extrusion lines to vectors of extrusion junctions, and add those to the binned insets.
+        for (const Arachne::ExtrusionLine &line : path) {
+            // Sort into the right bin, ...
+            size_t     bin_index;
+            const bool in_hole_region = line.region_id != outer_region_id && line.region_id != 0;
+            if (center_last && line.is_odd) {
+                bin_index = inset_index > 0;
+            } else if (pack_regions_by_inset) {
+                bin_index = std::min(inset_index, static_cast<size_t>(1)) + 2 * (in_hole_region ? line.region_id : 0) + center_last * 2;
+            } else {
+                bin_index = inset_index + (in_hole_region ? (max_inset_index + 1) : 0) + center_last * 2;
+            }
+            insets[bin_index].emplace_back(line.junctions.begin(), line.junctions.end());
+
+            // Collect all bins that have zero-inset indices in them, if needed:
+            if (inset_index == 0 && p_bins_with_index_zero_insets != nullptr)
+                p_bins_with_index_zero_insets->insert(bin_index);
+        }
+    }
+    return insets;
+}
+
+BinJunctions WallToolPaths::getBinJunctions(std::set<size_t> &bins_with_index_zero_insets)
+{
+    if (!toolpaths_generated)
+        generate();
+
+    return variableWidthPathToBinJunctions(toolpaths, true, true, &bins_with_index_zero_insets);
+}
+
+} // namespace Slic3r::Arachne

src/libslic3r/Arachne/WallToolPaths.hpp

@@ -0,0 +1,130 @@
+// Copyright (c) 2020 Ultimaker B.V.
+// CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#ifndef CURAENGINE_WALLTOOLPATHS_H
+#define CURAENGINE_WALLTOOLPATHS_H
+
+#include <memory>
+
+#include "BeadingStrategy/BeadingStrategyFactory.hpp"
+#include "utils/ExtrusionLine.hpp"
+#include "../Polygon.hpp"
+#include "../PrintConfig.hpp"
+
+namespace Slic3r::Arachne
+{
+
+constexpr bool    fill_outline_gaps                        = true;
+constexpr coord_t meshfix_maximum_resolution               = scaled<coord_t>(0.5);
+constexpr coord_t meshfix_maximum_deviation                = scaled<coord_t>(0.025);
+constexpr coord_t meshfix_maximum_extrusion_area_deviation = scaled<coord_t>(2.);
+
+class WallToolPaths
+{
+public:
+    /*!
+     * A class that creates the toolpaths given an outline, nominal bead width and maximum amount of walls
+     * \param outline An outline of the area in which the ToolPaths are to be generated
+     * \param nominal_bead_width The nominal bead width used in the generation of the toolpaths
+     * \param inset_count The maximum number of parallel extrusion lines that make up the wall
+     * \param wall_0_inset How far to inset the outer wall, to make it adhere better to other walls.
+     */
+    WallToolPaths(const Polygons& outline, const coord_t nominal_bead_width, const size_t inset_count, const coord_t wall_0_inset, const PrintConfig &print_config);
+
+    /*!
+     * A class that creates the toolpaths given an outline, nominal bead width and maximum amount of walls
+     * \param outline An outline of the area in which the ToolPaths are to be generated
+     * \param bead_width_0 The bead width of the first wall used in the generation of the toolpaths
+     * \param bead_width_x The bead width of the inner walls used in the generation of the toolpaths
+     * \param inset_count The maximum number of parallel extrusion lines that make up the wall
+     * \param wall_0_inset How far to inset the outer wall, to make it adhere better to other walls.
+     */
+    WallToolPaths(const Polygons& outline, const coord_t bead_width_0, const coord_t bead_width_x, const size_t inset_count, const coord_t wall_0_inset, const PrintConfig &print_config);
+
+    /*!
+     * Generates the Toolpaths
+     * \return A reference to the newly create  ToolPaths
+     */
+    const VariableWidthPaths& generate();
+
+    /*!
+     * Gets the toolpaths, if this called before \p generate() it will first generate the Toolpaths
+     * \return a reference to the toolpaths
+     */
+    const VariableWidthPaths& getToolPaths();
+
+    BinJunctions getBinJunctions(std::set<size_t> &bins_with_index_zero_insets);
+
+    /*!
+     * Compute the inner contour of the walls. This contour indicates where the walled area ends and its infill begins.
+     * The inside can then be filled, e.g. with skin/infill for the walls of a part, or with a pattern in the case of
+     * infill with extra infill walls.
+     */
+    void computeInnerContour();
+
+    /*!
+     * Gets the inner contour of the area which is inside of the generated tool
+     * paths.
+     *
+     * If the walls haven't been generated yet, this will lazily call the
+     * \p generate() function to generate the walls with variable width.
+     * The resulting polygon will snugly match the inside of the variable-width
+     * walls where the walls get limited by the LimitedBeadingStrategy to a
+     * maximum wall count.
+     * If there are no walls, the outline will be returned.
+     * \return The inner contour of the generated walls.
+     */
+    const Polygons& getInnerContour();
+
+    /*!
+     * Removes empty paths from the toolpaths
+     * \param toolpaths the VariableWidthPaths generated with \p generate()
+     * \return true if there are still paths left. If all toolpaths were removed it returns false
+     */
+    static bool removeEmptyToolPaths(VariableWidthPaths& toolpaths);
+
+    /*!
+     * Stitches toolpaths together to form contours.
+     *
+     * All toolpaths are used. Paths that are not closed will get closed in the
+     * output by virtue of becoming polygons. As such, the input is expected to
+     * consist of almost completely closed contours, which may be split up into
+     * different polylines.
+     * This function combines those polylines into the polygons they are
+     * probably intended to depict.
+     * \param input The paths to stitch together.
+     * \param stitch_distance Any endpoints closer than this distance can be
+     * stitched together. An additional line segment will bridge the gap.
+     * \param output Where to store the output polygons.
+     */
+    static void stitchContours(const VariableWidthPaths& input, const coord_t stitch_distance, Polygons& output) ;
+
+protected:
+    /*!
+     * Simplifies the variable-width toolpaths by calling the simplify on every line in the toolpath using the provided
+     * settings.
+     * \param settings The settings as provided by the user
+     * \return
+     */
+    static void simplifyToolPaths(VariableWidthPaths& toolpaths/*, const Settings& settings*/);
+
+private:
+    const Polygons& outline; //<! A reference to the outline polygon that is the designated area
+    coord_t bead_width_0; //<! The nominal or first extrusion line width with which libArachne generates its walls
+    coord_t bead_width_x; //<! The subsequently extrusion line width with which libArachne generates its walls if WallToolPaths was called with the nominal_bead_width Constructor this is the same as bead_width_0
+    size_t inset_count; //<! The maximum number of walls to generate
+    coord_t wall_0_inset; //<! How far to inset the outer wall. Should only be applied when printing the actual walls, not extra infill/skin/support walls.
+    BeadingStrategyType strategy_type; //<! The wall generating strategy
+    bool print_thin_walls; //<! Whether to enable the widening beading meta-strategy for thin features
+    coord_t min_feature_size; //<! The minimum size of the features that can be widened by the widening beading meta-strategy. Features thinner than that will not be printed
+    coord_t min_bead_width;  //<! The minimum bead size to use when widening thin model features with the widening beading meta-strategy
+    double small_area_length; //<! The length of the small features which are to be filtered out, this is squared into a surface
+    bool toolpaths_generated; //<! Are the toolpaths generated
+    VariableWidthPaths toolpaths; //<! The generated toolpaths
+    Polygons inner_contour;  //<! The inner contour of the generated toolpaths
+    const PrintConfig &print_config;
+};
+
+} // namespace Slic3r::Arachne
+
+#endif // CURAENGINE_WALLTOOLPATHS_H

src/libslic3r/Arachne/utils/ExtrusionJunction.cpp

@@ -0,0 +1,23 @@
+//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, const size_t region_id)
+    : p(p),
+      w(w),
+      perimeter_index(perimeter_index),
+      region_id(region_id)
+{}
+
+}

src/libslic3r/Arachne/utils/ExtrusionJunction.hpp

@@ -0,0 +1,68 @@
+//Copyright (c) 2020 Ultimaker B.V.
+//CuraEngine is released under the terms of the AGPLv3 or higher.
+
+
+#ifndef UTILS_EXTRUSION_JUNCTION_H
+#define UTILS_EXTRUSION_JUNCTION_H
+
+#include "../../Point.hpp"
+
+namespace Slic3r::Arachne
+{
+
+/*!
+ * This struct represents one vertex in an extruded path.
+ *
+ * It contains information on how wide the extruded path must be at this point,
+ * and which perimeter it represents.
+ */
+struct ExtrusionJunction
+{
+	/*!
+	 * The position of the centreline of the path when it reaches this junction.
+	 * This is the position that should end up in the g-code eventually.
+	 */
+    Point p;
+
+	/*!
+	 * The width of the extruded path at this junction.
+	 */
+    coord_t w;
+
+	/*!
+	 * Which perimeter this junction is part of.
+	 *
+	 * Perimeters are counted from the outside inwards. The outer wall has index
+	 * 0.
+	 */
+    size_t perimeter_index;
+
+    /*!
+     * Which region this juntion is part of. A solid polygon without holes has only one region.
+     * A polygon with holes has 2. Disconnected parts of the polygon are also separate regions.
+     * Will be 0 if no region was given.
+     */
+    size_t region_id;
+
+    ExtrusionJunction(const Point p, const coord_t w, const coord_t perimeter_index, const size_t region_id = 0);
+
+    bool operator==(const ExtrusionJunction& other) const;
+};
+
+inline Point operator-(const ExtrusionJunction& a, const ExtrusionJunction& b)
+{
+    return a.p - b.p;
+}
+
+// Identity function, used to be able to make templated algorithms that do their operations on 'point-like' input.
+inline const Point& make_point(const ExtrusionJunction& ej)
+{
+    return ej.p;
+}
+
+using LineJunctions = std::vector<ExtrusionJunction>; //<! Vector of Lines
+using PathJunctions = std::vector<LineJunctions>; //<! Vector of paths
+using BinJunctions = std::vector<PathJunctions>; //<! Vector of insets (bins)
+
+}
+#endif // UTILS_EXTRUSION_JUNCTION_H

src/libslic3r/Arachne/utils/ExtrusionLine.cpp

@@ -0,0 +1,241 @@
+//Copyright (c) 2020 Ultimaker B.V.
+//CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#include <algorithm>
+
+#include "ExtrusionLine.hpp"
+#include "linearAlg2D.hpp"
+
+namespace Slic3r::Arachne
+{
+
+ExtrusionLine::ExtrusionLine(const size_t inset_idx, const bool is_odd, const size_t region_id)
+: inset_idx(inset_idx)
+, is_odd(is_odd)
+, region_id(region_id)
+{}
+
+coord_t ExtrusionLine::getLength() const
+{
+    if (junctions.empty())
+    {
+        return 0;
+    }
+    coord_t len = 0;
+    ExtrusionJunction prev = junctions.front();
+    for (const ExtrusionJunction& next : junctions)
+    {
+        len += (next.p - prev.p).cast<int64_t>().norm();
+        prev = next;
+    }
+    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::appendJunctionsTo(LineJunctions& result) const
+{
+    result.insert(result.end(), junctions.begin(), junctions.end());
+}
+
+void ExtrusionLine::simplify(const int64_t smallest_line_segment_squared, const int64_t allowed_error_distance_squared, const int64_t maximum_extrusion_area_deviation)
+{
+    if (junctions.size() <= 3)
+    {
+        return;
+    }
+
+    /* ExtrusionLines are treated as (open) polylines, so in case an ExtrusionLine is actually a closed polygon, its
+     * starting and ending points will be equal (or almost equal). Therefore, the simplification of the ExtrusionLine
+     * should not touch the first and last points. As a result, start simplifying from point at index 1.
+     * */
+    std::vector<ExtrusionJunction> new_junctions;
+    // Starting junction should always exist in the simplified path
+    new_junctions.emplace_back(junctions.front());
+
+    /* Initially, previous_previous is always the same as previous because, for open ExtrusionLines the last junction
+     * cannot be taken into consideration when checking the points at index 1. For closed ExtrusionLines, the first and
+     * last junctions are anyway the same.
+     * */
+    ExtrusionJunction previous_previous = junctions.front();
+    ExtrusionJunction previous = junctions.front();
+
+    /* When removing a vertex, we check the height of the triangle of the area
+     being removed from the original polygon by the simplification. However,
+     when consecutively removing multiple vertices the height of the previously
+     removed vertices w.r.t. the shortcut path changes.
+     In order to not recompute the new height value of previously removed
+     vertices we compute the height of a representative triangle, which covers
+     the same amount of area as the area being cut off. We use the Shoelace
+     formula to accumulate the area under the removed segments. This works by
+     computing the area in a 'fan' where each of the blades of the fan go from
+     the origin to one of the segments. While removing vertices the area in
+     this fan accumulates. By subtracting the area of the blade connected to
+     the short-cutting segment we obtain the total area of the cutoff region.
+     From this area we compute the height of the representative triangle using
+     the standard formula for a triangle area: A = .5*b*h
+     */
+    const ExtrusionJunction& initial = junctions.at(1);
+    int64_t accumulated_area_removed = int64_t(previous.p.x()) * int64_t(initial.p.y()) - int64_t(previous.p.y()) * int64_t(initial.p.x()); // Twice the Shoelace formula for area of polygon per line segment.
+
+    for (size_t point_idx = 1; point_idx < junctions.size() - 1; point_idx++)
+    {
+        const ExtrusionJunction& current = junctions[point_idx];
+
+        // Spill over in case of overflow, unless the [next] vertex will then be equal to [previous].
+        const bool spill_over = point_idx + 1 == junctions.size() && new_junctions.size() > 1;
+        ExtrusionJunction& next = spill_over ? new_junctions[0] : junctions[point_idx + 1];
+
+        const int64_t removed_area_next = int64_t(current.p.x()) * int64_t(next.p.y()) - int64_t(current.p.y()) * int64_t(next.p.x()); // Twice the Shoelace formula for area of polygon per line segment.
+        const int64_t negative_area_closing = int64_t(next.p.x()) * int64_t(previous.p.y()) - int64_t(next.p.y()) * int64_t(previous.p.x()); // Area between the origin and the short-cutting segment
+        accumulated_area_removed += removed_area_next;
+
+        const int64_t length2 = (current - previous).cast<int64_t>().squaredNorm();
+        if (length2 < scaled<coord_t>(0.025))
+        {
+            // We're allowed to always delete segments of less than 5 micron. The width in this case doesn't matter that much.
+            continue;
+        }
+
+        const int64_t area_removed_so_far = accumulated_area_removed + negative_area_closing; // Close the shortcut area polygon
+        const int64_t base_length_2 = (next - previous).cast<int64_t>().squaredNorm();
+
+        if (base_length_2 == 0) // Two line segments form a line back and forth with no area.
+        {
+            continue; // Remove the junction (vertex).
+        }
+        //We want to check if the height of the triangle formed by previous, current and next vertices is less than allowed_error_distance_squared.
+        //1/2 L = A           [actual area is half of the computed shoelace value] // Shoelace formula is .5*(...) , but we simplify the computation and take out the .5
+        //A = 1/2 * b * h     [triangle area formula]
+        //L = b * h           [apply above two and take out the 1/2]
+        //h = L / b           [divide by b]
+        //h^2 = (L / b)^2     [square it]
+        //h^2 = L^2 / b^2     [factor the divisor]
+        const int64_t height_2 = int64_t(double(area_removed_so_far) * double(area_removed_so_far) / double(base_length_2));
+        coord_t weighted_average_width;
+        const int64_t extrusion_area_error = calculateExtrusionAreaDeviationError(previous, current, next, weighted_average_width);
+        if ((height_2 <= 1 //Almost exactly colinear (barring rounding errors).
+             && Line::distance_to_infinite(current.p, previous.p, next.p) <= 1.) // Make sure that height_2 is not small because of cancellation of positive and negative areas
+            // We shouldn't remove middle junctions of colinear segments if the area changed for the C-P segment is exceeding the maximum allowed
+             && extrusion_area_error <= maximum_extrusion_area_deviation)
+        {
+            // Adjust the width of the entire P-N line as a weighted average of the widths of the P-C and C-N lines and
+            // then remove the current junction (vertex).
+            next.w = weighted_average_width;
+            continue;
+        }
+
+        if (length2 < smallest_line_segment_squared
+            && height_2 <= allowed_error_distance_squared) // Removing the junction (vertex) doesn't introduce too much error.
+        {
+            const int64_t next_length2 = (current - next).cast<int64_t>().squaredNorm();
+            if (next_length2 > smallest_line_segment_squared)
+            {
+                // Special case; The next line is long. If we were to remove this, it could happen that we get quite noticeable artifacts.
+                // We should instead move this point to a location where both edges are kept and then remove the previous point that we wanted to keep.
+                // By taking the intersection of these two lines, we get a point that preserves the direction (so it makes the corner a bit more pointy).
+                // We just need to be sure that the intersection point does not introduce an artifact itself.
+                Point intersection_point;
+                bool has_intersection = Line(previous_previous.p, previous.p).intersection_infinite(Line(current.p, next.p), &intersection_point);
+                if (!has_intersection
+                    || Line::distance_to_infinite_squared(intersection_point, previous.p, current.p) > double(allowed_error_distance_squared)
+                    || (intersection_point - previous.p).cast<int64_t>().squaredNorm() > smallest_line_segment_squared  // The intersection point is way too far from the 'previous'
+                    || (intersection_point - next.p).cast<int64_t>().squaredNorm() > smallest_line_segment_squared)     // and 'next' points, so it shouldn't replace 'current'
+                {
+                    // We can't find a better spot for it, but the size of the line is more than 5 micron.
+                    // So the only thing we can do here is leave it in...
+                }
+                else
+                {
+                    // New point seems like a valid one.
+                    const ExtrusionJunction new_to_add = ExtrusionJunction(intersection_point, current.w, current.perimeter_index, current.region_id);
+                    // If there was a previous point added, remove it.
+                    if(!new_junctions.empty())
+                    {
+                        new_junctions.pop_back();
+                        previous = previous_previous;
+                    }
+
+                    // The junction (vertex) is replaced by the new one.
+                    accumulated_area_removed = removed_area_next; // So that in the next iteration it's the area between the origin, [previous] and [current]
+                    previous_previous = previous;
+                    previous = new_to_add; // Note that "previous" is only updated if we don't remove the junction (vertex).
+                    new_junctions.push_back(new_to_add);
+                    continue;
+                }
+            }
+            else
+            {
+                continue; // Remove the junction (vertex).
+            }
+        }
+        // The junction (vertex) isn't removed.
+        accumulated_area_removed = removed_area_next; // So that in the next iteration it's the area between the origin, [previous] and [current]
+        previous_previous = previous;
+        previous = current; // Note that "previous" is only updated if we don't remove the junction (vertex).
+        new_junctions.push_back(current);
+    }
+
+    // Ending junction (vertex) should always exist in the simplified path
+    new_junctions.emplace_back(junctions.back());
+
+    /* In case this is a closed polygon (instead of a poly-line-segments), the invariant that the first and last points are the same should be enforced.
+     * Since one of them didn't move, and the other can't have been moved further than the constraints, if originally equal, they can simply be equated.
+     */
+    if ((junctions.front().p - junctions.back().p).cast<int64_t>().squaredNorm() == 0)
+    {
+        new_junctions.back().p = junctions.front().p;
+    }
+
+    junctions = new_junctions;
+}
+
+int64_t ExtrusionLine::calculateExtrusionAreaDeviationError(ExtrusionJunction A, ExtrusionJunction B, ExtrusionJunction C, coord_t& weighted_average_width)
+{
+    /*
+     * A             B                          C              A                                        C
+     * ---------------                                         **************
+     * |             |                                         ------------------------------------------
+     * |             |--------------------------|  B removed   |            |***************************|
+     * |             |                          |  --------->  |            |                           |
+     * |             |--------------------------|              |            |***************************|
+     * |             |                                         ------------------------------------------
+     * ---------------             ^                           **************
+     *       ^                    C.w                                             ^
+     *      B.w                                                     new_width = weighted_average_width
+     *
+     *
+     * ******** denote the total extrusion area deviation error in the consecutive segments as a result of using the
+     * weighted-average width for the entire extrusion line.
+     *
+     * */
+    const int64_t ab_length = (B - A).cast<int64_t>().norm();
+    const int64_t bc_length = (C - B).cast<int64_t>().norm();
+    const int64_t width_diff = llabs(B.w - C.w);
+    if (width_diff > 1)
+    {
+        // Adjust the width only if there is a difference, or else the rounding errors may produce the wrong
+        // weighted average value.
+        assert(((int64_t(ab_length) * int64_t(B.w) + int64_t(bc_length) * int64_t(C.w)) / (C - A).cast<int64_t>().norm()) <= std::numeric_limits<int64_t>::max());
+        weighted_average_width = (ab_length * int64_t(B.w) + bc_length * int64_t(C.w)) / (C - A).cast<int64_t>().norm();
+        assert((double(llabs(B.w - weighted_average_width)) * double(ab_length) + double(llabs(C.w - weighted_average_width)) * double(bc_length)) <= double(std::numeric_limits<int64_t>::max()));
+        return int64_t(llabs(B.w - weighted_average_width)) * ab_length + int64_t(llabs(C.w - weighted_average_width)) * bc_length;
+    }
+    else
+    {
+        // If the width difference is very small, then select the width of the segment that is longer
+        weighted_average_width = ab_length > bc_length ? B.w : C.w;
+        assert((int64_t(width_diff) * int64_t(bc_length)) <= std::numeric_limits<coord_t>::max());
+        assert((int64_t(width_diff) * int64_t(ab_length)) <= std::numeric_limits<coord_t>::max());
+        return ab_length > bc_length ? width_diff * bc_length : width_diff * ab_length;
+    }
+}
+
+}

src/libslic3r/Arachne/utils/ExtrusionLine.hpp

@@ -0,0 +1,159 @@
+//Copyright (c) 2020 Ultimaker B.V.
+//CuraEngine is released under the terms of the AGPLv3 or higher.
+
+
+#ifndef UTILS_EXTRUSION_LINE_H
+#define UTILS_EXTRUSION_LINE_H
+
+#include "ExtrusionJunction.hpp"
+#include "../../Polyline.hpp"
+#include "../../Polygon.hpp"
+
+namespace Slic3r {
+class ThickPolyline;
+}
+
+namespace Slic3r::Arachne
+{
+
+/*!
+ * Represents a polyline (not just a line) that is to be extruded with variable
+ * line width.
+ *
+ * This polyline is a sequence of \ref ExtrusionJunction, with a bit of metadata
+ * about which inset it represents.
+ */
+struct ExtrusionLine
+{
+    /*!
+     * Which inset this path represents, counted from the outside inwards.
+     *
+     * The outer wall has index 0.
+     */
+    size_t inset_idx;
+
+    /*!
+     * If a thin piece needs to be printed with an odd number of walls (e.g. 5
+     * walls) then there will be one wall in the middle that is not a loop. This
+     * field indicates whether this path is such a line through the middle, that
+     * has no companion line going back on the other side and is not a closed
+     * loop.
+     */
+    bool is_odd;
+
+    /*!
+     * Which region this line is part of. A solid polygon without holes has only one region.
+     * A polygon with holes has 2. Disconnected parts of the polygon are also separate regions.
+     * Will be 0 if no region was given.
+     */
+    size_t region_id;
+
+    /*!
+     * The list of vertices along which this path runs.
+     *
+     * Each junction has a width, making this path a variable-width path.
+     */
+    std::vector<ExtrusionJunction> junctions;
+
+    ExtrusionLine(const size_t inset_idx, const bool is_odd, const size_t region_id = 0);
+
+    /*!
+     * Sum the total length of this path.
+     */
+    coord_t getLength() const;
+
+    /*!
+     * Get the minimal width of this path
+     */
+    coord_t getMinimalWidth() const;
+
+    /*!
+     * Export the included junctions as vector.
+     */
+    void appendJunctionsTo(LineJunctions& result) const;
+
+    /*!
+     * Removes vertices of the ExtrusionLines to make sure that they are not too high
+     * resolution.
+     *
+     * This removes junctions which are connected to line segments that are shorter
+     * than the `smallest_line_segment`, unless that would introduce a deviation
+     * in the contour of more than `allowed_error_distance`.
+     *
+     * Criteria:
+     * 1. Never remove a junction if either of the connected segments is larger than \p smallest_line_segment
+     * 2. Never remove a junction if the distance between that junction and the final resulting polygon would be higher
+     *    than \p allowed_error_distance
+     * 3. The direction of segments longer than \p smallest_line_segment always
+     *    remains unaltered (but their end points may change if it is connected to
+     *    a small segment)
+     * 4. Never remove a junction if it has a distinctively different width than the next junction, as this can
+     *    introduce unwanted irregularities on the wall widths.
+     *
+     * Simplify uses a heuristic and doesn't necessarily remove all removable
+     * vertices under the above criteria, but simplify may never violate these
+     * criteria. Unless the segments or the distance is smaller than the
+     * rounding error of 5 micron.
+     *
+     * Vertices which introduce an error of less than 5 microns are removed
+     * anyway, even if the segments are longer than the smallest line segment.
+     * This makes sure that (practically) co-linear line segments are joined into
+     * a single line segment.
+     * \param smallest_line_segment Maximal length of removed line segments.
+     * \param allowed_error_distance If removing a vertex introduces a deviation
+     *         from the original path that is more than this distance, the vertex may
+     *         not be removed.
+     * \param maximum_extrusion_area_deviation The maximum extrusion area deviation allowed when removing intermediate
+     *        junctions from a straight ExtrusionLine
+     */
+    void simplify(int64_t smallest_line_segment_squared, int64_t allowed_error_distance_squared, int64_t maximum_extrusion_area_deviation);
+
+    /*!
+     * Computes and returns the total area error (in μm²) of the AB and BC segments of an ABC straight ExtrusionLine
+     * when the junction B with a width B.w is removed from the ExtrusionLine. The area changes due to the fact that the
+     * new simplified line AC has a uniform width which equals to the weighted average of the width of the subsegments
+     * (based on their length).
+     *
+     * \param A Start point of the 3-point-straight line
+     * \param B Intermediate point of the 3-point-straight line
+     * \param C End point of the 3-point-straight line
+     * \param weighted_average_width The weighted average of the widths of the two colinear extrusion segments
+     * */
+    static int64_t calculateExtrusionAreaDeviationError(ExtrusionJunction A, ExtrusionJunction B, ExtrusionJunction C, coord_t& weighted_average_width);
+};
+
+using VariableWidthLines = std::vector<ExtrusionLine>; //<! The ExtrusionLines generated by libArachne for each Path
+using VariableWidthPaths = std::vector<VariableWidthLines>; //<! The toolpaths generated by libArachne
+
+static inline Slic3r::ThickPolyline to_thick_polyline(const Arachne::LineJunctions &line_junctions)
+{
+    assert(line_junctions.size() >= 2);
+    Slic3r::ThickPolyline out;
+    out.points.emplace_back(line_junctions.front().p);
+    out.width.emplace_back(line_junctions.front().w);
+    out.points.emplace_back(line_junctions[1].p);
+    out.width.emplace_back(line_junctions[1].w);
+
+    auto it_prev = line_junctions.begin() + 1;
+    for (auto it = line_junctions.begin() + 2; it != line_junctions.end(); ++it) {
+        out.points.emplace_back(it->p);
+        out.width.emplace_back(it_prev->w);
+        out.width.emplace_back(it->w);
+        it_prev = it;
+    }
+
+    return out;
+}
+
+static inline Polygon to_polygon(const ExtrusionLine& line) {
+    Polygon out;
+    assert(line.junctions.size() >= 3);
+    assert(line.junctions.front().p == line.junctions.back().p);
+    out.points.reserve(line.junctions.size() - 1);
+    for (auto it = line.junctions.begin(); it != line.junctions.end() - 1; ++it)
+        out.points.emplace_back(it->p);
+    return out;
+}
+
+} // namespace Slic3r::Arachne
+#endif // UTILS_EXTRUSION_LINE_H

src/libslic3r/Arachne/utils/HalfEdge.hpp

@@ -0,0 +1,39 @@
+//Copyright (c) 2020 Ultimaker B.V.
+//CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#ifndef UTILS_HALF_EDGE_H
+#define UTILS_HALF_EDGE_H
+
+#include <forward_list>
+#include <optional>
+
+namespace Slic3r::Arachne
+{
+
+template<typename node_data_t, typename edge_data_t, typename derived_node_t, typename derived_edge_t>
+class HalfEdgeNode;
+
+
+template<typename node_data_t, typename edge_data_t, typename derived_node_t, typename derived_edge_t>
+class HalfEdge
+{
+    using edge_t = derived_edge_t;
+    using node_t = derived_node_t;
+public:
+    edge_data_t data;
+    edge_t* twin = nullptr;
+    edge_t* next = nullptr;
+    edge_t* prev = nullptr;
+    node_t* from = nullptr;
+    node_t* to = nullptr;
+    HalfEdge(edge_data_t data)
+    : data(data)
+    {}
+    bool operator==(const edge_t& other)
+    {
+        return this == &other;
+    }
+};
+
+} // namespace Slic3r::Arachne
+#endif // UTILS_HALF_EDGE_H

src/libslic3r/Arachne/utils/HalfEdgeGraph.hpp

@@ -0,0 +1,29 @@
+//Copyright (c) 2020 Ultimaker B.V.
+//CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#ifndef UTILS_HALF_EDGE_GRAPH_H
+#define UTILS_HALF_EDGE_GRAPH_H
+
+
+#include <list>
+#include <cassert>
+
+
+
+#include "HalfEdge.hpp"
+#include "HalfEdgeNode.hpp"
+
+namespace Slic3r::Arachne
+{
+template<class node_data_t, class edge_data_t, class derived_node_t, class derived_edge_t> // types of data contained in nodes and edges
+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;
+};
+
+} // namespace Slic3r::Arachne
+#endif // UTILS_HALF_EDGE_GRAPH_H

src/libslic3r/Arachne/utils/HalfEdgeNode.hpp

@@ -0,0 +1,38 @@
+//Copyright (c) 2020 Ultimaker B.V.
+//CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#ifndef UTILS_HALF_EDGE_NODE_H
+#define UTILS_HALF_EDGE_NODE_H
+
+#include <list>
+
+#include "../../Point.hpp"
+
+namespace Slic3r::Arachne
+{
+
+template<typename node_data_t, typename edge_data_t, typename derived_node_t, typename derived_edge_t>
+class HalfEdge;
+
+template<typename node_data_t, typename edge_data_t, typename derived_node_t, typename derived_edge_t>
+class HalfEdgeNode
+{
+    using edge_t = derived_edge_t;
+    using node_t = derived_node_t;
+public:
+    node_data_t data;
+    Point p;
+    edge_t* incident_edge = nullptr;
+    HalfEdgeNode(node_data_t data, Point p)
+    : data(data)
+    , p(p)
+    {}
+
+    bool operator==(const node_t& other)
+    {
+        return this == &other;
+    }
+};
+
+} // namespace Slic3r::Arachne
+#endif // UTILS_HALF_EDGE_NODE_H

src/libslic3r/Arachne/utils/PolygonsPointIndex.hpp

@@ -0,0 +1,131 @@
+//Copyright (c) 2018 Ultimaker B.V.
+//CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#ifndef UTILS_POLYGONS_POINT_INDEX_H
+#define UTILS_POLYGONS_POINT_INDEX_H
+
+#include <vector>
+
+#include "../../Point.hpp"
+#include "../../Polygon.hpp"
+
+
+namespace Slic3r::Arachne
+{
+
+/*!
+ * A class for iterating over the points in one of the polygons in a \ref Polygons object
+ */
+class PolygonsPointIndex
+{
+public:
+    /*!
+     * The polygons into which this index is indexing.
+     */
+    const Polygons* polygons; // (pointer to const polygons)
+
+    unsigned int poly_idx; //!< The index of the polygon in \ref PolygonsPointIndex::polygons
+
+    unsigned int point_idx; //!< The index of the point in the polygon in \ref PolygonsPointIndex::polygons
+
+    /*!
+     * Constructs an empty point index to no polygon.
+     *
+     * This is used as a placeholder for when there is a zero-construction
+     * needed. Since the `polygons` field is const you can't ever make this
+     * initialisation useful.
+     */
+    PolygonsPointIndex() : polygons(nullptr), poly_idx(0), point_idx(0) {}
+
+    /*!
+     * Constructs a new point index to a vertex of a polygon.
+     * \param polygons The Polygons instance to which this index points.
+     * \param poly_idx The index of the sub-polygon to point to.
+     * \param point_idx The index of the vertex in the sub-polygon.
+     */
+    PolygonsPointIndex(const Polygons *polygons, unsigned int poly_idx, unsigned int point_idx)
+        : polygons(polygons), poly_idx(poly_idx), point_idx(point_idx) {}
+
+    /*!
+     * Copy constructor to copy these indices.
+     */
+    PolygonsPointIndex(const PolygonsPointIndex& original) = default;
+
+    Point p() const
+    {
+        if (!polygons)
+            return {0, 0};
+
+        return (*polygons)[poly_idx][point_idx];
+    }
+
+    /*!
+     * Test whether two iterators refer to the same polygon in the same polygon list.
+     * 
+     * \param other The PolygonsPointIndex to test for equality
+     * \return Wether the right argument refers to the same polygon in the same ListPolygon as the left argument.
+     */
+    bool operator==(const PolygonsPointIndex &other) const
+    {
+        return polygons == other.polygons && poly_idx == other.poly_idx && point_idx == other.point_idx;
+    }
+    bool                operator!=(const PolygonsPointIndex &other) const { return !(*this == other); }
+    bool                operator<(const PolygonsPointIndex &other) const { return this->p() < other.p(); }
+    PolygonsPointIndex &operator=(const PolygonsPointIndex &other)
+    {
+        polygons  = other.polygons;
+        poly_idx  = other.poly_idx;
+        point_idx = other.point_idx;
+        return *this;
+    }
+    //! move the iterator forward (and wrap around at the end)
+    PolygonsPointIndex &operator++()
+    {
+        point_idx = (point_idx + 1) % (*polygons)[poly_idx].size();
+        return *this;
+    }
+    //! move the iterator backward (and wrap around at the beginning)
+    PolygonsPointIndex &operator--()
+    {
+        if (point_idx == 0)
+            point_idx = (*polygons)[poly_idx].size();
+        point_idx--;
+        return *this;
+    }
+    //! move the iterator forward (and wrap around at the end)
+    PolygonsPointIndex next() const
+    {
+        PolygonsPointIndex ret(*this);
+        ++ret;
+        return ret;
+    }
+    //! move the iterator backward (and wrap around at the beginning)
+    PolygonsPointIndex prev() const
+    {
+        PolygonsPointIndex ret(*this);
+        --ret;
+        return ret;
+    }
+};
+
+
+}//namespace Slic3r::Arachne
+
+namespace std
+{
+/*!
+ * Hash function for \ref PolygonsPointIndex
+ */
+template <>
+struct hash<Slic3r::Arachne::PolygonsPointIndex>
+{
+    size_t operator()(const Slic3r::Arachne::PolygonsPointIndex& lpi) const
+    {
+        return Slic3r::PointHash{}(lpi.p());
+    }
+};
+}//namespace std
+
+
+
+#endif//UTILS_POLYGONS_POINT_INDEX_H

src/libslic3r/Arachne/utils/PolygonsSegmentIndex.hpp

@@ -0,0 +1,31 @@
+//Copyright (c) 2020 Ultimaker B.V.
+//CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#ifndef UTILS_POLYGONS_SEGMENT_INDEX_H
+#define UTILS_POLYGONS_SEGMENT_INDEX_H
+
+#include <vector>
+
+#include "PolygonsPointIndex.hpp"
+
+namespace Slic3r::Arachne
+{
+
+/*!
+ * A class for iterating over the points in one of the polygons in a \ref Polygons object
+ */
+class PolygonsSegmentIndex : public PolygonsPointIndex
+{
+public:
+    PolygonsSegmentIndex() : PolygonsPointIndex(){};
+    PolygonsSegmentIndex(const Polygons *polygons, unsigned int poly_idx, unsigned int point_idx) : PolygonsPointIndex(polygons, poly_idx, point_idx){};
+
+    Point from() const { return PolygonsPointIndex::p(); }
+
+    Point to() const { return PolygonsSegmentIndex::next().p(); }
+};
+
+} // namespace Slic3r::Arachne
+
+
+#endif//UTILS_POLYGONS_SEGMENT_INDEX_H

src/libslic3r/Arachne/utils/SparseGrid.hpp

@@ -0,0 +1,133 @@
+//Copyright (c) 2016 Scott Lenser
+//Copyright (c) 2018 Ultimaker B.V.
+//CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#ifndef UTILS_SPARSE_GRID_H
+#define UTILS_SPARSE_GRID_H
+
+#include <cassert>
+#include <unordered_map>
+#include <vector>
+#include <functional>
+
+#include "../../Point.hpp"
+#include "SquareGrid.hpp"
+
+namespace Slic3r::Arachne {
+
+/*! \brief Sparse grid which can locate spatially nearby elements efficiently.
+ * 
+ * \note This is an abstract template class which doesn't have any functions to insert elements.
+ * \see SparsePointGrid
+ *
+ * \tparam ElemT The element type to store.
+ */
+template<class ElemT> class SparseGrid : public SquareGrid
+{
+public:
+    using Elem = ElemT;
+
+    using GridPoint    = SquareGrid::GridPoint;
+    using grid_coord_t = SquareGrid::grid_coord_t;
+    using GridMap       = std::unordered_multimap<GridPoint, Elem, PointHash>;
+
+    using iterator       = typename GridMap::iterator;
+    using const_iterator = typename GridMap::const_iterator;
+
+    /*! \brief Constructs a sparse grid with the specified cell size.
+     *
+     * \param[in] cell_size The size to use for a cell (square) in the grid.
+     *    Typical values would be around 0.5-2x of expected query radius.
+     * \param[in] elem_reserve Number of elements to research space for.
+     * \param[in] max_load_factor Maximum average load factor before rehashing.
+     */
+    SparseGrid(coord_t cell_size, size_t elem_reserve=0U, float max_load_factor=1.0f);
+
+    iterator begin() { return m_grid.begin(); }
+    iterator end() { return m_grid.end(); }
+    const_iterator begin() const { return m_grid.begin(); }
+    const_iterator end() const { return m_grid.end(); }
+
+    /*! \brief Returns all data within radius of query_pt.
+     *
+     * Finds all elements with location within radius of \p query_pt.  May
+     * return additional elements that are beyond radius.
+     *
+     * Average running time is a*(1 + 2 * radius / cell_size)**2 +
+     * b*cnt where a and b are proportionality constance and cnt is
+     * the number of returned items.  The search will return items in
+     * an area of (2*radius + cell_size)**2 on average.  The max range
+     * of an item from the query_point is radius + cell_size.
+     *
+     * \param[in] query_pt The point to search around.
+     * \param[in] radius The search radius.
+     * \return Vector of elements found
+     */
+    std::vector<Elem> getNearby(const Point &query_pt, coord_t radius) const;
+
+    /*! \brief Process elements from cells that might contain sought after points.
+     *
+     * Processes elements from cell that might have elements within \p
+     * radius of \p query_pt.  Processes all elements that are within
+     * radius of query_pt.  May process elements that are up to radius +
+     * cell_size from query_pt.
+     *
+     * \param[in] query_pt The point to search around.
+     * \param[in] radius The search radius.
+     * \param[in] process_func Processes each element.  process_func(elem) is
+     *    called for each element in the cell. Processing stops if function returns false.
+     * \return Whether we need to continue processing after this function
+     */
+    bool processNearby(const Point &query_pt, coord_t radius, const std::function<bool(const ElemT &)> &process_func) const;
+
+protected:
+    /*! \brief Process elements from the cell indicated by \p grid_pt.
+     *
+     * \param[in] grid_pt The grid coordinates of the cell.
+     * \param[in] process_func Processes each element.  process_func(elem) is
+     *    called for each element in the cell. Processing stops if function returns false.
+     * \return Whether we need to continue processing a next cell.
+     */
+    bool processFromCell(const GridPoint &grid_pt, const std::function<bool(const Elem &)> &process_func) const;
+
+    /*! \brief Map from grid locations (GridPoint) to elements (Elem). */
+    GridMap m_grid;
+};
+
+template<class ElemT> SparseGrid<ElemT>::SparseGrid(coord_t cell_size, size_t elem_reserve, float max_load_factor) : SquareGrid(cell_size)
+{
+    // Must be before the reserve call.
+    m_grid.max_load_factor(max_load_factor);
+    if (elem_reserve != 0U)
+        m_grid.reserve(elem_reserve);
+}
+
+template<class ElemT> bool SparseGrid<ElemT>::processFromCell(const GridPoint &grid_pt, const std::function<bool(const Elem &)> &process_func) const
+{
+    auto grid_range = m_grid.equal_range(grid_pt);
+    for (auto iter = grid_range.first; iter != grid_range.second; ++iter)
+        if (!process_func(iter->second))
+            return false;
+    return true;
+}
+
+template<class ElemT>
+bool SparseGrid<ElemT>::processNearby(const Point &query_pt, coord_t radius, const std::function<bool(const Elem &)> &process_func) const
+{
+    return SquareGrid::processNearby(query_pt, radius, [&process_func, this](const GridPoint &grid_pt) { return processFromCell(grid_pt, process_func); });
+}
+
+template<class ElemT> std::vector<typename SparseGrid<ElemT>::Elem> SparseGrid<ElemT>::getNearby(const Point &query_pt, coord_t radius) const
+{
+    std::vector<Elem>                       ret;
+    const std::function<bool(const Elem &)> process_func = [&ret](const Elem &elem) {
+        ret.push_back(elem);
+        return true;
+    };
+    processNearby(query_pt, radius, process_func);
+    return ret;
+}
+
+} // namespace Slic3r::Arachne
+
+#endif // UTILS_SPARSE_GRID_H

src/libslic3r/Arachne/utils/SparseLineGrid.hpp

@@ -0,0 +1,77 @@
+//Copyright (c) 2018 Ultimaker B.V.
+//CuraEngine is released under the terms of the AGPLv3 or higher.
+
+
+#ifndef UTILS_SPARSE_LINE_GRID_H
+#define UTILS_SPARSE_LINE_GRID_H
+
+#include <cassert>
+#include <unordered_map>
+#include <vector>
+#include <functional>
+
+#include "SparseGrid.hpp"
+
+namespace Slic3r::Arachne {
+
+/*! \brief Sparse grid which can locate spatially nearby elements efficiently.
+ *
+ * \tparam ElemT The element type to store.
+ * \tparam Locator The functor to get the start and end locations from ElemT.
+ *    must have: std::pair<Point, Point> operator()(const ElemT &elem) const
+ *    which returns the location associated with val.
+ */
+template<class ElemT, class Locator> class SparseLineGrid : public SparseGrid<ElemT>
+{
+public:
+    using Elem = ElemT;
+
+    /*! \brief Constructs a sparse grid with the specified cell size.
+     *
+     * \param[in] cell_size The size to use for a cell (square) in the grid.
+     *    Typical values would be around 0.5-2x of expected query radius.
+     * \param[in] elem_reserve Number of elements to research space for.
+     * \param[in] max_load_factor Maximum average load factor before rehashing.
+     */
+    SparseLineGrid(coord_t cell_size, size_t elem_reserve = 0U, float max_load_factor = 1.0f);
+
+    /*! \brief Inserts elem into the sparse grid.
+     *
+     * \param[in] elem The element to be inserted.
+     */
+    void insert(const Elem &elem);
+
+protected:
+    using GridPoint = typename SparseGrid<ElemT>::GridPoint;
+
+    /*! \brief Accessor for getting locations from elements. */
+    Locator m_locator;
+};
+
+template<class ElemT, class Locator>
+SparseLineGrid<ElemT, Locator>::SparseLineGrid(coord_t cell_size, size_t elem_reserve, float max_load_factor)
+    : SparseGrid<ElemT>(cell_size, elem_reserve, max_load_factor) {}
+
+template<class ElemT, class Locator> void SparseLineGrid<ElemT, Locator>::insert(const Elem &elem)
+{
+    const std::pair<Point, Point> line = m_locator(elem);
+    using GridMap                      = std::unordered_multimap<GridPoint, Elem, PointHash>;
+    // below is a workaround for the fact that lambda functions cannot access private or protected members
+    // first we define a lambda which works on any GridMap and then we bind it to the actual protected GridMap of the parent class
+    std::function<bool(GridMap *, const GridPoint)> process_cell_func_ = [&elem](GridMap *m_grid, const GridPoint grid_loc) {
+        m_grid->emplace(grid_loc, elem);
+        return true;
+    };
+    using namespace std::placeholders; // for _1, _2, _3...
+    GridMap                             *m_grid = &(this->m_grid);
+    std::function<bool(const GridPoint)> process_cell_func(std::bind(process_cell_func_, m_grid, _1));
+
+    SparseGrid<ElemT>::processLineCells(line, process_cell_func);
+}
+
+#undef SGI_TEMPLATE
+#undef SGI_THIS
+
+} // namespace Slic3r::Arachne
+
+#endif // UTILS_SPARSE_LINE_GRID_H

src/libslic3r/Arachne/utils/SquareGrid.cpp

@@ -0,0 +1,147 @@
+//Copyright (c) 2021 Ultimaker B.V.
+//CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#include "SquareGrid.hpp"
+#include "../../Point.hpp"
+
+using namespace Slic3r::Arachne;
+
+
+SquareGrid::SquareGrid(coord_t cell_size) : cell_size(cell_size)
+{
+    assert(cell_size > 0U);
+}
+
+
+SquareGrid::GridPoint SquareGrid::toGridPoint(const Vec2i64 &point)  const
+{
+    return Point(toGridCoord(point.x()), toGridCoord(point.y()));
+}
+
+
+SquareGrid::grid_coord_t SquareGrid::toGridCoord(const int64_t &coord)  const
+{
+    // This mapping via truncation results in the cells with
+    // GridPoint.x==0 being twice as large and similarly for
+    // GridPoint.y==0.  This doesn't cause any incorrect behavior,
+    // just changes the running time slightly.  The change in running
+    // time from this is probably not worth doing a proper floor
+    // operation.
+    return coord / cell_size;
+}
+
+coord_t SquareGrid::toLowerCoord(const grid_coord_t& grid_coord)  const
+{
+    // This mapping via truncation results in the cells with
+    // GridPoint.x==0 being twice as large and similarly for
+    // GridPoint.y==0.  This doesn't cause any incorrect behavior,
+    // just changes the running time slightly.  The change in running
+    // time from this is probably not worth doing a proper floor
+    // operation.
+    return grid_coord * cell_size;
+}
+
+
+bool SquareGrid::processLineCells(const std::pair<Point, Point> line, const std::function<bool (GridPoint)>& process_cell_func)
+{
+    return static_cast<const SquareGrid*>(this)->processLineCells(line, process_cell_func);
+}
+
+
+bool SquareGrid::processLineCells(const std::pair<Point, Point> line, const std::function<bool (GridPoint)>& process_cell_func) const
+{
+    Point start = line.first;
+    Point end = line.second;
+    if (end.x() < start.x())
+    { // make sure X increases between start and end
+        std::swap(start, end);
+    }
+
+    const GridPoint start_cell = toGridPoint(start.cast<int64_t>());
+    const GridPoint end_cell = toGridPoint(end.cast<int64_t>());
+    const int64_t y_diff = int64_t(end.y() - start.y());
+    const grid_coord_t y_dir = nonzeroSign(y_diff);
+
+    /* This line drawing algorithm iterates over the range of Y coordinates, and
+    for each Y coordinate computes the range of X coordinates crossed in one
+    unit of Y. These ranges are rounded to be inclusive, so effectively this
+    creates a "fat" line, marking more cells than a strict one-cell-wide path.*/
+    grid_coord_t x_cell_start = start_cell.x();
+    for (grid_coord_t cell_y = start_cell.y(); cell_y * y_dir <= end_cell.y() * y_dir; cell_y += y_dir)
+    { // for all Y from start to end
+        // nearest y coordinate of the cells in the next row
+        const coord_t nearest_next_y = toLowerCoord(cell_y + ((nonzeroSign(cell_y) == y_dir || cell_y == 0) ? y_dir : coord_t(0)));
+        grid_coord_t x_cell_end; // the X coord of the last cell to include from this row
+        if (y_diff == 0)
+        {
+            x_cell_end = end_cell.x();
+        }
+        else
+        {
+            const int64_t area = int64_t(end.x() - start.x()) * int64_t(nearest_next_y - start.y());
+            // corresponding_x: the x coordinate corresponding to nearest_next_y
+            int64_t corresponding_x = int64_t(start.x()) + area / y_diff;
+            x_cell_end = toGridCoord(corresponding_x + ((corresponding_x < 0) && ((area % y_diff) != 0)));
+            if (x_cell_end < start_cell.x())
+            { // process at least one cell!
+                x_cell_end = x_cell_start;
+            }
+        }
+
+        for (grid_coord_t cell_x = x_cell_start; cell_x <= x_cell_end; ++cell_x)
+        {
+            GridPoint grid_loc(cell_x, cell_y);
+            if (! process_cell_func(grid_loc))
+            {
+                return false;
+            }
+            if (grid_loc == end_cell)
+            {
+                return true;
+            }
+        }
+        // TODO: this causes at least a one cell overlap for each row, which
+        // includes extra cells when crossing precisely on the corners
+        // where positive slope where x > 0 and negative slope where x < 0
+        x_cell_start = x_cell_end;
+    }
+    assert(false && "We should have returned already before here!");
+    return false;
+}
+
+bool SquareGrid::processNearby
+(
+    const Point &query_pt,
+    coord_t radius,
+    const std::function<bool (const GridPoint&)>& process_func
+) const
+{
+    const Point min_loc(query_pt.x() - radius, query_pt.y() - radius);
+    const Point max_loc(query_pt.x() + radius, query_pt.y() + radius);
+
+    GridPoint min_grid = toGridPoint(min_loc.cast<int64_t>());
+    GridPoint max_grid = toGridPoint(max_loc.cast<int64_t>());
+
+    for (coord_t grid_y = min_grid.y(); grid_y <= max_grid.y(); ++grid_y)
+    {
+        for (coord_t grid_x = min_grid.x(); grid_x <= max_grid.x(); ++grid_x)
+        {
+            GridPoint grid_pt(grid_x,grid_y);
+            if (!process_func(grid_pt))
+            {
+                return false;
+            }
+        }
+    }
+    return true;
+}
+
+SquareGrid::grid_coord_t SquareGrid::nonzeroSign(const grid_coord_t z) const
+{
+    return (z >= 0) - (z < 0);
+}
+
+coord_t SquareGrid::getCellSize() const
+{
+    return cell_size;
+}

src/libslic3r/Arachne/utils/SquareGrid.hpp

@@ -0,0 +1,110 @@
+//Copyright (c) 2021 Ultimaker B.V.
+//CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#ifndef UTILS_SQUARE_GRID_H
+#define UTILS_SQUARE_GRID_H
+
+#include "../../Point.hpp"
+
+#include <cassert>
+#include <unordered_map>
+#include <vector>
+#include <functional>
+
+namespace Slic3r::Arachne {
+
+/*!
+ * Helper class to calculate coordinates on a square grid, and providing some
+ * utility functions to process grids.
+ *
+ * Doesn't contain any data, except cell size. The purpose is only to
+ * automatically generate coordinates on a grid, and automatically feed them to
+ * functions.
+ * The grid is theoretically infinite (bar integer limits).
+ */
+class SquareGrid
+{
+public:
+    /*! \brief Constructs a grid with the specified cell size.
+     * \param[in] cell_size The size to use for a cell (square) in the grid.
+     */
+    SquareGrid(const coord_t cell_size);
+
+    /*!
+     * Get the cell size this grid was created for.
+     */
+    coord_t getCellSize() const;
+
+    using GridPoint = Point;
+    using grid_coord_t = coord_t;
+
+    /*! \brief Process cells along a line indicated by \p line.
+     *
+     * \param line The line along which to process cells.
+     * \param process_func Processes each cell. ``process_func(elem)`` is called
+     * for each cell. Processing stops if function returns false.
+     * \return Whether we need to continue processing after this function.
+     */
+    bool processLineCells(const std::pair<Point, Point> line, const std::function<bool (GridPoint)>& process_cell_func);
+
+    /*! \brief Process cells along a line indicated by \p line.
+     *
+     * \param line The line along which to process cells
+     * \param process_func Processes each cell. ``process_func(elem)`` is called
+     * for each cell. Processing stops if function returns false.
+     * \return Whether we need to continue processing after this function.
+     */
+    bool processLineCells(const std::pair<Point, Point> line, const std::function<bool (GridPoint)>& process_cell_func) const;
+
+    /*! \brief Process cells that might contain sought after points.
+     *
+     * Processes cells that might be within a square with twice \p radius as
+     * width, centered around \p query_pt.
+     * May process elements that are up to radius + cell_size from query_pt.
+     * \param query_pt The point to search around.
+     * \param radius The search radius.
+     * \param process_func Processes each cell. ``process_func(loc)`` is called
+     * for each cell coord within range. Processing stops if function returns
+     * ``false``.
+     * \return Whether we need to continue processing after this function.
+     */
+    bool processNearby(const Point &query_pt, coord_t radius, const std::function<bool(const GridPoint &)> &process_func) const;
+
+    /*! \brief Compute the grid coordinates of a point.
+     * \param point The actual location.
+     * \return The grid coordinates that correspond to \p point.
+     */
+    GridPoint toGridPoint(const Vec2i64 &point) const;
+
+    /*! \brief Compute the grid coordinate of a real space coordinate.
+     * \param coord The actual location.
+     * \return The grid coordinate that corresponds to \p coord.
+     */
+    grid_coord_t toGridCoord(const int64_t &coord) const;
+
+    /*! \brief Compute the lowest coord in a grid cell.
+     * The lowest point is the point in the grid cell closest to the origin.
+     *
+     * \param grid_coord The grid coordinate.
+     * \return The print space coordinate that corresponds to \p grid_coord.
+     */
+    coord_t toLowerCoord(const grid_coord_t &grid_coord) const;
+
+protected:
+    /*! \brief The cell (square) size. */
+    coord_t cell_size;
+
+    /*!
+     * Compute the sign of a number.
+     *
+     * The number 0 will result in a positive sign (1).
+     * \param z The number to find the sign of.
+     * \return 1 if the number is positive or 0, or -1 if the number is
+     * negative.
+     */
+    grid_coord_t nonzeroSign(grid_coord_t z) const;
+};
+
+} // namespace Slic3r::Arachne
+
+#endif //UTILS_SQUARE_GRID_H

src/libslic3r/Arachne/utils/VoronoiUtils.cpp

@@ -0,0 +1,250 @@
+//Copyright (c) 2021 Ultimaker B.V.
+//CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#include <stack>
+#include <optional>
+#include <boost/log/trivial.hpp>
+
+#include "linearAlg2D.hpp"
+#include "VoronoiUtils.hpp"
+
+namespace Slic3r::Arachne
+{
+
+Vec2i64 VoronoiUtils::p(const vd_t::vertex_type *node)
+{
+    const double x = node->x();
+    const double y = node->y();
+    assert(x <= double(std::numeric_limits<int64_t>::max()) && x >= std::numeric_limits<int64_t>::lowest());
+    assert(y <= double(std::numeric_limits<int64_t>::max()) && y >= std::numeric_limits<int64_t>::lowest());
+    return Vec2i64(int64_t(x + 0.5 - (x < 0)), int64_t(y + 0.5 - (y < 0))); // Round to the nearest integer coordinates.
+}
+
+Point VoronoiUtils::getSourcePoint(const vd_t::cell_type& cell, const std::vector<Segment>& segments)
+{
+    assert(cell.contains_point());
+    if(!cell.contains_point())
+        BOOST_LOG_TRIVIAL(debug) << "Voronoi cell doesn't contain a source point!";
+
+    switch (cell.source_category()) {
+        case boost::polygon::SOURCE_CATEGORY_SINGLE_POINT:
+            assert(false && "Voronoi diagram is always constructed using segments, so cell.source_category() shouldn't be SOURCE_CATEGORY_SINGLE_POINT!\n");
+            BOOST_LOG_TRIVIAL(error) << "Voronoi diagram is always constructed using segments, so cell.source_category() shouldn't be SOURCE_CATEGORY_SINGLE_POINT!";
+            break;
+        case boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT:
+            assert(cell.source_index() < segments.size());
+            return segments[cell.source_index()].to();
+            break;
+        case boost::polygon::SOURCE_CATEGORY_SEGMENT_END_POINT:
+            assert(cell.source_index() < segments.size());
+            return segments[cell.source_index()].from();
+            break;
+        default:
+            assert(false && "getSourcePoint should only be called on point cells!\n");
+            break;
+    }
+
+    assert(false && "cell.source_category() is equal to an invalid value!\n");
+    BOOST_LOG_TRIVIAL(error) << "cell.source_category() is equal to an invalid value!";
+    return {};
+}
+
+PolygonsPointIndex VoronoiUtils::getSourcePointIndex(const vd_t::cell_type& cell, const std::vector<Segment>& segments)
+{
+    assert(cell.contains_point());
+    if(!cell.contains_point())
+        BOOST_LOG_TRIVIAL(debug) << "Voronoi cell doesn't contain a source point!";
+
+    assert(cell.source_category() != boost::polygon::SOURCE_CATEGORY_SINGLE_POINT);
+    switch (cell.source_category()) {
+        case boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT: {
+            assert(cell.source_index() < segments.size());
+            PolygonsPointIndex ret = segments[cell.source_index()];
+            ++ret;
+            return ret;
+            break;
+        }
+        case boost::polygon::SOURCE_CATEGORY_SEGMENT_END_POINT: {
+            assert(cell.source_index() < segments.size());
+            return segments[cell.source_index()];
+            break;
+        }
+        default:
+            assert(false && "getSourcePoint should only be called on point cells!\n");
+            break;
+    }
+    PolygonsPointIndex ret = segments[cell.source_index()];
+    return ++ret;
+}
+
+const VoronoiUtils::Segment &VoronoiUtils::getSourceSegment(const vd_t::cell_type &cell, const std::vector<Segment> &segments)
+{
+    assert(cell.contains_segment());
+    if (!cell.contains_segment())
+        BOOST_LOG_TRIVIAL(debug) << "Voronoi cell doesn't contain a source segment!";
+
+    return segments[cell.source_index()];
+}
+
+class PointMatrix
+{
+public:
+    double matrix[4];
+
+    PointMatrix()
+    {
+        matrix[0] = 1;
+        matrix[1] = 0;
+        matrix[2] = 0;
+        matrix[3] = 1;
+    }
+
+    PointMatrix(double rotation)
+    {
+        rotation = rotation / 180 * M_PI;
+        matrix[0] = cos(rotation);
+        matrix[1] = -sin(rotation);
+        matrix[2] = -matrix[1];
+        matrix[3] = matrix[0];
+    }
+
+    PointMatrix(const Point p)
+    {
+        matrix[0] = p.x();
+        matrix[1] = p.y();
+        double f = sqrt((matrix[0] * matrix[0]) + (matrix[1] * matrix[1]));
+        matrix[0] /= f;
+        matrix[1] /= f;
+        matrix[2] = -matrix[1];
+        matrix[3] = matrix[0];
+    }
+
+    static PointMatrix scale(double s)
+    {
+        PointMatrix ret;
+        ret.matrix[0] = s;
+        ret.matrix[3] = s;
+        return ret;
+    }
+
+    Point apply(const Point p) const
+    {
+        return Point(coord_t(p.x() * matrix[0] + p.y() * matrix[1]), coord_t(p.x() * matrix[2] + p.y() * matrix[3]));
+    }
+
+    Point unapply(const Point p) const
+    {
+        return Point(coord_t(p.x() * matrix[0] + p.y() * matrix[2]), coord_t(p.x() * matrix[1] + p.y() * matrix[3]));
+    }
+};
+std::vector<Point> VoronoiUtils::discretizeParabola(const Point& p, const Segment& segment, Point s, Point e, coord_t approximate_step_size, float transitioning_angle)
+{
+    std::vector<Point> discretized;
+    // x is distance of point projected on the segment ab
+    // xx is point projected on the segment ab
+    const Point a = segment.from();
+    const Point b = segment.to();
+    const Point ab = b - a;
+    const Point as = s - a;
+    const Point ae = e - a;
+    const coord_t ab_size = ab.cast<int64_t>().norm();
+    const coord_t sx = as.cast<int64_t>().dot(ab.cast<int64_t>()) / ab_size;
+    const coord_t ex = ae.cast<int64_t>().dot(ab.cast<int64_t>()) / ab_size;
+    const coord_t sxex = ex - sx;
+
+    assert((as.cast<int64_t>().dot(ab.cast<int64_t>()) / int64_t(ab_size)) <= std::numeric_limits<coord_t>::max());
+    assert((ae.cast<int64_t>().dot(ab.cast<int64_t>()) / int64_t(ab_size)) <= std::numeric_limits<coord_t>::max());
+
+    const Point ap = p - a;
+    const coord_t px = ap.cast<int64_t>().dot(ab.cast<int64_t>()) / ab_size;
+
+    assert((ap.cast<int64_t>().dot(ab.cast<int64_t>()) / int64_t(ab_size)) <= std::numeric_limits<coord_t>::max());
+
+    Point pxx;
+    Line(a, b).distance_to_infinite_squared(p, &pxx);
+    const Point ppxx = pxx - p;
+    const coord_t d = ppxx.cast<int64_t>().norm();
+    const PointMatrix rot = PointMatrix(ppxx.rotate_90_degree_ccw());
+
+    if (d == 0)
+    {
+        discretized.emplace_back(s);
+        discretized.emplace_back(e);
+        return discretized;
+    }
+    
+    const float marking_bound = atan(transitioning_angle * 0.5);
+    int64_t msx = - marking_bound * int64_t(d); // projected marking_start
+    int64_t mex = marking_bound * int64_t(d); // projected marking_end
+
+    assert(msx <= std::numeric_limits<coord_t>::max());
+    assert(double(msx) * double(msx) <= double(std::numeric_limits<int64_t>::max()));
+    assert(mex <= std::numeric_limits<coord_t>::max());
+    assert(double(msx) * double(msx) / double(2 * d) + double(d / 2) <= std::numeric_limits<coord_t>::max());
+
+    const coord_t marking_start_end_h = msx * msx / (2 * d) + d / 2;
+    Point marking_start = rot.unapply(Point(coord_t(msx), marking_start_end_h)) + pxx;
+    Point marking_end = rot.unapply(Point(coord_t(mex), marking_start_end_h)) + pxx;
+    const int dir = (sx > ex) ? -1 : 1;
+    if (dir < 0)
+    {
+        std::swap(marking_start, marking_end);
+        std::swap(msx, mex);
+    }
+    
+    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 = rot.unapply(Point(0, d / 2)) + 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);
+    if(add_marking_start && add_marking_end && !add_apex)
+    {
+        BOOST_LOG_TRIVIAL(warning) << "Failing to discretize parabola! Must add an apex or one of the endpoints.";
+    }
+    
+    const coord_t step_count = static_cast<coord_t>(static_cast<float>(std::abs(ex - sx)) / approximate_step_size + 0.5);
+    
+    discretized.emplace_back(s);
+    for (coord_t step = 1; step < step_count; step++)
+    {
+        assert(double(sxex) * double(step) <= double(std::numeric_limits<int64_t>::max()));
+        const int64_t x = int64_t(sx) + int64_t(sxex) * int64_t(step) / int64_t(step_count) - int64_t(px);
+        assert(double(x) * double(x) <= double(std::numeric_limits<int64_t>::max()));
+        assert(double(x) * double(x) / double(2 * d) + double(d / 2) <= double(std::numeric_limits<int64_t>::max()));
+        const int64_t y = int64_t(x) * int64_t(x) / int64_t(2 * d) + int64_t(d / 2);
+        
+        if (add_marking_start && msx * int64_t(dir) < int64_t(x) * int64_t(dir))
+        {
+            discretized.emplace_back(marking_start);
+            add_marking_start = false;
+        }
+        if (add_apex && int64_t(x) * int64_t(dir) > 0)
+        {
+            discretized.emplace_back(apex);
+            add_apex = false; // only add the apex just before the 
+        }
+        if (add_marking_end && mex * int64_t(dir) < int64_t(x) * int64_t(dir))
+        {
+            discretized.emplace_back(marking_end);
+            add_marking_end = false;
+        }
+        assert(x <= std::numeric_limits<coord_t>::max() && x >= std::numeric_limits<coord_t>::lowest());
+        assert(y <= std::numeric_limits<coord_t>::max() && y >= std::numeric_limits<coord_t>::lowest());
+        const Point result = rot.unapply(Point(x, y)) + pxx;
+        discretized.emplace_back(result);
+    }
+    if (add_apex)
+    {
+        discretized.emplace_back(apex);
+    }
+    if (add_marking_end)
+    {
+        discretized.emplace_back(marking_end);
+    }
+    discretized.emplace_back(e);
+    return discretized;
+}
+
+}//namespace Slic3r::Arachne

src/libslic3r/Arachne/utils/VoronoiUtils.hpp

@@ -0,0 +1,42 @@
+//Copyright (c) 2020 Ultimaker B.V.
+//CuraEngine is released under the terms of the AGPLv3 or higher.
+
+
+#ifndef UTILS_VORONOI_UTILS_H
+#define UTILS_VORONOI_UTILS_H
+
+#include <vector>
+
+
+#include <boost/polygon/voronoi.hpp>
+
+#include "PolygonsSegmentIndex.hpp"
+
+namespace Slic3r::Arachne
+{
+
+/*!
+ */
+class VoronoiUtils
+{
+public:
+    using Segment        = PolygonsSegmentIndex;
+    using voronoi_data_t = double;
+    using vd_t           = boost::polygon::voronoi_diagram<voronoi_data_t>;
+
+    static Point              getSourcePoint(const vd_t::cell_type &cell, const std::vector<Segment> &segments);
+    static const Segment     &getSourceSegment(const vd_t::cell_type &cell, const std::vector<Segment> &segments);
+    static PolygonsPointIndex getSourcePointIndex(const vd_t::cell_type &cell, const std::vector<Segment> &segments);
+
+    static Vec2i64 p(const vd_t::vertex_type *node);
+
+    /*!
+     * Discretize a parabola based on (approximate) step size.
+     * The \p approximate_step_size is measured parallel to the \p source_segment, not along the parabola.
+     */
+    static std::vector<Point> discretizeParabola(const Point &source_point, const Segment &source_segment, Point start, Point end, coord_t approximate_step_size, float transitioning_angle);
+};
+
+} // namespace Slic3r::Arachne
+
+#endif // UTILS_VORONOI_UTILS_H

src/libslic3r/Arachne/utils/linearAlg2D.hpp

@@ -0,0 +1,131 @@
+//Copyright (c) 2020 Ultimaker B.V.
+//CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#ifndef UTILS_LINEAR_ALG_2D_H
+#define UTILS_LINEAR_ALG_2D_H
+
+#include "../../Point.hpp"
+
+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.cast<int64_t>().norm();
+        if (_len < 1)
+            return Point(len, 0);
+        return (p0.cast<int64_t>() * int64_t(len) / _len).cast<coord_t>();
+    };
+
+    auto rotate_90_degree_ccw = [](const Vec2i64 &p) -> Vec2i64 {
+        return Vec2i64(-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 Vec2i64 bq = query_point - b.cast<int64_t>();
+    const Vec2i64 perpendicular = rotate_90_degree_ccw(bq); //The query projects to this perpendicular to coordinate 0.
+
+    assert(ba.cast<double>().dot(perpendicular.cast<double>()) <= double(std::numeric_limits<int64_t>::max()) && ba.cast<double>().dot(perpendicular.cast<double>()) >= double(std::numeric_limits<int64_t>::lowest()));
+    assert(bc.cast<double>().dot(perpendicular.cast<double>()) <= double(std::numeric_limits<int64_t>::max()) && bc.cast<double>().dot(perpendicular.cast<double>()) >= double(std::numeric_limits<int64_t>::lowest()));
+
+    const int64_t project_a_perpendicular = ba.cast<int64_t>().dot(perpendicular); //Project vertex A on the perpendicular line.
+    const int64_t project_c_perpendicular = bc.cast<int64_t>().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.
+    {
+        assert(ba.cast<double>().dot(bq.cast<double>()) <= double(std::numeric_limits<int64_t>::max()) && ba.cast<double>().dot(bq.cast<double>()) >= double(std::numeric_limits<int64_t>::lowest()));
+        assert(bc.cast<double>().dot(bq.cast<double>()) <= double(std::numeric_limits<int64_t>::max()) && bc.cast<double>().dot(bq.cast<double>()) >= double(std::numeric_limits<int64_t>::lowest()));
+
+        const int64_t project_a_parallel = ba.cast<int64_t>().dot(bq); //Project not on the perpendicular, but on the original.
+        const int64_t project_c_parallel = bc.cast<int64_t>().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.
+     * 
+     * The returned value is zero for \p p lying (approximately) on the line going through \p a and \p b
+     * The value is positive for values lying to the left and negative for values lying to the right when looking from \p a to \p b.
+     * 
+     * \param p the point to check
+     * \param a the from point of the line
+     * \param b the to point of the line
+     * \return a positive value when \p p lies to the left of the line from \p a to \p b
+ */
+static inline int64_t pointIsLeftOfLine(const Point& p, const Point& a, const Point& b)
+{
+    return int64_t(b.x() - a.x()) * int64_t(p.y() - a.y()) - int64_t(b.y() - a.y()) * int64_t(p.x() - a.x());
+}
+
+/*!
+     * Compute the angle between two consecutive line segments.
+     *
+     * The angle is computed from the left side of b when looking from a.
+     *
+     *   c
+     *    \                     .
+     *     \ b
+     * angle|
+     *      |
+     *      a
+     *
+     * \param a start of first line segment
+     * \param b end of first segment and start of second line segment
+     * \param c end of second line segment
+     * \return the angle in radians between 0 and 2 * pi of the corner in \p b
+ */
+static inline float getAngleLeft(const Point& a, const Point& b, const Point& c)
+{
+    const Vec2i64 ba   = (a - b).cast<int64_t>();
+    const Vec2i64 bc   = (c - b).cast<int64_t>();
+    const int64_t dott = ba.dot(bc);      // dot product
+    const int64_t det  = cross2(ba, bc); // determinant
+    if (det == 0) {
+        if ((ba.x() != 0 && (ba.x() > 0) == (bc.x() > 0)) || (ba.x() == 0 && (ba.y() > 0) == (bc.y() > 0)))
+            return 0; // pointy bit
+        else
+            return float(M_PI); // straight bit
+    }
+    const float angle = -atan2(double(det), double(dott)); // from -pi to pi
+    if (angle >= 0)
+        return angle;
+    else
+        return M_PI * 2 + angle;
+}
+
+}//namespace Slic3r::Arachne
+#endif//UTILS_LINEAR_ALG_2D_H

src/libslic3r/CMakeLists.txt

@@ -291,6 +291,45 @@ add_library(libslic3r STATIC
     SLA/Clustering.hpp
     SLA/Clustering.cpp
     SLA/ReprojectPointsOnMesh.hpp
+
+    Arachne/BeadingStrategy/BeadingStrategy.hpp
+    Arachne/BeadingStrategy/BeadingStrategy.cpp
+    Arachne/BeadingStrategy/BeadingStrategyFactory.hpp
+    Arachne/BeadingStrategy/BeadingStrategyFactory.cpp
+    Arachne/BeadingStrategy/CenterDeviationBeadingStrategy.hpp
+    Arachne/BeadingStrategy/CenterDeviationBeadingStrategy.cpp
+    Arachne/BeadingStrategy/DistributedBeadingStrategy.hpp
+    Arachne/BeadingStrategy/DistributedBeadingStrategy.cpp
+    Arachne/BeadingStrategy/LimitedBeadingStrategy.hpp
+    Arachne/BeadingStrategy/LimitedBeadingStrategy.cpp
+    Arachne/BeadingStrategy/OuterWallInsetBeadingStrategy.hpp
+    Arachne/BeadingStrategy/OuterWallInsetBeadingStrategy.cpp
+    Arachne/BeadingStrategy/RedistributeBeadingStrategy.hpp
+    Arachne/BeadingStrategy/RedistributeBeadingStrategy.cpp
+    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
+    Arachne/utils/HalfEdgeGraph.hpp
+    Arachne/utils/HalfEdgeNode.hpp
+    Arachne/utils/SparseGrid.hpp
+    Arachne/utils/SquareGrid.hpp
+    Arachne/utils/SquareGrid.cpp
+    Arachne/utils/PolygonsPointIndex.hpp
+    Arachne/utils/PolygonsSegmentIndex.hpp
+    Arachne/utils/VoronoiUtils.hpp
+    Arachne/utils/VoronoiUtils.cpp
+    Arachne/SkeletalTrapezoidation.hpp
+    Arachne/SkeletalTrapezoidation.cpp
+    Arachne/SkeletalTrapezoidationEdge.hpp
+    Arachne/SkeletalTrapezoidationGraph.hpp
+    Arachne/SkeletalTrapezoidationGraph.cpp
+    Arachne/SkeletalTrapezoidationJoint.hpp
+    Arachne/WallToolPaths.hpp
+    Arachne/WallToolPaths.cpp
 )
 
 if (SLIC3R_STATIC)

src/libslic3r/LayerRegion.cpp

@@ -100,8 +100,11 @@ void LayerRegion::make_perimeters(const SurfaceCollection &slices, SurfaceCollec
     g.ext_perimeter_flow    = this->flow(frExternalPerimeter);
     g.overhang_flow         = this->bridging_flow(frPerimeter);
     g.solid_infill_flow     = this->flow(frSolidInfill);
-    
-    g.process();
+
+    if (print_config.slicing_engine.value == SlicingEngine::Arachne)
+        g.process_arachne();
+    else
+        g.process_classic();
 }
 
 //#define EXTERNAL_SURFACES_OFFSET_PARAMETERS ClipperLib::jtMiter, 3.

src/libslic3r/Line.hpp

@@ -82,6 +82,44 @@ double distance_to(const L &line, const Vec<Dim<L>, Scalar<L>> &point)
     return std::sqrt(distance_to_squared(line, point));
 }
 
+// Returns a squared distance to the closest point on the infinite.
+// Returned nearest_point (and returned squared distance to this point) could be beyond the 'a' and 'b' ends of the segment.
+template<class L>
+double distance_to_infinite_squared(const L &line, const Vec<Dim<L>, Scalar<L>> &point, Vec<Dim<L>, Scalar<L>> *closest_point)
+{
+    const Vec<Dim<L>, double> v  = (get_b(line) - get_a(line)).template cast<double>();
+    const Vec<Dim<L>, double> va = (point - get_a(line)).template cast<double>();
+    const double              l2 = v.squaredNorm(); // avoid a sqrt
+    if (l2 == 0.) {
+        // a == b case
+        *closest_point = get_a(line);
+        return va.squaredNorm();
+    }
+    // Consider the line extending the segment, parameterized as a + t (b - a).
+    // We find projection of this point onto the line.
+    // It falls where t = [(this-a) . (b-a)] / |b-a|^2
+    const double t = va.dot(v) / l2;
+    *closest_point = (get_a(line).template cast<double>() + t * v).template cast<Scalar<L>>();
+    return (t * v - va).squaredNorm();
+}
+
+// Returns a squared distance to the closest point on the infinite.
+// Closest point (and returned squared distance to this point) could be beyond the 'a' and 'b' ends of the segment.
+template<class L>
+double distance_to_infinite_squared(const L &line, const Vec<Dim<L>, Scalar<L>> &point)
+{
+    Vec<Dim<L>, Scalar<L>> nearest_point;
+    return distance_to_infinite_squared<L>(line, point, &nearest_point);
+}
+
+// Returns a distance to the closest point on the infinite.
+// Closest point (and returned squared distance to this point) could be beyond the 'a' and 'b' ends of the segment.
+template<class L>
+double distance_to_infinite(const L &line, const Vec<Dim<L>, Scalar<L>> &point)
+{
+    return std::sqrt(distance_to_infinite_squared(line, point));
+}
+
 } // namespace line_alg
 
 class Line
@@ -102,6 +140,7 @@ public:
     double distance_to_squared(const Point &point) const { return distance_to_squared(point, this->a, this->b); }
     double distance_to_squared(const Point &point, Point *closest_point) const { return line_alg::distance_to_squared(*this, point, closest_point); }
     double distance_to(const Point &point) const { return distance_to(point, this->a, this->b); }
+    double distance_to_infinite_squared(const Point &point, Point *closest_point) const { return line_alg::distance_to_infinite_squared(*this, point, closest_point); }
     double perp_distance_to(const Point &point) const;
     bool   parallel_to(double angle) const;
     bool   parallel_to(const Line& line) const;
@@ -122,6 +161,11 @@ public:
     static inline double distance_to_squared(const Point &point, const Point &a, const Point &b) { return line_alg::distance_to_squared(Line{a, b}, Vec<2, coord_t>{point}); }
     static double distance_to(const Point &point, const Point &a, const Point &b) { return sqrt(distance_to_squared(point, a, b)); }
 
+    // Returns a distance to the closest point on the infinite.
+    // Closest point (and returned squared distance to this point) could be beyond the 'a' and 'b' ends of the segment.
+    static inline double distance_to_infinite_squared(const Point &point, const Point &a, const Point &b) { return line_alg::distance_to_infinite_squared(Line{a, b}, Vec<2, coord_t>{point}); }
+    static double distance_to_infinite(const Point &point, const Point &a, const Point &b) { return sqrt(distance_to_infinite_squared(point, a, b)); }
+
     Point a;
     Point b;
 

src/libslic3r/PerimeterGenerator.cpp

@@ -2,6 +2,7 @@
 #include "ClipperUtils.hpp"
 #include "ExtrusionEntityCollection.hpp"
 #include "ShortestPath.hpp"
+#include "Arachne/WallToolPaths.hpp"
 
 #include <cmath>
 #include <cassert>
@@ -275,7 +276,106 @@ static ExtrusionEntityCollection traverse_loops(const PerimeterGenerator &perime
     return out;
 }
 
-void PerimeterGenerator::process()
+// Thanks, Cura developers, for implementing an algorithm for generating perimeters with variable width (Arachne) that is based on the paper
+// "A framework for adaptive width control of dense contour-parallel toolpaths in fused deposition modeling"
+void PerimeterGenerator::process_arachne()
+{
+    // other perimeters
+    m_mm3_per_mm               		= this->perimeter_flow.mm3_per_mm();
+    coord_t perimeter_width         = this->perimeter_flow.scaled_width();
+
+    // external perimeters
+    m_ext_mm3_per_mm           		= this->ext_perimeter_flow.mm3_per_mm();
+    coord_t ext_perimeter_width     = this->ext_perimeter_flow.scaled_width();
+
+    // overhang perimeters
+    m_mm3_per_mm_overhang      		= this->overhang_flow.mm3_per_mm();
+
+    // solid infill
+    coord_t solid_infill_spacing    = this->solid_infill_flow.scaled_spacing();
+
+    // prepare grown lower layer slices for overhang detection
+    if (this->lower_slices != NULL && this->config->overhangs) {
+        // We consider overhang any part where the entire nozzle diameter is not supported by the
+        // lower layer, so we take lower slices and offset them by half the nozzle diameter used
+        // in the current layer
+        double nozzle_diameter = this->print_config->nozzle_diameter.get_at(this->config->perimeter_extruder-1);
+        m_lower_slices_polygons = offset(*this->lower_slices, float(scale_(+nozzle_diameter/2)));
+    }
+
+    // we need to process each island separately because we might have different
+    // extra perimeters for each one
+    for (const Surface &surface : this->slices->surfaces) {
+        // detect how many perimeters must be generated for this island
+        int        loop_number = this->config->perimeters + surface.extra_perimeters - 1; // 0-indexed loops
+        ExPolygons last        = union_ex(surface.expolygon.simplify_p(m_scaled_resolution));
+        Polygons   last_p      = to_polygons(last);
+
+        coord_t bead_width_0 = ext_perimeter_width;
+        coord_t bead_width_x = perimeter_width;
+        coord_t wall_0_inset = 0;
+
+        Arachne::WallToolPaths wallToolPaths(last_p, bead_width_0, bead_width_x, coord_t(loop_number + 1), wall_0_inset, *this->print_config);
+        wallToolPaths.generate();
+
+        std::set<size_t>      bins_with_index_zero_perimeters;
+        Arachne::BinJunctions perimeters = wallToolPaths.getBinJunctions(bins_with_index_zero_perimeters);
+
+        int start_perimeter = int(perimeters.size()) - 1;
+        int end_perimeter   = -1;
+        int direction       = -1;
+
+        if (this->config->external_perimeters_first) {
+            start_perimeter = 0;
+            end_perimeter   = int(perimeters.size());
+            direction       = 1;
+        }
+
+        for (int perimeter_idx = start_perimeter; perimeter_idx != end_perimeter; perimeter_idx += direction) {
+            if (perimeters[perimeter_idx].empty())
+                continue;
+
+            ThickPolylines thick_polylines;
+            for (const Arachne::LineJunctions &ej : perimeters[perimeter_idx])
+                thick_polylines.emplace_back(Arachne::to_thick_polyline(ej));
+            ExtrusionEntityCollection entities_coll;
+            if (bins_with_index_zero_perimeters.count(perimeter_idx) > 0) // Print using outer wall config.
+                variable_width(thick_polylines, erExternalPerimeter, this->ext_perimeter_flow, entities_coll.entities);
+            else
+                variable_width(thick_polylines, erPerimeter, this->perimeter_flow, entities_coll.entities);
+            this->loops->append(entities_coll);
+        }
+
+        ExPolygons infill_contour = union_ex(wallToolPaths.getInnerContour());
+        // create one more offset to be used as boundary for fill
+        // we offset by half the perimeter spacing (to get to the actual infill boundary)
+        // and then we offset back and forth by half the infill spacing to only consider the
+        // non-collapsing regions
+        coord_t inset =
+            (loop_number < 0) ? 0 :
+            (loop_number == 0) ?
+                                // one loop
+                ext_perimeter_width:
+                // two or more loops?
+                perimeter_width;
+
+        inset = coord_t(scale_(this->config->get_abs_value("infill_overlap", unscale<double>(inset))));
+        Polygons pp;
+        for (ExPolygon &ex : infill_contour)
+            ex.simplify_p(m_scaled_resolution, &pp);
+        // collapse too narrow infill areas
+        coord_t min_perimeter_infill_spacing = coord_t(solid_infill_spacing * (1. - INSET_OVERLAP_TOLERANCE));
+        // append infill areas to fill_surfaces
+        this->fill_surfaces->append(
+            offset_ex(offset2_ex(
+                          union_ex(pp),
+                          float(-min_perimeter_infill_spacing / 2.),
+                          float(min_perimeter_infill_spacing / 2.)), inset),
+            stInternal);
+    }
+}
+
+void PerimeterGenerator::process_classic()
 {
     // other perimeters
     m_mm3_per_mm               		= this->perimeter_flow.mm3_per_mm();

src/libslic3r/PerimeterGenerator.hpp

@@ -55,7 +55,8 @@ public:
             m_ext_mm3_per_mm(-1), m_mm3_per_mm(-1), m_mm3_per_mm_overhang(-1)
         {}
 
-    void        process();
+    void        process_classic();
+    void        process_arachne();
 
     double      ext_mm3_per_mm()        const { return m_ext_mm3_per_mm; }
     double      mm3_per_mm()            const { return m_mm3_per_mm; }

src/libslic3r/Point.hpp

@@ -161,6 +161,7 @@ public:
     Point  rotated(double angle) const { Point res(*this); res.rotate(angle); return res; }
     Point  rotated(double cos_a, double sin_a) const { Point res(*this); res.rotate(cos_a, sin_a); return res; }
     Point  rotated(double angle, const Point &center) const { Point res(*this); res.rotate(angle, center); return res; }
+    Point  rotate_90_degree_ccw() const { return Point(-this->y(), this->x()); }
     int    nearest_point_index(const Points &points) const;
     int    nearest_point_index(const PointConstPtrs &points) const;
     int    nearest_point_index(const PointPtrs &points) const;
@@ -248,6 +249,15 @@ inline bool has_duplicate_successive_points_closed(const std::vector<Point> &pts
     return has_duplicate_successive_points(pts) || (pts.size() >= 2 && pts.front() == pts.back());
 }
 
+inline bool shorter_then(const Point& p0, const coord_t len)
+{
+    if (p0.x() > len || p0.x() < -len)
+        return false;
+    if (p0.y() > len || p0.y() < -len)
+        return false;
+    return p0.cast<int64_t>().squaredNorm() <= Slic3r::sqr(int64_t(len));
+}
+
 namespace int128 {
     // Exact orientation predicate,
     // returns +1: CCW, 0: collinear, -1: CW.

src/libslic3r/Preset.cpp

@@ -448,7 +448,9 @@ static std::vector<std::string> s_Preset_print_options {
     "top_infill_extrusion_width", "support_material_extrusion_width", "infill_overlap", "infill_anchor", "infill_anchor_max", "bridge_flow_ratio", "clip_multipart_objects",
     "elefant_foot_compensation", "xy_size_compensation", "threads", "resolution", "gcode_resolution", "wipe_tower", "wipe_tower_x", "wipe_tower_y",
     "wipe_tower_width", "wipe_tower_rotation_angle", "wipe_tower_brim_width", "wipe_tower_bridging", "single_extruder_multi_material_priming", "mmu_segmented_region_max_width",
-    "wipe_tower_no_sparse_layers", "compatible_printers", "compatible_printers_condition", "inherits"
+    "wipe_tower_no_sparse_layers", "compatible_printers", "compatible_printers_condition", "inherits",
+    "slicing_engine", "beading_strategy_type", "wall_transition_length", "wall_transition_filter_distance", "wall_transition_angle",
+    "wall_distribution_count", "wall_split_middle_threshold", "wall_add_middle_threshold", "min_feature_size", "min_bead_width"
 };
 
 static std::vector<std::string> s_Preset_filament_options {

src/libslic3r/Print.cpp

@@ -222,6 +222,18 @@ bool Print::invalidate_state_by_config_options(const ConfigOptionResolver & /* n
             osteps.emplace_back(posInfill);
             osteps.emplace_back(posSupportMaterial);
             steps.emplace_back(psSkirtBrim);
+        } else if (
+               opt_key == "slicing_engine"
+            || opt_key == "beading_strategy_type"
+            || opt_key == "wall_transition_length"
+            || opt_key == "wall_transition_filter_distance"
+            || opt_key == "wall_transition_angle"
+            || opt_key == "wall_distribution_count"
+            || opt_key == "wall_split_middle_threshold"
+            || opt_key == "wall_add_middle_threshold"
+            || opt_key == "min_feature_size"
+            || opt_key == "min_bead_width") {
+            osteps.emplace_back(posSlice);
         } else {
             // for legacy, if we can't handle this option let's invalidate all steps
             //FIXME invalidate all steps of all objects as well?

src/libslic3r/PrintConfig.cpp

@@ -195,6 +195,19 @@ static const t_config_enum_values s_keys_map_ForwardCompatibilitySubstitutionRul
 };
 CONFIG_OPTION_ENUM_DEFINE_STATIC_MAPS(ForwardCompatibilitySubstitutionRule)
 
+static t_config_enum_values s_keys_map_SlicingEngine {
+    { "classic", int(SlicingEngine::Classic) },
+    { "arachne", int(SlicingEngine::Arachne) }
+};
+CONFIG_OPTION_ENUM_DEFINE_STATIC_MAPS(SlicingEngine)
+
+static t_config_enum_values s_keys_map_BeadingStrategyType {
+    { "center_deviation",   int(BeadingStrategyType::Center) },
+    { "distributed",        int(BeadingStrategyType::Distributed) },
+    { "inward_distributed", int(BeadingStrategyType::InwardDistributed) }
+};
+CONFIG_OPTION_ENUM_DEFINE_STATIC_MAPS(BeadingStrategyType)
+
 static void assign_printer_technology_to_unknown(t_optiondef_map &options, PrinterTechnology printer_technology)
 {
     for (std::pair<const t_config_option_key, ConfigOptionDef> &kvp : options)
@@ -3037,6 +3050,136 @@ void PrintConfigDef::init_fff_params()
     def->mode = comAdvanced;
     def->set_default_value(new ConfigOptionFloat(0));
 
+    def = this->add("slicing_engine", coEnum);
+    def->label = L("Slicing engine");
+    def->category = L("Advanced");
+    def->tooltip = L("Classic slicing engine produces perimeters with constant extrusion width and for"
+                      " very thing areas is used gap-fill."
+                      "Arachne produces perimeters with variable extrusion width.");
+    def->enum_keys_map = &ConfigOptionEnum<SlicingEngine>::get_enum_values();
+    def->enum_values.push_back("classic");
+    def->enum_values.push_back("arachne");
+    def->enum_labels.push_back(L("Classic"));
+    def->enum_labels.push_back(L("Arachne"));
+    def->mode = comExpert;
+    def->set_default_value(new ConfigOptionEnum<SlicingEngine>(SlicingEngine::Classic));
+
+    def = this->add("beading_strategy_type", coEnum);
+    def->label = L("Variable Line Strategy");
+    def->category = L("Advanced");
+    def->tooltip = L("Strategy to use to print the width of a part with a number of walls. This determines "
+                                 "how many walls it will use for a certain total width, and how wide each of"
+                                 " these lines are. \"Center Deviation\" will print all walls at the nominal"
+                                 " line width except the central one(s), causing big variations in the center"
+                                 " but very consistent outsides. \"Distributed\" distributes the width equally"
+                                 " over all walls. \"Inward Distributed\" is a balance between the other two, "
+                                 "distributing the changes in width over all walls but keeping the walls on the"
+                                 " outside slightly more consistent.");
+    def->enum_keys_map = &ConfigOptionEnum<BeadingStrategyType>::get_enum_values();
+    def->enum_values.push_back("center_deviation");
+    def->enum_values.push_back("distributed");
+    def->enum_values.push_back("inward_distributed");
+    def->enum_labels.push_back(L("Center Deviation"));
+    def->enum_labels.push_back(L("Distributed"));
+    def->enum_labels.push_back(L("Inward Distributed"));
+    def->mode = comExpert;
+    def->set_default_value(new ConfigOptionEnum<BeadingStrategyType>(BeadingStrategyType::InwardDistributed));
+
+    def = this->add("wall_transition_length", coFloat);
+    def->label = L("Wall Transition Length");
+    def->category = L("Advanced");
+    def->tooltip  = L("When transitioning between different numbers of walls as the part becomes"
+                       "thinner, a certain amount of space is allotted to split or join the wall lines.");
+    def->sidetext = L("mm");
+    def->mode = comExpert;
+    def->min = 0;
+    def->set_default_value(new ConfigOptionFloat(0.4));
+
+    def = this->add("wall_transition_filter_distance", coFloat);
+    def->label = L("Wall Transition Distance Filter");
+    def->category = L("Advanced");
+    def->tooltip  = L("If it would be transitioning back and forth between different numbers of walls in "
+                       "quick succession, don't transition at all. Remove transitions if they are closer "
+                       "together than this distance.");
+    def->sidetext = L("mm");
+    def->mode = comExpert;
+    def->min = 0;
+    def->set_default_value(new ConfigOptionFloat(1.4));
+
+    def = this->add("wall_transition_angle", coFloat);
+    def->label = L("Wall Transition Angle");
+    def->category = L("Advanced");
+    def->tooltip  = L("When transitioning between different numbers of walls as the part becomes thinner, "
+                       "two adjacent walls will join together at this angle. This can make the walls come "
+                       "together faster than what the Wall Transition Length indicates, filling the space "
+                       "better.");
+    def->sidetext = L("°");
+    def->mode = comExpert;
+    def->min = 1.;
+    def->max = 59.;
+    def->set_default_value(new ConfigOptionFloat(10.));
+
+    def = this->add("wall_distribution_count", coInt);
+    def->label = L("Wall Distribution Count");
+    def->category = L("Advanced");
+    def->tooltip  = L("The number of walls, counted from the center, over which the variation needs to be "
+                       "spread. Lower values mean that the outer walls don't change in width.");
+    def->mode = comExpert;
+    def->min = 1;
+    def->set_default_value(new ConfigOptionInt(1));
+
+    def = this->add("wall_split_middle_threshold", coPercent);
+    def->label = L("Split Middle Line Threshold");
+    def->category = L("Advanced");
+    def->tooltip  = L("The smallest line width, as a factor of the normal line width, above which the middle "
+                       "line (if there is one) will be split into two. Reduce this setting to use more, thinner "
+                       "lines. Increase to use fewer, wider lines. Note that this applies -as if- the entire "
+                       "shape should be filled with wall, so the middle here refers to the middle of the object "
+                       "between two outer edges of the shape, even if there actually is fill or (other) skin in "
+                       "the print instead of wall.");
+    def->sidetext = L("%");
+    def->mode = comExpert;
+    def->min = 1;
+    def->max = 99;
+    def->set_default_value(new ConfigOptionPercent(90));
+
+    def = this->add("wall_add_middle_threshold", coPercent);
+    def->label = L("Add Middle Line Threshold");
+    def->category = L("Advanced");
+    def->tooltip  = L("The smallest line width, as a factor of the normal line width, above which a middle "
+                       "line (if there wasn't one already) will be added. Reduce this setting to use more, "
+                       "thinner lines. Increase to use fewer, wider lines. Note that this applies -as if- the "
+                       "entire shape should be filled with wall, so the middle here refers to the middle of the "
+                       "object between two outer edges of the shape, even if there actually is fill or (other) "
+                       "skin in the print instead of wall.");
+    def->sidetext = L("%");
+    def->mode = comExpert;
+    def->min = 1;
+    def->max = 99;
+    def->set_default_value(new ConfigOptionPercent(80));
+
+    def = this->add("min_feature_size", coFloat);
+    def->label = L("Minimum Feature Size");
+    def->category = L("Advanced");
+    def->tooltip  = L("Minimum thickness of thin features. Model features that are thinner than this value will "
+                       "not be printed, while features thicker than the Minimum Feature Size will be widened to "
+                       "the Minimum Wall Line Width.");
+    def->sidetext = L("mm");
+    def->mode = comExpert;
+    def->min = 0;
+    def->set_default_value(new ConfigOptionFloat(0.1));
+
+    def = this->add("min_bead_width", coFloat);
+    def->label = L("Minimum Wall Line Width");
+    def->category = L("Advanced");
+    def->tooltip  = L("Width of the wall that will replace thin features (according to the Minimum Feature Size) "
+                       "of the model. If the Minimum Wall Line Width is thinner than the thickness of the feature,"
+                       " the wall will become as thick as the feature itself.");
+    def->sidetext = L("mm");
+    def->mode = comExpert;
+    def->min = 0;
+    def->set_default_value(new ConfigOptionFloat(0.2));
+
     // Declare retract values for filament profile, overriding the printer's extruder profile.
     for (const char *opt_key : {
         // floats
@@ -3968,6 +4111,13 @@ void DynamicPrintConfig::normalize_fdm()
     if (auto *opt_gcode_resolution = this->opt<ConfigOptionFloat>("gcode_resolution", false); opt_gcode_resolution)
         // Resolution will be above 1um.
         opt_gcode_resolution->value = std::max(opt_gcode_resolution->value, 0.001);
+
+    if (auto *opt_min_bead_width = this->opt<ConfigOptionFloat>("min_bead_width", false); opt_min_bead_width)
+        opt_min_bead_width->value = std::max(opt_min_bead_width->value, 0.001);
+    if (auto *opt_wall_transition_length = this->opt<ConfigOptionFloat>("wall_transition_length", false); opt_wall_transition_length)
+        opt_wall_transition_length->value = std::max(opt_wall_transition_length->value, 0.001);
+    if (auto *opt_wall_transition_filter_distance = this->opt<ConfigOptionFloat>("wall_transition_filter_distance", false); opt_wall_transition_filter_distance)
+        opt_wall_transition_filter_distance->value = std::max(opt_wall_transition_filter_distance->value, 0.001);
 }
 
 void  handle_legacy_sla(DynamicPrintConfig &config)

src/libslic3r/PrintConfig.hpp

@@ -127,6 +127,24 @@ enum DraftShield {
     dsDisabled, dsLimited, dsEnabled
 };
 
+enum class SlicingEngine
+{
+    // Classic perimeter generator using Clipper offsets with constant extrusion width.
+    Classic,
+    // Perimeter generator with variable extrusion width based on the paper
+    // "A framework for adaptive width control of dense contour-parallel toolpaths in fused deposition modeling" ported from Cura.
+    Arachne
+};
+
+enum class BeadingStrategyType
+{
+    Center,
+    Distributed,
+    InwardDistributed,
+    None,
+    Count
+};
+
 #define CONFIG_OPTION_ENUM_DECLARE_STATIC_MAPS(NAME) \
     template<> const t_config_enum_names& ConfigOptionEnum<NAME>::get_enum_names(); \
     template<> const t_config_enum_values& ConfigOptionEnum<NAME>::get_enum_values();
@@ -149,6 +167,8 @@ CONFIG_OPTION_ENUM_DECLARE_STATIC_MAPS(SLAPillarConnectionMode)
 CONFIG_OPTION_ENUM_DECLARE_STATIC_MAPS(BrimType)
 CONFIG_OPTION_ENUM_DECLARE_STATIC_MAPS(DraftShield)
 CONFIG_OPTION_ENUM_DECLARE_STATIC_MAPS(ForwardCompatibilitySubstitutionRule)
+CONFIG_OPTION_ENUM_DECLARE_STATIC_MAPS(SlicingEngine)
+CONFIG_OPTION_ENUM_DECLARE_STATIC_MAPS(BeadingStrategyType)
 
 #undef CONFIG_OPTION_ENUM_DECLARE_STATIC_MAPS
 
@@ -748,6 +768,16 @@ PRINT_CONFIG_CLASS_DERIVED_DEFINE(
     ((ConfigOptionFloat,              skirt_distance))
     ((ConfigOptionInt,                skirt_height))
     ((ConfigOptionInt,                skirts))
+    ((ConfigOptionEnum<SlicingEngine>, slicing_engine))
+    ((ConfigOptionEnum<BeadingStrategyType>, beading_strategy_type))
+    ((ConfigOptionFloat,              wall_transition_length))
+    ((ConfigOptionFloat,              wall_transition_filter_distance))
+    ((ConfigOptionFloat,              wall_transition_angle))
+    ((ConfigOptionInt,                wall_distribution_count))
+    ((ConfigOptionPercent,            wall_split_middle_threshold))
+    ((ConfigOptionPercent,            wall_add_middle_threshold))
+    ((ConfigOptionFloat,              min_feature_size))
+    ((ConfigOptionFloat,              min_bead_width))
     ((ConfigOptionInts,               slowdown_below_layer_time))
     ((ConfigOptionBool,               spiral_vase))
     ((ConfigOptionInt,                standby_temperature_delta))

src/slic3r/GUI/Tab.cpp

@@ -1670,6 +1670,18 @@ void TabPrint::build()
         optgroup = page->new_optgroup(L("Other"));
         optgroup->append_single_option_line("clip_multipart_objects");
 
+        optgroup = page->new_optgroup(L("Experimental"));
+        optgroup->append_single_option_line("slicing_engine");
+        optgroup->append_single_option_line("beading_strategy_type");
+        optgroup->append_single_option_line("wall_transition_length");
+        optgroup->append_single_option_line("wall_transition_filter_distance");
+        optgroup->append_single_option_line("wall_transition_angle");
+        optgroup->append_single_option_line("wall_distribution_count");
+        optgroup->append_single_option_line("wall_split_middle_threshold");
+        optgroup->append_single_option_line("wall_add_middle_threshold");
+        optgroup->append_single_option_line("min_feature_size");
+        optgroup->append_single_option_line("min_bead_width");
+
     page = add_options_page(L("Output options"), "output+page_white");
         optgroup = page->new_optgroup(L("Sequential printing"));
         optgroup->append_single_option_line("complete_objects", "sequential-printing_124589");