Merge prusa 2.6.1

src/libslic3r/Arrange/Core/NFP/CircularEdgeIterator.hpp

@@ -0,0 +1,111 @@
+
+#ifndef CIRCULAR_EDGEITERATOR_HPP
+#define CIRCULAR_EDGEITERATOR_HPP
+
+#include <libslic3r/Polygon.hpp>
+#include <libslic3r/Line.hpp>
+
+namespace Slic3r {
+
+// Circular iterator over a polygon yielding individual edges as Line objects
+// if flip_lines is true, the orientation of each line is flipped (not the
+// direction of traversal)
+template<bool flip_lines = false>
+class CircularEdgeIterator_ {
+    const Polygon *m_poly = nullptr;
+    size_t m_i = 0;
+    size_t m_c = 0; // counting how many times the iterator has circled over
+
+public:
+
+    // i: vertex position of first line's starting vertex
+    // poly: target polygon
+    CircularEdgeIterator_(size_t i, const Polygon &poly)
+        : m_poly{&poly}
+        , m_i{!poly.empty() ? i % poly.size() : 0}
+        , m_c{!poly.empty() ? i / poly.size() : 0}
+    {}
+
+    explicit CircularEdgeIterator_ (const Polygon &poly)
+        : CircularEdgeIterator_(0, poly) {}
+
+    using iterator_category = std::forward_iterator_tag;
+    using difference_type   = std::ptrdiff_t;
+    using value_type        = Line;
+    using pointer           = Line*;
+    using reference         = Line&;
+
+    CircularEdgeIterator_ & operator++()
+    {
+        assert (m_poly);
+        ++m_i;
+        if (m_i == m_poly->size()) { // faster than modulo (?)
+            m_i = 0;
+            ++m_c;
+        }
+
+        return *this;
+    }
+
+    CircularEdgeIterator_ operator++(int)
+    {
+        auto cpy = *this; ++(*this); return cpy;
+    }
+
+    Line operator*() const
+    {
+        size_t nx = m_i == m_poly->size() - 1 ? 0 : m_i + 1;
+        Line ret;
+        if constexpr (flip_lines)
+            ret = Line((*m_poly)[nx], (*m_poly)[m_i]);
+        else
+            ret = Line((*m_poly)[m_i], (*m_poly)[nx]);
+
+        return ret;
+    }
+
+    Line operator->() const { return *(*this); }
+
+    bool operator==(const CircularEdgeIterator_& other) const
+    {
+        return m_i == other.m_i && m_c == other.m_c;
+    }
+
+    bool operator!=(const CircularEdgeIterator_& other) const
+    {
+        return !(*this == other);
+    }
+
+    CircularEdgeIterator_& operator +=(size_t dist)
+    {
+        m_i = (m_i + dist) % m_poly->size();
+        m_c = (m_i + (m_c * m_poly->size()) + dist) / m_poly->size();
+
+        return *this;
+    }
+
+    CircularEdgeIterator_ operator +(size_t dist)
+    {
+        auto cpy = *this;
+        cpy += dist;
+
+        return cpy;
+    }
+};
+
+using CircularEdgeIterator = CircularEdgeIterator_<>;
+using CircularReverseEdgeIterator = CircularEdgeIterator_<true>;
+
+inline Range<CircularEdgeIterator> line_range(const Polygon &poly)
+{
+    return Range{CircularEdgeIterator{0, poly}, CircularEdgeIterator{poly.size(), poly}};
+}
+
+inline Range<CircularReverseEdgeIterator> line_range_flp(const Polygon &poly)
+{
+    return Range{CircularReverseEdgeIterator{0, poly}, CircularReverseEdgeIterator{poly.size(), poly}};
+}
+
+} // namespace Slic3r
+
+#endif // CIRCULAR_EDGEITERATOR_HPP

src/libslic3r/Arrange/Core/NFP/EdgeCache.cpp

@@ -0,0 +1,100 @@
+
+#include "EdgeCache.hpp"
+#include "CircularEdgeIterator.hpp"
+
+namespace Slic3r { namespace arr2 {
+
+void EdgeCache::create_cache(const ExPolygon &sh)
+{
+    m_contour.distances.reserve(sh.contour.size());
+    m_holes.reserve(sh.holes.size());
+
+    m_contour.poly = &sh.contour;
+
+    fill_distances(sh.contour, m_contour.distances);
+
+    for (const Polygon &hole : sh.holes) {
+        auto &hc = m_holes.emplace_back();
+        hc.poly = &hole;
+        fill_distances(hole, hc.distances);
+    }
+}
+
+Vec2crd EdgeCache::coords(const ContourCache &cache, double distance) const
+{
+    assert(cache.poly);
+    return arr2::coords(*cache.poly, cache.distances, distance);
+}
+
+void EdgeCache::sample_contour(double accuracy, std::vector<ContourLocation> &samples)
+{
+    const auto N = m_contour.distances.size();
+    const auto S = stride(N, accuracy);
+
+    if (N == 0 || S == 0)
+        return;
+
+    samples.reserve(N / S + 1);
+    for(size_t i = 0; i < N; i += S) {
+        samples.emplace_back(
+            ContourLocation{0, m_contour.distances[i] / m_contour.distances.back()});
+    }
+
+    for (size_t hidx = 1; hidx <= m_holes.size(); ++hidx) {
+        auto& hc = m_holes[hidx - 1];
+
+        const auto NH = hc.distances.size();
+        const auto SH = stride(NH, accuracy);
+
+        if (NH == 0 || SH == 0)
+            continue;
+
+        samples.reserve(samples.size() + NH / SH + 1);
+        for (size_t i = 0; i < NH; i += SH) {
+            samples.emplace_back(
+                ContourLocation{hidx, hc.distances[i] / hc.distances.back()});
+        }
+    }
+}
+
+Vec2crd coords(const Polygon &poly, const std::vector<double> &distances, double distance)
+{
+    assert(poly.size() > 1 && distance >= .0 && distance <= 1.0);
+
+    // distance is from 0.0 to 1.0, we scale it up to the full length of
+    // the circumference
+    double d = distance * distances.back();
+
+    // Magic: we find the right edge in log time
+    auto it = std::lower_bound(distances.begin(), distances.end(), d);
+
+    assert(it != distances.end());
+
+    auto idx  = it - distances.begin(); // get the index of the edge
+    auto &pts = poly.points;
+    auto edge = idx == long(pts.size() - 1) ? Line(pts.back(), pts.front()) :
+                                              Line(pts[idx], pts[idx + 1]);
+
+    // Get the remaining distance on the target edge
+    auto ed = d - (idx > 0 ? *std::prev(it) : 0 );
+
+    double t = ed / edge.length();
+    Vec2d n {double(edge.b.x()) - edge.a.x(), double(edge.b.y()) - edge.a.y()};
+    Vec2crd ret = (edge.a.cast<double>() + t * n).cast<coord_t>();
+
+    return ret;
+}
+
+void fill_distances(const Polygon &poly, std::vector<double> &distances)
+{
+    distances.reserve(poly.size());
+
+    double dist = 0.;
+    auto lrange = line_range(poly);
+    for (const Line &l : lrange) {
+        dist += l.length();
+        distances.emplace_back(dist);
+    }
+}
+
+}} // namespace Slic3r::arr2

src/libslic3r/Arrange/Core/NFP/EdgeCache.hpp

@@ -0,0 +1,72 @@
+
+#ifndef EDGECACHE_HPP
+#define EDGECACHE_HPP
+
+#include <vector>
+
+#include <libslic3r/ExPolygon.hpp>
+
+namespace Slic3r { namespace arr2 {
+
+// Position on the circumference of an ExPolygon.
+// countour_id:  0th is contour, 1..N are holes
+// dist: position given as a floating point number within <0., 1.>
+struct ContourLocation { size_t contour_id; double dist; };
+
+void fill_distances(const Polygon &poly, std::vector<double> &distances);
+
+Vec2crd coords(const Polygon &poly, const std::vector<double>& distances, double distance);
+
+// A class for getting a point on the circumference of the polygon (in log time)
+//
+// This is a transformation of the provided polygon to be able to pinpoint
+// locations on the circumference. The optimizer will pass a floating point
+// value e.g. within <0,1> and we have to transform this value quickly into a
+// coordinate on the circumference. By definition 0 should yield the first
+// vertex and 1.0 would be the last (which should coincide with first).
+//
+// We also have to make this work for the holes of the captured polygon.
+class EdgeCache {
+    struct ContourCache {
+        const Polygon *poly;
+        std::vector<double> distances;
+    } m_contour;
+
+    std::vector<ContourCache> m_holes;
+
+    void create_cache(const ExPolygon& sh);
+
+    Vec2crd coords(const ContourCache& cache, double distance) const;
+
+public:
+
+    explicit EdgeCache(const ExPolygon *sh)
+    {
+        create_cache(*sh);
+    }
+
+    // Given coeff for accuracy <0., 1.>, return the number of vertices to skip
+    // when fetching corners.
+    static inline size_t stride(const size_t N, double accuracy)
+    {
+        size_t n = std::max(size_t{1}, N);
+        return static_cast<coord_t>(
+            std::round(N / std::pow(n, std::pow(accuracy, 1./3.)))
+            );
+    }
+
+    void sample_contour(double accuracy, std::vector<ContourLocation> &samples);
+
+    Vec2crd coords(const ContourLocation &loc) const
+    {
+        assert(loc.contour_id <= m_holes.size());
+
+        return loc.contour_id > 0 ?
+                   coords(m_holes[loc.contour_id - 1], loc.dist) :
+                   coords(m_contour, loc.dist);
+    }
+};
+
+}} // namespace Slic3r::arr2
+
+#endif // EDGECACHE_HPP

src/libslic3r/Arrange/Core/NFP/Kernels/CompactifyKernel.hpp

@@ -0,0 +1,62 @@
+
+#ifndef COMPACTIFYKERNEL_HPP
+#define COMPACTIFYKERNEL_HPP
+
+#include <numeric>
+
+#include "libslic3r/Arrange/Core/NFP/NFPArrangeItemTraits.hpp"
+#include "libslic3r/Arrange/Core/Beds.hpp"
+
+#include <libslic3r/Geometry/ConvexHull.hpp>
+#include <libslic3r/ClipperUtils.hpp>
+
+#include "KernelUtils.hpp"
+
+namespace Slic3r { namespace arr2 {
+
+struct CompactifyKernel {
+    ExPolygons merged_pile;
+
+    template<class ArrItem>
+    double placement_fitness(const ArrItem &itm, const Vec2crd &transl) const
+    {
+        auto pile = merged_pile;
+
+        ExPolygons itm_tr = to_expolygons(envelope_outline(itm));
+        for (auto &p : itm_tr)
+            p.translate(transl);
+
+        append(pile, std::move(itm_tr));
+
+        pile = union_ex(pile);
+
+        Polygon chull = Geometry::convex_hull(pile);
+
+        return -(chull.area());
+    }
+
+    template<class ArrItem, class Bed, class Context, class RemIt>
+    bool on_start_packing(ArrItem &itm,
+                          const Bed &bed,
+                          const Context &packing_context,
+                          const Range<RemIt> & /*remaining_items*/)
+    {
+        bool ret = find_initial_position(itm, bounding_box(bed).center(), bed,
+                                         packing_context);
+
+        merged_pile.clear();
+        for (const auto &gitm : all_items_range(packing_context)) {
+            append(merged_pile, to_expolygons(fixed_outline(gitm)));
+        }
+        merged_pile = union_ex(merged_pile);
+
+        return ret;
+    }
+
+    template<class ArrItem>
+    bool on_item_packed(ArrItem &itm) { return true; }
+};
+
+}} // namespace Slic3r::arr2
+
+#endif // COMPACTIFYKERNEL_HPP

src/libslic3r/Arrange/Core/NFP/Kernels/GravityKernel.hpp

@@ -0,0 +1,59 @@
+
+#ifndef GRAVITYKERNEL_HPP
+#define GRAVITYKERNEL_HPP
+
+#include "libslic3r/Arrange/Core/NFP/NFPArrangeItemTraits.hpp"
+#include "libslic3r/Arrange/Core/Beds.hpp"
+
+#include "KernelUtils.hpp"
+
+namespace Slic3r { namespace arr2 {
+
+struct GravityKernel {
+    std::optional<Vec2crd> sink;
+    std::optional<Vec2crd> item_sink;
+    Vec2d active_sink;
+
+    GravityKernel(Vec2crd gravity_center) : sink{gravity_center} {}
+    GravityKernel() = default;
+
+    template<class ArrItem>
+    double placement_fitness(const ArrItem &itm, const Vec2crd &transl) const
+    {
+        Vec2d center = unscaled(envelope_centroid(itm));
+
+        center += unscaled(transl);
+
+        return - (center - active_sink).squaredNorm();
+    }
+
+    template<class ArrItem, class Bed, class Ctx, class RemIt>
+    bool on_start_packing(ArrItem &itm,
+                          const Bed &bed,
+                          const Ctx &packing_context,
+                          const Range<RemIt> & /*remaining_items*/)
+    {
+        bool ret = false;
+
+        item_sink = get_gravity_sink(itm);
+
+        if (!sink) {
+            sink = bounding_box(bed).center();
+        }
+
+        if (item_sink)
+            active_sink = unscaled(*item_sink);
+        else
+            active_sink = unscaled(*sink);
+
+        ret = find_initial_position(itm, scaled(active_sink), bed, packing_context);
+
+        return ret;
+    }
+
+    template<class ArrItem> bool on_item_packed(ArrItem &itm) { return true; }
+};
+
+}} // namespace Slic3r::arr2
+
+#endif // GRAVITYKERNEL_HPP

src/libslic3r/Arrange/Core/NFP/Kernels/KernelTraits.hpp

@@ -0,0 +1,58 @@
+
+#ifndef KERNELTRAITS_HPP
+#define KERNELTRAITS_HPP
+
+#include "libslic3r/Arrange/Core/ArrangeItemTraits.hpp"
+
+namespace Slic3r { namespace arr2 {
+
+// An arrangement kernel that specifies the object function to the arrangement
+// optimizer and additional callback functions to be able to track the state
+// of the arranged pile during arrangement.
+template<class Kernel, class En = void> struct KernelTraits_
+{
+    // Has to return a score value marking the quality of the arrangement. The
+    // higher this value is, the better a particular placement of the item is.
+    // parameter transl is the translation needed for the item to be moved to
+    // the candidate position.
+    // To discard the item, return NaN as score for every translation.
+    template<class ArrItem>
+    static double placement_fitness(const Kernel  &k,
+                                    const ArrItem &itm,
+                                    const Vec2crd &transl)
+    {
+        return k.placement_fitness(itm, transl);
+    }
+
+    // Called whenever a new item is about to be processed by the optimizer.
+    // The current state of the arrangement can be saved by the kernel: the
+    // already placed items and the remaining items that need to fit into a
+    // particular bed.
+    // Returns true if the item is can be packed immediately, false if it
+    // should be processed further. This way, a kernel have the power to
+    // choose an initial position for the item that is not on the NFP.
+    template<class ArrItem, class Bed, class Ctx, class RemIt>
+    static bool on_start_packing(Kernel &k,
+                                 ArrItem &itm,
+                                 const Bed &bed,
+                                 const Ctx &packing_context,
+                                 const Range<RemIt> &remaining_items)
+    {
+        return k.on_start_packing(itm, bed, packing_context, remaining_items);
+    }
+
+    // Called when an item has been succesfully packed. itm should have the
+    // final translation and rotation already set.
+    // Can return false to discard the item after the optimization.
+    template<class ArrItem>
+    static bool on_item_packed(Kernel &k, ArrItem &itm)
+    {
+        return k.on_item_packed(itm);
+    }
+};
+
+template<class K> using KernelTraits = KernelTraits_<StripCVRef<K>>;
+
+}} // namespace Slic3r::arr2
+
+#endif // KERNELTRAITS_HPP

src/libslic3r/Arrange/Core/NFP/Kernels/KernelUtils.hpp

@@ -0,0 +1,76 @@
+
+#ifndef ARRANGEKERNELUTILS_HPP
+#define ARRANGEKERNELUTILS_HPP
+
+#include <type_traits>
+
+#include "libslic3r/Arrange/Core/NFP/NFPArrangeItemTraits.hpp"
+#include "libslic3r/Arrange/Core/Beds.hpp"
+#include "libslic3r/Arrange/Core/DataStoreTraits.hpp"
+
+namespace Slic3r { namespace arr2 {
+
+template<class Itm, class Bed, class Context>
+bool find_initial_position(Itm &itm,
+                           const Vec2crd &sink,
+                           const Bed &bed,
+                           const Context &packing_context)
+{
+    bool ret = false;
+
+    if constexpr (std::is_convertible_v<Bed, RectangleBed> ||
+                  std::is_convertible_v<Bed, InfiniteBed> ||
+                  std::is_convertible_v<Bed, CircleBed>)
+    {
+        if (all_items_range(packing_context).empty()) {
+            auto rotations = allowed_rotations(itm);
+            auto chull     = envelope_convex_hull(itm);
+
+            for (double rot : rotations) {
+                auto chullcpy = chull;
+                chullcpy.rotate(rot);
+                auto bbitm = bounding_box(chullcpy);
+
+                Vec2crd cb = sink;
+                Vec2crd ci = bbitm.center();
+
+                Vec2crd d = cb - ci;
+                bbitm.translate(d);
+
+                if (bounding_box(bed).contains(bbitm)) {
+                    rotate(itm, rot);
+                    translate(itm, d);
+                    ret = true;
+                    break;
+                }
+            }
+        }
+    }
+
+    return ret;
+}
+
+template<class ArrItem> std::optional<Vec2crd> get_gravity_sink(const ArrItem &itm)
+{
+    constexpr const char * SinkKey = "sink";
+
+    std::optional<Vec2crd> ret;
+
+    auto ptr = get_data<Vec2crd>(itm, SinkKey);
+
+    if (ptr)
+        ret = *ptr;
+
+    return ret;
+}
+
+template<class ArrItem> bool is_wipe_tower(const ArrItem &itm)
+{
+    constexpr const char * Key = "is_wipe_tower";
+
+    return has_key(itm, Key);
+}
+
+}} // namespace Slic3r::arr2
+
+#endif // ARRANGEKERNELUTILS_HPP

src/libslic3r/Arrange/Core/NFP/Kernels/RectangleOverfitKernelWrapper.hpp

@@ -0,0 +1,95 @@
+
+#ifndef RECTANGLEOVERFITKERNELWRAPPER_HPP
+#define RECTANGLEOVERFITKERNELWRAPPER_HPP
+
+#include "KernelTraits.hpp"
+
+#include "libslic3r/Arrange/Core/NFP/NFPArrangeItemTraits.hpp"
+#include "libslic3r/Arrange/Core/Beds.hpp"
+
+namespace Slic3r { namespace arr2 {
+
+// This is a kernel wrapper that will apply a penality to the object function
+// if the result cannot fit into the given rectangular bounds. This can be used
+// to arrange into rectangular boundaries without calculating the IFP of the
+// rectangle bed. Note that after the arrangement, what is garanteed is that
+// the resulting pile will fit into the rectangular boundaries, but it will not
+// be within the given rectangle. The items need to be moved afterwards manually.
+// Use RectangeOverfitPackingStrategy to automate this post process step.
+template<class Kernel>
+struct RectangleOverfitKernelWrapper {
+    Kernel &k;
+    BoundingBox binbb;
+    BoundingBox pilebb;
+
+    RectangleOverfitKernelWrapper(Kernel &kern, const BoundingBox &limits)
+        : k{kern}
+          , binbb{limits}
+    {}
+
+    double overfit(const BoundingBox &itmbb) const
+    {
+        auto fullbb = pilebb;
+        fullbb.merge(itmbb);
+        auto fullbbsz = fullbb.size();
+        auto binbbsz  = binbb.size();
+
+        auto wdiff = fullbbsz.x() - binbbsz.x() - SCALED_EPSILON;
+        auto hdiff = fullbbsz.y() - binbbsz.y() - SCALED_EPSILON;
+        double miss = .0;
+        if (wdiff > 0)
+            miss += double(wdiff);
+        if (hdiff > 0)
+            miss += double(hdiff);
+
+        miss = miss > 0? miss : 0;
+
+        return miss;
+    }
+
+    template<class ArrItem>
+    double placement_fitness(const ArrItem &item, const Vec2crd &transl) const
+    {
+        double score = KernelTraits<Kernel>::placement_fitness(k, item, transl);
+
+        auto itmbb = envelope_bounding_box(item);
+        itmbb.translate(transl);
+        double miss = overfit(itmbb);
+        score -= miss * miss;
+
+        return score;
+    }
+
+    template<class ArrItem, class Bed, class Ctx, class RemIt>
+    bool on_start_packing(ArrItem &itm,
+                          const Bed &bed,
+                          const Ctx &packing_context,
+                          const Range<RemIt> &remaining_items)
+    {
+        pilebb = BoundingBox{};
+
+        for (auto &fitm : all_items_range(packing_context))
+            pilebb.merge(fixed_bounding_box(fitm));
+
+        return KernelTraits<Kernel>::on_start_packing(k, itm, RectangleBed{binbb},
+                                                      packing_context,
+                                                      remaining_items);
+    }
+
+    template<class ArrItem>
+    bool on_item_packed(ArrItem &itm)
+    {
+        bool ret = KernelTraits<Kernel>::on_item_packed(k, itm);
+
+        double miss = overfit(envelope_bounding_box(itm));
+
+        if (miss > 0.)
+            ret = false;
+
+        return ret;
+    }
+};
+
+}} // namespace Slic3r::arr2
+
+#endif // RECTANGLEOVERFITKERNELWRAPPER_H

src/libslic3r/Arrange/Core/NFP/Kernels/SVGDebugOutputKernelWrapper.hpp

@@ -0,0 +1,97 @@
+
+#ifndef SVGDEBUGOUTPUTKERNELWRAPPER_HPP
+#define SVGDEBUGOUTPUTKERNELWRAPPER_HPP
+
+#include <memory>
+
+#include "KernelTraits.hpp"
+
+#include "libslic3r/Arrange/Core/PackingContext.hpp"
+#include "libslic3r/Arrange/Core/NFP/NFPArrangeItemTraits.hpp"
+#include "libslic3r/Arrange/Core/Beds.hpp"
+
+#include <libslic3r/SVG.hpp>
+
+namespace Slic3r { namespace arr2 {
+
+template<class Kernel>
+struct SVGDebugOutputKernelWrapper {
+    Kernel &k;
+    std::unique_ptr<Slic3r::SVG> svg;
+    BoundingBox drawbounds;
+
+    template<class... Args>
+    SVGDebugOutputKernelWrapper(const BoundingBox &bounds, Kernel &kern)
+        : k{kern}, drawbounds{bounds}
+    {}
+
+    template<class ArrItem, class Bed, class Context, class RemIt>
+    bool on_start_packing(ArrItem &itm,
+                          const Bed &bed,
+                          const Context &packing_context,
+                          const Range<RemIt> &rem)
+    {
+        using namespace Slic3r;
+
+        bool ret = KernelTraits<Kernel>::on_start_packing(k, itm, bed,
+                                                          packing_context,
+                                                          rem);
+
+        if (arr2::get_bed_index(itm) < 0)
+            return ret;
+
+        svg.reset();
+        auto bounds = drawbounds;
+        auto fixed = all_items_range(packing_context);
+        svg = std::make_unique<SVG>(std::string("arrange_bed") +
+                                        std::to_string(
+                                            arr2::get_bed_index(itm)) +
+                                        "_" + std::to_string(fixed.size()) +
+                                        ".svg",
+                                    bounds, 0, false);
+
+        svg->draw(ExPolygon{arr2::to_rectangle(drawbounds)}, "blue", .2f);
+
+        auto nfp = calculate_nfp(itm, packing_context, bed);
+        svg->draw_outline(nfp);
+        svg->draw(nfp, "green", 0.2f);
+
+        for (const auto &fixeditm : fixed) {
+            ExPolygons fixeditm_outline = to_expolygons(fixed_outline(fixeditm));
+            svg->draw_outline(fixeditm_outline);
+            svg->draw(fixeditm_outline, "yellow", 0.5f);
+        }
+
+        return ret;
+    }
+
+    template<class ArrItem>
+    double placement_fitness(const ArrItem &item, const Vec2crd &transl) const
+    {
+        return KernelTraits<Kernel>::placement_fitness(k, item, transl);
+    }
+
+    template<class ArrItem>
+    bool on_item_packed(ArrItem &itm)
+    {
+        using namespace Slic3r;
+        using namespace Slic3r::arr2;
+
+        bool ret = KernelTraits<Kernel>::on_item_packed(k, itm);
+
+        if (svg) {
+            ExPolygons itm_outline = to_expolygons(fixed_outline(itm));
+
+            svg->draw_outline(itm_outline);
+            svg->draw(itm_outline, "grey");
+
+            svg->Close();
+        }
+
+        return ret;
+    }
+};
+
+}} // namespace Slic3r::arr2
+
+#endif // SVGDEBUGOUTPUTKERNELWRAPPER_HPP

src/libslic3r/Arrange/Core/NFP/Kernels/TMArrangeKernel.hpp

@@ -0,0 +1,271 @@
+
+#ifndef TMARRANGEKERNEL_HPP
+#define TMARRANGEKERNEL_HPP
+
+#include "libslic3r/Arrange/Core/NFP/NFPArrangeItemTraits.hpp"
+#include "libslic3r/Arrange/Core/Beds.hpp"
+
+#include "KernelUtils.hpp"
+
+#include <boost/geometry/index/rtree.hpp>
+#include <libslic3r/BoostAdapter.hpp>
+
+namespace Slic3r { namespace arr2 {
+
+// Summon the spatial indexing facilities from boost
+namespace bgi = boost::geometry::index;
+using SpatElement = std::pair<BoundingBox, unsigned>;
+using SpatIndex = bgi::rtree<SpatElement, bgi::rstar<16, 4> >;
+
+class TMArrangeKernel {
+    SpatIndex   m_rtree;        // spatial index for the normal (bigger) objects
+    SpatIndex   m_smallsrtree;  // spatial index for only the smaller items
+    BoundingBox m_pilebb;
+    double      m_bin_area = NaNd;
+    double      m_norm;
+    size_t      m_rem_cnt = 0;
+    size_t      m_item_cnt = 0;
+
+
+    struct ItemStats { double area = 0.; BoundingBox bb; };
+    std::vector<ItemStats> m_itemstats;
+
+    // A coefficient used in separating bigger items and smaller items.
+    static constexpr double BigItemTreshold = 0.02;
+
+    template<class T> ArithmeticOnly<T, double> norm(T val) const
+    {
+        return double(val) / m_norm;
+    }
+
+    // Treat big items (compared to the print bed) differently
+    bool is_big(double a) const { return a / m_bin_area > BigItemTreshold; }
+
+protected:
+    std::optional<Point> sink;
+    std::optional<Point> item_sink;
+    Point                active_sink;
+
+    const BoundingBox & pilebb() const { return m_pilebb; }
+
+public:
+    TMArrangeKernel() = default;
+    TMArrangeKernel(Vec2crd gravity_center, size_t itm_cnt, double bedarea = NaNd)
+        : sink{gravity_center}
+        , m_bin_area(bedarea)
+        , m_item_cnt{itm_cnt}
+    {}
+
+    TMArrangeKernel(size_t itm_cnt, double bedarea = NaNd)
+        : m_bin_area(bedarea), m_item_cnt{itm_cnt}
+    {}
+
+    template<class ArrItem>
+    double placement_fitness(const ArrItem &item, const Vec2crd &transl) const
+    {
+        // Candidate item bounding box
+        auto ibb = envelope_bounding_box(item);
+        ibb.translate(transl);
+        auto itmcntr = envelope_centroid(item);
+        itmcntr += transl;
+
+        // Calculate the full bounding box of the pile with the candidate item
+        auto fullbb = m_pilebb;
+        fullbb.merge(ibb);
+
+        // The bounding box of the big items (they will accumulate in the center
+        // of the pile
+        BoundingBox bigbb;
+        if(m_rtree.empty()) {
+            bigbb = fullbb;
+        }
+        else {
+            auto boostbb = m_rtree.bounds();
+            boost::geometry::convert(boostbb, bigbb);
+        }
+
+        // Will hold the resulting score
+        double score = 0;
+
+        // Density is the pack density: how big is the arranged pile
+        double density = 0;
+
+        // Distinction of cases for the arrangement scene
+        enum e_cases {
+            // This branch is for big items in a mixed (big and small) scene
+            // OR for all items in a small-only scene.
+            BIG_ITEM,
+
+            // This branch is for the last big item in a mixed scene
+            LAST_BIG_ITEM,
+
+            // For small items in a mixed scene.
+            SMALL_ITEM,
+
+            WIPE_TOWER,
+        } compute_case;
+
+        bool is_wt = is_wipe_tower(item);
+        bool bigitems = is_big(envelope_area(item)) || m_rtree.empty();
+        if (is_wt)
+            compute_case = WIPE_TOWER;
+        else if (bigitems && m_rem_cnt > 0)
+            compute_case = BIG_ITEM;
+        else if (bigitems && m_rem_cnt == 0)
+            compute_case = LAST_BIG_ITEM;
+        else
+            compute_case = SMALL_ITEM;
+
+        switch (compute_case) {
+        case WIPE_TOWER: {
+            score = (unscaled(itmcntr) - unscaled(active_sink)).squaredNorm();
+            break;
+        }
+        case BIG_ITEM: {
+            const Point& minc = ibb.min; // bottom left corner
+            const Point& maxc = ibb.max; // top right corner
+
+            // top left and bottom right corners
+            Point top_left{minc.x(), maxc.y()};
+            Point bottom_right{maxc.x(), minc.y()};
+
+            // Now the distance of the gravity center will be calculated to the
+            // five anchor points and the smallest will be chosen.
+            std::array<double, 5> dists;
+            auto cc = fullbb.center(); // The gravity center
+            dists[0] = (minc - cc).cast<double>().norm();
+            dists[1] = (maxc - cc).cast<double>().norm();
+            dists[2] = (itmcntr - cc).template cast<double>().norm();
+            dists[3] = (top_left - cc).cast<double>().norm();
+            dists[4] = (bottom_right - cc).cast<double>().norm();
+
+            // The smalles distance from the arranged pile center:
+            double dist = norm(*(std::min_element(dists.begin(), dists.end())));
+            double bindist = norm((ibb.center() - active_sink).template cast<double>().norm());
+            dist = 0.8 * dist + 0.2 * bindist;
+
+            // Prepare a variable for the alignment score.
+            // This will indicate: how well is the candidate item
+            // aligned with its neighbors. We will check the alignment
+            // with all neighbors and return the score for the best
+            // alignment. So it is enough for the candidate to be
+            // aligned with only one item.
+            auto alignment_score = 1.0;
+
+            auto query = bgi::intersects(ibb);
+            auto& index = is_big(envelope_area(item)) ? m_rtree : m_smallsrtree;
+
+            // Query the spatial index for the neighbors
+            std::vector<SpatElement> result;
+            result.reserve(index.size());
+
+            index.query(query, std::back_inserter(result));
+
+            // now get the score for the best alignment
+            for(auto& e : result) {
+                auto idx = e.second;
+                const ItemStats& p = m_itemstats[idx];
+                auto parea = p.area;
+                if(std::abs(1.0 - parea / fixed_area(item)) < 1e-6) {
+                    auto bb = p.bb;
+                    bb.merge(ibb);
+                    auto bbarea = area(bb);
+                    auto ascore = 1.0 - (fixed_area(item) + parea) / bbarea;
+
+                    if(ascore < alignment_score)
+                        alignment_score = ascore;
+                }
+            }
+
+            auto fullbbsz = fullbb.size();
+            density = std::sqrt(norm(fullbbsz.x()) * norm(fullbbsz.y()));
+            double R = double(m_rem_cnt) / (m_item_cnt);
+
+            // The final mix of the score is the balance between the
+            // distance from the full pile center, the pack density and
+            // the alignment with the neighbors
+            if (result.empty())
+                score = 0.50 * dist + 0.50 * density;
+            else
+                // Let the density matter more when fewer objects remain
+                score = 0.50 * dist + (1.0 - R) * 0.20 * density +
+                        0.30 * alignment_score;
+
+            break;
+        }
+        case LAST_BIG_ITEM: {
+            score = norm((itmcntr - m_pilebb.center()).template cast<double>().norm());
+            break;
+        }
+        case SMALL_ITEM: {
+            // Here there are the small items that should be placed around the
+            // already processed bigger items.
+            // No need to play around with the anchor points, the center will be
+            // just fine for small items
+            score = norm((itmcntr - bigbb.center()).template cast<double>().norm());
+            break;
+        }
+        }
+
+        return -score;
+    }
+
+    template<class ArrItem, class Bed, class Context, class RemIt>
+    bool on_start_packing(ArrItem &itm,
+                          const Bed &bed,
+                          const Context &packing_context,
+                          const Range<RemIt> &remaining_items)
+    {
+        item_sink = get_gravity_sink(itm);
+
+        if (!sink) {
+            sink = bounding_box(bed).center();
+        }
+
+        if (item_sink)
+            active_sink = *item_sink;
+        else
+            active_sink = *sink;
+
+        auto fixed = all_items_range(packing_context);
+
+        bool ret = find_initial_position(itm, active_sink, bed, packing_context);
+
+        m_rem_cnt = remaining_items.size();
+
+        if (m_item_cnt == 0)
+            m_item_cnt = m_rem_cnt + fixed.size() + 1;
+
+        if (std::isnan(m_bin_area))
+            m_bin_area = area(bed);
+
+        m_norm = std::sqrt(m_bin_area);
+
+        m_itemstats.clear();
+        m_itemstats.reserve(fixed.size());
+        m_rtree.clear();
+        m_smallsrtree.clear();
+        m_pilebb = {};
+        unsigned idx = 0;
+        for (auto &fixitem : fixed) {
+            auto fixitmbb = fixed_bounding_box(fixitem);
+            m_itemstats.emplace_back(ItemStats{fixed_area(fixitem), fixitmbb});
+            m_pilebb.merge(fixitmbb);
+
+            if(is_big(fixed_area(fixitem)))
+                m_rtree.insert({fixitmbb, idx});
+
+            m_smallsrtree.insert({fixitmbb, idx});
+            idx++;
+        }
+
+        return ret;
+    }
+
+    template<class ArrItem>
+    bool on_item_packed(ArrItem &itm) { return true; }
+};
+
+}} // namespace Slic3r::arr2
+
+#endif // TMARRANGEKERNEL_HPP

src/libslic3r/Arrange/Core/NFP/NFP.cpp

@@ -0,0 +1,419 @@
+
+#ifndef NFP_CPP
+#define NFP_CPP
+
+#include "NFP.hpp"
+#include "CircularEdgeIterator.hpp"
+
+#include "NFPConcave_Tesselate.hpp"
+
+#if !defined(_MSC_VER) && defined(__SIZEOF_INT128__) && !defined(__APPLE__)
+namespace Slic3r { using LargeInt = __int128; }
+#else
+#include <boost/multiprecision/integer.hpp>
+namespace Slic3r { using LargeInt = boost::multiprecision::int128_t; }
+#endif
+
+#include <boost/rational.hpp>
+
+namespace Slic3r {
+
+static bool line_cmp(const Line& e1, const Line& e2)
+{
+    using Ratio = boost::rational<LargeInt>;
+
+    const Vec<2, int64_t> ax(1, 0); // Unit vector for the X axis
+
+    Vec<2, int64_t> p1 = (e1.b - e1.a).cast<int64_t>();
+    Vec<2, int64_t> p2 = (e2.b - e2.a).cast<int64_t>();
+
+    // Quadrant mapping array. The quadrant of a vector can be determined
+    // from the dot product of the vector and its perpendicular pair
+    // with the unit vector X axis. The products will carry the values
+    // lcos = dot(p, ax) = l * cos(phi) and
+    // lsin = -dotperp(p, ax) = l * sin(phi) where
+    // l is the length of vector p. From the signs of these values we can
+    // construct an index which has the sign of lcos as MSB and the
+    // sign of lsin as LSB. This index can be used to retrieve the actual
+    // quadrant where vector p resides using the following map:
+    // (+ is 0, - is 1)
+    // cos | sin | decimal | quadrant
+    //  +  |  +  |    0    |    0
+    //  +  |  -  |    1    |    3
+    //  -  |  +  |    2    |    1
+    //  -  |  -  |    3    |    2
+    std::array<int, 4> quadrants {0, 3, 1, 2 };
+
+    std::array<int, 2> q {0, 0}; // Quadrant indices for p1 and p2
+
+    using TDots = std::array<int64_t, 2>;
+    TDots lcos { p1.dot(ax), p2.dot(ax) };
+    TDots lsin { -dotperp(p1, ax), -dotperp(p2, ax) };
+
+    // Construct the quadrant indices for p1 and p2
+    for(size_t i = 0; i < 2; ++i) {
+        if (lcos[i] == 0)
+            q[i] = lsin[i] > 0 ? 1 : 3;
+        else if (lsin[i] == 0)
+            q[i] = lcos[i] > 0 ? 0 : 2;
+        else
+            q[i] = quadrants[((lcos[i] < 0) << 1) + (lsin[i] < 0)];
+    }
+
+    if (q[0] == q[1]) { // only bother if p1 and p2 are in the same quadrant
+        auto lsq1 = p1.squaredNorm();     // squared magnitudes, avoid sqrt
+        auto lsq2 = p2.squaredNorm();     // squared magnitudes, avoid sqrt
+
+        // We will actually compare l^2 * cos^2(phi) which saturates the
+        // cos function. But with the quadrant info we can get the sign back
+        int sign = q[0] == 1 || q[0] == 2 ? -1 : 1;
+
+        // If Ratio is an actual rational type, there is no precision loss
+        auto pcos1 = Ratio(lcos[0]) / lsq1 * sign * lcos[0];
+        auto pcos2 = Ratio(lcos[1]) / lsq2 * sign * lcos[1];
+
+        return q[0] < 2 ? pcos1 > pcos2 : pcos1 < pcos2;
+    }
+
+    // If in different quadrants, compare the quadrant indices only.
+    return q[0] < q[1];
+}
+
+static inline bool vsort(const Vec2crd& v1, const Vec2crd& v2)
+{
+    return v1.y() == v2.y() ? v1.x() < v2.x() : v1.y() < v2.y();
+}
+
+ExPolygons ifp_convex(const arr2::RectangleBed &obed, const Polygon &convexpoly)
+{
+    ExPolygon ret;
+
+    auto sbox = bounding_box(convexpoly);
+    auto sboxsize = sbox.size();
+    coord_t sheight = sboxsize.y();
+    coord_t swidth  = sboxsize.x();
+    Point   sliding_top = reference_vertex(convexpoly);
+    auto    leftOffset   = sliding_top.x() - sbox.min.x();
+    auto    rightOffset  = sliding_top.x() - sbox.max.x();
+    coord_t topOffset    = 0;
+    auto    bottomOffset = sheight;
+
+    auto bedbb = obed.bb;
+//    bedbb.offset(1);
+    auto bedsz = bedbb.size();
+    auto boxWidth  = bedsz.x();
+    auto boxHeight = bedsz.y();
+
+    auto bedMinx = bedbb.min.x();
+    auto bedMiny = bedbb.min.y();
+    auto bedMaxx = bedbb.max.x();
+    auto bedMaxy = bedbb.max.y();
+
+    Polygon innerNfp{ Point{bedMinx + leftOffset, bedMaxy + topOffset},
+                      Point{bedMaxx + rightOffset, bedMaxy + topOffset},
+                      Point{bedMaxx + rightOffset, bedMiny + bottomOffset},
+                      Point{bedMinx + leftOffset, bedMiny + bottomOffset},
+                      Point{bedMinx + leftOffset, bedMaxy + topOffset} };
+
+    if (sheight <= boxHeight && swidth <= boxWidth)
+        ret.contour = std::move(innerNfp);
+
+    return {ret};
+}
+
+Polygon ifp_convex_convex(const Polygon &fixed, const Polygon &movable)
+{
+    auto subnfps = reserve_polygons(fixed.size());
+
+    // For each edge of the bed polygon, determine the nfp of convexpoly and
+    // the zero area polygon formed by the edge. The union of all these sub-nfps
+    // will contain a hole that is the actual ifp.
+    auto lrange = line_range(fixed);
+    for (const Line &l : lrange) { // Older mac compilers generate warnging if line_range is called in-place
+        Polygon fixed = {l.a, l.b};
+        subnfps.emplace_back(nfp_convex_convex_legacy(fixed, movable));
+    }
+
+    // Do the union and then keep only the holes (should be only one or zero, if
+    // the convexpoly cannot fit into the bed)
+    Polygons ifp = union_(subnfps);
+    Polygon ret;
+
+    // find the first hole
+    auto it = std::find_if(ifp.begin(), ifp.end(), [](const Polygon &subifp){
+        return subifp.is_clockwise();
+    });
+
+    if (it != ifp.end()) {
+        ret = std::move(*it);
+        std::reverse(ret.begin(), ret.end());
+    }
+
+    return ret;
+}
+
+ExPolygons ifp_convex(const arr2::CircleBed &bed, const Polygon &convexpoly)
+{
+    Polygon circle = approximate_circle_with_polygon(bed);
+
+    return {ExPolygon{ifp_convex_convex(circle, convexpoly)}};
+}
+
+ExPolygons ifp_convex(const arr2::IrregularBed &bed, const Polygon &convexpoly)
+{
+    auto bb = get_extents(bed.poly);
+    bb.offset(scaled(1.));
+
+    Polygon rect = arr2::to_rectangle(bb);
+
+    ExPolygons blueprint = diff_ex(rect, bed.poly);
+    Polygons ifp;
+    for (const ExPolygon &part : blueprint) {
+        Polygons triangles = Slic3r::convex_decomposition_tess(part);
+        for (const Polygon &tr : triangles) {
+            Polygon subifp = nfp_convex_convex_legacy(tr, convexpoly);
+            ifp.emplace_back(std::move(subifp));
+        }
+    }
+
+    ifp = union_(ifp);
+
+    Polygons ret;
+
+    std::copy_if(ifp.begin(), ifp.end(), std::back_inserter(ret),
+                 [](const Polygon &p) { return p.is_clockwise(); });
+
+    for (Polygon &p : ret)
+        std::reverse(p.begin(), p.end());
+
+    return to_expolygons(ret);
+}
+
+Vec2crd reference_vertex(const Polygon &poly)
+{
+    Vec2crd ret{std::numeric_limits<coord_t>::min(),
+                std::numeric_limits<coord_t>::min()};
+
+    auto it = std::max_element(poly.points.begin(), poly.points.end(), vsort);
+    if (it != poly.points.end())
+        ret = std::max(ret, static_cast<const Vec2crd &>(*it), vsort);
+
+    return ret;
+}
+
+Vec2crd reference_vertex(const ExPolygon &expoly)
+{
+    return reference_vertex(expoly.contour);
+}
+
+Vec2crd reference_vertex(const Polygons &outline)
+{
+    Vec2crd ret{std::numeric_limits<coord_t>::min(),
+                std::numeric_limits<coord_t>::min()};
+
+    for (const Polygon &poly : outline)
+        ret = std::max(ret, reference_vertex(poly), vsort);
+
+    return ret;
+}
+
+Vec2crd reference_vertex(const ExPolygons &outline)
+{
+    Vec2crd ret{std::numeric_limits<coord_t>::min(),
+                std::numeric_limits<coord_t>::min()};
+
+    for (const ExPolygon &expoly : outline)
+        ret = std::max(ret, reference_vertex(expoly), vsort);
+
+    return ret;
+}
+
+Vec2crd min_vertex(const Polygon &poly)
+{
+    Vec2crd ret{std::numeric_limits<coord_t>::max(),
+                std::numeric_limits<coord_t>::max()};
+
+    auto it = std::min_element(poly.points.begin(), poly.points.end(), vsort);
+    if (it != poly.points.end())
+        ret = std::min(ret, static_cast<const Vec2crd&>(*it), vsort);
+
+    return ret;
+}
+
+// Find the vertex corresponding to the edge with minimum angle to X axis.
+// Only usable with CircularEdgeIterator<> template.
+template<class It> It find_min_anglex_edge(It from)
+{
+    bool found = false;
+    auto it = from;
+    while (!found ) {
+        found = !line_cmp(*it, *std::next(it));
+        ++it;
+    }
+
+    return it;
+}
+
+// Only usable if both fixed and movable polygon is convex. In that case,
+// their edges are already sorted by angle to X axis, only the starting
+// (lowest X axis) edge needs to be found first.
+void nfp_convex_convex(const Polygon &fixed, const Polygon &movable, Polygon &poly)
+{
+    if (fixed.empty() || movable.empty())
+        return;
+
+    // Clear poly and adjust its capacity. Nothing happens if poly is
+    // already sufficiently large and and empty.
+    poly.clear();
+    poly.points.reserve(fixed.size() + movable.size());
+
+    // Find starting positions on the fixed and moving polygons
+    auto it_fx = find_min_anglex_edge(CircularEdgeIterator{fixed});
+    auto it_mv = find_min_anglex_edge(CircularReverseEdgeIterator{movable});
+
+    // End positions are at the same vertex after completing one full circle
+    auto end_fx = it_fx + fixed.size();
+    auto end_mv = it_mv + movable.size();
+
+    // Pos zero is just fine as starting point:
+    poly.points.emplace_back(0, 0);
+
+    // Output iterator adapter for std::merge
+    struct OutItAdaptor {
+        using value_type [[maybe_unused]] = Line;
+        using difference_type [[maybe_unused]] = std::ptrdiff_t;
+        using pointer [[maybe_unused]] = Line*;
+        using reference [[maybe_unused]] = Line& ;
+        using iterator_category [[maybe_unused]] = std::output_iterator_tag;
+
+        Polygon *outpoly;
+        OutItAdaptor(Polygon &out) : outpoly{&out} {}
+
+        OutItAdaptor &operator *() { return *this; }
+        void operator=(const Line &l)
+        {
+            // Yielding l.b, offsetted so that l.a touches the last vertex in
+            // in outpoly
+            outpoly->points.emplace_back(l.b + outpoly->back() - l.a);
+        }
+
+        OutItAdaptor& operator++() { return *this; };
+    };
+
+    // Use std algo to merge the edges from the two polygons
+    std::merge(it_fx, end_fx, it_mv, end_mv, OutItAdaptor{poly}, line_cmp);
+}
+
+Polygon nfp_convex_convex(const Polygon &fixed, const Polygon &movable)
+{
+    Polygon ret;
+    nfp_convex_convex(fixed, movable, ret);
+
+    return ret;
+}
+
+static void buildPolygon(const std::vector<Line>& edgelist,
+                         Polygon& rpoly,
+                         Point& top_nfp)
+{
+    auto& rsh = rpoly.points;
+
+    rsh.reserve(2 * edgelist.size());
+
+    // Add the two vertices from the first edge into the final polygon.
+    rsh.emplace_back(edgelist.front().a);
+    rsh.emplace_back(edgelist.front().b);
+
+    // Sorting function for the nfp reference vertex search
+
+    // the reference (rightmost top) vertex so far
+    top_nfp = *std::max_element(std::cbegin(rsh), std::cend(rsh), vsort);
+
+    auto tmp = std::next(std::begin(rsh));
+
+    // Construct final nfp by placing each edge to the end of the previous
+    for(auto eit = std::next(edgelist.begin()); eit != edgelist.end(); ++eit) {
+        auto d = *tmp - eit->a;
+        Vec2crd p = eit->b + d;
+
+        rsh.emplace_back(p);
+
+        // Set the new reference vertex
+        if (vsort(top_nfp, p))
+            top_nfp = p;
+
+        tmp = std::next(tmp);
+    }
+}
+
+Polygon nfp_convex_convex_legacy(const Polygon &fixed, const Polygon &movable)
+{
+    assert (!fixed.empty());
+    assert (!movable.empty());
+
+    Polygon rsh;   // Final nfp placeholder
+    Point   max_nfp;
+    std::vector<Line> edgelist;
+
+    auto cap = fixed.points.size() + movable.points.size();
+
+    // Reserve the needed memory
+    edgelist.reserve(cap);
+    rsh.points.reserve(cap);
+
+    auto add_edge = [&edgelist](const Point &v1, const Point &v2) {
+        Line e{v1, v2};
+        if ((e.b - e.a).cast<int64_t>().squaredNorm() > 0)
+            edgelist.emplace_back(e);
+    };
+
+    Point max_fixed = fixed.points.front();
+    { // place all edges from fixed into edgelist
+        auto first = std::cbegin(fixed);
+        auto next = std::next(first);
+
+        while(next != std::cend(fixed)) {
+            add_edge(*(first), *(next));
+            max_fixed = std::max(max_fixed, *first, vsort);
+
+            ++first; ++next;
+        }
+
+        add_edge(*std::crbegin(fixed), *std::cbegin(fixed));
+        max_fixed = std::max(max_fixed, *std::crbegin(fixed), vsort);
+    }
+
+    Point max_movable = movable.points.front();
+    Point min_movable = movable.points.front();
+    { // place all edges from movable into edgelist
+        auto first = std::cbegin(movable);
+        auto next = std::next(first);
+
+        while(next != std::cend(movable)) {
+            add_edge(*(next), *(first));
+            min_movable = std::min(min_movable, *first, vsort);
+            max_movable = std::max(max_movable, *first, vsort);
+
+            ++first; ++next;
+        }
+
+        add_edge(*std::cbegin(movable), *std::crbegin(movable));
+        min_movable = std::min(min_movable, *std::crbegin(movable), vsort);
+        max_movable = std::max(max_movable, *std::crbegin(movable), vsort);
+    }
+
+    std::sort(edgelist.begin(), edgelist.end(), line_cmp);
+
+    buildPolygon(edgelist, rsh, max_nfp);
+
+    auto dtouch = max_fixed - min_movable;
+    auto top_other = max_movable + dtouch;
+    auto dnfp = top_other - max_nfp;
+    rsh.translate(dnfp);
+
+    return rsh;
+}
+
+} // namespace Slic3r
+
+#endif // NFP_CPP

src/libslic3r/Arrange/Core/NFP/NFP.hpp

@@ -0,0 +1,51 @@
+
+#ifndef NFP_HPP
+#define NFP_HPP
+
+#include <libslic3r/ExPolygon.hpp>
+#include <libslic3r/Arrange/Core/Beds.hpp>
+
+namespace Slic3r {
+
+template<class Unit = int64_t, class T>
+Unit dotperp(const Vec<2, T> &a, const Vec<2, T> &b)
+{
+    return Unit(a.x()) * Unit(b.y()) - Unit(a.y()) * Unit(b.x());
+}
+
+// Convex-Convex nfp in linear time (fixed.size() + movable.size()),
+// no memory allocations (if out param is used).
+// FIXME: Currently broken for very sharp triangles.
+Polygon nfp_convex_convex(const Polygon &fixed, const Polygon &movable);
+void nfp_convex_convex(const Polygon &fixed, const Polygon &movable, Polygon &out);
+Polygon nfp_convex_convex_legacy(const Polygon &fixed, const Polygon &movable);
+
+Polygon ifp_convex_convex(const Polygon &fixed, const Polygon &movable);
+
+ExPolygons ifp_convex(const arr2::RectangleBed &bed, const Polygon &convexpoly);
+ExPolygons ifp_convex(const arr2::CircleBed &bed, const Polygon &convexpoly);
+ExPolygons ifp_convex(const arr2::IrregularBed &bed, const Polygon &convexpoly);
+inline ExPolygons ifp_convex(const arr2::InfiniteBed &bed, const Polygon &convexpoly)
+{
+    return {};
+}
+
+inline ExPolygons ifp_convex(const arr2::ArrangeBed &bed, const Polygon &convexpoly)
+{
+    ExPolygons ret;
+    auto visitor = [&ret, &convexpoly](const auto &b) { ret = ifp_convex(b, convexpoly); };
+    boost::apply_visitor(visitor, bed);
+
+    return ret;
+}
+
+Vec2crd reference_vertex(const Polygon &outline);
+Vec2crd reference_vertex(const ExPolygon &outline);
+Vec2crd reference_vertex(const Polygons &outline);
+Vec2crd reference_vertex(const ExPolygons &outline);
+
+Vec2crd min_vertex(const Polygon &outline);
+
+} // namespace Slic3r
+
+#endif // NFP_HPP

src/libslic3r/Arrange/Core/NFP/NFPArrangeItemTraits.hpp

@@ -0,0 +1,197 @@
+
+#ifndef NFPARRANGEITEMTRAITS_HPP
+#define NFPARRANGEITEMTRAITS_HPP
+
+#include <numeric>
+
+#include "libslic3r/Arrange/Core/ArrangeBase.hpp"
+
+#include "libslic3r/ExPolygon.hpp"
+#include "libslic3r/BoundingBox.hpp"
+
+namespace Slic3r { namespace arr2 {
+
+// Additional methods that an ArrangeItem object has to implement in order
+// to be usable with PackStrategyNFP.
+template<class ArrItem, class En = void> struct NFPArrangeItemTraits_
+{
+    template<class Context, class Bed, class StopCond = DefaultStopCondition>
+    static ExPolygons calculate_nfp(const ArrItem &item,
+                                    const Context &packing_context,
+                                    const Bed &bed,
+                                    StopCond stop_condition = {})
+    {
+        static_assert(always_false<ArrItem>::value,
+                      "NFP unimplemented for this item type.");
+        return {};
+    }
+
+    static Vec2crd reference_vertex(const ArrItem &item)
+    {
+        return item.reference_vertex();
+    }
+
+    static BoundingBox envelope_bounding_box(const ArrItem &itm)
+    {
+        return itm.envelope_bounding_box();
+    }
+
+    static BoundingBox fixed_bounding_box(const ArrItem &itm)
+    {
+        return itm.fixed_bounding_box();
+    }
+
+    static const Polygons & envelope_outline(const ArrItem &itm)
+    {
+        return itm.envelope_outline();
+    }
+
+    static const Polygons & fixed_outline(const ArrItem &itm)
+    {
+        return itm.fixed_outline();
+    }
+
+    static const Polygon & envelope_convex_hull(const ArrItem &itm)
+    {
+        return itm.envelope_convex_hull();
+    }
+
+    static const Polygon & fixed_convex_hull(const ArrItem &itm)
+    {
+        return itm.fixed_convex_hull();
+    }
+
+    static double envelope_area(const ArrItem &itm)
+    {
+        return itm.envelope_area();
+    }
+
+    static double fixed_area(const ArrItem &itm)
+    {
+        return itm.fixed_area();
+    }
+
+    static auto allowed_rotations(const ArrItem &)
+    {
+        return std::array{0.};
+    }
+
+    static Vec2crd fixed_centroid(const ArrItem &itm)
+    {
+        return fixed_bounding_box(itm).center();
+    }
+
+    static Vec2crd envelope_centroid(const ArrItem &itm)
+    {
+        return envelope_bounding_box(itm).center();
+    }
+};
+
+template<class T>
+using NFPArrangeItemTraits = NFPArrangeItemTraits_<StripCVRef<T>>;
+
+template<class ArrItem,
+         class Context,
+         class Bed,
+         class StopCond = DefaultStopCondition>
+ExPolygons calculate_nfp(const ArrItem &itm,
+                         const Context &context,
+                         const Bed &bed,
+                         StopCond stopcond = {})
+{
+    return NFPArrangeItemTraits<ArrItem>::calculate_nfp(itm, context, bed,
+                                                        std::move(stopcond));
+}
+
+template<class ArrItem> Vec2crd reference_vertex(const ArrItem &itm)
+{
+    return NFPArrangeItemTraits<ArrItem>::reference_vertex(itm);
+}
+
+template<class ArrItem> BoundingBox envelope_bounding_box(const ArrItem &itm)
+{
+    return NFPArrangeItemTraits<ArrItem>::envelope_bounding_box(itm);
+}
+
+template<class ArrItem> BoundingBox fixed_bounding_box(const ArrItem &itm)
+{
+    return NFPArrangeItemTraits<ArrItem>::fixed_bounding_box(itm);
+}
+
+template<class ArrItem> decltype(auto) envelope_convex_hull(const ArrItem &itm)
+{
+    return NFPArrangeItemTraits<ArrItem>::envelope_convex_hull(itm);
+}
+
+template<class ArrItem> decltype(auto) fixed_convex_hull(const ArrItem &itm)
+{
+    return NFPArrangeItemTraits<ArrItem>::fixed_convex_hull(itm);
+}
+
+template<class ArrItem> decltype(auto) envelope_outline(const ArrItem &itm)
+{
+    return NFPArrangeItemTraits<ArrItem>::envelope_outline(itm);
+}
+
+template<class ArrItem> decltype(auto) fixed_outline(const ArrItem &itm)
+{
+    return NFPArrangeItemTraits<ArrItem>::fixed_outline(itm);
+}
+
+template<class ArrItem> double envelope_area(const ArrItem &itm)
+{
+    return NFPArrangeItemTraits<ArrItem>::envelope_area(itm);
+}
+
+template<class ArrItem> double fixed_area(const ArrItem &itm)
+{
+    return NFPArrangeItemTraits<ArrItem>::fixed_area(itm);
+}
+
+template<class ArrItem> Vec2crd fixed_centroid(const ArrItem &itm)
+{
+    return NFPArrangeItemTraits<ArrItem>::fixed_centroid(itm);
+}
+
+template<class ArrItem> Vec2crd envelope_centroid(const ArrItem &itm)
+{
+    return NFPArrangeItemTraits<ArrItem>::envelope_centroid(itm);
+}
+
+template<class ArrItem>
+auto allowed_rotations(const ArrItem &itm)
+{
+    return NFPArrangeItemTraits<ArrItem>::allowed_rotations(itm);
+}
+
+template<class It>
+BoundingBox bounding_box(const Range<It> &itms) noexcept
+{
+    auto pilebb =
+        std::accumulate(itms.begin(), itms.end(), BoundingBox{},
+                        [](BoundingBox bb, const auto &itm) {
+                            bb.merge(fixed_bounding_box(itm));
+                            return bb;
+                        });
+
+    return pilebb;
+}
+
+template<class It>
+BoundingBox bounding_box_on_bedidx(const Range<It> &itms, int bed_index) noexcept
+{
+    auto pilebb =
+        std::accumulate(itms.begin(), itms.end(), BoundingBox{},
+                        [bed_index](BoundingBox bb, const auto &itm) {
+                            if (bed_index == get_bed_index(itm))
+                                bb.merge(fixed_bounding_box(itm));
+
+                            return bb;
+                        });
+
+    return pilebb;
+}
+
+}} // namespace Slic3r::arr2
+
+#endif // ARRANGEITEMTRAITSNFP_HPP

src/libslic3r/Arrange/Core/NFP/NFPConcave_CGAL.cpp

@@ -0,0 +1,112 @@
+
+#include "NFP.hpp"
+#include "NFPConcave_CGAL.hpp"
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/partition_2.h>
+#include <CGAL/Partition_traits_2.h>
+#include <CGAL/property_map.h>
+#include <CGAL/Polygon_vertical_decomposition_2.h>
+
+#include "libslic3r/ClipperUtils.hpp"
+
+namespace Slic3r {
+
+using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+using Partition_traits_2 = CGAL::Partition_traits_2<K, CGAL::Pointer_property_map<K::Point_2>::type >;
+using Point_2 = Partition_traits_2::Point_2;
+using Polygon_2 = Partition_traits_2::Polygon_2;  // a polygon of indices
+
+ExPolygons nfp_concave_concave_cgal(const ExPolygon &fixed, const ExPolygon &movable)
+{
+    Polygons fixed_decomp = convex_decomposition_cgal(fixed);
+    Polygons movable_decomp = convex_decomposition_cgal(movable);
+
+    auto refs_mv = reserve_vector<Vec2crd>(movable_decomp.size());
+
+    for (const Polygon &p : movable_decomp)
+        refs_mv.emplace_back(reference_vertex(p));
+
+    auto nfps = reserve_polygons(fixed_decomp.size() *movable_decomp.size());
+
+    Vec2crd ref_whole = reference_vertex(movable);
+    for (const Polygon &fixed_part : fixed_decomp) {
+        size_t mvi = 0;
+        for (const Polygon &movable_part : movable_decomp) {
+            Polygon subnfp = nfp_convex_convex(fixed_part, movable_part);
+            const Vec2crd &ref_mp = refs_mv[mvi];
+            auto d = ref_whole - ref_mp;
+            subnfp.translate(d);
+            nfps.emplace_back(subnfp);
+            mvi++;
+        }
+    }
+
+    return union_ex(nfps);
+}
+
+// TODO: holes
+Polygons convex_decomposition_cgal(const ExPolygon &expoly)
+{
+    CGAL::Polygon_vertical_decomposition_2<K> decomp;
+
+    CGAL::Polygon_2<K> contour;
+    for (auto &p : expoly.contour.points)
+        contour.push_back({unscaled(p.x()), unscaled(p.y())});
+
+    CGAL::Polygon_with_holes_2<K> cgalpoly{contour};
+    for (const Polygon &h : expoly.holes) {
+        CGAL::Polygon_2<K> hole;
+        for (auto &p : h.points)
+            hole.push_back({unscaled(p.x()), unscaled(p.y())});
+
+        cgalpoly.add_hole(hole);
+    }
+
+    std::vector<CGAL::Polygon_2<K>> out;
+    decomp(cgalpoly, std::back_inserter(out));
+
+    Polygons ret;
+    for (auto &pwh : out) {
+        Polygon poly;
+        for (auto &p : pwh)
+            poly.points.emplace_back(scaled(p.x()), scaled(p.y()));
+        ret.emplace_back(std::move(poly));
+    }
+
+    return ret; //convex_decomposition_cgal(expoly.contour);
+}
+
+Polygons convex_decomposition_cgal(const Polygon &poly)
+{
+    auto pts = reserve_vector<K::Point_2>(poly.size());
+
+    for (const Point &p : poly.points)
+        pts.emplace_back(unscaled(p.x()), unscaled(p.y()));
+
+    Partition_traits_2 traits(CGAL::make_property_map(pts));
+
+    Polygon_2 polyidx;
+    for (size_t i = 0; i < pts.size(); ++i)
+        polyidx.push_back(i);
+
+    std::vector<Polygon_2> outp;
+
+    CGAL::optimal_convex_partition_2(polyidx.vertices_begin(),
+                                     polyidx.vertices_end(),
+                                     std::back_inserter(outp),
+                                     traits);
+
+    Polygons ret;
+    for (const Polygon_2& poly : outp){
+        Polygon r;
+        for(Point_2 p : poly.container())
+            r.points.emplace_back(scaled(pts[p].x()), scaled(pts[p].y()));
+
+        ret.emplace_back(std::move(r));
+    }
+
+    return ret;
+}
+
+} // namespace Slic3r

src/libslic3r/Arrange/Core/NFP/NFPConcave_CGAL.hpp

@@ -0,0 +1,15 @@
+
+#ifndef NFPCONCAVE_CGAL_HPP
+#define NFPCONCAVE_CGAL_HPP
+
+#include <libslic3r/ExPolygon.hpp>
+
+namespace Slic3r {
+
+Polygons convex_decomposition_cgal(const Polygon &expoly);
+Polygons convex_decomposition_cgal(const ExPolygon &expoly);
+ExPolygons nfp_concave_concave_cgal(const ExPolygon &fixed, const ExPolygon &movable);
+
+} // namespace Slic3r
+
+#endif // NFPCONCAVE_CGAL_HPP

src/libslic3r/Arrange/Core/NFP/NFPConcave_Tesselate.cpp

@@ -0,0 +1,71 @@
+
+#include "NFPConcave_Tesselate.hpp"
+
+#include <libslic3r/ClipperUtils.hpp>
+#include <libslic3r/Tesselate.hpp>
+
+#include "NFP.hpp"
+
+namespace Slic3r {
+
+Polygons convex_decomposition_tess(const Polygon &expoly)
+{
+    return convex_decomposition_tess(ExPolygon{expoly});
+}
+
+Polygons convex_decomposition_tess(const ExPolygon &expoly)
+{
+    std::vector<Vec2d> tr = Slic3r::triangulate_expolygon_2d(expoly);
+
+    auto ret = Slic3r::reserve_polygons(tr.size() / 3);
+    for (size_t i = 0; i < tr.size(); i += 3) {
+        ret.emplace_back(
+            Polygon{scaled(tr[i]), scaled(tr[i + 1]), scaled(tr[i + 2])});
+    }
+
+    return ret;
+}
+
+Polygons convex_decomposition_tess(const ExPolygons &expolys)
+{
+    constexpr size_t AvgTriangleCountGuess = 50;
+
+    auto ret = reserve_polygons(AvgTriangleCountGuess * expolys.size());
+    for (const ExPolygon &expoly : expolys) {
+        Polygons convparts = convex_decomposition_tess(expoly);
+        std::move(convparts.begin(), convparts.end(), std::back_inserter(ret));
+    }
+
+    return ret;
+}
+
+ExPolygons nfp_concave_concave_tess(const ExPolygon &fixed,
+                                    const ExPolygon &movable)
+{
+    Polygons fixed_decomp = convex_decomposition_tess(fixed);
+    Polygons movable_decomp = convex_decomposition_tess(movable);
+
+    auto refs_mv = reserve_vector<Vec2crd>(movable_decomp.size());
+
+    for (const Polygon &p : movable_decomp)
+        refs_mv.emplace_back(reference_vertex(p));
+
+    auto nfps = reserve_polygons(fixed_decomp.size() * movable_decomp.size());
+
+    Vec2crd ref_whole = reference_vertex(movable);
+    for (const Polygon &fixed_part : fixed_decomp) {
+        size_t mvi = 0;
+        for (const Polygon &movable_part : movable_decomp) {
+            Polygon subnfp = nfp_convex_convex(fixed_part, movable_part);
+            const Vec2crd &ref_mp = refs_mv[mvi];
+            auto d = ref_whole - ref_mp;
+            subnfp.translate(d);
+            nfps.emplace_back(subnfp);
+            mvi++;
+        }
+    }
+
+    return union_ex(nfps);
+}
+
+} // namespace Slic3r

src/libslic3r/Arrange/Core/NFP/NFPConcave_Tesselate.hpp

@@ -0,0 +1,16 @@
+
+#ifndef NFPCONCAVE_TESSELATE_HPP
+#define NFPCONCAVE_TESSELATE_HPP
+
+#include <libslic3r/ExPolygon.hpp>
+
+namespace Slic3r {
+
+Polygons convex_decomposition_tess(const Polygon &expoly);
+Polygons convex_decomposition_tess(const ExPolygon &expoly);
+Polygons convex_decomposition_tess(const ExPolygons &expolys);
+ExPolygons nfp_concave_concave_tess(const ExPolygon &fixed, const ExPolygon &movable);
+
+} // namespace Slic3r
+
+#endif // NFPCONCAVE_TESSELATE_HPP

src/libslic3r/Arrange/Core/NFP/PackStrategyNFP.hpp

@@ -0,0 +1,286 @@
+
+#ifndef PACKSTRATEGYNFP_HPP
+#define PACKSTRATEGYNFP_HPP
+
+#include "libslic3r/Arrange/Core/ArrangeBase.hpp"
+
+#include "EdgeCache.hpp"
+#include "Kernels/KernelTraits.hpp"
+
+#include "NFPArrangeItemTraits.hpp"
+
+#include "libslic3r/Optimize/NLoptOptimizer.hpp"
+#include "libslic3r/Execution/ExecutionSeq.hpp"
+
+namespace Slic3r { namespace arr2 {
+
+struct NFPPackingTag{};
+
+struct DummyArrangeKernel
+{
+    template<class ArrItem>
+    double placement_fitness(const ArrItem &itm, const Vec2crd &dest_pos) const
+    {
+        return NaNd;
+    }
+
+    template<class ArrItem, class Bed, class Context, class RemIt>
+    bool on_start_packing(ArrItem            &itm,
+                          const Bed          &bed,
+                          const Context      &packing_context,
+                          const Range<RemIt> &remaining_items)
+    {
+        return true;
+    }
+
+    template<class ArrItem> bool on_item_packed(ArrItem &itm) { return true; }
+};
+
+template<class Strategy> using OptAlg = typename Strategy::OptAlg;
+
+template<class ArrangeKernel = DummyArrangeKernel,
+         class ExecPolicy = ExecutionSeq,
+         class OptMethod = opt::AlgNLoptSubplex,
+         class StopCond  = DefaultStopCondition>
+struct PackStrategyNFP {
+    using OptAlg = OptMethod;
+
+    ArrangeKernel kernel;
+    ExecPolicy ep;
+    double accuracy = 1.;
+    opt::Optimizer<OptMethod> solver;
+    StopCond stop_condition;
+
+    PackStrategyNFP(opt::Optimizer<OptMethod> slv,
+                    ArrangeKernel k = {},
+                    ExecPolicy execpolicy = {},
+                    double accur = 1.,
+                    StopCond stop_cond = {})
+        : kernel{std::move(k)},
+          ep{std::move(execpolicy)},
+          accuracy{accur},
+          solver{std::move(slv)},
+          stop_condition{std::move(stop_cond)}
+    {}
+
+    PackStrategyNFP(ArrangeKernel k = {},
+                    ExecPolicy execpolicy = {},
+                    double accur = 1.,
+                    StopCond stop_cond = {})
+        : PackStrategyNFP{opt::Optimizer<OptMethod>{}, std::move(k),
+                          std::move(execpolicy), accur, std::move(stop_cond)}
+    {
+        // Defaults for AlgNLoptSubplex
+        auto iters = static_cast<unsigned>(std::floor(1000 * accuracy));
+        auto optparams =
+            opt::StopCriteria{}.max_iterations(iters).rel_score_diff(
+                1e-20) /*.abs_score_diff(1e-20)*/;
+
+        solver.set_criteria(optparams);
+    }
+};
+
+template<class...Args>
+struct PackStrategyTag_<PackStrategyNFP<Args...>>
+{
+    using Tag = NFPPackingTag;
+};
+
+
+template<class ArrItem, class Bed, class PStrategy>
+double pick_best_spot_on_nfp_verts_only(ArrItem            &item,
+                                        const ExPolygons   &nfp,
+                                        const Bed          &bed,
+                                        const PStrategy    &strategy)
+{
+    using KernelT = KernelTraits<decltype(strategy.kernel)>;
+
+    auto    score   = -std::numeric_limits<double>::infinity();
+    Vec2crd orig_tr = get_translation(item);
+    Vec2crd translation{0, 0};
+
+    auto eval_fitness = [&score, &strategy, &item, &translation,
+                         &orig_tr](const Vec2crd &p) {
+        set_translation(item, orig_tr);
+        Vec2crd ref_v = reference_vertex(item);
+        Vec2crd tr    = p - ref_v;
+        double fitness = KernelT::placement_fitness(strategy.kernel, item, tr);
+        if (fitness > score) {
+            score       = fitness;
+            translation = tr;
+        }
+    };
+
+    for (const ExPolygon &expoly : nfp) {
+        for (const Point &p : expoly.contour) {
+            eval_fitness(p);
+        }
+
+        for (const Polygon &h : expoly.holes)
+            for (const Point &p : h.points)
+                eval_fitness(p);
+    }
+
+    set_translation(item, orig_tr + translation);
+
+    return score;
+}
+
+struct CornerResult
+{
+    size_t         contour_id;
+    opt::Result<1> oresult;
+};
+
+template<class ArrItem, class Bed, class... Args>
+double pick_best_spot_on_nfp(ArrItem                        &item,
+                             const ExPolygons               &nfp,
+                             const Bed                      &bed,
+                             const PackStrategyNFP<Args...> &strategy)
+{
+    auto &ex_policy = strategy.ep;
+    using KernelT = KernelTraits<decltype(strategy.kernel)>;
+
+    auto score = -std::numeric_limits<double>::infinity();
+    Vec2crd orig_tr = get_translation(item);
+    Vec2crd translation{0, 0};
+    Vec2crd ref_v = reference_vertex(item);
+
+    auto edge_caches = reserve_vector<EdgeCache>(nfp.size());
+    auto sample_sets = reserve_vector<std::vector<ContourLocation>>(
+        nfp.size());
+
+    for (const ExPolygon &expoly : nfp) {
+        edge_caches.emplace_back(EdgeCache{&expoly});
+        edge_caches.back().sample_contour(strategy.accuracy,
+                                          sample_sets.emplace_back());
+    }
+
+    auto nthreads = execution::max_concurrency(ex_policy);
+
+    std::vector<CornerResult> gresults(edge_caches.size());
+
+    auto resultcmp = [](auto &a, auto &b) {
+        return a.oresult.score < b.oresult.score;
+    };
+
+    execution::for_each(
+        ex_policy, size_t(0), edge_caches.size(),
+        [&](size_t edge_cache_idx) {
+            auto &ec_contour = edge_caches[edge_cache_idx];
+            auto &corners = sample_sets[edge_cache_idx];
+            std::vector<CornerResult> results(corners.size());
+
+            auto cornerfn = [&](size_t i) {
+                ContourLocation cr = corners[i];
+                auto objfn = [&](opt::Input<1> &in) {
+                    Vec2crd p = ec_contour.coords(ContourLocation{cr.contour_id, in[0]});
+                    Vec2crd tr = p - ref_v;
+
+                    return KernelT::placement_fitness(strategy.kernel, item, tr);
+                };
+
+                // Assuming that solver is a lightweight object
+                auto solver = strategy.solver;
+                solver.to_max();
+                auto oresult = solver.optimize(objfn,
+                                               opt::initvals({cr.dist}),
+                                               opt::bounds({{0., 1.}}));
+
+                results[i] = CornerResult{cr.contour_id, oresult};
+            };
+
+            execution::for_each(ex_policy, size_t(0), results.size(),
+                                cornerfn, nthreads);
+
+            auto it = std::max_element(results.begin(), results.end(),
+                                       resultcmp);
+
+            if (it != results.end())
+                gresults[edge_cache_idx] = *it;
+        },
+        nthreads);
+
+    auto it = std::max_element(gresults.begin(), gresults.end(), resultcmp);
+    if (it != gresults.end()) {
+        score = it->oresult.score;
+        size_t path_id = std::distance(gresults.begin(), it);
+        size_t contour_id = it->contour_id;
+        double dist = it->oresult.optimum[0];
+
+        Vec2crd pos = edge_caches[path_id].coords(ContourLocation{contour_id, dist});
+        Vec2crd tr = pos - ref_v;
+
+        set_translation(item, orig_tr + tr);
+    }
+
+    return score;
+}
+
+template<class Strategy, class ArrItem, class Bed, class RemIt>
+bool pack(Strategy &strategy,
+          const Bed &bed,
+          ArrItem &item,
+          const PackStrategyContext<Strategy, ArrItem> &packing_context,
+          const Range<RemIt> &remaining_items,
+          const NFPPackingTag &)
+{
+    using KernelT = KernelTraits<decltype(strategy.kernel)>;
+
+    // The kernel might pack the item immediately
+    bool packed = KernelT::on_start_packing(strategy.kernel, item, bed,
+                                            packing_context, remaining_items);
+
+    double  orig_rot    = get_rotation(item);
+    double  final_rot   = 0.;
+    double  final_score = -std::numeric_limits<double>::infinity();
+    Vec2crd orig_tr     = get_translation(item);
+    Vec2crd final_tr    = orig_tr;
+
+    bool cancelled = strategy.stop_condition();
+    const auto & rotations = allowed_rotations(item);
+
+    // Check all rotations but only if item is not already packed
+    for (auto rot_it = rotations.begin();
+         !cancelled && !packed && rot_it != rotations.end(); ++rot_it) {
+
+        double rot = *rot_it;
+
+        set_rotation(item, orig_rot + rot);
+        set_translation(item, orig_tr);
+
+        auto nfp = calculate_nfp(item, packing_context, bed,
+                                 strategy.stop_condition);
+        double score = NaNd;
+        if (!nfp.empty()) {
+            score = pick_best_spot_on_nfp(item, nfp, bed, strategy);
+
+            cancelled = strategy.stop_condition();
+            if (score > final_score) {
+                final_score = score;
+                final_rot   = rot;
+                final_tr    = get_translation(item);
+            }
+        }
+    }
+
+    // If the score is not valid, and the item is not already packed, or
+    // the packing was cancelled asynchronously by stop condition, then
+    // discard the packing
+    bool is_score_valid = !std::isnan(final_score) && !std::isinf(final_score);
+    packed = !cancelled && (packed || is_score_valid);
+
+    if (packed) {
+        set_translation(item, final_tr);
+        set_rotation(item, orig_rot + final_rot);
+
+        // Finally, consult the kernel if the packing is sane
+        packed = KernelT::on_item_packed(strategy.kernel, item);
+    }
+
+    return packed;
+}
+
+}} // namespace Slic3r::arr2
+
+#endif // PACKSTRATEGYNFP_HPP

src/libslic3r/Arrange/Core/NFP/RectangleOverfitPackingStrategy.hpp

@@ -0,0 +1,142 @@
+
+#ifndef RECTANGLEOVERFITPACKINGSTRATEGY_HPP
+#define RECTANGLEOVERFITPACKINGSTRATEGY_HPP
+
+#include "Kernels/RectangleOverfitKernelWrapper.hpp"
+
+#include "libslic3r/Arrange/Core/NFP/PackStrategyNFP.hpp"
+#include "libslic3r/Arrange/Core/Beds.hpp"
+
+namespace Slic3r { namespace arr2 {
+
+using PostAlignmentFn = std::function<Vec2crd(const BoundingBox &bedbb,
+                                              const BoundingBox &pilebb)>;
+
+struct CenterAlignmentFn {
+    Vec2crd operator() (const BoundingBox &bedbb,
+                        const BoundingBox &pilebb)
+    {
+        return bedbb.center() - pilebb.center();
+    }
+};
+
+template<class ArrItem>
+struct RectangleOverfitPackingContext : public DefaultPackingContext<ArrItem>
+{
+    BoundingBox limits;
+    int bed_index;
+    PostAlignmentFn post_alignment_fn;
+
+    explicit RectangleOverfitPackingContext(const BoundingBox limits,
+                     int bedidx,
+                     PostAlignmentFn alignfn = CenterAlignmentFn{})
+        : limits{limits}, bed_index{bedidx}, post_alignment_fn{alignfn}
+    {}
+
+    void align_pile()
+    {
+        // Here, the post alignment can be safely done. No throwing
+        // functions are called!
+        if (fixed_items_range(*this).empty()) {
+            auto itms = packed_items_range(*this);
+            auto pilebb = bounding_box(itms);
+
+            for (auto &itm : itms) {
+                translate(itm, post_alignment_fn(limits, pilebb));
+            }
+        }
+    }
+
+    ~RectangleOverfitPackingContext() { align_pile(); }
+};
+
+// With rectange bed, and no fixed items, an infinite bed with
+// RectangleOverfitKernelWrapper can produce better results than a pure
+// RectangleBed with inner-fit polygon calculation.
+template<class ...Args>
+struct RectangleOverfitPackingStrategy {
+    PackStrategyNFP<Args...> base_strategy;
+
+    PostAlignmentFn post_alignment_fn = CenterAlignmentFn{};
+
+    template<class ArrItem>
+    using Context = RectangleOverfitPackingContext<ArrItem>;
+
+    RectangleOverfitPackingStrategy(PackStrategyNFP<Args...> s,
+                                    PostAlignmentFn post_align_fn)
+        : base_strategy{std::move(s)}, post_alignment_fn{post_align_fn}
+    {}
+
+    RectangleOverfitPackingStrategy(PackStrategyNFP<Args...> s)
+        : base_strategy{std::move(s)}
+    {}
+};
+
+struct RectangleOverfitPackingStrategyTag {};
+
+template<class... Args>
+struct PackStrategyTag_<RectangleOverfitPackingStrategy<Args...>> {
+    using Tag = RectangleOverfitPackingStrategyTag;
+};
+
+template<class... Args>
+struct PackStrategyTraits_<RectangleOverfitPackingStrategy<Args...>> {
+    template<class ArrItem>
+    using Context = typename RectangleOverfitPackingStrategy<
+        Args...>::template Context<StripCVRef<ArrItem>>;
+
+    template<class ArrItem, class Bed>
+    static Context<ArrItem> create_context(
+        RectangleOverfitPackingStrategy<Args...> &ps,
+        const Bed &bed,
+        int bed_index)
+    {
+        return Context<ArrItem>{bounding_box(bed), bed_index,
+                                ps.post_alignment_fn};
+    }
+};
+
+template<class ArrItem>
+struct PackingContextTraits_<RectangleOverfitPackingContext<ArrItem>>
+    : public PackingContextTraits_<DefaultPackingContext<ArrItem>>
+{
+    static void add_packed_item(RectangleOverfitPackingContext<ArrItem> &ctx, ArrItem &itm)
+    {
+        ctx.add_packed_item(itm);
+
+        // to prevent coords going out of range
+        ctx.align_pile();
+    }
+};
+
+template<class Strategy, class ArrItem, class Bed, class RemIt>
+bool pack(Strategy &strategy,
+          const Bed &bed,
+          ArrItem &item,
+          const PackStrategyContext<Strategy, ArrItem> &packing_context,
+          const Range<RemIt> &remaining_items,
+          const RectangleOverfitPackingStrategyTag &)
+{
+    bool ret = false;
+
+    if (fixed_items_range(packing_context).empty()) {
+        auto &base = strategy.base_strategy;
+        PackStrategyNFP modded_strategy{
+            base.solver,
+            RectangleOverfitKernelWrapper{base.kernel, packing_context.limits},
+            base.ep, base.accuracy};
+
+        ret = pack(modded_strategy,
+                   InfiniteBed{packing_context.limits.center()}, item,
+                   packing_context, remaining_items, NFPPackingTag{});
+    } else {
+        ret = pack(strategy.base_strategy, bed, item, packing_context,
+                   remaining_items, NFPPackingTag{});
+    }
+
+    return ret;
+}
+
+}} // namespace Slic3r::arr2
+
+#endif // RECTANGLEOVERFITPACKINGSTRATEGY_HPP