WIP "ensure verticall wall thickness" rework:

1) New region expansion code to propagate wave from a boundary
   of a region inside of it.
2) get_extents() extended with a template attribute to work with
   zero area data sets.
3) ClipperZUtils.hpp for handling Clipper operation with Z coordinate
   (for source contour identification)

src/clipper/clipper.cpp

@@ -4093,19 +4093,40 @@ void AddPolyNodeToPaths(const PolyNode& polynode, NodeType nodetype, Paths& path
   for (int i = 0; i < polynode.ChildCount(); ++i)
     AddPolyNodeToPaths(*polynode.Childs[i], nodetype, paths);
 }
+
+void AddPolyNodeToPaths(PolyNode&& polynode, NodeType nodetype, Paths& paths)
+{
+  bool match = true;
+  if (nodetype == ntClosed) match = !polynode.IsOpen();
+  else if (nodetype == ntOpen) return;
+
+  if (!polynode.Contour.empty() && match)
+    paths.push_back(std::move(polynode.Contour));
+  for (int i = 0; i < polynode.ChildCount(); ++i)
+    AddPolyNodeToPaths(std::move(*polynode.Childs[i]), nodetype, paths);
+}
+
 //------------------------------------------------------------------------------
 
 void PolyTreeToPaths(const PolyTree& polytree, Paths& paths)
 {
-  paths.resize(0); 
+  paths.clear();
   paths.reserve(polytree.Total());
   AddPolyNodeToPaths(polytree, ntAny, paths);
 }
+
+void PolyTreeToPaths(PolyTree&& polytree, Paths& paths)
+{
+  paths.clear();
+  paths.reserve(polytree.Total());
+  AddPolyNodeToPaths(std::move(polytree), ntAny, paths);
+}
+
 //------------------------------------------------------------------------------
 
 void ClosedPathsFromPolyTree(const PolyTree& polytree, Paths& paths)
 {
-  paths.resize(0); 
+  paths.clear();
   paths.reserve(polytree.Total());
   AddPolyNodeToPaths(polytree, ntClosed, paths);
 }
@@ -4113,7 +4134,7 @@ void ClosedPathsFromPolyTree(const PolyTree& polytree, Paths& paths)
 
 void OpenPathsFromPolyTree(PolyTree& polytree, Paths& paths)
 {
-  paths.resize(0); 
+  paths.clear();
   paths.reserve(polytree.Total());
   //Open paths are top level only, so ...
   for (int i = 0; i < polytree.ChildCount(); ++i)

src/clipper/clipper.hpp

@@ -206,6 +206,7 @@ void MinkowskiSum(const Path& pattern, const Paths& paths, Paths& solution, bool
 void MinkowskiDiff(const Path& poly1, const Path& poly2, Paths& solution);
 
 void PolyTreeToPaths(const PolyTree& polytree, Paths& paths);
+void PolyTreeToPaths(PolyTree&& polytree, Paths& paths);
 void ClosedPathsFromPolyTree(const PolyTree& polytree, Paths& paths);
 void OpenPathsFromPolyTree(PolyTree& polytree, Paths& paths);
 

src/libslic3r/Algorithm/RegionExpansion.cpp

@@ -0,0 +1,351 @@
+#include "RegionExpansion.hpp"
+
+#include <libslic3r/ClipperZUtils.hpp>
+#include <libslic3r/ClipperUtils.hpp>
+
+namespace Slic3r {
+namespace Algorithm {
+
+// similar to expolygons_to_zpaths(), but each contour is expanded before converted to zpath.
+// The expanded contours are then opened (the first point is repeated at the end).
+static ClipperLib_Z::Paths expolygons_to_zpaths_expanded_opened(
+    const ExPolygons &src, const float expansion, coord_t base_idx)
+{
+    ClipperLib_Z::Paths out;
+    out.reserve(2 * std::accumulate(src.begin(), src.end(), size_t(0),
+        [](const size_t acc, const ExPolygon &expoly) { return acc + expoly.num_contours(); }));
+    coord_t z = base_idx;
+    ClipperLib::ClipperOffset offsetter;
+    offsetter.ShortestEdgeLength = expansion * ClipperOffsetShortestEdgeFactor;
+    ClipperLib::Paths expansion_cache;
+    for (const ExPolygon &expoly : src) {
+        for (size_t icontour = 0; icontour < expoly.num_contours(); ++ icontour) {
+            // Execute reorients the contours so that the outer most contour has a positive area. Thus the output
+            // contours will be CCW oriented even though the input paths are CW oriented.
+            // Offset is applied after contour reorientation, thus the signum of the offset value is reversed.
+            offsetter.Clear();
+            offsetter.AddPath(expoly.contour_or_hole(icontour).points, ClipperLib::jtSquare, ClipperLib::etClosedPolygon);
+            expansion_cache.clear();
+            offsetter.Execute(expansion_cache, icontour == 0 ? expansion : -expansion);
+            append(out, ClipperZUtils::to_zpaths<true>(expansion_cache, z));
+        }
+        ++ z;
+    }
+    return out;
+}
+
+// Paths were created by splitting closed polygons into open paths and then by clipping them.
+// Thus some pieces of the clipped polygons may now become split at the ends of the source polygons.
+// Those ends are sorted lexicographically in "splits".
+// Reconnect those split pieces.
+static inline void merge_splits(ClipperLib_Z::Paths &paths, std::vector<std::pair<ClipperLib_Z::IntPoint, int>> &splits)
+{
+    for (auto it_path = paths.begin(); it_path != paths.end(); ) {
+        ClipperLib_Z::Path &path = *it_path;
+        assert(path.size() >= 2);
+        bool merged = false;
+        if (path.size() >= 2) {
+            const ClipperLib_Z::IntPoint &front = path.front();
+            const ClipperLib_Z::IntPoint &back  = path.back();
+            // The path before clipping was supposed to cross the clipping boundary or be fully out of it.
+            // Thus the clipped contour is supposed to become open.
+            assert(front.x() != back.x() || front.y() != back.y());
+            if (front.x() != back.x() || front.y() != back.y()) {
+                // Look up the ends in "splits", possibly join the contours.
+                // "splits" maps into the other piece connected to the same end point.
+                auto find_end = [&splits](const ClipperLib_Z::IntPoint &pt) -> std::pair<ClipperLib_Z::IntPoint, int>* {
+                    auto it = std::lower_bound(splits.begin(), splits.end(), pt,
+                        [](const auto &l, const auto &r){ return ClipperZUtils::zpoint_lower(l.first, r); });
+                    return it != splits.end() && it->first == pt ? &(*it) : nullptr;
+                };
+                auto *end = find_end(front);
+                bool  end_front = true;
+                if (! end) {
+                    end_front = false;
+                    end = find_end(back);
+                }
+                if (end) {
+                    // This segment ends at a split point of the source closed contour before clipping.
+                    if (end->second == -1) {
+                        // Open end was found, not matched yet.
+                        end->second = int(it_path - paths.begin());
+                    } else {
+                        // Open end was found and matched with end->second
+                        ClipperLib_Z::Path &other_path = paths[end->second];
+                        polylines_merge(other_path, other_path.front() == end->first, std::move(path), end_front);
+                        if (std::next(it_path) == paths.end()) {
+                            paths.pop_back();
+                            break;
+                        }
+                        path = std::move(paths.back());
+                        paths.pop_back();
+                        merged = true;
+                    }
+                }
+            }
+        }
+        if (! merged)
+            ++ it_path;
+    }
+}
+
+static ClipperLib::Paths wavefront_initial(ClipperLib::ClipperOffset &co, const ClipperLib::Paths &polylines, float offset)
+{
+    ClipperLib::Paths out;
+    out.reserve(polylines.size());
+    ClipperLib::Paths out_this;
+    for (const ClipperLib::Path &path : polylines) {
+        co.Clear();
+        co.AddPath(path, jtRound, ClipperLib::etOpenRound);
+        co.Execute(out_this, offset);
+        append(out, std::move(out_this));
+    }
+    return out;
+}
+
+// Input polygons may consist of multiple expolygons, even nested expolygons.
+// After inflation some polygons may thus overlap, however the overlap is being resolved during the successive
+// clipping operation, thus it is not being done here.
+static ClipperLib::Paths wavefront_step(ClipperLib::ClipperOffset &co, const ClipperLib::Paths &polygons, float offset)
+{
+    ClipperLib::Paths out;
+    out.reserve(polygons.size());
+    ClipperLib::Paths out_this;
+    for (const ClipperLib::Path &polygon : polygons) {
+        co.Clear();
+        // Execute reorients the contours so that the outer most contour has a positive area. Thus the output
+        // contours will be CCW oriented even though the input paths are CW oriented.
+        // Offset is applied after contour reorientation, thus the signum of the offset value is reversed.
+        co.AddPath(polygon, jtRound, ClipperLib::etClosedPolygon);
+        bool ccw = ClipperLib::Orientation(polygon);
+        co.Execute(out_this, ccw ? offset : - offset);
+        if (! ccw) {
+            // Reverse the resulting contours.
+            for (ClipperLib::Path &path : out_this)
+                std::reverse(path.begin(), path.end());
+        }
+        append(out, std::move(out_this));
+    }
+    return out;
+}
+
+static ClipperLib::Paths wavefront_clip(const ClipperLib::Paths &wavefront, const Polygons &clipping)
+{
+    ClipperLib::Clipper clipper;
+    clipper.AddPaths(wavefront, ClipperLib::ptSubject, true);
+    clipper.AddPaths(ClipperUtils::PolygonsProvider(clipping),  ClipperLib::ptClip, true);
+    ClipperLib::Paths out;
+    clipper.Execute(ClipperLib::ctIntersection, out, ClipperLib::pftPositive, ClipperLib::pftPositive);
+    return out;
+}
+
+static Polygons propagate_wave_from_boundary(
+    ClipperLib::ClipperOffset   &co,
+    // Seed of the wave: Open polylines very close to the boundary.
+    const ClipperLib::Paths     &seed,
+    // Boundary inside which the waveform will propagate.
+    const ExPolygon             &boundary,
+    // How much to inflate the seed lines to produce the first wave area.
+    const float                  initial_step,
+    // How much to inflate the first wave area and the successive wave areas in each step.
+    const float                  other_step,
+    // Number of inflate steps after the initial step.
+    const size_t                 num_other_steps,
+    // Maximum inflation of seed contours over the boundary. Used to trim boundary to speed up
+    // clipping during wave propagation.
+    const float                  max_inflation)
+{
+    assert(! seed.empty() && seed.front().size() >= 2);
+    Polygons clipping = ClipperUtils::clip_clipper_polygons_with_subject_bbox(boundary, get_extents<true>(seed).inflated(max_inflation));
+    ClipperLib::Paths polygons = wavefront_clip(wavefront_initial(co, seed, initial_step), clipping);
+    // Now offset the remaining 
+    for (size_t ioffset = 0; ioffset < num_other_steps; ++ ioffset)
+        polygons = wavefront_clip(wavefront_step(co, polygons, other_step), clipping);
+    return to_polygons(polygons);
+}
+
+// Calculating radius discretization according to ClipperLib offsetter code, see void ClipperOffset::DoOffset(double delta)
+inline double clipper_round_offset_error(double offset, double arc_tolerance)
+{
+    static constexpr const double def_arc_tolerance = 0.25;
+    const double y =
+        arc_tolerance <= 0 ?
+            def_arc_tolerance :
+            arc_tolerance > offset * def_arc_tolerance ?
+                offset * def_arc_tolerance :
+                arc_tolerance;
+    double steps = std::min(M_PI / std::acos(1. - y / offset), offset * M_PI);
+    return offset * (1. - cos(M_PI / steps));
+}
+
+// Returns regions per source ExPolygon expanded into boundary.
+std::vector<Polygons> expand_expolygons(
+    // Source regions that are supposed to touch the boundary.
+    // Boundaries of source regions touching the "boundary" regions will be expanded into the "boundary" region.
+    const ExPolygons    &src,
+    const ExPolygons    &boundary,
+    // Scaled expansion value
+    float                full_expansion,
+    // Expand by waves of expansion_step size (expansion_step is scaled).
+    float                expansion_step,
+    // Don't take more than max_nr_steps for small expansion_step.
+    size_t               max_nr_expansion_steps)
+{
+    assert(full_expansion > 0);
+    assert(expansion_step > 0);
+    assert(max_nr_expansion_steps > 0);
+
+    // Initial expansion of src to make the source regions intersect with boundary regions just a bit.
+    float                  tiny_expansion;
+    // How much to inflate the seed lines to produce the first wave area.
+    float                  initial_step;
+    // How much to inflate the first wave area and the successive wave areas in each step.
+    float                  other_step;
+    // Number of inflate steps after the initial step.
+    size_t                 num_other_steps;
+    // Maximum inflation of seed contours over the boundary. Used to trim boundary to speed up
+    // clipping during wave propagation.
+    float                  max_inflation;
+
+    // Offsetter to be applied for all inflation waves. Its accuracy is set with the block below.
+    ClipperLib::ClipperOffset co;
+
+    {
+        // Initial expansion of src to make the source regions intersect with boundary regions just a bit.
+        // The expansion should not be too tiny, but also small enough, so the following expansion will
+        // compensate for tiny_expansion and bring the wave back to the boundary without producing
+        // ugly cusps where it touches the boundary.
+        tiny_expansion = std::min(0.25f * full_expansion, scaled<float>(0.05f));
+        size_t nsteps = size_t(ceil((full_expansion - tiny_expansion) / expansion_step));
+        if (max_nr_expansion_steps > 0)
+            nsteps = std::min(nsteps, max_nr_expansion_steps);
+        assert(nsteps > 0);
+        initial_step = (full_expansion - tiny_expansion) / nsteps;
+        if (nsteps > 1 && 0.25 * initial_step < tiny_expansion) {
+            // Decrease the step size by lowering number of steps.
+            nsteps       = std::max<size_t>(1, (floor((full_expansion - tiny_expansion) / (4. * tiny_expansion))));
+            initial_step = (full_expansion - tiny_expansion) / nsteps;
+        }
+        if (0.25 * initial_step < tiny_expansion || nsteps == 1) {
+            tiny_expansion = 0.2f * full_expansion;
+            initial_step   = 0.8f * full_expansion;
+        }
+        other_step = initial_step;
+        num_other_steps = nsteps - 1;
+
+        // Accuracy of the offsetter for wave propagation.
+        co.ArcTolerance       = float(scale_(0.1));
+        co.ShortestEdgeLength = std::abs(initial_step * ClipperOffsetShortestEdgeFactor);
+
+        // Maximum inflation of seed contours over the boundary. Used to trim boundary to speed up
+        // clipping during wave propagation. Needs to be in sync with the offsetter accuracy.
+        // Clipper positive round offset should rather offset less than more.
+        // Still a little bit of additional offset was added.
+        max_inflation = (tiny_expansion + nsteps * initial_step) * 1.1;
+//                (clipper_round_offset_error(tiny_expansion, co.ArcTolerance) + nsteps * clipper_round_offset_error(initial_step, co.ArcTolerance) * 1.5; // Account for uncertainty 
+    }
+
+    using Intersection  = ClipperZUtils::ClipperZIntersectionVisitor::Intersection;
+    using Intersections = ClipperZUtils::ClipperZIntersectionVisitor::Intersections;
+
+    ClipperLib_Z::Paths expansion_seeds;
+    Intersections       intersections;
+
+    coord_t             idx_boundary_begin = 1;
+    coord_t             idx_boundary_end;
+    coord_t             idx_src_end;
+
+    {
+        ClipperLib_Z::Clipper zclipper;
+        ClipperZUtils::ClipperZIntersectionVisitor visitor(intersections);
+        zclipper.ZFillFunction(visitor.clipper_callback());
+        // as closed contours
+        {
+            ClipperLib_Z::Paths zboundary = ClipperZUtils::expolygons_to_zpaths(boundary, idx_boundary_begin);
+            idx_boundary_end = idx_boundary_begin + coord_t(zboundary.size());
+            zclipper.AddPaths(zboundary, ClipperLib_Z::ptClip, true);
+        }
+        // as open contours
+        std::vector<std::pair<ClipperLib_Z::IntPoint, int>> zsrc_splits;
+        {
+            ClipperLib_Z::Paths zsrc = expolygons_to_zpaths_expanded_opened(src, tiny_expansion, idx_boundary_end);
+            zclipper.AddPaths(zsrc, ClipperLib_Z::ptSubject, false);
+            idx_src_end = idx_boundary_end + coord_t(zsrc.size());
+            zsrc_splits.reserve(zsrc.size());
+            for (const ClipperLib_Z::Path &path : zsrc) {
+                assert(path.size() >= 2);
+                assert(path.front() == path.back());
+                zsrc_splits.emplace_back(path.front(), -1);
+            }
+            std::sort(zsrc_splits.begin(), zsrc_splits.end(), [](const auto &l, const auto &r){ return ClipperZUtils::zpoint_lower(l.first, r.first); });
+        }
+        ClipperLib_Z::PolyTree polytree;
+        zclipper.Execute(ClipperLib_Z::ctIntersection, polytree, ClipperLib_Z::pftNonZero, ClipperLib_Z::pftNonZero);
+        ClipperLib_Z::PolyTreeToPaths(std::move(polytree), expansion_seeds);
+        merge_splits(expansion_seeds, zsrc_splits);
+    }
+
+    // Sort paths into their respective islands.
+    // Each src x boundary will be processed (wave expanded) independently.
+    // Multiple pieces of a single src may intersect the same boundary.
+    struct SeedOrigin {
+        int     src;
+        int     boundary;
+        int     seed;
+    };
+    std::vector<SeedOrigin> map_seeds;
+    map_seeds.reserve(expansion_seeds.size());
+    int iseed = 0;
+    for (const ClipperLib_Z::Path &path : expansion_seeds) {
+        assert(path.size() >= 2);
+        const ClipperLib_Z::IntPoint &front = path.front();
+        const ClipperLib_Z::IntPoint &back  = path.back();
+        // Both ends of a seed segment are supposed to be inside a single boundary expolygon.
+        assert(front.z() < 0);
+        assert(back.z() < 0);
+        const Intersection *intersection = nullptr;
+        auto intersection_point_valid = [idx_boundary_end, idx_src_end](const Intersection &is) {
+            return is.first >= 1 && is.first < idx_boundary_end &&
+                   is.second >= idx_boundary_end && is.second < idx_src_end;
+        };
+        if (front.z() < 0) {
+            const Intersection &is = intersections[- front.z() - 1];
+            assert(intersection_point_valid(is));
+            if (intersection_point_valid(is))
+                intersection = &is;
+        }
+        if (! intersection && back.z() < 0) {
+            const Intersection &is = intersections[- back.z() - 1];
+            assert(intersection_point_valid(is));
+            if (intersection_point_valid(is))
+                intersection = &is;
+        }
+        if (intersection) {
+            // The path intersects the boundary contour at least at one side. 
+            map_seeds.push_back({ intersection->second - idx_boundary_end, intersection->first - 1, iseed });
+        }
+        ++ iseed;
+    }
+    // Sort the seeds by their intersection boundary and source contour.
+    std::sort(map_seeds.begin(), map_seeds.end(), [](const auto &l, const auto &r){ 
+        return l.boundary < r.boundary || (l.boundary == r.boundary && l.src < r.src); 
+    });
+
+    std::vector<Polygons> out(src.size(), Polygons{});
+    ClipperLib::Paths     paths;
+    for (auto it_seed = map_seeds.begin(); it_seed != map_seeds.end();) {
+        auto it = it_seed;
+        paths.clear();
+        for (; it != map_seeds.end() && it->boundary == it_seed->boundary && it->src == it_seed->src; ++ it)
+            paths.emplace_back(ClipperZUtils::from_zpath(expansion_seeds[it->seed]));
+        // Propagate the wavefront while clipping it with the trimmed boundary.
+        // Collect the expanded polygons, merge them with the source polygons.
+        append(out[it_seed->src], propagate_wave_from_boundary(co, paths, boundary[it_seed->boundary], initial_step, other_step, num_other_steps, max_inflation));
+        it_seed = it;
+    }
+
+    return out;
+}
+
+} // Algorithm
+} // Slic3r

src/libslic3r/Algorithm/RegionExpansion.hpp

@@ -0,0 +1,26 @@
+#ifndef SRC_LIBSLIC3R_ALGORITHM_REGION_EXPANSION_HPP_
+#define SRC_LIBSLIC3R_ALGORITHM_REGION_EXPANSION_HPP_
+
+#include <cstdint>
+
+namespace Slic3r {
+
+class Polygon;
+using Polygons = std::vector<Polygon>;
+class ExPolygon;
+using ExPolygons = std::vector<ExPolygon>;
+
+namespace Algorithm {
+
+std::vector<Polygons> expand_expolygons(const ExPolygons &src, const ExPolygons &boundary,
+    // Scaled expansion value
+    float expansion, 
+    // Expand by waves of expansion_step size (expansion_step is scaled).
+    float expansion_step,
+    // Don't take more than max_nr_steps for small expansion_step.
+    size_t max_nr_steps);
+
+} // Algorithm
+} // Slic3r
+
+#endif /* SRC_LIBSLIC3R_ALGORITHM_REGION_EXPANSION_HPP_ */

src/libslic3r/BoundingBox.hpp

@@ -23,23 +23,8 @@ public:
         min(p1), max(p1), defined(false) { merge(p2); merge(p3); }
 
     template<class It, class = IteratorOnly<It>>
-    BoundingBoxBase(It from, It to) : min(PointClass::Zero()), max(PointClass::Zero())
+    BoundingBoxBase(It from, It to)
-    {
+        { construct(*this, from, to); }
-        if (from == to) {
-            this->defined = false;
-            // throw Slic3r::InvalidArgument("Empty point set supplied to BoundingBoxBase constructor");
-        } else {
-            auto it = from;
-            this->min = it->template cast<typename PointClass::Scalar>();
-            this->max = this->min;
-            for (++ it; it != to; ++ it) {
-                auto vec = it->template cast<typename PointClass::Scalar>();
-                this->min = this->min.cwiseMin(vec);
-                this->max = this->max.cwiseMax(vec);
-            }
-            this->defined = (this->min.x() < this->max.x()) && (this->min.y() < this->max.y());
-        }
-    }
 
     BoundingBoxBase(const std::vector<PointClass> &points)
         : BoundingBoxBase(points.begin(), points.end())
@@ -70,6 +55,30 @@ public:
     }
     bool operator==(const BoundingBoxBase<PointClass> &rhs) { return this->min == rhs.min && this->max == rhs.max; }
     bool operator!=(const BoundingBoxBase<PointClass> &rhs) { return ! (*this == rhs); }
+
+private:
+    // to access construct()
+    friend BoundingBox get_extents<false>(const Points &pts);
+    friend BoundingBox get_extents<true>(const Points &pts);
+
+    // if IncludeBoundary, then a bounding box is defined even for a single point.
+    // otherwise a bounding box is only defined if it has a positive area.
+    // The output bounding box is expected to be set to "undefined" initially.
+    template<bool IncludeBoundary = false, class BoundingBoxType, class It, class = IteratorOnly<It>>
+    static void construct(BoundingBoxType &out, It from, It to)
+    {
+        if (from != to) {
+            auto it = from;
+            out.min = it->template cast<typename PointClass::Scalar>();
+            out.max = out.min;
+            for (++ it; it != to; ++ it) {
+                auto vec = it->template cast<typename PointClass::Scalar>();
+                out.min = out.min.cwiseMin(vec);
+                out.max = out.max.cwiseMax(vec);
+            }
+            out.defined = IncludeBoundary || (out.min.x() < out.max.x() && out.min.y() < out.max.y());
+        }
+    }
 };
 
 template <class PointClass>

src/libslic3r/Brim.cpp

@@ -641,7 +641,7 @@ ExtrusionEntityCollection make_brim(const Print &print, PrintTryCancel try_cance
 			// perform operation
 			ClipperLib_Z::PolyTree loops_trimmed_tree;
 			clipper.Execute(ClipperLib_Z::ctDifference, loops_trimmed_tree, ClipperLib_Z::pftNonZero, ClipperLib_Z::pftNonZero);
-			ClipperLib_Z::PolyTreeToPaths(loops_trimmed_tree, loops_trimmed);
+			ClipperLib_Z::PolyTreeToPaths(std::move(loops_trimmed_tree), loops_trimmed);
 		}
 
 		// Third, produce the extrusions, sorted by the source loop indices.

src/libslic3r/CMakeLists.txt

@@ -22,6 +22,8 @@ set(SLIC3R_SOURCES
     AABBTreeLines.hpp
     AABBMesh.hpp
     AABBMesh.cpp
+    Algorithm/RegionExpansion.cpp
+    Algorithm/RegionExpansion.hpp
     BoundingBox.cpp
     BoundingBox.hpp
     BridgeDetector.cpp
@@ -35,6 +37,7 @@ set(SLIC3R_SOURCES
     clipper.hpp
     ClipperUtils.cpp
     ClipperUtils.hpp
+    ClipperZUtils.hpp
     Color.cpp
     Color.hpp
     Config.cpp

src/libslic3r/ClipperUtils.cpp

@@ -8,8 +8,6 @@
 #include "SVG.hpp"
 #endif /* CLIPPER_UTILS_DEBUG */
 
-#define CLIPPER_OFFSET_SHORTEST_EDGE_FACTOR (0.005f)
-
 namespace Slic3r {
 
 #ifdef CLIPPER_UTILS_DEBUG
@@ -267,7 +265,7 @@ static ClipperLib::Paths raw_offset(PathsProvider &&paths, float offset, Clipper
         co.ArcTolerance = miterLimit;
     else
         co.MiterLimit = miterLimit;
-    co.ShortestEdgeLength = double(std::abs(offset * CLIPPER_OFFSET_SHORTEST_EDGE_FACTOR));
+    co.ShortestEdgeLength = std::abs(offset * ClipperOffsetShortestEdgeFactor);
     for (const ClipperLib::Path &path : paths) {
         co.Clear();
         // Execute reorients the contours so that the outer most contour has a positive area. Thus the output
@@ -414,7 +412,7 @@ static int offset_expolygon_inner(const Slic3r::ExPolygon &expoly, const float d
             co.ArcTolerance = miterLimit;
         else
             co.MiterLimit = miterLimit;
-        co.ShortestEdgeLength = double(std::abs(delta * CLIPPER_OFFSET_SHORTEST_EDGE_FACTOR));
+        co.ShortestEdgeLength = std::abs(delta * ClipperOffsetShortestEdgeFactor);
         co.AddPath(expoly.contour.points, joinType, ClipperLib::etClosedPolygon);
         co.Execute(contours, delta);
     }
@@ -435,7 +433,7 @@ static int offset_expolygon_inner(const Slic3r::ExPolygon &expoly, const float d
                     co.ArcTolerance = miterLimit;
                 else
                     co.MiterLimit = miterLimit;
-                co.ShortestEdgeLength = double(std::abs(delta * CLIPPER_OFFSET_SHORTEST_EDGE_FACTOR));
+                co.ShortestEdgeLength = std::abs(delta * ClipperOffsetShortestEdgeFactor);
                 co.AddPath(hole.points, joinType, ClipperLib::etClosedPolygon);
                 ClipperLib::Paths out2;
                 // Execute reorients the contours so that the outer most contour has a positive area. Thus the output
@@ -1055,7 +1053,7 @@ ClipperLib::Path mittered_offset_path_scaled(const Points &contour, const std::v
 		};
 
 		// Minimum edge length, squared.
-		double lmin  = *std::max_element(deltas.begin(), deltas.end()) * CLIPPER_OFFSET_SHORTEST_EDGE_FACTOR;
+		double lmin  = *std::max_element(deltas.begin(), deltas.end()) * ClipperOffsetShortestEdgeFactor;
 		double l2min = lmin * lmin;
 		// Minimum angle to consider two edges to be parallel.
 		// Vojtech's estimate.

src/libslic3r/ClipperUtils.hpp

@@ -39,6 +39,9 @@ static constexpr const Slic3r::ClipperLib::JoinType DefaultLineJoinType     = Sl
 // Miter limit is ignored for jtSquare.
 static constexpr const double                       DefaultLineMiterLimit   = 0.;
 
+// Decimation factor applied on input contour when doing offset, multiplied by the offset distance.
+static constexpr const double                       ClipperOffsetShortestEdgeFactor = 0.005;
+
 enum class ApplySafetyOffset {
     No,
     Yes
@@ -599,6 +602,6 @@ Polygons  variable_offset_outer(const ExPolygon &expoly, const std::vector<std::
 ExPolygons variable_offset_outer_ex(const ExPolygon &expoly, const std::vector<std::vector<float>> &deltas, double miter_limit = 2.);
 ExPolygons variable_offset_inner_ex(const ExPolygon &expoly, const std::vector<std::vector<float>> &deltas, double miter_limit = 2.);
 
-}
+} // namespace Slic3r
 
-#endif
+#endif // slic3r_ClipperUtils_hpp_

src/libslic3r/ClipperZUtils.hpp

@@ -0,0 +1,143 @@
+#ifndef slic3r_ClipperZUtils_hpp_
+#define slic3r_ClipperZUtils_hpp_
+
+#include <vector>
+
+#include <clipper/clipper_z.hpp>
+#include <libslic3r/Point.hpp>
+
+namespace Slic3r {
+
+namespace ClipperZUtils {
+
+using ZPoint  = ClipperLib_Z::IntPoint;
+using ZPoints = ClipperLib_Z::Path;
+using ZPath   = ClipperLib_Z::Path;
+using ZPaths  = ClipperLib_Z::Paths;
+
+inline bool zpoint_lower(const ZPoint &l, const ZPoint &r)
+{
+    return l.x() < r.x() || (l.x() == r.x() && (l.y() < r.y() || (l.y() == r.y() && l.z() < r.z())));
+}
+
+// Convert a single path to path with a given Z coordinate.
+// If Open, then duplicate the first point at the end.
+template<bool Open = false>
+inline ZPath to_zpath(const Points &path, coord_t z)
+{
+    ZPath out;
+    if (! path.empty()) {
+        out.reserve(path.size() + Open ? 1 : 0);
+        for (const Point &p : path)
+            out.emplace_back(p.x(), p.y(), z);
+        if (Open)
+            out.emplace_back(out.front());
+    }
+    return out;
+}
+
+// Convert multiple paths to paths with a given Z coordinate.
+// If Open, then duplicate the first point of each path at its end.
+template<bool Open = false>
+inline ZPaths to_zpaths(const std::vector<Points> &paths, coord_t z)
+{
+    ZPaths out;
+    out.reserve(paths.size());
+    for (const Points &path : paths)
+        out.emplace_back(to_zpath<Open>(path, z));
+    return out;
+}
+
+// Convert multiple expolygons into z-paths with Z specified by an index of the source expolygon
+// offsetted by base_index.
+// If Open, then duplicate the first point of each path at its end.
+template<bool Open = false>
+inline ZPaths expolygons_to_zpaths(const ExPolygons &src, coord_t base_idx)
+{
+    ZPaths out;
+    out.reserve(std::accumulate(src.begin(), src.end(), size_t(0),
+        [](const size_t acc, const ExPolygon &expoly) { return acc + expoly.num_contours(); }));
+    coord_t z = base_idx;
+    for (const ExPolygon &expoly : src) {
+        out.emplace_back(to_zpath<Open>(expoly.contour.points, z));
+        for (const Polygon &hole : expoly.holes)
+            out.emplace_back(to_zpath<Open>(hole.points, z));
+        ++ z;
+    }
+    return out;
+}
+
+// Convert a single path to path with a given Z coordinate.
+// If Open, then duplicate the first point at the end.
+template<bool Open = false>
+inline Points from_zpath(const ZPoints &path)
+{
+    Points out;
+    if (! path.empty()) {
+        out.reserve(path.size() + Open ? 1 : 0);
+        for (const ZPoint &p : path)
+            out.emplace_back(p.x(), p.y());
+        if (Open)
+            out.emplace_back(out.front());
+    }
+    return out;
+}
+
+// Convert multiple paths to paths with a given Z coordinate.
+// If Open, then duplicate the first point of each path at its end.
+template<bool Open = false>
+inline void from_zpaths(const ZPaths &paths, std::vector<Points> &out)
+{
+    out.reserve(out.size() + paths.size());
+    for (const ZPoints &path : paths)
+        out.emplace_back(from_zpath<Open>(path));
+}
+template<bool Open = false>
+inline std::vector<Points> from_zpaths(const ZPaths &paths)
+{
+    std::vector<Points> out;
+    from_zpaths(paths, out);
+    return out;
+}
+
+class ClipperZIntersectionVisitor {
+public:
+    using Intersection  = std::pair<coord_t, coord_t>;
+    using Intersections = std::vector<Intersection>;
+    ClipperZIntersectionVisitor(Intersections &intersections) : m_intersections(intersections) {}
+    void reset() { m_intersections.clear(); }
+    void operator()(const ZPoint &e1bot, const ZPoint &e1top, const ZPoint &e2bot, const ZPoint &e2top, ZPoint &pt) {
+        coord_t srcs[4]{ e1bot.z(), e1top.z(), e2bot.z(), e2top.z() };
+        coord_t *begin = srcs;
+        coord_t *end = srcs + 4;
+        //FIXME bubble sort manually?
+        std::sort(begin, end);
+        end = std::unique(begin, end);
+        if (begin + 1 == end) {
+            // Self intersection may happen on source contour. Just copy the Z value.
+            pt.z() = *begin;
+        } else {
+            assert(begin + 2 == end);
+            if (begin + 2 <= end) {
+                // store a -1 based negative index into the "intersections" vector here.
+                m_intersections.emplace_back(srcs[0], srcs[1]);
+                pt.z() = -coord_t(m_intersections.size());
+            }
+        }
+    }
+    ClipperLib_Z::ZFillCallback clipper_callback() {
+        return [this](const ZPoint &e1bot, const ZPoint &e1top, 
+                 const ZPoint &e2bot, const ZPoint &e2top, ZPoint &pt)
+        { return (*this)(e1bot, e1top, e2bot, e2top, pt); };
+    }
+
+    const std::vector<std::pair<coord_t, coord_t>>& intersections() const { return m_intersections; }
+
+private:
+    std::vector<std::pair<coord_t, coord_t>> &m_intersections;
+};
+
+} // namespace ClipperZUtils
+} // namespace Slic3r
+
+#endif // slic3r_ClipperZUtils_hpp_

src/libslic3r/Layer.cpp

@@ -1,5 +1,5 @@
 #include "Layer.hpp"
-#include <clipper/clipper_z.hpp>
+#include "ClipperZUtils.hpp"
 #include "ClipperUtils.hpp"
 #include "Print.hpp"
 #include "Fill/Fill.hpp"
@@ -304,48 +304,16 @@ void Layer::build_up_down_graph(Layer& below, Layer& above)
     coord_t             paths_end = paths_above_offset + coord_t(above.lslices.size());
 #endif // NDEBUG
 
-    class ZFill {
-    public:
-        ZFill() = default;
-        void reset() { m_intersections.clear(); }
-        void operator()(
-            const ClipperLib_Z::IntPoint& e1bot, const ClipperLib_Z::IntPoint& e1top,
-            const ClipperLib_Z::IntPoint& e2bot, const ClipperLib_Z::IntPoint& e2top,
-            ClipperLib_Z::IntPoint& pt) {
-            coord_t srcs[4]{ e1bot.z(), e1top.z(), e2bot.z(), e2top.z() };
-            coord_t* begin = srcs;
-            coord_t* end = srcs + 4;
-            std::sort(begin, end);
-            end = std::unique(begin, end);
-            if (begin + 1 == end) {
-                // Self intersection may happen on source contour. Just copy the Z value.
-                pt.z() = *begin;
-            } else {
-                assert(begin + 2 == end);
-                if (begin + 2 <= end) {
-                    // store a -1 based negative index into the "intersections" vector here.
-                    m_intersections.emplace_back(srcs[0], srcs[1]);
-                    pt.z() = -coord_t(m_intersections.size());
-                }
-            }
-        }
-        const std::vector<std::pair<coord_t, coord_t>>& intersections() const { return m_intersections; }
-
-    private:
-        std::vector<std::pair<coord_t, coord_t>> m_intersections;
-    } zfill;
-
     ClipperLib_Z::Clipper  clipper;
     ClipperLib_Z::PolyTree result;
-    clipper.ZFillFunction(
+    ClipperZUtils::ClipperZIntersectionVisitor::Intersections intersections;
-        [&zfill](const ClipperLib_Z::IntPoint &e1bot, const ClipperLib_Z::IntPoint &e1top, 
+    ClipperZUtils::ClipperZIntersectionVisitor visitor(intersections);
-                 const ClipperLib_Z::IntPoint &e2bot, const ClipperLib_Z::IntPoint &e2top, ClipperLib_Z::IntPoint &pt)
+    clipper.ZFillFunction(visitor.clipper_callback());
-        { return zfill(e1bot, e1top, e2bot, e2top, pt); });
     clipper.AddPaths(paths_below, ClipperLib_Z::ptSubject, true);
     clipper.AddPaths(paths_above, ClipperLib_Z::ptClip, true);
     clipper.Execute(ClipperLib_Z::ctIntersection, result, ClipperLib_Z::pftNonZero, ClipperLib_Z::pftNonZero);
 
-    connect_layer_slices(below, above, result, zfill.intersections(), paths_below_offset, paths_above_offset
+    connect_layer_slices(below, above, result, intersections, paths_below_offset, paths_above_offset
 #ifndef NDEBUG
         , paths_end
 #endif // NDEBUG

src/libslic3r/PerimeterGenerator.cpp

@@ -433,10 +433,12 @@ static ClipperLib_Z::Paths clip_extrusion(const ClipperLib_Z::Path &subject, con
     clipper.AddPath(subject, ClipperLib_Z::ptSubject, false);
     clipper.AddPaths(clip, ClipperLib_Z::ptClip, true);
 
-    ClipperLib_Z::PolyTree clipped_polytree;
     ClipperLib_Z::Paths    clipped_paths;
+    {
+        ClipperLib_Z::PolyTree clipped_polytree;
         clipper.Execute(clipType, clipped_polytree, ClipperLib_Z::pftNonZero, ClipperLib_Z::pftNonZero);
-    ClipperLib_Z::PolyTreeToPaths(clipped_polytree, clipped_paths);
+        ClipperLib_Z::PolyTreeToPaths(std::move(clipped_polytree), clipped_paths);
+    }
 
     // Clipped path could contain vertices from the clip with a Z coordinate equal to zero.
     // For those vertices, we must assign value based on the subject.

src/libslic3r/Point.cpp

@@ -84,18 +84,15 @@ Points collect_duplicates(Points pts /* Copy */)
     return duplicits;
 }
 
+template<bool IncludeBoundary>
 BoundingBox get_extents(const Points &pts)
 { 
-    return BoundingBox(pts);
+    BoundingBox out;
-}
+    BoundingBox::construct<IncludeBoundary>(out, pts.begin(), pts.end());
-
+    return out;
-BoundingBox get_extents(const std::vector<Points> &pts)
-{
-    BoundingBox bbox;
-    for (const Points &p : pts)
-        bbox.merge(get_extents(p));
-    return bbox;
 }
+template BoundingBox get_extents<false>(const Points &pts);
+template BoundingBox get_extents<true>(const Points &pts);
 
 BoundingBoxf get_extents(const std::vector<Vec2d> &pts)
 {

src/libslic3r/Point.hpp

@@ -229,8 +229,24 @@ inline Point lerp(const Point &a, const Point &b, double t)
     return ((1. - t) * a.cast<double>() + t * b.cast<double>()).cast<coord_t>();
 }
 
+// if IncludeBoundary, then a bounding box is defined even for a single point.
+// otherwise a bounding box is only defined if it has a positive area.
+template<bool IncludeBoundary = false>
 BoundingBox get_extents(const Points &pts);
-BoundingBox get_extents(const std::vector<Points> &pts);
+extern template BoundingBox get_extents<false>(const Points& pts);
+extern template BoundingBox get_extents<true>(const Points& pts);
+
+// if IncludeBoundary, then a bounding box is defined even for a single point.
+// otherwise a bounding box is only defined if it has a positive area.
+template<bool IncludeBoundary = false>
+BoundingBox get_extents(const std::vector<Points> &pts)
+{
+    BoundingBox bbox;
+    for (const Points &p : pts)
+        bbox.merge(get_extents<IncludeBoundary>(p));
+    return bbox;
+}
+
 BoundingBoxf get_extents(const std::vector<Vec2d> &pts);
 
 int nearest_point_index(const Points &points, const Point &pt);

src/libslic3r/Polyline.hpp

@@ -153,6 +153,25 @@ inline void polylines_append(Polylines &dst, Polylines &&src)
     }
 }
 
+// Merge polylines at their respective end points.
+// dst_first: the merge point is at dst.begin() or dst.end()?
+// src_first: the merge point is at src.begin() or src.end()?
+// The orientation of the resulting polyline is unknown, the output polyline may start
+// either with src piece or dst piece.
+template<typename PointType>
+inline void polylines_merge(std::vector<PointType> &dst, bool dst_first, std::vector<PointType> &&src, bool src_first)
+{
+    if (dst_first) {
+        if (src_first)
+            std::reverse(dst.begin(), dst.end());
+        else
+            std::swap(dst, src);
+    } else if (! src_first)
+        std::reverse(src.begin(), src.end());
+    // Merge src into dst.
+    append(dst, std::move(src));
+}
+
 const Point& leftmost_point(const Polylines &polylines);
 
 bool remove_degenerate(Polylines &polylines);

src/libslic3r/libslic3r.h

@@ -124,8 +124,7 @@ inline void append(std::vector<T>& dest, std::vector<T>&& src)
         dest.insert(dest.end(),
             std::make_move_iterator(src.begin()),
             std::make_move_iterator(src.end()));
-
+        // Release memory of the source contour now.
-        // Vojta wants back compatibility
         src.clear();
         src.shrink_to_fit();
     }
@@ -161,8 +160,7 @@ inline void append_reversed(std::vector<T>& dest, std::vector<T>&& src)
         dest.insert(dest.end(), 
             std::make_move_iterator(src.rbegin()),
             std::make_move_iterator(src.rend()));
-
+    // Release memory of the source contour now.
-    // Vojta wants back compatibility
     src.clear();
     src.shrink_to_fit();
 }

tests/fff_print/test_clipper.cpp

@@ -1,7 +1,7 @@
 #include <catch2/catch.hpp>
 
 #include "test_data.hpp"
-#include "clipper/clipper_z.hpp"
+#include "libslic3r/ClipperZUtils.hpp"
 #include "libslic3r/clipper.hpp"
 
 using namespace Slic3r;
@@ -132,3 +132,72 @@ SCENARIO("Clipper Z", "[ClipperZ]")
         REQUIRE(pt.z() == 1);
 }
 
+SCENARIO("Intersection with multiple polylines", "[ClipperZ]")
+{
+    // 1000x1000 CCQ square
+    ClipperLib_Z::Path clip { { 0, 0, 1 }, { 1000, 0, 1 }, { 1000, 1000, 1 }, { 0, 1000, 1 } };
+    // Two lines interseting inside the square above, crossing the bottom edge of the square.
+    ClipperLib_Z::Path line1 { { +100, -100, 2 }, { +900, +900, 2 } };
+    ClipperLib_Z::Path line2 { { +100, +900, 3 }, { +900, -100, 3 } };
+
+    ClipperLib_Z::Clipper clipper;
+    ClipperZUtils::ClipperZIntersectionVisitor::Intersections intersections;
+    ClipperZUtils::ClipperZIntersectionVisitor visitor(intersections);
+    clipper.ZFillFunction(visitor.clipper_callback());
+    clipper.AddPath(line1, ClipperLib_Z::ptSubject, false);
+    clipper.AddPath(line2, ClipperLib_Z::ptSubject, false);
+    clipper.AddPath(clip,  ClipperLib_Z::ptClip, true);
+
+    ClipperLib_Z::PolyTree polytree;
+    ClipperLib_Z::Paths    paths;
+    clipper.Execute(ClipperLib_Z::ctIntersection, polytree, ClipperLib_Z::pftNonZero, ClipperLib_Z::pftNonZero);
+    ClipperLib_Z::PolyTreeToPaths(polytree, paths);
+
+    REQUIRE(paths.size() == 2);
+
+    THEN("First output polyline is a trimmed 2nd line") {
+        // Intermediate point (intersection) was removed)
+        REQUIRE(paths.front().size() == 2);
+        REQUIRE(paths.front().front().z() == 3);
+        REQUIRE(paths.front().back().z() < 0);
+        REQUIRE(intersections[- paths.front().back().z() - 1] == std::pair<coord_t, coord_t>(1, 3));
+    }
+
+    THEN("Second output polyline is a trimmed 1st line") {
+        // Intermediate point (intersection) was removed)
+        REQUIRE(paths[1].size() == 2);
+        REQUIRE(paths[1].front().z() < 0);
+        REQUIRE(paths[1].back().z() == 2);
+        REQUIRE(intersections[- paths[1].front().z() - 1] == std::pair<coord_t, coord_t>(1, 2));
+    }
+}
+
+SCENARIO("Interseting a closed loop as an open polyline", "[ClipperZ]")
+{
+    // 1000x1000 CCQ square
+    ClipperLib_Z::Path clip{ { 0, 0, 1 }, { 1000, 0, 1 }, { 1000, 1000, 1 }, { 0, 1000, 1 } };
+    // Two lines interseting inside the square above, crossing the bottom edge of the square.
+    ClipperLib_Z::Path rect{ { 500, 500, 2}, { 500, 1500, 2 }, { 1500, 1500, 2}, { 500, 1500, 2}, { 500, 500, 2 } };
+
+    ClipperLib_Z::Clipper clipper;
+    clipper.AddPath(rect, ClipperLib_Z::ptSubject, false);
+    clipper.AddPath(clip, ClipperLib_Z::ptClip, true);
+
+    ClipperLib_Z::PolyTree polytree;
+    ClipperLib_Z::Paths    paths;
+    ClipperZUtils::ClipperZIntersectionVisitor::Intersections intersections;
+    ClipperZUtils::ClipperZIntersectionVisitor visitor(intersections);
+    clipper.ZFillFunction(visitor.clipper_callback());
+    clipper.Execute(ClipperLib_Z::ctIntersection, polytree, ClipperLib_Z::pftNonZero, ClipperLib_Z::pftNonZero);
+    ClipperLib_Z::PolyTreeToPaths(std::move(polytree), paths);
+
+    THEN("Open polyline is clipped into two pieces") {
+        REQUIRE(paths.size() == 2);
+        REQUIRE(paths.front().size() == 2);
+        REQUIRE(paths.back().size() == 2);
+        REQUIRE(paths.front().front().z() == 2);
+        REQUIRE(paths.back().back().z() == 2);
+        REQUIRE(paths.front().front().x() == paths.back().back().x());
+        REQUIRE(paths.front().front().y() == paths.back().back().y());
+    }
+}

tests/libslic3r/CMakeLists.txt

@@ -23,6 +23,7 @@ add_executable(${_TEST_NAME}_tests
 	test_stl.cpp
 	test_meshboolean.cpp
 	test_marchingsquares.cpp
+	test_region_expansion.cpp
 	test_timeutils.cpp
 	test_utils.cpp
 	test_voronoi.cpp

tests/libslic3r/test_region_expansion.cpp

@@ -0,0 +1,254 @@
+#include <catch2/catch.hpp>
+
+#include <libslic3r/libslic3r.h>
+#include <libslic3r/Algorithm/RegionExpansion.hpp>
+#include <libslic3r/ClipperUtils.hpp>
+#include <libslic3r/ExPolygon.hpp>
+#include <libslic3r/Polygon.hpp>
+#include <libslic3r/SVG.cpp>
+
+using namespace Slic3r;
+
+//#define DEBUG_TEMP_DIR "d:\\temp\\"
+
+SCENARIO("Region expansion basics", "[RegionExpansion]") {
+    static constexpr const coord_t ten = scaled<coord_t>(10.);
+    GIVEN("two touching squares") {
+        Polygon square1{ { 1 * ten, 1 * ten }, { 2 * ten, 1 * ten }, { 2 * ten, 2 * ten }, { 1 * ten, 2 * ten } };
+        Polygon square2{ { 2 * ten, 1 * ten }, { 3 * ten, 1 * ten }, { 3 * ten, 2 * ten }, { 2 * ten, 2 * ten } };
+        Polygon square3{ { 1 * ten, 2 * ten }, { 2 * ten, 2 * ten }, { 2 * ten, 3 * ten }, { 1 * ten, 3 * ten } };
+        static constexpr const float expansion = scaled<float>(1.);
+        auto test_expansion = [](const Polygon &src, const Polygon &boundary) {
+            std::vector<Polygons> expanded = Algorithm::expand_expolygons({ ExPolygon{src} }, { ExPolygon{boundary} },
+                expansion,
+                scaled<float>(0.3), // expansion step
+                5);                 // max num steps
+            THEN("Single anchor is produced") {
+                REQUIRE(expanded.size() == 1);
+            }
+            THEN("The area of the anchor is 10mm2") {
+                REQUIRE(area(expanded.front()) == Approx(expansion * ten));
+            }
+        };
+
+        WHEN("second square expanded into the first square (to left)") {
+            test_expansion(square2, square1);
+        }
+        WHEN("first square expanded into the second square (to right)") {
+            test_expansion(square1, square2);
+        }
+        WHEN("third square expanded into the first square (down)") {
+            test_expansion(square3, square1);
+        }
+        WHEN("first square expanded into the third square (up)") {
+            test_expansion(square1, square3);
+        }
+    }
+
+    GIVEN("simple bridge") {
+        Polygon square1{ { 1 * ten, 1 * ten }, { 2 * ten, 1 * ten }, { 2 * ten, 2 * ten }, { 1 * ten, 2 * ten } };
+        Polygon square2{ { 2 * ten, 1 * ten }, { 3 * ten, 1 * ten }, { 3 * ten, 2 * ten }, { 2 * ten, 2 * ten } };
+        Polygon square3{ { 3 * ten, 1 * ten }, { 4 * ten, 1 * ten }, { 4 * ten, 2 * ten }, { 3 * ten, 2 * ten } };
+
+        WHEN("expanded") {
+            static constexpr const float expansion = scaled<float>(1.);
+            std::vector<Polygons> expanded = Algorithm::expand_expolygons({ ExPolygon{square2} }, { ExPolygon{square1}, ExPolygon{square3} },
+                expansion,
+                scaled<float>(0.3), // expansion step
+                5);                 // max num steps
+            THEN("Two anchors are produced") {
+                REQUIRE(expanded.size() == 1);
+                REQUIRE(expanded.front().size() == 2);
+            }
+            THEN("The area of each anchor is 10mm2") {
+                REQUIRE(area(expanded.front().front()) == Approx(expansion * ten));
+                REQUIRE(area(expanded.front().back()) == Approx(expansion * ten));
+            }
+        }
+
+        WHEN("fully expanded") {
+            static constexpr const float expansion = scaled<float>(10.1);
+            std::vector<Polygons> expanded = Algorithm::expand_expolygons({ ExPolygon{square2} }, { ExPolygon{square1}, ExPolygon{square3} },
+                expansion,
+                scaled<float>(2.3), // expansion step
+                5);                 // max num steps
+            THEN("Two anchors are produced") {
+                REQUIRE(expanded.size() == 1);
+                REQUIRE(expanded.front().size() == 2);
+            }
+            THEN("The area of each anchor is 100mm2") {
+                REQUIRE(area(expanded.front().front()) == Approx(sqr<double>(ten)));
+                REQUIRE(area(expanded.front().back()) == Approx(sqr<double>(ten)));
+            }
+        }
+    }
+
+    GIVEN("two bridges") {
+        Polygon left_support  { { 1 * ten, 1 * ten }, { 2 * ten, 1 * ten }, { 2 * ten, 4 * ten }, { 1 * ten, 4 * ten } };
+        Polygon right_support { { 3 * ten, 1 * ten }, { 4 * ten, 1 * ten }, { 4 * ten, 4 * ten }, { 3 * ten, 4 * ten } };
+        Polygon bottom_bridge { { 2 * ten, 1 * ten }, { 3 * ten, 1 * ten }, { 3 * ten, 2 * ten }, { 2 * ten, 2 * ten } };
+        Polygon top_bridge    { { 2 * ten, 3 * ten }, { 3 * ten, 3 * ten }, { 3 * ten, 4 * ten }, { 2 * ten, 4 * ten } };
+
+        WHEN("expanded") {
+            static constexpr const float expansion = scaled<float>(1.);
+            std::vector<Polygons> expanded = Algorithm::expand_expolygons({ ExPolygon{bottom_bridge}, ExPolygon{top_bridge} }, { ExPolygon{left_support}, ExPolygon{right_support} },
+                expansion,
+                scaled<float>(0.3), // expansion step
+                5);                 // max num steps
+#if 0
+            SVG::export_expolygons(DEBUG_TEMP_DIR "two_bridges-out.svg",
+                { { { { ExPolygon{left_support}, ExPolygon{right_support} } }, { "supports", "orange", 0.5f } },
+                  { { { ExPolygon{bottom_bridge}, ExPolygon{top_bridge} } },   { "bridges",  "blue",   0.5f } },
+                  { { union_ex(union_(expanded.front(), expanded.back())) },   { "expanded", "red",    "black", "", scaled<coord_t>(0.1f), 0.5f } } });
+#endif
+            THEN("Two anchors are produced for each bridge") {
+                REQUIRE(expanded.size() == 2);
+                REQUIRE(expanded.front().size() == 2);
+                REQUIRE(expanded.back().size() == 2);
+            }
+            THEN("The area of each anchor is 10mm2") {
+                double a = expansion * ten + M_PI * sqr(expansion) / 4;
+                double eps = sqr(scaled<double>(0.1));
+                REQUIRE(is_approx(area(expanded.front().front()), a, eps));
+                REQUIRE(is_approx(area(expanded.front().back()), a, eps));
+                REQUIRE(is_approx(area(expanded.back().front()), a, eps));
+                REQUIRE(is_approx(area(expanded.back().back()), a, eps));
+            }
+        }
+    }
+
+    GIVEN("rectangle with rhombic cut-out") {
+        double  diag = 1 * ten * sqrt(2.) / 4.;
+        Polygon square_with_rhombic_cutout{ { 0, 0 }, { 1 * ten, 0 }, { ten / 2, ten / 2 }, { 1 * ten, 1 * ten }, { 0, 1 * ten } };
+        Polygon rhombic { { ten / 2, ten / 2 }, { 3 * ten / 4, ten / 4 }, { 1 * ten, ten / 2 }, { 3 * ten / 4, 3 * ten / 4 } };
+
+        WHEN("expanded") {
+            static constexpr const float expansion = scaled<float>(1.);
+            std::vector<Polygons> expanded = Algorithm::expand_expolygons({ ExPolygon{rhombic} }, { ExPolygon{square_with_rhombic_cutout} },
+                expansion,
+                scaled<float>(0.1), // expansion step
+                11);                // max num steps
+#if 0
+            SVG::export_expolygons(DEBUG_TEMP_DIR "rectangle_with_rhombic_cut-out.svg",
+                { { { { ExPolygon{square_with_rhombic_cutout} } }, { "square_with_rhombic_cutout", "orange", 0.5f } },
+                  { { { ExPolygon{rhombic} } },                    { "rhombic",                    "blue",   0.5f } },
+                  { { union_ex(expanded.front()) },                { "bridges",                    "red",    "black", "", scaled<coord_t>(0.1f), 0.5f } } });
+#endif
+            THEN("Single anchor is produced") {
+                REQUIRE(expanded.size() == 1);
+            }
+            THEN("The area of anchor is correct") {
+                double area_calculated = area(expanded.front());
+                double area_expected = 2. * diag * expansion + M_PI * sqr(expansion) * 0.75;
+                REQUIRE(is_approx(area_expected, area_calculated, sqr(scaled<double>(0.2))));
+            }
+        }
+
+        WHEN("extra expanded") {
+            static constexpr const float expansion = scaled<float>(2.5);
+            std::vector<Polygons> expanded = Algorithm::expand_expolygons({ ExPolygon{rhombic} }, { ExPolygon{square_with_rhombic_cutout} },
+                expansion,
+                scaled<float>(0.25), // expansion step
+                11);                 // max num steps
+#if 0
+            SVG::export_expolygons(DEBUG_TEMP_DIR "rectangle_with_rhombic_cut-out2.svg",
+                { { { { ExPolygon{square_with_rhombic_cutout} } }, { "square_with_rhombic_cutout", "orange", 0.5f } },
+                  { { { ExPolygon{rhombic} } },                    { "rhombic",                    "blue",   0.5f } },
+                  { { union_ex(expanded.front()) },                { "bridges",                    "red",    "black", "", scaled<coord_t>(0.1f), 0.5f } } });
+#endif
+            THEN("Single anchor is produced") {
+                REQUIRE(expanded.size() == 1);
+            }
+            THEN("The area of anchor is correct") {
+                double area_calculated = area(expanded.front());
+                double area_expected = 2. * diag * expansion + M_PI * sqr(expansion) * 0.75;
+                REQUIRE(is_approx(area_expected, area_calculated, sqr(scaled<double>(0.3))));
+            }
+        }
+    }
+
+    GIVEN("square with two holes") {
+        Polygon outer{ { 0, 0 }, { 3 * ten, 0 }, { 3 * ten, 5 * ten }, { 0, 5 * ten } };
+        Polygon hole1{ { 1 * ten, 1 * ten }, { 1 * ten, 2 * ten }, { 2 * ten, 2 * ten }, { 2 * ten, 1 * ten } };
+        Polygon hole2{ { 1 * ten, 3 * ten }, { 1 * ten, 4 * ten }, { 2 * ten, 4 * ten }, { 2 * ten, 3 * ten } };
+        ExPolygon boundary(outer);
+        boundary.holes = { hole1, hole2 };
+
+        Polygon anchor{ { -1 * ten, coord_t(1.5 * ten) }, { 0 * ten, coord_t(1.5 * ten) }, { 0, coord_t(3.5 * ten) }, { -1 * ten, coord_t(3.5 * ten) } };
+
+        WHEN("expanded") {
+            static constexpr const float expansion = scaled<float>(5.);
+            std::vector<Polygons> expanded = Algorithm::expand_expolygons({ ExPolygon{anchor} }, { boundary },
+                expansion,
+                scaled<float>(0.4), // expansion step
+                15);                // max num steps
+#if 0
+            SVG::export_expolygons(DEBUG_TEMP_DIR "square_with_two_holes-out.svg",
+                { { { { ExPolygon{anchor} } }, { "anchor", "orange", 0.5f } },
+                  { { { boundary } },          { "boundary",  "blue", 0.5f } },
+                  { { union_ex(expanded.front()) }, { "expanded", "red", "black", "", scaled<coord_t>(0.1f), 0.5f } } });
+#endif
+            THEN("The anchor expands into a single region") {
+                REQUIRE(expanded.size() == 1);
+                REQUIRE(expanded.front().size() == 1);
+            }
+            THEN("The area of anchor is correct") {
+                double area_calculated = area(expanded.front());
+                double area_expected = double(expansion) * 2. * double(ten) + M_PI * sqr(expansion) * 0.5;
+                REQUIRE(is_approx(area_expected, area_calculated, sqr(scaled<double>(0.45))));
+            }
+        }
+        WHEN("expanded even more") {
+            static constexpr const float expansion = scaled<float>(25.);
+            std::vector<Polygons> expanded = Algorithm::expand_expolygons({ ExPolygon{anchor} }, { boundary },
+                expansion,
+                scaled<float>(2.), // expansion step
+                15);               // max num steps
+#if 0
+            SVG::export_expolygons(DEBUG_TEMP_DIR "square_with_two_holes-expanded2-out.svg",
+                { { { { ExPolygon{anchor} } }, { "anchor", "orange", 0.5f } },
+                  { { { boundary } },          { "boundary",  "blue", 0.5f } },
+                  { { union_ex(expanded.front()) }, { "expanded", "red", "black", "", scaled<coord_t>(0.1f), 0.5f } } });
+#endif
+            THEN("The anchor expands into a single region") {
+                REQUIRE(expanded.size() == 1);
+                REQUIRE(expanded.front().size() == 1);
+            }
+        }
+        WHEN("expanded yet even more") {
+            static constexpr const float expansion = scaled<float>(28.);
+            std::vector<Polygons> expanded = Algorithm::expand_expolygons({ ExPolygon{anchor} }, { boundary },
+                expansion,
+                scaled<float>(2.), // expansion step
+                20);               // max num steps
+#if 0
+            SVG::export_expolygons(DEBUG_TEMP_DIR "square_with_two_holes-expanded3-out.svg",
+                { { { { ExPolygon{anchor} } }, { "anchor", "orange", 0.5f } },
+                  { { { boundary } },          { "boundary",  "blue", 0.5f } },
+                  { { union_ex(expanded.front()) }, { "expanded", "red", "black", "", scaled<coord_t>(0.1f), 0.5f } } });
+#endif
+            THEN("The anchor expands into a single region with two holes") {
+                REQUIRE(expanded.size() == 1);
+                REQUIRE(expanded.front().size() == 3);
+            }
+        }
+        WHEN("expanded fully") {
+            static constexpr const float expansion = scaled<float>(35.);
+            std::vector<Polygons> expanded = Algorithm::expand_expolygons({ ExPolygon{anchor} }, { boundary },
+                expansion,
+                scaled<float>(2.), // expansion step
+                25);               // max num steps
+#if 0
+            SVG::export_expolygons(DEBUG_TEMP_DIR "square_with_two_holes-expanded_fully-out.svg",
+                { { { { ExPolygon{anchor} } }, { "anchor", "orange", 0.5f } },
+                  { { { boundary } },          { "boundary",  "blue", 0.5f } },
+                  { { union_ex(expanded.front()) }, { "expanded", "red", "black", "", scaled<coord_t>(0.1f), 0.5f } } });
+#endif
+            THEN("The anchor expands into a single region with two holes, fully covering the boundary") {
+                REQUIRE(expanded.size() == 1);
+                REQUIRE(expanded.front().size() == 3);
+                REQUIRE(area(expanded.front()) == Approx(area(boundary)));
+            }
+        }
+    }
+}