QIDISlicer1.0.0

src/libslic3r/Arachne/utils/ExtrusionJunction.cpp

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

src/libslic3r/Arachne/utils/ExtrusionJunction.hpp

@@ -0,0 +1,59 @@
+//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;
+
+    ExtrusionJunction(const Point p, const coord_t w, const coord_t perimeter_index);
+
+    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>; //<! The junctions along a line without further information. See \ref ExtrusionLine for a more extensive class.
+
+}
+#endif // UTILS_EXTRUSION_JUNCTION_H

src/libslic3r/Arachne/utils/ExtrusionLine.cpp

@@ -0,0 +1,270 @@
+//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"
+#include "../../PerimeterGenerator.hpp"
+
+namespace Slic3r::Arachne
+{
+
+ExtrusionLine::ExtrusionLine(const size_t inset_idx, const bool is_odd) : inset_idx(inset_idx), is_odd(is_odd), is_closed(false) {}
+
+int64_t ExtrusionLine::getLength() const
+{
+    if (junctions.empty())
+        return 0;
+
+    int64_t           len  = 0;
+    ExtrusionJunction prev = junctions.front();
+    for (const ExtrusionJunction &next : junctions) {
+        len += (next.p - prev.p).cast<int64_t>().norm();
+        prev = next;
+    }
+    if (is_closed)
+        len += (front().p - back().p).cast<int64_t>().norm();
+
+    return len;
+}
+
+coord_t ExtrusionLine::getMinimalWidth() const
+{
+    return std::min_element(junctions.cbegin(), junctions.cend(),
+                            [](const ExtrusionJunction& l, const ExtrusionJunction& r)
+                            {
+                                return l.w < r.w;
+                            })->w;
+}
+
+void ExtrusionLine::simplify(const int64_t smallest_line_segment_squared, const int64_t allowed_error_distance_squared, const int64_t maximum_extrusion_area_deviation)
+{
+    const size_t min_path_size = is_closed ? 3 : 2;
+    if (junctions.size() <= min_path_size)
+        return;
+
+    // TODO: allow for the first point to be removed in case of simplifying closed Extrusionlines.
+
+    /* 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 auto    height_2 = int64_t(double(area_removed_so_far) * double(area_removed_so_far) / double(base_length_2));
+        const int64_t extrusion_area_error = calculateExtrusionAreaDeviationError(previous, current, next);
+        if ((height_2 <= scaled<coord_t>(0.001) //Almost exactly colinear (barring rounding errors).
+             && Line::distance_to_infinite(current.p, previous.p, next.p) <= scaled<double>(0.001)) // 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)
+        {
+            // Remove the current junction (vertex).
+            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 > 4 * 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);
+                    // 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) {
+    /*
+     * A             B                          C              A                                        C
+     * ---------------                                         **************
+     * |             |                                         ------------------------------------------
+     * |             |--------------------------|  B removed   |            |***************************|
+     * |             |                          |  --------->  |            |                           |
+     * |             |--------------------------|              |            |***************************|
+     * |             |                                         ------------------------------------------
+     * ---------------             ^                           **************
+     *       ^                B.w + C.w / 2                                       ^
+     *  A.w + B.w / 2                                               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.p - A.p).cast<int64_t>().norm();
+    const int64_t bc_length = (C.p - B.p).cast<int64_t>().norm();
+    if (const coord_t width_diff = std::max(std::abs(B.w - A.w), std::abs(C.w - B.w)); width_diff > 1) {
+        // Adjust the width only if there is a difference, or else the rounding errors may produce the wrong
+        // weighted average value.
+        const int64_t ab_weight              = (A.w + B.w) / 2;
+        const int64_t bc_weight              = (B.w + C.w) / 2;
+        const int64_t weighted_average_width = (ab_length * ab_weight + bc_length * bc_weight) / (ab_length + bc_length);
+        const int64_t ac_length              = (C.p - A.p).cast<int64_t>().norm();
+        return std::abs((ab_weight * ab_length + bc_weight * bc_length) - (weighted_average_width * ac_length));
+    } else {
+        // If the width difference is very small, then select the width of the segment that is longer
+        return ab_length > bc_length ? int64_t(width_diff) * bc_length : int64_t(width_diff) * ab_length;
+    }
+}
+
+bool ExtrusionLine::is_contour() const
+{
+    if (!this->is_closed)
+        return false;
+
+    Polygon poly;
+    poly.points.reserve(this->junctions.size());
+    for (const ExtrusionJunction &junction : this->junctions)
+        poly.points.emplace_back(junction.p);
+
+    // Arachne produces contour with clockwise orientation and holes with counterclockwise orientation.
+    return poly.is_clockwise();
+}
+
+double ExtrusionLine::area() const
+{
+    assert(this->is_closed);
+    double a = 0.;
+    if (this->junctions.size() >= 3) {
+        Vec2d p1 = this->junctions.back().p.cast<double>();
+        for (const ExtrusionJunction &junction : this->junctions) {
+            Vec2d p2 = junction.p.cast<double>();
+            a += cross2(p1, p2);
+            p1 = p2;
+        }
+    }
+    return 0.5 * a;
+}
+
+} // namespace Slic3r::Arachne
+
+namespace Slic3r {
+void extrusion_paths_append(ExtrusionPaths &dst, const ClipperLib_Z::Paths &extrusion_paths, const ExtrusionRole role, const Flow &flow)
+{
+    for (const ClipperLib_Z::Path &extrusion_path : extrusion_paths) {
+        ThickPolyline thick_polyline = Arachne::to_thick_polyline(extrusion_path);
+        Slic3r::append(dst, PerimeterGenerator::thick_polyline_to_multi_path(thick_polyline, role, flow, scaled<float>(0.05), float(SCALED_EPSILON)).paths);
+    }
+}
+
+void extrusion_paths_append(ExtrusionPaths &dst, const Arachne::ExtrusionLine &extrusion, const ExtrusionRole role, const Flow &flow)
+{
+    ThickPolyline thick_polyline = Arachne::to_thick_polyline(extrusion);
+    Slic3r::append(dst, PerimeterGenerator::thick_polyline_to_multi_path(thick_polyline, role, flow, scaled<float>(0.05), float(SCALED_EPSILON)).paths);
+}
+} // namespace Slic3r

src/libslic3r/Arachne/utils/ExtrusionLine.hpp

@@ -0,0 +1,306 @@
+//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"
+#include "../../BoundingBox.hpp"
+#include "../../ExtrusionEntity.hpp"
+#include "../../Flow.hpp"
+#include "../../../clipper/clipper_z.hpp"
+
+namespace Slic3r {
+struct 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;
+
+    /*!
+     * Whether this is a closed polygonal path
+     */
+    bool is_closed;
+
+    /*!
+     * Gets the number of vertices in this polygon.
+     * \return The number of vertices in this polygon.
+     */
+    size_t size() const { return junctions.size(); }
+
+    /*!
+     * Whether there are no junctions.
+     */
+    bool empty() const { return junctions.empty(); }
+
+    /*!
+     * 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);
+    ExtrusionLine() : inset_idx(-1), is_odd(true), is_closed(false) {}
+    ExtrusionLine(const ExtrusionLine &other) : inset_idx(other.inset_idx), is_odd(other.is_odd), is_closed(other.is_closed), junctions(other.junctions) {}
+
+    ExtrusionLine &operator=(ExtrusionLine &&other)
+    {
+        junctions = std::move(other.junctions);
+        inset_idx = other.inset_idx;
+        is_odd    = other.is_odd;
+        is_closed = other.is_closed;
+        return *this;
+    }
+
+    ExtrusionLine &operator=(const ExtrusionLine &other)
+    {
+        junctions = other.junctions;
+        inset_idx = other.inset_idx;
+        is_odd    = other.is_odd;
+        is_closed = other.is_closed;
+        return *this;
+    }
+
+    std::vector<ExtrusionJunction>::const_iterator begin() const { return junctions.begin(); }
+    std::vector<ExtrusionJunction>::const_iterator end() const { return junctions.end(); }
+    std::vector<ExtrusionJunction>::const_reverse_iterator rbegin() const { return junctions.rbegin(); }
+    std::vector<ExtrusionJunction>::const_reverse_iterator rend() const { return junctions.rend(); }
+    std::vector<ExtrusionJunction>::const_reference front() const { return junctions.front(); }
+    std::vector<ExtrusionJunction>::const_reference back() const { return junctions.back(); }
+    const ExtrusionJunction &operator[](unsigned int index) const { return junctions[index]; }
+    ExtrusionJunction &operator[](unsigned int index) { return junctions[index]; }
+    std::vector<ExtrusionJunction>::iterator begin() { return junctions.begin(); }
+    std::vector<ExtrusionJunction>::iterator end() { return junctions.end(); }
+    std::vector<ExtrusionJunction>::reference front() { return junctions.front(); }
+    std::vector<ExtrusionJunction>::reference back() { return junctions.back(); }
+
+    template<typename... Args> void emplace_back(Args &&...args) { junctions.emplace_back(args...); }
+    void remove(unsigned int index) { junctions.erase(junctions.begin() + index); }
+    void insert(size_t index, const ExtrusionJunction &p) { junctions.insert(junctions.begin() + index, p); }
+
+    template<class iterator>
+    std::vector<ExtrusionJunction>::iterator insert(std::vector<ExtrusionJunction>::const_iterator pos, iterator first, iterator last)
+    {
+        return junctions.insert(pos, first, last);
+    }
+
+    void clear() { junctions.clear(); }
+    void reverse() { std::reverse(junctions.begin(), junctions.end()); }
+
+    /*!
+     * Sum the total length of this path.
+     */
+    int64_t getLength() const;
+    int64_t polylineLength() const { return getLength(); }
+
+    /*!
+     * Put all junction locations into a polygon object.
+     *
+     * When this path is not closed the returned Polygon should be handled as a polyline, rather than a polygon.
+     */
+    Polygon toPolygon() const
+    {
+        Polygon ret;
+        for (const ExtrusionJunction &j : junctions)
+            ret.points.emplace_back(j.p);
+
+        return ret;
+    }
+
+    /*!
+     * Get the minimal width of this path
+     */
+    coord_t getMinimalWidth() const;
+
+    /*!
+     * Removes vertices of the ExtrusionLines to make sure that they are not too high
+     * resolution.
+     *
+     * 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
+     * */
+    static int64_t calculateExtrusionAreaDeviationError(ExtrusionJunction A, ExtrusionJunction B, ExtrusionJunction C);
+
+    bool is_contour() const;
+
+    double area() const;
+};
+
+static inline Slic3r::ThickPolyline to_thick_polyline(const Arachne::ExtrusionLine &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 Slic3r::ThickPolyline to_thick_polyline(const ClipperLib_Z::Path &path)
+{
+    assert(path.size() >= 2);
+    Slic3r::ThickPolyline out;
+    out.points.emplace_back(path.front().x(), path.front().y());
+    out.width.emplace_back(path.front().z());
+    out.points.emplace_back(path[1].x(), path[1].y());
+    out.width.emplace_back(path[1].z());
+
+    auto it_prev = path.begin() + 1;
+    for (auto it = path.begin() + 2; it != path.end(); ++it) {
+        out.points.emplace_back(it->x(), it->y());
+        out.width.emplace_back(it_prev->z());
+        out.width.emplace_back(it->z());
+        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;
+}
+
+#if 0
+static BoundingBox get_extents(const ExtrusionLine &extrusion_line)
+{
+    BoundingBox bbox;
+    for (const ExtrusionJunction &junction : extrusion_line.junctions)
+        bbox.merge(junction.p);
+    return bbox;
+}
+
+static BoundingBox get_extents(const std::vector<ExtrusionLine> &extrusion_lines)
+{
+    BoundingBox bbox;
+    for (const ExtrusionLine &extrusion_line : extrusion_lines)
+        bbox.merge(get_extents(extrusion_line));
+    return bbox;
+}
+
+static BoundingBox get_extents(const std::vector<const ExtrusionLine *> &extrusion_lines)
+{
+    BoundingBox bbox;
+    for (const ExtrusionLine *extrusion_line : extrusion_lines) {
+        assert(extrusion_line != nullptr);
+        bbox.merge(get_extents(*extrusion_line));
+    }
+    return bbox;
+}
+
+static Points to_points(const ExtrusionLine &extrusion_line)
+{
+    Points points;
+    points.reserve(extrusion_line.junctions.size());
+    for (const ExtrusionJunction &junction : extrusion_line.junctions)
+        points.emplace_back(junction.p);
+    return points;
+}
+
+static std::vector<Points> to_points(const std::vector<const ExtrusionLine *> &extrusion_lines)
+{
+    std::vector<Points> points;
+    for (const ExtrusionLine *extrusion_line : extrusion_lines) {
+        assert(extrusion_line != nullptr);
+        points.emplace_back(to_points(*extrusion_line));
+    }
+    return points;
+}
+#endif
+
+using VariableWidthLines = std::vector<ExtrusionLine>; //<! The ExtrusionLines generated by libArachne
+
+} // namespace Slic3r::Arachne
+
+namespace Slic3r {
+
+void extrusion_paths_append(ExtrusionPaths &dst, const ClipperLib_Z::Paths &extrusion_paths, const ExtrusionRole role, const Flow &flow);
+void extrusion_paths_append(ExtrusionPaths &dst, const Arachne::ExtrusionLine &extrusion, const ExtrusionRole role, const Flow &flow);
+
+} // namespace Slic3r
+
+#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,31 @@
+//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;
+    using Edges = std::list<edge_t>;
+    using Nodes = std::list<node_t>;
+    Edges edges;
+    Nodes 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,180 @@
+//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
+{
+
+// Identity function, used to be able to make templated algorithms where the input is sometimes points, sometimes things that contain or can be converted to points.
+inline const Point &make_point(const Point &p) { return p; }
+
+/*!
+ * A class for iterating over the points in one of the polygons in a \ref Polygons object
+ */
+template<typename Paths>
+class PathsPointIndex
+{
+public:
+    /*!
+     * The polygons into which this index is indexing.
+     */
+    const Paths* 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.
+     */
+    PathsPointIndex() : 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.
+     */
+    PathsPointIndex(const Paths *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.
+     */
+    PathsPointIndex(const PathsPointIndex& original) = default;
+
+    Point p() const
+    {
+        if (!polygons)
+            return {0, 0};
+
+        return make_point((*polygons)[poly_idx][point_idx]);
+    }
+
+    /*!
+     * \brief Returns whether this point is initialised.
+     */
+    bool initialized() const { return polygons; }
+
+    /*!
+     * Get the polygon to which this PolygonsPointIndex refers
+     */
+    const Polygon &getPolygon() const { return (*polygons)[poly_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 PathsPointIndex &other) const
+    {
+        return polygons == other.polygons && poly_idx == other.poly_idx && point_idx == other.point_idx;
+    }
+    bool operator!=(const PathsPointIndex &other) const
+    {
+        return !(*this == other);
+    }
+    bool operator<(const PathsPointIndex &other) const
+    {
+        return this->p() < other.p();
+    }
+    PathsPointIndex &operator=(const PathsPointIndex &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)
+    PathsPointIndex &operator++()
+    {
+        point_idx = (point_idx + 1) % (*polygons)[poly_idx].size();
+        return *this;
+    }
+    //! move the iterator backward (and wrap around at the beginning)
+    PathsPointIndex &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)
+    PathsPointIndex next() const
+    {
+        PathsPointIndex ret(*this);
+        ++ret;
+        return ret;
+    }
+    //! move the iterator backward (and wrap around at the beginning)
+    PathsPointIndex prev() const
+    {
+        PathsPointIndex ret(*this);
+        --ret;
+        return ret;
+    }
+};
+
+using PolygonsPointIndex = PathsPointIndex<Polygons>;
+
+/*!
+ * 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];
+        Point          start          = poly[val.point_idx];
+        unsigned int   next_point_idx = (val.point_idx + 1) % poly.size();
+        Point          end            = poly[next_point_idx];
+        return std::pair<Point, Point>(start, end);
+    }
+};
+
+/*!
+ * Locator of a \ref PolygonsPointIndex
+ */
+template<typename Paths>
+struct PathsPointIndexLocator
+{
+    Point operator()(const PathsPointIndex<Paths>& val) const
+    {
+        return make_point(val.p());
+    }
+};
+
+using PolygonsPointIndexLocator = PathsPointIndexLocator<Polygons>;
+
+}//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/PolylineStitcher.cpp

@@ -0,0 +1,42 @@
+//Copyright (c) 2022 Ultimaker B.V.
+//CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#include "PolylineStitcher.hpp"
+#include "ExtrusionLine.hpp"
+
+namespace Slic3r::Arachne {
+
+template<> bool PolylineStitcher<VariableWidthLines, ExtrusionLine, ExtrusionJunction>::canReverse(const PathsPointIndex<VariableWidthLines> &ppi)
+{
+    if ((*ppi.polygons)[ppi.poly_idx].is_odd)
+        return true;
+    else
+        return false;
+}
+
+template<> bool PolylineStitcher<Polygons, Polygon, Point>::canReverse(const PathsPointIndex<Polygons> &)
+{
+    return true;
+}
+
+template<> bool PolylineStitcher<VariableWidthLines, ExtrusionLine, ExtrusionJunction>::canConnect(const ExtrusionLine &a, const ExtrusionLine &b)
+{
+    return a.is_odd == b.is_odd;
+}
+
+template<> bool PolylineStitcher<Polygons, Polygon, Point>::canConnect(const Polygon &, const Polygon &)
+{
+    return true;
+}
+
+template<> bool PolylineStitcher<VariableWidthLines, ExtrusionLine, ExtrusionJunction>::isOdd(const ExtrusionLine &line)
+{
+    return line.is_odd;
+}
+
+template<> bool PolylineStitcher<Polygons, Polygon, Point>::isOdd(const Polygon &)
+{
+    return false;
+}
+
+} // namespace Slic3r::Arachne

src/libslic3r/Arachne/utils/PolylineStitcher.hpp

@@ -0,0 +1,233 @@
+//Copyright (c) 2022 Ultimaker B.V.
+//CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#ifndef UTILS_POLYLINE_STITCHER_H
+#define UTILS_POLYLINE_STITCHER_H
+
+#include "SparsePointGrid.hpp"
+#include "PolygonsPointIndex.hpp"
+#include "../../Polygon.hpp"
+#include <cassert>
+
+namespace Slic3r::Arachne
+{
+
+/*!
+ * Class for stitching polylines into longer polylines or into polygons
+ */
+template<typename Paths, typename Path, typename Junction>
+class PolylineStitcher
+{
+public:
+    /*!
+     * Stitch together the separate \p lines into \p result_lines and if they
+     * can be closed into \p result_polygons.
+     *
+     * Only introduce new segments shorter than \p max_stitch_distance, and
+     * larger than \p snap_distance but always try to take the shortest
+     * connection possible.
+     *
+     * Only stitch polylines into closed polygons if they are larger than 3 *
+     * \p max_stitch_distance, in order to prevent small segments to
+     * accidentally get closed into a polygon.
+     *
+     * \warning Tiny polylines (smaller than 3 * max_stitch_distance) will not
+     * be closed into polygons.
+     *
+     * \note Resulting polylines and polygons are added onto the existing
+     * containers, so you can directly output onto a polygons container with
+     * existing polygons in it. However, you shouldn't call this function with
+     * the same parameter in \p lines as \p result_lines, because that would
+     * duplicate (some of) the polylines.
+     * \param lines The lines to stitch together.
+     * \param result_lines[out] The stitched parts that are not closed polygons
+     * will be stored in here.
+     * \param result_polygons[out] The stitched parts that were closed as
+     * polygons will be stored in here.
+     * \param max_stitch_distance The maximum distance that will be bridged to
+     * connect two lines.
+     * \param snap_distance Points closer than this distance are considered to
+     * be the same point.
+     */
+    static void stitch(const Paths& lines, Paths& result_lines, Paths& result_polygons, coord_t max_stitch_distance = scaled<coord_t>(0.1), coord_t snap_distance = scaled<coord_t>(0.01))
+    {
+        if (lines.empty())
+            return;
+
+        SparsePointGrid<PathsPointIndex<Paths>, PathsPointIndexLocator<Paths>> grid(max_stitch_distance, lines.size() * 2);
+
+        // populate grid
+        for (size_t line_idx = 0; line_idx < lines.size(); line_idx++)
+        {
+            const auto line = lines[line_idx];
+            grid.insert(PathsPointIndex<Paths>(&lines, line_idx, 0));
+            grid.insert(PathsPointIndex<Paths>(&lines, line_idx, line.size() - 1));
+        }
+
+        std::vector<bool> processed(lines.size(), false);
+
+        for (size_t line_idx = 0; line_idx < lines.size(); line_idx++)
+        {
+            if (processed[line_idx])
+            {
+                continue;
+            }
+            processed[line_idx] = true;
+            const auto line = lines[line_idx];
+            bool should_close = isOdd(line);
+
+            Path chain = line;
+            bool closest_is_closing_polygon = false;
+            for (bool go_in_reverse_direction : { false, true }) // first go in the unreversed direction, to try to prevent the chain.reverse() operation.
+            { // NOTE: Implementation only works for this order; we currently only re-reverse the chain when it's closed.
+                if (go_in_reverse_direction)
+                { // try extending chain in the other direction
+                    chain.reverse();
+                }
+                int64_t chain_length = chain.polylineLength();
+
+                while (true)
+                {
+                    Point from = make_point(chain.back());
+
+                    PathsPointIndex<Paths> closest;
+                    coord_t closest_distance = std::numeric_limits<coord_t>::max();
+                    grid.processNearby(from, max_stitch_distance,
+                                       std::function<bool (const PathsPointIndex<Paths>&)> (
+                                           [from, &chain, &closest, &closest_is_closing_polygon, &closest_distance, &processed, &chain_length, go_in_reverse_direction, max_stitch_distance, snap_distance, should_close]
+                                           (const PathsPointIndex<Paths>& nearby)->bool
+                                           {
+                                               bool is_closing_segment = false;
+                                               coord_t dist = (nearby.p().template cast<int64_t>() - from.template cast<int64_t>()).norm();
+                                               if (dist > max_stitch_distance)
+                                               {
+                                                   return true; // keep looking
+                                               }
+                                               if ((nearby.p().template cast<int64_t>() - make_point(chain.front()).template cast<int64_t>()).squaredNorm() < snap_distance * snap_distance)
+                                               {
+                                                   if (chain_length + dist < 3 * max_stitch_distance // prevent closing of small poly, cause it might be able to continue making a larger polyline
+                                                       || chain.size() <= 2) // don't make 2 vert polygons
+                                                   {
+                                                       return true; // look for a better next line
+                                                   }
+                                                   is_closing_segment = true;
+                                                   if (!should_close)
+                                                   {
+                                                       dist += scaled<coord_t>(0.01); // prefer continuing polyline over closing a polygon; avoids closed zigzags from being printed separately
+                                                       // continue to see if closing segment is also the closest
+                                                       // there might be a segment smaller than [max_stitch_distance] which closes the polygon better
+                                                   }
+                                                   else
+                                                   {
+                                                       dist -= scaled<coord_t>(0.01); //Prefer closing the polygon if it's 100% even lines. Used to create closed contours.
+                                                       //Continue to see if closing segment is also the closest.
+                                                   }
+                                               }
+                                               else if (processed[nearby.poly_idx])
+                                               { // it was already moved to output
+                                                   return true; // keep looking for a connection
+                                               }
+                                               bool nearby_would_be_reversed = nearby.point_idx != 0;
+                                               nearby_would_be_reversed = nearby_would_be_reversed != go_in_reverse_direction; // flip nearby_would_be_reversed when searching in the reverse direction
+                                               if (!canReverse(nearby) && nearby_would_be_reversed)
+                                               { // connecting the segment would reverse the polygon direction
+                                                   return true; // keep looking for a connection
+                                               }
+                                               if (!canConnect(chain, (*nearby.polygons)[nearby.poly_idx]))
+                                               {
+                                                   return true; // keep looking for a connection
+                                               }
+                                               if (dist < closest_distance)
+                                               {
+                                                   closest_distance = dist;
+                                                   closest = nearby;
+                                                   closest_is_closing_polygon = is_closing_segment;
+                                               }
+                                               if (dist < snap_distance)
+                                               { // we have found a good enough next line
+                                                   return false; // stop looking for alternatives
+                                               }
+                                               return true; // keep processing elements
+                                           })
+                    );
+
+                    if (!closest.initialized()          // we couldn't find any next line
+                        || closest_is_closing_polygon   // we closed the polygon
+                    )
+                    {
+                        break;
+                    }
+
+                    coord_t segment_dist = (make_point(chain.back()).template cast<int64_t>() - closest.p().template cast<int64_t>()).norm();
+                    assert(segment_dist <= max_stitch_distance + scaled<coord_t>(0.01));
+                    const size_t old_size = chain.size();
+                    if (closest.point_idx == 0)
+                    {
+                        auto start_pos = (*closest.polygons)[closest.poly_idx].begin();
+                        if (segment_dist < snap_distance)
+                        {
+                            ++start_pos;
+                        }
+                        chain.insert(chain.end(), start_pos, (*closest.polygons)[closest.poly_idx].end());
+                    }
+                    else
+                    {
+                        auto start_pos = (*closest.polygons)[closest.poly_idx].rbegin();
+                        if (segment_dist < snap_distance)
+                        {
+                            ++start_pos;
+                        }
+                        chain.insert(chain.end(), start_pos, (*closest.polygons)[closest.poly_idx].rend());
+                    }
+                    for(size_t i = old_size; i < chain.size(); ++i) //Update chain length.
+                    {
+                        chain_length += (make_point(chain[i]).template cast<int64_t>() - make_point(chain[i - 1]).template cast<int64_t>()).norm();
+                    }
+                    should_close = should_close & !isOdd((*closest.polygons)[closest.poly_idx]); //If we connect an even to an odd line, we should no longer try to close it.
+                    assert( ! processed[closest.poly_idx]);
+                    processed[closest.poly_idx] = true;
+                }
+                if (closest_is_closing_polygon)
+                {
+                    if (go_in_reverse_direction)
+                    { // re-reverse chain to retain original direction
+                        // NOTE: not sure if this code could ever be reached, since if a polygon can be closed that should be already possible in the forward direction
+                        chain.reverse();
+                    }
+
+                    break; // don't consider reverse direction
+                }
+            }
+            if (closest_is_closing_polygon)
+            {
+                result_polygons.emplace_back(chain);
+            }
+            else
+            {
+                PathsPointIndex<Paths> ppi_here(&lines, line_idx, 0);
+                if ( ! canReverse(ppi_here))
+                { // Since closest_is_closing_polygon is false we went through the second iterations of the for-loop, where go_in_reverse_direction is true
+                    // the polyline isn't allowed to be reversed, so we re-reverse it.
+                    chain.reverse();
+                }
+                result_lines.emplace_back(chain);
+            }
+        }
+    }
+
+    /*!
+     * Whether a polyline is allowed to be reversed. (Not true for wall polylines which are not odd)
+     */
+    static bool canReverse(const PathsPointIndex<Paths> &polyline);
+
+    /*!
+     * Whether two paths are allowed to be connected.
+     * (Not true for an odd and an even wall.)
+     */
+    static bool canConnect(const Path &a, const Path &b);
+
+    static bool isOdd(const Path &line);
+};
+
+} // namespace Slic3r::Arachne
+#endif // UTILS_POLYLINE_STITCHER_H

src/libslic3r/Arachne/utils/SparseGrid.hpp

@@ -0,0 +1,132 @@
+//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 <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,76 @@
+//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 <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/SparsePointGrid.hpp

@@ -0,0 +1,89 @@
+// Copyright (c) 2016 Scott Lenser
+// Copyright (c) 2020 Ultimaker B.V.
+// CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#ifndef UTILS_SPARSE_POINT_GRID_H
+#define UTILS_SPARSE_POINT_GRID_H
+
+#include <cassert>
+#include <vector>
+
+#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 location from ElemT.  Locator
+ *    must have: Point operator()(const ElemT &elem) const
+ *    which returns the location associated with val.
+ */
+template<class ElemT, class Locator> class SparsePointGrid : 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.
+     */
+    SparsePointGrid(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);
+
+    /*!
+     * Get just any element that's within a certain radius of a point.
+     *
+     * Rather than giving a vector of nearby elements, this function just gives
+     * a single element, any element, in no particular order.
+     * \param query_pt The point to query for an object nearby.
+     * \param radius The radius of what is considered "nearby".
+     */
+    const ElemT *getAnyNearby(const Point &query_pt, coord_t radius);
+
+protected:
+    using GridPoint = typename SparseGrid<ElemT>::GridPoint;
+
+    /*! \brief Accessor for getting locations from elements. */
+    Locator m_locator;
+};
+
+template<class ElemT, class Locator>
+SparsePointGrid<ElemT, Locator>::SparsePointGrid(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 SparsePointGrid<ElemT, Locator>::insert(const Elem &elem)
+{
+    Point     loc      = m_locator(elem);
+    GridPoint grid_loc = SparseGrid<ElemT>::toGridPoint(loc.template cast<int64_t>());
+
+    SparseGrid<ElemT>::m_grid.emplace(grid_loc, elem);
+}
+
+template<class ElemT, class Locator>
+const ElemT *SparsePointGrid<ElemT, Locator>::getAnyNearby(const Point &query_pt, coord_t radius)
+{
+    const ElemT                              *ret          = nullptr;
+    const std::function<bool(const ElemT &)> &process_func = [&ret, query_pt, radius, this](const ElemT &maybe_nearby) {
+        if (shorter_then(m_locator(maybe_nearby) - query_pt, radius)) {
+            ret = &maybe_nearby;
+            return false;
+        }
+        return true;
+    };
+    SparseGrid<ElemT>::processNearby(query_pt, radius, process_func);
+
+    return ret;
+}
+
+} // namespace Slic3r::Arachne
+
+#endif // UTILS_SPARSE_POINT_GRID_H

src/libslic3r/Arachne/utils/SquareGrid.cpp

@@ -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,109 @@
+//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 <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,251 @@
+//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(std::isfinite(x) && std::isfinite(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 {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]));
+    }
+};
+Points VoronoiUtils::discretizeParabola(const Point& p, const Segment& segment, Point s, Point e, coord_t approximate_step_size, float transitioning_angle)
+{
+    Points 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(perp(ppxx));
+
+    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,47 @@
+//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 Points discretizeParabola(const Point &source_point, const Segment &source_segment, Point start, Point end, coord_t approximate_step_size, float transitioning_angle);
+
+    static inline bool is_finite(const VoronoiUtils::vd_t::vertex_type &vertex)
+    {
+        return std::isfinite(vertex.x()) && std::isfinite(vertex.y());
+    }
+};
+
+} // namespace Slic3r::Arachne
+
+#endif // UTILS_VORONOI_UTILS_H

src/libslic3r/Arachne/utils/linearAlg2D.hpp

@@ -0,0 +1,118 @@
+//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 {len, 0};
+        return (p0.cast<int64_t>() * int64_t(len) / _len).cast<coord_t>();
+    };
+
+    constexpr coord_t normal_length = 10000; //Create a normal vector of reasonable length in order to reduce rounding error.
+    const Point ba = normal(a - b, normal_length);
+    const Point bc = normal(c - b, normal_length);
+    const Vec2d bq = query_point.cast<double>() - b.cast<double>();
+    const Vec2d perpendicular = perp(bq); //The query projects to this perpendicular to coordinate 0.
+
+    const double project_a_perpendicular = ba.cast<double>().dot(perpendicular); //Project vertex A on the perpendicular line.
+    const double project_c_perpendicular = bc.cast<double>().dot(perpendicular); //Project vertex C on the perpendicular line.
+    if ((project_a_perpendicular > 0.) != (project_c_perpendicular > 0.)) //Query is between A and C on the projection.
+    {
+        return project_a_perpendicular > 0.; //Due to the winding order of corner ABC, this means that the query is inside.
+    }
+    else //Beyond either A or C, but it could still be inside of the polygon.
+    {
+        const double project_a_parallel = ba.cast<double>().dot(bq); //Project not on the perpendicular, but on the original.
+        const double project_c_parallel = bc.cast<double>().dot(bq);
+
+        //Either:
+        // * A is to the right of B (project_a_perpendicular > 0) and C is below A (project_c_parallel < project_a_parallel), or
+        // * A is to the left of B (project_a_perpendicular < 0) and C is above A (project_c_parallel > project_a_parallel).
+        return (project_c_parallel < project_a_parallel) == (project_a_perpendicular > 0.);
+    }
+}
+
+/*!
+     * Returns the determinant of the 2D matrix defined by the the vectors ab and ap as rows.
+     * 
+     * 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