diff --git a/src/libnest2d/CMakeLists.txt b/src/libnest2d/CMakeLists.txt
index 1faf542dd..7853a48d5 100644
--- a/src/libnest2d/CMakeLists.txt
+++ b/src/libnest2d/CMakeLists.txt
@@ -70,7 +70,10 @@ if(TBB_FOUND)
     # The Intel TBB library will use the std::exception_ptr feature of C++11.
     target_compile_definitions(libnest2d INTERFACE -DTBB_USE_CAPTURED_EXCEPTION=0)
 
-    target_link_libraries(libnest2d INTERFACE tbb)
+    find_package(Threads REQUIRED)
+    target_link_libraries(libnest2d INTERFACE ${TBB_LIBRARIES} ${CMAKE_DL_LIBS}
+        Threads::Threads
+        )
 else()
    find_package(OpenMP QUIET)
 
@@ -88,7 +91,7 @@ endif()
 add_subdirectory(${SRC_DIR}/libnest2d/backends/${LIBNEST2D_GEOMETRIES})
 add_subdirectory(${SRC_DIR}/libnest2d/optimizers/${LIBNEST2D_OPTIMIZER})
 
-#target_sources(libnest2d INTERFACE ${LIBNEST2D_SRCFILES})
+target_sources(libnest2d INTERFACE ${LIBNEST2D_SRCFILES})
 target_include_directories(libnest2d INTERFACE ${SRC_DIR})
 
 if(NOT LIBNEST2D_HEADER_ONLY)
diff --git a/src/libnest2d/include/libnest2d/backends/clipper/CMakeLists.txt b/src/libnest2d/include/libnest2d/backends/clipper/CMakeLists.txt
index aa53f957e..995afcc76 100644
--- a/src/libnest2d/include/libnest2d/backends/clipper/CMakeLists.txt
+++ b/src/libnest2d/include/libnest2d/backends/clipper/CMakeLists.txt
@@ -62,9 +62,9 @@ if(NOT Boost_INCLUDE_DIRS_FOUND)
 endif()
 
 target_include_directories(ClipperBackend INTERFACE ${Boost_INCLUDE_DIRS} )
-#target_sources(ClipperBackend INTERFACE
-#    ${CMAKE_CURRENT_SOURCE_DIR}/geometries.hpp
-#    ${SRC_DIR}/libnest2d/utils/boost_alg.hpp )
+target_sources(ClipperBackend INTERFACE
+    ${CMAKE_CURRENT_SOURCE_DIR}/geometries.hpp
+    ${SRC_DIR}/libnest2d/utils/boost_alg.hpp )
 
 target_compile_definitions(ClipperBackend INTERFACE LIBNEST2D_BACKEND_CLIPPER)
 
diff --git a/src/libnest2d/include/libnest2d/backends/clipper/geometries.hpp b/src/libnest2d/include/libnest2d/backends/clipper/geometries.hpp
index c05d08d0d..9f881e7e0 100644
--- a/src/libnest2d/include/libnest2d/backends/clipper/geometries.hpp
+++ b/src/libnest2d/include/libnest2d/backends/clipper/geometries.hpp
@@ -113,6 +113,7 @@ template<> struct CountourType<PolygonImpl> {
 
 template<> struct ShapeTag<PolygonImpl> { using Type = PolygonTag; };
 template<> struct ShapeTag<PathImpl> { using Type = PathTag; };
+template<> struct ShapeTag<PointImpl> { using Type = PointTag; };
 
 template<> struct ShapeTag<TMultiShape<PolygonImpl>> {
     using Type = MultiPolygonTag;
diff --git a/src/libnest2d/include/libnest2d/geometry_traits.hpp b/src/libnest2d/include/libnest2d/geometry_traits.hpp
index 828044afe..917f5280d 100644
--- a/src/libnest2d/include/libnest2d/geometry_traits.hpp
+++ b/src/libnest2d/include/libnest2d/geometry_traits.hpp
@@ -69,12 +69,14 @@ struct PointPair {
     RawPoint p2;
 };
 
+struct PointTag {};
 struct PolygonTag {};
 struct PathTag {};
 struct MultiPolygonTag {};
 struct BoxTag {};
 struct CircleTag {};
 
+/// Meta-functions to derive the tags
 template<class Shape> struct ShapeTag { using Type = typename Shape::Tag; };
 template<class S> using Tag = typename ShapeTag<remove_cvref_t<S>>::Type;
 
@@ -131,7 +133,7 @@ public:
     _Circle(const RawPoint& center, double r): center_(center), radius_(r) {}
 
     inline const RawPoint& center() const BP2D_NOEXCEPT { return center_; }
-    inline const void center(const RawPoint& c) { center_ = c; }
+    inline void center(const RawPoint& c) { center_ = c; }
 
     inline double radius() const BP2D_NOEXCEPT { return radius_; }
     inline void radius(double r) { radius_ = r; }
@@ -518,21 +520,19 @@ inline bool intersects(const RawShape& /*sh*/, const RawShape& /*sh*/)
     return false;
 }
 
-template<class RawShape>
-inline bool isInside(const TPoint<RawShape>& /*point*/,
-                     const RawShape& /*shape*/)
-{
-    static_assert(always_false<RawShape>::value,
-                  "shapelike::isInside(point, shape) unimplemented!");
+template<class TGuest, class THost>
+inline bool isInside(const TGuest&, const THost&,
+                     const PointTag&, const PolygonTag&) {
+    static_assert(always_false<THost>::value,
+                      "shapelike::isInside(point, path) unimplemented!");
     return false;
 }
 
-template<class RawShape>
-inline bool isInside(const RawShape& /*shape*/,
-                     const RawShape& /*shape*/)
-{
-    static_assert(always_false<RawShape>::value,
-                  "shapelike::isInside(shape, shape) unimplemented!");
+template<class TGuest, class THost>
+inline bool isInside(const TGuest&, const THost&,
+                     const PolygonTag&, const PolygonTag&) {
+    static_assert(always_false<THost>::value,
+                      "shapelike::isInside(shape, shape) unimplemented!");
     return false;
 }
 
@@ -651,7 +651,7 @@ template<class RawPath> inline bool isConvex(const RawPath& sh, const PathTag&)
 
 template<class RawShape>
 inline typename TContour<RawShape>::iterator
-begin(RawShape& sh, const PolygonTag& t)
+begin(RawShape& sh, const PolygonTag&)
 {
     return begin(contour(sh), PathTag());
 }
@@ -818,16 +818,16 @@ inline auto convexHull(const RawShape& sh)
     return convexHull(sh, Tag<RawShape>());
 }
 
-template<class RawShape>
-inline bool isInside(const TPoint<RawShape>& point,
-                     const _Circle<TPoint<RawShape>>& circ)
+template<class TP, class TC>
+inline bool isInside(const TP& point, const TC& circ,
+                     const PointTag&, const CircleTag&)
 {
     return pointlike::distance(point, circ.center()) < circ.radius();
 }
 
-template<class RawShape>
-inline bool isInside(const TPoint<RawShape>& point,
-                     const _Box<TPoint<RawShape>>& box)
+template<class TP, class TB>
+inline bool isInside(const TP& point, const TB& box,
+                     const PointTag&, const BoxTag&)
 {
     auto px = getX(point);
     auto py = getY(point);
@@ -839,27 +839,27 @@ inline bool isInside(const TPoint<RawShape>& point,
     return px > minx && px < maxx && py > miny && py < maxy;
 }
 
-template<class RawShape>
-inline bool isInside(const RawShape& sh,
-                     const _Circle<TPoint<RawShape>>& circ)
+template<class RawShape, class TC>
+inline bool isInside(const RawShape& sh, const TC& circ,
+                     const PolygonTag&, const CircleTag&)
 {
-    return std::all_of(cbegin(sh), cend(sh),
-                       [&circ](const TPoint<RawShape>& p){
-        return isInside<RawShape>(p, circ);
+    return std::all_of(cbegin(sh), cend(sh), [&circ](const TPoint<RawShape>& p)
+    {
+        return isInside(p, circ, PointTag(), CircleTag());
     });
 }
 
-template<class RawShape>
-inline bool isInside(const _Box<TPoint<RawShape>>& box,
-                     const _Circle<TPoint<RawShape>>& circ)
+template<class TB, class TC>
+inline bool isInside(const TB& box, const TC& circ,
+                     const BoxTag&, const CircleTag&)
 {
-    return isInside<RawShape>(box.minCorner(), circ) &&
-            isInside<RawShape>(box.maxCorner(), circ);
+    return isInside(box.minCorner(), circ, BoxTag(), CircleTag()) &&
+           isInside(box.maxCorner(), circ, BoxTag(), CircleTag());
 }
 
-template<class RawShape>
-inline bool isInside(const _Box<TPoint<RawShape>>& ibb,
-                     const _Box<TPoint<RawShape>>& box)
+template<class TBGuest, class TBHost>
+inline bool isInside(const TBGuest& ibb, const TBHost& box,
+                     const BoxTag&, const BoxTag&)
 {
     auto iminX = getX(ibb.minCorner());
     auto imaxX = getX(ibb.maxCorner());
@@ -874,6 +874,18 @@ inline bool isInside(const _Box<TPoint<RawShape>>& ibb,
     return iminX > minX && imaxX < maxX && iminY > minY && imaxY < maxY;
 }
 
+template<class RawShape, class TB>
+inline bool isInside(const RawShape& poly, const TB& box,
+                     const PolygonTag&, const BoxTag&)
+{
+    return isInside(boundingBox(poly), box, BoxTag(), BoxTag());
+}
+
+template<class TGuest, class THost>
+inline bool isInside(const TGuest& guest, const THost& host) {
+    return isInside(guest, host, Tag<TGuest>(), Tag<THost>());
+}
+
 template<class RawShape> // Potential O(1) implementation may exist
 inline TPoint<RawShape>& vertex(RawShape& sh, unsigned long idx,
                                 const PolygonTag&)
diff --git a/src/libnest2d/include/libnest2d/geometry_traits_nfp.hpp b/src/libnest2d/include/libnest2d/geometry_traits_nfp.hpp
index b57b8dc53..cb0580ef4 100644
--- a/src/libnest2d/include/libnest2d/geometry_traits_nfp.hpp
+++ b/src/libnest2d/include/libnest2d/geometry_traits_nfp.hpp
@@ -251,6 +251,460 @@ inline NfpResult<RawShape> nfpConvexOnly(const RawShape& sh,
     return {rsh, top_nfp};
 }
 
+template<class RawShape>
+NfpResult<RawShape> nfpSimpleSimple(const RawShape& cstationary,
+                                           const RawShape& cother)
+{
+
+    // Algorithms are from the original algorithm proposed in paper:
+    // https://eprints.soton.ac.uk/36850/1/CORMSIS-05-05.pdf
+
+    // /////////////////////////////////////////////////////////////////////////
+    // Algorithm 1: Obtaining the minkowski sum
+    // /////////////////////////////////////////////////////////////////////////
+
+    // I guess this is not a full minkowski sum of the two input polygons by
+    // definition. This yields a subset that is compatible with the next 2
+    // algorithms.
+
+    using Result = NfpResult<RawShape>;
+    using Vertex = TPoint<RawShape>;
+    using Coord = TCoord<Vertex>;
+    using Edge = _Segment<Vertex>;
+    namespace sl = shapelike;
+    using std::signbit;
+    using std::sort;
+    using std::vector;
+    using std::ref;
+    using std::reference_wrapper;
+
+    // TODO The original algorithms expects the stationary polygon in
+    // counter clockwise and the orbiter in clockwise order.
+    // So for preventing any further complication, I will make the input
+    // the way it should be, than make my way around the orientations.
+
+    // Reverse the stationary contour to counter clockwise
+    auto stcont = sl::contour(cstationary);
+    {
+        std::reverse(sl::begin(stcont), sl::end(stcont));
+        stcont.pop_back();
+        auto it = std::min_element(sl::begin(stcont), sl::end(stcont),
+                               [](const Vertex& v1, const Vertex& v2) {
+            return getY(v1) < getY(v2);
+        });
+        std::rotate(sl::begin(stcont), it, sl::end(stcont));
+        sl::addVertex(stcont, sl::front(stcont));
+    }
+    RawShape stationary;
+    sl::contour(stationary) = stcont;
+
+    // Reverse the orbiter contour to counter clockwise
+    auto orbcont = sl::contour(cother);
+    {
+        std::reverse(orbcont.begin(), orbcont.end());
+
+        // Step 1: Make the orbiter reverse oriented
+
+        orbcont.pop_back();
+        auto it = std::min_element(orbcont.begin(), orbcont.end(),
+                              [](const Vertex& v1, const Vertex& v2) {
+            return getY(v1) < getY(v2);
+        });
+
+        std::rotate(orbcont.begin(), it, orbcont.end());
+        orbcont.emplace_back(orbcont.front());
+
+        for(auto &v : orbcont) v = -v;
+
+    }
+
+    // Copy the orbiter (contour only), we will have to work on it
+    RawShape orbiter;
+    sl::contour(orbiter) = orbcont;
+
+    // An edge with additional data for marking it
+    struct MarkedEdge {
+        Edge e; Radians turn_angle = 0; bool is_turning_point = false;
+        MarkedEdge() = default;
+        MarkedEdge(const Edge& ed, Radians ta, bool tp):
+            e(ed), turn_angle(ta), is_turning_point(tp) {}
+
+        // debug
+        std::string label;
+    };
+
+    // Container for marked edges
+    using EdgeList = vector<MarkedEdge>;
+
+    EdgeList A, B;
+
+    // This is how an edge list is created from the polygons
+    auto fillEdgeList = [](EdgeList& L, const RawShape& ppoly, int dir) {
+        auto& poly = sl::contour(ppoly);
+
+        L.reserve(sl::contourVertexCount(poly));
+
+        if(dir > 0) {
+            auto it = poly.begin();
+            auto nextit = std::next(it);
+
+            double turn_angle = 0;
+            bool is_turn_point = false;
+
+            while(nextit != poly.end()) {
+                L.emplace_back(Edge(*it, *nextit), turn_angle, is_turn_point);
+                it++; nextit++;
+            }
+        } else {
+            auto it = sl::rbegin(poly);
+            auto nextit = std::next(it);
+
+            double turn_angle = 0;
+            bool is_turn_point = false;
+
+            while(nextit != sl::rend(poly)) {
+                L.emplace_back(Edge(*it, *nextit), turn_angle, is_turn_point);
+                it++; nextit++;
+            }
+        }
+
+        auto getTurnAngle = [](const Edge& e1, const Edge& e2) {
+            auto phi = e1.angleToXaxis();
+            auto phi_prev = e2.angleToXaxis();
+            auto turn_angle = phi-phi_prev;
+            if(turn_angle > Pi) turn_angle -= TwoPi;
+            if(turn_angle < -Pi) turn_angle += TwoPi;
+            return turn_angle;
+        };
+
+        auto eit = L.begin();
+        auto enext = std::next(eit);
+
+        eit->turn_angle = getTurnAngle(L.front().e, L.back().e);
+
+        while(enext != L.end()) {
+            enext->turn_angle = getTurnAngle( enext->e, eit->e);
+            eit->is_turning_point =
+                    signbit(enext->turn_angle) != signbit(eit->turn_angle);
+            ++eit; ++enext;
+        }
+
+        L.back().is_turning_point = signbit(L.back().turn_angle) !=
+                                    signbit(L.front().turn_angle);
+
+    };
+
+    // Step 2: Fill the edgelists
+    fillEdgeList(A, stationary, 1);
+    fillEdgeList(B, orbiter, 1);
+
+    int i = 1;
+    for(MarkedEdge& me : A) {
+        std::cout << "a" << i << ":\n\t"
+                  << getX(me.e.first()) << " " << getY(me.e.first()) << "\n\t"
+                  << getX(me.e.second()) << " " << getY(me.e.second()) << "\n\t"
+                  << "Turning point: " << (me.is_turning_point ? "yes" : "no")
+                  << std::endl;
+
+        me.label = "a"; me.label += std::to_string(i);
+        i++;
+    }
+
+    i = 1;
+    for(MarkedEdge& me : B) {
+        std::cout << "b" << i << ":\n\t"
+                  << getX(me.e.first()) << " " << getY(me.e.first()) << "\n\t"
+                  << getX(me.e.second()) << " " << getY(me.e.second()) << "\n\t"
+                  << "Turning point: " << (me.is_turning_point ? "yes" : "no")
+                  << std::endl;
+        me.label = "b"; me.label += std::to_string(i);
+        i++;
+    }
+
+    // A reference to a marked edge that also knows its container
+    struct MarkedEdgeRef {
+        reference_wrapper<MarkedEdge> eref;
+        reference_wrapper<vector<MarkedEdgeRef>> container;
+        Coord dir = 1;  // Direction modifier
+
+        inline Radians angleX() const { return eref.get().e.angleToXaxis(); }
+        inline const Edge& edge() const { return eref.get().e; }
+        inline Edge& edge() { return eref.get().e; }
+        inline bool isTurningPoint() const {
+            return eref.get().is_turning_point;
+        }
+        inline bool isFrom(const vector<MarkedEdgeRef>& cont ) {
+            return &(container.get()) == &cont;
+        }
+        inline bool eq(const MarkedEdgeRef& mr) {
+            return &(eref.get()) == &(mr.eref.get());
+        }
+
+        MarkedEdgeRef(reference_wrapper<MarkedEdge> er,
+                      reference_wrapper<vector<MarkedEdgeRef>> ec):
+            eref(er), container(ec), dir(1) {}
+
+        MarkedEdgeRef(reference_wrapper<MarkedEdge> er,
+                      reference_wrapper<vector<MarkedEdgeRef>> ec,
+                      Coord d):
+            eref(er), container(ec), dir(d) {}
+    };
+
+    using EdgeRefList = vector<MarkedEdgeRef>;
+
+    // Comparing two marked edges
+    auto sortfn = [](const MarkedEdgeRef& e1, const MarkedEdgeRef& e2) {
+        return e1.angleX() < e2.angleX();
+    };
+
+    EdgeRefList Aref, Bref;     // We create containers for the references
+    Aref.reserve(A.size()); Bref.reserve(B.size());
+
+    // Fill reference container for the stationary polygon
+    std::for_each(A.begin(), A.end(), [&Aref](MarkedEdge& me) {
+        Aref.emplace_back( ref(me), ref(Aref) );
+    });
+
+    // Fill reference container for the orbiting polygon
+    std::for_each(B.begin(), B.end(), [&Bref](MarkedEdge& me) {
+        Bref.emplace_back( ref(me), ref(Bref) );
+    });
+
+    auto mink = [sortfn] // the Mink(Q, R, direction) sub-procedure
+            (const EdgeRefList& Q, const EdgeRefList& R, bool positive)
+    {
+
+        // Step 1 "merge sort_list(Q) and sort_list(R) to form merge_list(Q,R)"
+        // Sort the containers of edge references and merge them.
+        // Q could be sorted only once and be reused here but we would still
+        // need to merge it with sorted(R).
+
+        EdgeRefList merged;
+        EdgeRefList S, seq;
+        merged.reserve(Q.size() + R.size());
+
+        merged.insert(merged.end(), R.begin(), R.end());
+        std::stable_sort(merged.begin(), merged.end(), sortfn);
+        merged.insert(merged.end(), Q.begin(), Q.end());
+        std::stable_sort(merged.begin(), merged.end(), sortfn);
+
+        // Step 2 "set i = 1, k = 1, direction = 1, s1 = q1"
+        // we don't use i, instead, q is an iterator into Q. k would be an index
+        // into the merged sequence but we use "it" as an iterator for that
+
+        // here we obtain references for the containers for later comparisons
+        const auto& Rcont = R.begin()->container.get();
+        const auto& Qcont = Q.begin()->container.get();
+
+        // Set the initial direction
+        Coord dir = 1;
+
+        // roughly i = 1 (so q = Q.begin()) and s1 = q1 so S[0] = q;
+        if(positive) {
+            auto q = Q.begin();
+            S.emplace_back(*q);
+
+            // Roughly step 3
+
+            std::cout << "merged size: " << merged.size() << std::endl;
+            auto mit = merged.begin();
+            for(bool finish = false; !finish && q != Q.end();) {
+                ++q; // "Set i = i + 1"
+
+                while(!finish && mit != merged.end()) {
+                    if(mit->isFrom(Rcont)) {
+                        auto s = *mit;
+                        s.dir = dir;
+                        S.emplace_back(s);
+                    }
+
+                    if(mit->eq(*q)) {
+                        S.emplace_back(*q);
+                        if(mit->isTurningPoint()) dir = -dir;
+                        if(q == Q.begin()) finish = true;
+                        break;
+                    }
+
+                    mit += dir;
+    //                __nfp::advance(mit, merged, dir > 0);
+                }
+            }
+        } else {
+            auto q = Q.rbegin();
+            S.emplace_back(*q);
+
+            // Roughly step 3
+
+            std::cout << "merged size: " << merged.size() << std::endl;
+            auto mit = merged.begin();
+            for(bool finish = false; !finish && q != Q.rend();) {
+                ++q; // "Set i = i + 1"
+
+                while(!finish && mit != merged.end()) {
+                    if(mit->isFrom(Rcont)) {
+                        auto s = *mit;
+                        s.dir = dir;
+                        S.emplace_back(s);
+                    }
+
+                    if(mit->eq(*q)) {
+                        S.emplace_back(*q);
+                        S.back().dir = -1;
+                        if(mit->isTurningPoint()) dir = -dir;
+                        if(q == Q.rbegin()) finish = true;
+                        break;
+                    }
+
+                    mit += dir;
+            //                __nfp::advance(mit, merged, dir > 0);
+                }
+            }
+        }
+
+
+        // Step 4:
+
+        // "Let starting edge r1 be in position si in sequence"
+        // whaaat? I guess this means the following:
+        auto it = S.begin();
+        while(!it->eq(*R.begin())) ++it;
+
+        // "Set j = 1, next = 2, direction = 1, seq1 = si"
+        // we don't use j, seq is expanded dynamically.
+        dir = 1;
+        auto next = std::next(R.begin()); seq.emplace_back(*it);
+
+        // Step 5:
+        // "If all si edges have been allocated to seqj" should mean that
+        // we loop until seq has equal size with S
+        auto send = it; //it == S.begin() ? it : std::prev(it);
+        while(it != S.end()) {
+            ++it; if(it == S.end()) it = S.begin();
+            if(it == send) break;
+
+            if(it->isFrom(Qcont)) {
+                seq.emplace_back(*it); // "If si is from Q, j = j + 1, seqj = si"
+
+                // "If si is a turning point in Q,
+                // direction = - direction, next = next + direction"
+                if(it->isTurningPoint()) {
+                    dir = -dir;
+                    next += dir;
+//                    __nfp::advance(next, R, dir > 0);
+                }
+            }
+
+            if(it->eq(*next) /*&& dir == next->dir*/) { // "If si = direction.rnext"
+                // "j = j + 1, seqj = si, next = next + direction"
+                seq.emplace_back(*it);
+                next += dir;
+//                __nfp::advance(next, R, dir > 0);
+            }
+        }
+
+        return seq;
+    };
+
+    std::vector<EdgeRefList> seqlist;
+    seqlist.reserve(Bref.size());
+
+    EdgeRefList Bslope = Bref;  // copy Bref, we will make a slope diagram
+
+    // make the slope diagram of B
+    std::sort(Bslope.begin(), Bslope.end(), sortfn);
+
+    auto slopeit = Bslope.begin(); // search for the first turning point
+    while(!slopeit->isTurningPoint() && slopeit != Bslope.end()) slopeit++;
+
+    if(slopeit == Bslope.end()) {
+        // no turning point means convex polygon.
+        seqlist.emplace_back(mink(Aref, Bref, true));
+    } else {
+        int dir = 1;
+
+        auto firstturn = Bref.begin();
+        while(!firstturn->eq(*slopeit)) ++firstturn;
+
+        assert(firstturn != Bref.end());
+
+        EdgeRefList bgroup; bgroup.reserve(Bref.size());
+        bgroup.emplace_back(*slopeit);
+
+        auto b_it = std::next(firstturn);
+        while(b_it != firstturn) {
+            if(b_it == Bref.end()) b_it = Bref.begin();
+
+            while(!slopeit->eq(*b_it)) {
+                __nfp::advance(slopeit, Bslope, dir > 0);
+            }
+
+            if(!slopeit->isTurningPoint()) {
+                bgroup.emplace_back(*slopeit);
+            } else {
+                if(!bgroup.empty()) {
+                    if(dir > 0) bgroup.emplace_back(*slopeit);
+                    for(auto& me : bgroup) {
+                        std::cout << me.eref.get().label << ", ";
+                    }
+                    std::cout << std::endl;
+                    seqlist.emplace_back(mink(Aref, bgroup, dir == 1 ? true : false));
+                    bgroup.clear();
+                    if(dir < 0) bgroup.emplace_back(*slopeit);
+                } else {
+                    bgroup.emplace_back(*slopeit);
+                }
+
+                dir *= -1;
+            }
+            ++b_it;
+        }
+    }
+
+//    while(it != Bref.end()) // This is step 3 and step 4 in one loop
+//        if(it->isTurningPoint()) {
+//            R = {R.last, it++};
+//            auto seq = mink(Q, R, orientation);
+
+//            // TODO step 6 (should be 5 shouldn't it?): linking edges from A
+//            // I don't get this step
+
+//            seqlist.insert(seqlist.end(), seq.begin(), seq.end());
+//            orientation = !orientation;
+//        } else ++it;
+
+//    if(seqlist.empty()) seqlist = mink(Q, {Bref.begin(), Bref.end()}, true);
+
+    // /////////////////////////////////////////////////////////////////////////
+    // Algorithm 2: breaking Minkowski sums into track line trips
+    // /////////////////////////////////////////////////////////////////////////
+
+
+    // /////////////////////////////////////////////////////////////////////////
+    // Algorithm 3: finding the boundary of the NFP from track line trips
+    // /////////////////////////////////////////////////////////////////////////
+
+
+    for(auto& seq : seqlist) {
+        std::cout << "seqlist size: " << seq.size() << std::endl;
+        for(auto& s : seq) {
+            std::cout << (s.dir > 0 ? "" : "-") << s.eref.get().label << ", ";
+        }
+        std::cout << std::endl;
+    }
+
+    auto& seq = seqlist.front();
+    RawShape rsh;
+    Vertex top_nfp;
+    std::vector<Edge> edgelist; edgelist.reserve(seq.size());
+    for(auto& s : seq) {
+        edgelist.emplace_back(s.eref.get().e);
+    }
+
+    __nfp::buildPolygon(edgelist, rsh, top_nfp);
+
+    return Result(rsh, top_nfp);
+}
+
 // Specializable NFP implementation class. Specialize it if you have a faster
 // or better NFP implementation
 template<class RawShape, NfpLevel nfptype>
diff --git a/src/libnest2d/include/libnest2d/libnest2d.hpp b/src/libnest2d/include/libnest2d/libnest2d.hpp
index aac62e094..49baa65f2 100644
--- a/src/libnest2d/include/libnest2d/libnest2d.hpp
+++ b/src/libnest2d/include/libnest2d/libnest2d.hpp
@@ -482,17 +482,40 @@ public:
 
 template<class RawShape>
 inline bool _Item<RawShape>::isInside(const _Box<TPoint<RawShape>>& box) const {
-    return sl::isInside<RawShape>(boundingBox(), box);
+    return sl::isInside(boundingBox(), box);
 }
 
 template<class RawShape> inline bool
 _Item<RawShape>::isInside(const _Circle<TPoint<RawShape>>& circ) const {
-    return sl::isInside<RawShape>(transformedShape(), circ);
+    return sl::isInside(transformedShape(), circ);
 }
 
+template<class RawShape> using _ItemRef = std::reference_wrapper<_Item<RawShape>>;
+template<class RawShape> using _ItemGroup = std::vector<_ItemRef<RawShape>>;
 
-template<class I> using _ItemRef = std::reference_wrapper<I>;
-template<class I> using _ItemGroup = std::vector<_ItemRef<I>>;
+/**
+ * \brief A list of packed item vectors. Each vector represents a bin.
+ */
+template<class RawShape>
+using _PackGroup = std::vector<std::vector<_ItemRef<RawShape>>>;
+
+/**
+ * \brief A list of packed (index, item) pair vectors. Each vector represents a
+ * bin.
+ *
+ * The index is points to the position of the item in the original input
+ * sequence. This way the caller can use the items as a transformation data
+ * carrier and transform the original objects manually.
+ */
+template<class RawShape>
+using _IndexedPackGroup = std::vector<
+                               std::vector<
+                                   std::pair<
+                                       unsigned,
+                                       _ItemRef<RawShape>
+                                   >
+                               >
+                          >;
 
 template<class Iterator>
 struct ConstItemRange {
@@ -524,8 +547,10 @@ class PlacementStrategyLike {
     PlacementStrategy impl_;
 public:
 
+    using RawShape = typename PlacementStrategy::ShapeType;
+
     /// The item type that the placer works with.
-    using Item = typename PlacementStrategy::Item;
+    using Item = _Item<RawShape>;
 
     /// The placer's config type. Should be a simple struct but can be anything.
     using Config = typename PlacementStrategy::Config;
@@ -544,8 +569,7 @@ public:
      */
     using PackResult = typename PlacementStrategy::PackResult;
 
-    using ItemRef = _ItemRef<Item>;
-    using ItemGroup = _ItemGroup<Item>;
+    using ItemGroup = _ItemGroup<RawShape>;
     using DefaultIterator = typename ItemGroup::const_iterator;
 
     /**
@@ -619,6 +643,15 @@ public:
         return impl_.pack(item, remaining);
     }
 
+    /**
+     * This method makes possible to "preload" some items into the placer. It
+     * will not move these items but will consider them as already packed.
+     */
+    inline void preload(const ItemGroup& packeditems)
+    {
+        impl_.preload(packeditems);
+    }
+
     /// Unpack the last element (remove it from the list of packed items).
     inline void unpackLast() { impl_.unpackLast(); }
 
@@ -649,11 +682,11 @@ template<class SelectionStrategy>
 class SelectionStrategyLike {
     SelectionStrategy impl_;
 public:
-    using Item = typename SelectionStrategy::Item;
+    using RawShape = typename SelectionStrategy::ShapeType;
+    using Item = _Item<RawShape>;
+    using PackGroup = _PackGroup<RawShape>;
     using Config = typename SelectionStrategy::Config;
 
-    using ItemRef = std::reference_wrapper<Item>;
-    using ItemGroup = std::vector<ItemRef>;
 
     /**
      * @brief Provide a different configuration for the selection strategy.
@@ -703,60 +736,29 @@ public:
                                  std::forward<PConfig>(config));
     }
 
-    /**
-     * \brief Get the number of bins opened by the selection algorithm.
-     *
-     * Initially it is zero and after the call to packItems it will return
-     * the number of bins opened by the packing procedure.
-     *
-     * \return The number of bins opened.
-     */
-    inline size_t binCount() const { return impl_.binCount(); }
-
     /**
      * @brief Get the items for a particular bin.
      * @param binIndex The index of the requested bin.
      * @return Returns a list of all items packed into the requested bin.
      */
-    inline ItemGroup itemsForBin(size_t binIndex) {
-        return impl_.itemsForBin(binIndex);
+    inline const PackGroup& getResult() const {
+        return impl_.getResult();
     }
 
-    /// Same as itemsForBin but for a const context.
-    inline const ItemGroup itemsForBin(size_t binIndex) const {
-        return impl_.itemsForBin(binIndex);
-    }
+    /**
+     * @brief Loading a group of already packed bins. It is best to use a result
+     * from a previous packing. The algorithm will consider this input as if the
+     * objects are already packed and not move them. If any of these items are
+     * outside the bin, it is up to the placer algorithm what will happen.
+     * Packing additional items can fail for the bottom-left and nfp placers.
+     * @param pckgrp A packgroup which is a vector of item vectors. Each item
+     * vector corresponds to a packed bin.
+     */
+    inline void preload(const PackGroup& pckgrp) { impl_.preload(pckgrp); }
+
+    void clear() { impl_.clear(); }
 };
 
-
-/**
- * \brief A list of packed item vectors. Each vector represents a bin.
- */
-template<class RawShape>
-using _PackGroup = std::vector<
-                        std::vector<
-                            std::reference_wrapper<_Item<RawShape>>
-                        >
-                   >;
-
-/**
- * \brief A list of packed (index, item) pair vectors. Each vector represents a
- * bin.
- *
- * The index is points to the position of the item in the original input
- * sequence. This way the caller can use the items as a transformation data
- * carrier and transform the original objects manually.
- */
-template<class RawShape>
-using _IndexedPackGroup = std::vector<
-                               std::vector<
-                                   std::pair<
-                                       unsigned,
-                                       std::reference_wrapper<_Item<RawShape>>
-                                   >
-                               >
-                          >;
-
 /**
  * The Arranger is the front-end class for the libnest2d library. It takes the
  * input items and outputs the items with the proper transformations to be
@@ -868,17 +870,29 @@ public:
     }
 
     /// Set a predicate to tell when to abort nesting.
-    inline Nester& stopCondition(StopCondition fn) {
+    inline Nester& stopCondition(StopCondition fn)
+    {
         selector_.stopCondition(fn); return *this;
     }
 
-    inline PackGroup lastResult() {
-        PackGroup ret;
-        for(size_t i = 0; i < selector_.binCount(); i++) {
-            auto items = selector_.itemsForBin(i);
-            ret.push_back(items);
+    inline const PackGroup& lastResult() const
+    {
+        return selector_.getResult();
+    }
+
+    inline void preload(const PackGroup& pgrp)
+    {
+        selector_.preload(pgrp);
+    }
+
+    inline void preload(const IndexedPackGroup& ipgrp)
+    {
+        PackGroup pgrp; pgrp.reserve(ipgrp.size());
+        for(auto& ig : ipgrp) {
+            pgrp.emplace_back(); pgrp.back().reserve(ig.size());
+            for(auto& r : ig) pgrp.back().emplace_back(r.second);
         }
-        return ret;
+        preload(pgrp);
     }
 
 private:
@@ -892,7 +906,7 @@ private:
              // have to exist for the lifetime of this call.
              class T = enable_if_t< std::is_convertible<IT, TPItem>::value, IT>
              >
-    inline PackGroup _execute(TIterator from, TIterator to, bool = false)
+    inline const PackGroup& _execute(TIterator from, TIterator to, bool = false)
     {
         __execute(from, to);
         return lastResult();
@@ -902,7 +916,7 @@ private:
              class IT = remove_cvref_t<typename TIterator::value_type>,
              class T = enable_if_t<!std::is_convertible<IT, TPItem>::value, IT>
              >
-    inline PackGroup _execute(TIterator from, TIterator to, int = false)
+    inline const PackGroup& _execute(TIterator from, TIterator to, int = false)
     {
         item_cache_ = {from, to};
 
@@ -946,10 +960,12 @@ private:
                                                    TSel& selector)
     {
         IndexedPackGroup pg;
-        pg.reserve(selector.binCount());
+        pg.reserve(selector.getResult().size());
 
-        for(size_t i = 0; i < selector.binCount(); i++) {
-            auto items = selector.itemsForBin(i);
+        const PackGroup& pckgrp = selector.getResult();
+
+        for(size_t i = 0; i < pckgrp.size(); i++) {
+            auto items = pckgrp[i];
             pg.push_back({});
             pg[i].reserve(items.size());
 
diff --git a/src/libnest2d/include/libnest2d/optimizers/nlopt/CMakeLists.txt b/src/libnest2d/include/libnest2d/optimizers/nlopt/CMakeLists.txt
index 2a32019f4..5559ad645 100644
--- a/src/libnest2d/include/libnest2d/optimizers/nlopt/CMakeLists.txt
+++ b/src/libnest2d/include/libnest2d/optimizers/nlopt/CMakeLists.txt
@@ -48,12 +48,12 @@ else()
     target_link_libraries(NloptOptimizer INTERFACE Nlopt::Nlopt)
 endif()
 
-#target_sources( NloptOptimizer INTERFACE
-#${CMAKE_CURRENT_SOURCE_DIR}/simplex.hpp
-#${CMAKE_CURRENT_SOURCE_DIR}/subplex.hpp
-#${CMAKE_CURRENT_SOURCE_DIR}/genetic.hpp
-#${CMAKE_CURRENT_SOURCE_DIR}/nlopt_boilerplate.hpp
-#)
+target_sources( NloptOptimizer INTERFACE
+${CMAKE_CURRENT_SOURCE_DIR}/simplex.hpp
+${CMAKE_CURRENT_SOURCE_DIR}/subplex.hpp
+${CMAKE_CURRENT_SOURCE_DIR}/genetic.hpp
+${CMAKE_CURRENT_SOURCE_DIR}/nlopt_boilerplate.hpp
+)
 
 target_compile_definitions(NloptOptimizer INTERFACE LIBNEST2D_OPTIMIZER_NLOPT)
 
diff --git a/src/libnest2d/include/libnest2d/placers/nfpplacer.hpp b/src/libnest2d/include/libnest2d/placers/nfpplacer.hpp
index 28659c512..6fb717a7a 100644
--- a/src/libnest2d/include/libnest2d/placers/nfpplacer.hpp
+++ b/src/libnest2d/include/libnest2d/placers/nfpplacer.hpp
@@ -130,7 +130,7 @@ namespace placers {
 template<class RawShape>
 struct NfpPConfig {
 
-    using ItemGroup = _ItemGroup<_Item<RawShape>>;
+    using ItemGroup = _ItemGroup<RawShape>;
 
     enum class Alignment {
         CENTER,
@@ -138,6 +138,8 @@ struct NfpPConfig {
         BOTTOM_RIGHT,
         TOP_LEFT,
         TOP_RIGHT,
+        DONT_ALIGN      //!> Warning: parts may end up outside the bin with the
+                        //! default object function.
     };
 
     /// Which angles to try out for better results.
@@ -545,8 +547,8 @@ public:
     _NofitPolyPlacer& operator=(const _NofitPolyPlacer&) = default;
 
 #ifndef BP2D_COMPILER_MSVC12 // MSVC2013 does not support default move ctors
-    _NofitPolyPlacer(_NofitPolyPlacer&&) /*BP2D_NOEXCEPT*/ = default;
-    _NofitPolyPlacer& operator=(_NofitPolyPlacer&&) /*BP2D_NOEXCEPT*/ = default;
+    _NofitPolyPlacer(_NofitPolyPlacer&&) = default;
+    _NofitPolyPlacer& operator=(_NofitPolyPlacer&&) = default;
 #endif
 
     static inline double overfit(const Box& bb, const RawShape& bin) {
@@ -905,26 +907,44 @@ private:
 
                 // This is the kernel part of the object function that is
                 // customizable by the library client
-                auto _objfunc = config_.object_function?
-                            config_.object_function :
-                            [norm, bin, binbb, pbb](const Item& item)
-                {
-                    auto ibb = item.boundingBox();
-                    auto fullbb = boundingBox(pbb, ibb);
+                std::function<double(const Item&)> _objfunc;
+                if(config_.object_function) _objfunc = config_.object_function;
+                else {
 
-                    double score = pl::distance(ibb.center(), binbb.center());
-                    score /= norm;
+                    // Inside check has to be strict if no alignment was enabled
+                    std::function<double(const Box&)> ins_check;
+                    if(config_.alignment == Config::Alignment::DONT_ALIGN)
+                        ins_check = [&binbb, norm](const Box& fullbb) {
+                            double ret = 0;
+                            if(!sl::isInside(fullbb, binbb))
+                                ret += norm;
+                            return ret;
+                        };
+                    else
+                        ins_check = [&bin](const Box& fullbb) {
+                            double miss = overfit(fullbb, bin);
+                            miss = miss > 0? miss : 0;
+                            return std::pow(miss, 2);
+                        };
 
-                    double miss = overfit(fullbb, bin);
-                    miss = miss > 0? miss : 0;
-                    score += std::pow(miss, 2);
+                    _objfunc = [norm, binbb, pbb, ins_check](const Item& item)
+                    {
+                        auto ibb = item.boundingBox();
+                        auto fullbb = boundingBox(pbb, ibb);
 
-                    return score;
-                };
+                        double score = pl::distance(ibb.center(),
+                                                    binbb.center());
+                        score /= norm;
+
+                        score += ins_check(fullbb);
+
+                        return score;
+                    };
+                }
 
                 // Our object function for placement
-                auto rawobjfunc =
-                        [_objfunc, iv, startpos] (Vertex v, Item& itm)
+                auto rawobjfunc = [_objfunc, iv, startpos]
+                        (Vertex v, Item& itm)
                 {
                     auto d = v - iv;
                     d += startpos;
@@ -938,9 +958,10 @@ private:
                             ecache[opt.nfpidx].coords(opt.hidx, opt.relpos);
                 };
 
-                auto boundaryCheck =
-                    [&merged_pile, &getNfpPoint, &item, &bin, &iv, &startpos]
-                    (const Optimum& o)
+                auto alignment = config_.alignment;
+
+                auto boundaryCheck = [alignment, &merged_pile, &getNfpPoint,
+                        &item, &bin, &iv, &startpos] (const Optimum& o)
                 {
                     auto v = getNfpPoint(o);
                     auto d = v - iv;
@@ -951,7 +972,12 @@ private:
                     auto chull = sl::convexHull(merged_pile);
                     merged_pile.pop_back();
 
-                    return overfit(chull, bin);
+                    double miss = 0;
+                    if(alignment == Config::Alignment::DONT_ALIGN)
+                       miss = sl::isInside(chull, bin) ? -1.0 : 1.0;
+                    else miss = overfit(chull, bin);
+
+                    return miss;
                 };
 
                 Optimum optimum(0, 0);
@@ -1101,7 +1127,9 @@ private:
     }
 
     inline void finalAlign(_Circle<TPoint<RawShape>> cbin) {
-        if(items_.empty()) return;
+        if(items_.empty() ||
+                config_.alignment == Config::Alignment::DONT_ALIGN) return;
+
         nfp::Shapes<RawShape> m;
         m.reserve(items_.size());
         for(Item& item : items_) m.emplace_back(item.transformedShape());
@@ -1113,7 +1141,9 @@ private:
     }
 
     inline void finalAlign(Box bbin) {
-        if(items_.empty()) return;
+        if(items_.empty() ||
+                config_.alignment == Config::Alignment::DONT_ALIGN) return;
+
         nfp::Shapes<RawShape> m;
         m.reserve(items_.size());
         for(Item& item : items_) m.emplace_back(item.transformedShape());
@@ -1147,6 +1177,7 @@ private:
             cb = bbin.maxCorner();
             break;
         }
+        default: ; // DONT_ALIGN
         }
 
         auto d = cb - ci;
@@ -1184,6 +1215,7 @@ private:
             cb = bbin.maxCorner();
             break;
         }
+        default:;
         }
 
         auto d = cb - ci;
diff --git a/src/libnest2d/include/libnest2d/placers/placer_boilerplate.hpp b/src/libnest2d/include/libnest2d/placers/placer_boilerplate.hpp
index 9f940af4d..309a5007d 100644
--- a/src/libnest2d/include/libnest2d/placers/placer_boilerplate.hpp
+++ b/src/libnest2d/include/libnest2d/placers/placer_boilerplate.hpp
@@ -12,6 +12,7 @@ class PlacerBoilerplate {
     mutable bool farea_valid_ = false;
     mutable double farea_ = 0.0;
 public:
+    using ShapeType = RawShape;
     using Item = _Item<RawShape>;
     using Vertex = TPoint<RawShape>;
     using Segment = _Segment<Vertex>;
@@ -19,7 +20,7 @@ public:
     using Coord = TCoord<Vertex>;
     using Unit = Coord;
     using Config = Cfg;
-    using ItemGroup = _ItemGroup<Item>;
+    using ItemGroup = _ItemGroup<RawShape>;
     using DefaultIter = typename ItemGroup::const_iterator;
 
     class PackResult {
@@ -59,8 +60,7 @@ public:
     }
 
     template<class Range = ConstItemRange<DefaultIter>>
-    bool pack(Item& item,
-              const Range& rem = Range()) {
+    bool pack(Item& item, const Range& rem = Range()) {
         auto&& r = static_cast<Subclass*>(this)->trypack(item, rem);
         if(r) {
             items_.push_back(*(r.item_ptr_));
@@ -69,6 +69,11 @@ public:
         return r;
     }
 
+    void preload(const ItemGroup& packeditems) {
+        items_.insert(items_.end(), packeditems.begin(), packeditems.end());
+        farea_valid_ = false;
+    }
+
     void accept(PackResult& r) {
         if(r) {
             r.item_ptr_->translation(r.move_);
@@ -117,6 +122,7 @@ using Base::bin_;                 \
 using Base::items_;               \
 using Base::config_;              \
 public:                           \
+using typename Base::ShapeType;   \
 using typename Base::Item;        \
 using typename Base::ItemGroup;   \
 using typename Base::BinType;     \
diff --git a/src/libnest2d/include/libnest2d/selections/djd_heuristic.hpp b/src/libnest2d/include/libnest2d/selections/djd_heuristic.hpp
index 39761f557..b03534dc4 100644
--- a/src/libnest2d/include/libnest2d/selections/djd_heuristic.hpp
+++ b/src/libnest2d/include/libnest2d/selections/djd_heuristic.hpp
@@ -33,7 +33,7 @@ class _DJDHeuristic: public SelectionBoilerplate<RawShape> {
 
 public:
     using typename Base::Item;
-    using typename Base::ItemRef;
+    using ItemRef = std::reference_wrapper<Item>;
 
     /**
      * @brief The Config for DJD heuristic.
@@ -126,6 +126,8 @@ public:
 
         store_.clear();
         store_.reserve(last-first);
+
+        // TODO: support preloading
         packed_bins_.clear();
 
         std::copy(first, last, std::back_inserter(store_));
diff --git a/src/libnest2d/include/libnest2d/selections/filler.hpp b/src/libnest2d/include/libnest2d/selections/filler.hpp
index 5f95a6eff..19c44bfaa 100644
--- a/src/libnest2d/include/libnest2d/selections/filler.hpp
+++ b/src/libnest2d/include/libnest2d/selections/filler.hpp
@@ -34,6 +34,10 @@ public:
         store_.clear();
         auto total = last-first;
         store_.reserve(total);
+
+        // TODO: support preloading
+        packed_bins_.clear();
+
         packed_bins_.emplace_back();
 
         auto makeProgress = [this, &total](
diff --git a/src/libnest2d/include/libnest2d/selections/firstfit.hpp b/src/libnest2d/include/libnest2d/selections/firstfit.hpp
index d25487d6b..d521673b4 100644
--- a/src/libnest2d/include/libnest2d/selections/firstfit.hpp
+++ b/src/libnest2d/include/libnest2d/selections/firstfit.hpp
@@ -36,11 +36,19 @@ public:
 
         store_.clear();
         store_.reserve(last-first);
-        packed_bins_.clear();
 
         std::vector<Placer> placers;
         placers.reserve(last-first);
 
+        // If the packed_items array is not empty we have to create as many
+        // placers as there are elements in packed bins and preload each item
+        // into the appropriate placer
+        for(ItemGroup& ig : packed_bins_) {
+            placers.emplace_back(bin);
+            placers.back().configure(pconfig);
+            placers.back().preload(ig);
+        }
+
         std::copy(first, last, std::back_inserter(store_));
 
         auto sortfunc = [](Item& i1, Item& i2) {
diff --git a/src/libnest2d/include/libnest2d/selections/selection_boilerplate.hpp b/src/libnest2d/include/libnest2d/selections/selection_boilerplate.hpp
index 8351a99e7..fd6577d97 100644
--- a/src/libnest2d/include/libnest2d/selections/selection_boilerplate.hpp
+++ b/src/libnest2d/include/libnest2d/selections/selection_boilerplate.hpp
@@ -9,27 +9,23 @@ namespace libnest2d { namespace selections {
 template<class RawShape>
 class SelectionBoilerplate {
 public:
+    using ShapeType = RawShape;
     using Item = _Item<RawShape>;
-    using ItemRef = std::reference_wrapper<Item>;
-    using ItemGroup = std::vector<ItemRef>;
-    using PackGroup = std::vector<ItemGroup>;
+    using ItemGroup = _ItemGroup<RawShape>;
+    using PackGroup = _PackGroup<RawShape>;
 
-    size_t binCount() const { return packed_bins_.size(); }
-
-    ItemGroup itemsForBin(size_t binIndex) {
-        assert(binIndex < packed_bins_.size());
-        return packed_bins_[binIndex];
-    }
-
-    inline const ItemGroup itemsForBin(size_t binIndex) const {
-        assert(binIndex < packed_bins_.size());
-        return packed_bins_[binIndex];
+    inline const PackGroup& getResult() const {
+        return packed_bins_;
     }
 
     inline void progressIndicator(ProgressFunction fn) { progress_ = fn; }
 
     inline void stopCondition(StopCondition cond) { stopcond_ = cond; }
 
+    inline void preload(const PackGroup& pckgrp) { packed_bins_ = pckgrp; }
+
+    inline void clear() { packed_bins_.clear(); }
+
 protected:
 
     PackGroup packed_bins_;
diff --git a/src/libnest2d/include/libnest2d/utils/boost_alg.hpp b/src/libnest2d/include/libnest2d/utils/boost_alg.hpp
index c573edb47..a6988ca00 100644
--- a/src/libnest2d/include/libnest2d/utils/boost_alg.hpp
+++ b/src/libnest2d/include/libnest2d/utils/boost_alg.hpp
@@ -356,13 +356,15 @@ inline double area(const PolygonImpl& shape, const PolygonTag&)
 #endif
 
 template<>
-inline bool isInside(const PointImpl& point, const PolygonImpl& shape)
+inline bool isInside(const PointImpl& point, const PolygonImpl& shape,
+                     const PointTag&, const PolygonTag&)
 {
     return boost::geometry::within(point, shape);
 }
 
 template<>
-inline bool isInside(const PolygonImpl& sh1, const PolygonImpl& sh2)
+inline bool isInside(const PolygonImpl& sh1, const PolygonImpl& sh2,
+                     const PolygonTag&, const PolygonTag&)
 {
     return boost::geometry::within(sh1, sh2);
 }
diff --git a/src/libslic3r/ModelArrange.cpp b/src/libslic3r/ModelArrange.cpp
index 1f517375c..d527db9da 100644
--- a/src/libslic3r/ModelArrange.cpp
+++ b/src/libslic3r/ModelArrange.cpp
@@ -358,6 +358,29 @@ public:
         m_rtree.clear();
         return m_pck.executeIndexed(std::forward<Args>(args)...);
     }
+
+    inline void preload(const PackGroup& pg) {
+        m_pconf.alignment = PConfig::Alignment::DONT_ALIGN;
+        m_pconf.object_function = nullptr; // drop the special objectfunction
+        m_pck.preload(pg);
+
+        // Build the rtree for queries to work
+        for(const ItemGroup& grp : pg)
+        for(unsigned idx = 0; idx < grp.size(); ++idx) {
+            Item& itm = grp[idx];
+            m_rtree.insert({itm.boundingBox(), idx});
+        }
+
+        m_pck.configure(m_pconf);
+    }
+
+    bool is_colliding(const Item& item) {
+        if(m_rtree.empty()) return false;
+        std::vector<SpatElement> result;
+        m_rtree.query(bgi::intersects(item.boundingBox()),
+                      std::back_inserter(result));
+        return !result.empty();
+    }
 };
 
 // Arranger specialization for a Box shaped bin.
@@ -365,8 +388,8 @@ template<> class AutoArranger<Box>: public _ArrBase<Box> {
 public:
 
     AutoArranger(const Box& bin, Distance dist,
-                 std::function<void(unsigned)> progressind,
-                 std::function<bool(void)> stopcond):
+                 std::function<void(unsigned)> progressind = [](unsigned){},
+                 std::function<bool(void)> stopcond = [](){return false;}):
         _ArrBase<Box>(bin, dist, progressind, stopcond)
     {
 
@@ -411,8 +434,8 @@ template<> class AutoArranger<lnCircle>: public _ArrBase<lnCircle> {
 public:
 
     AutoArranger(const lnCircle& bin, Distance dist,
-                 std::function<void(unsigned)> progressind,
-                 std::function<bool(void)> stopcond):
+                 std::function<void(unsigned)> progressind = [](unsigned){},
+                 std::function<bool(void)> stopcond = [](){return false;}):
         _ArrBase<lnCircle>(bin, dist, progressind, stopcond) {
 
         // As with the box, only the inside check is different.
@@ -456,8 +479,8 @@ public:
 template<> class AutoArranger<PolygonImpl>: public _ArrBase<PolygonImpl> {
 public:
     AutoArranger(const PolygonImpl& bin, Distance dist,
-                 std::function<void(unsigned)> progressind,
-                 std::function<bool(void)> stopcond):
+                 std::function<void(unsigned)> progressind = [](unsigned){},
+                 std::function<bool(void)> stopcond = [](){return false;}):
         _ArrBase<PolygonImpl>(bin, dist, progressind, stopcond)
     {
         m_pconf.object_function = [this, &bin] (const Item &item) {
@@ -791,5 +814,174 @@ bool arrange(Model &model,              // The model with the geometries
     return ret && result.size() == 1;
 }
 
+void find_new_position(const Model &model,
+                       ModelInstancePtrs toadd,
+                       coord_t min_obj_distance,
+                       const Polyline &bed)
+{
+    // Get the 2D projected shapes with their 3D model instance pointers
+    auto shapemap = arr::projectModelFromTop(model);
+
+    // Copy the references for the shapes only as the arranger expects a
+    // sequence of objects convertible to Item or ClipperPolygon
+    PackGroup preshapes; preshapes.emplace_back();
+    ItemGroup shapes;
+    preshapes.front().reserve(shapemap.size());
+
+    std::vector<ModelInstance*> shapes_ptr; shapes_ptr.reserve(toadd.size());
+    IndexedPackGroup result;
+
+    // If there is no hint about the shape, we will try to guess
+    BedShapeHint bedhint = bedShape(bed);
+
+    BoundingBox bbb(bed);
+
+    auto binbb = Box({
+                         static_cast<libnest2d::Coord>(bbb.min(0)),
+                         static_cast<libnest2d::Coord>(bbb.min(1))
+                     },
+                     {
+                         static_cast<libnest2d::Coord>(bbb.max(0)),
+                         static_cast<libnest2d::Coord>(bbb.max(1))
+                     });
+
+    for(auto it = shapemap.begin(); it != shapemap.end(); ++it) {
+        if(std::find(toadd.begin(), toadd.end(), it->first) == toadd.end()) {
+           if(it->second.isInside(binbb)) // just ignore items which are outside
+               preshapes.front().emplace_back(std::ref(it->second));
+        }
+        else {
+            shapes_ptr.emplace_back(it->first);
+            shapes.emplace_back(std::ref(it->second));
+        }
+    }
+
+    auto try_first_to_center = [&shapes, &shapes_ptr, &binbb]
+            (std::function<bool(const Item&)> is_colliding,
+             std::function<void(Item&)> preload)
+    {
+        // Try to put the first item to the center, as the arranger will not
+        // do this for us.
+        auto shptrit = shapes_ptr.begin();
+        for(auto shit = shapes.begin(); shit != shapes.end(); ++shit, ++shptrit)
+        {
+            // Try to place items to the center
+            Item& itm = *shit;
+            auto ibb = itm.boundingBox();
+            auto d = binbb.center() - ibb.center();
+            itm.translate(d);
+            if(!is_colliding(itm)) {
+                preload(itm);
+
+                auto offset = itm.translation();
+                Radians rot = itm.rotation();
+                ModelInstance *minst = *shptrit;
+                Vec3d foffset(offset.X*SCALING_FACTOR,
+                              offset.Y*SCALING_FACTOR,
+                              minst->get_offset()(Z));
+
+                // write the transformation data into the model instance
+                minst->set_rotation(Z, rot);
+                minst->set_offset(foffset);
+
+                shit = shapes.erase(shit);
+                shptrit = shapes_ptr.erase(shptrit);
+                break;
+            }
+        }
+    };
+
+    switch(bedhint.type) {
+    case BedShapeType::BOX: {
+
+        // Create the arranger for the box shaped bed
+        AutoArranger<Box> arrange(binbb, min_obj_distance);
+
+        if(!preshapes.front().empty()) { // If there is something on the plate
+            arrange.preload(preshapes);
+            try_first_to_center(
+                [&arrange](const Item& itm) {return arrange.is_colliding(itm);},
+                [&arrange](Item& itm) { arrange.preload({{itm}}); }
+            );
+        }
+
+        // Arrange and return the items with their respective indices within the
+        // input sequence.
+        result = arrange(shapes.begin(), shapes.end());
+        break;
+    }
+    case BedShapeType::CIRCLE: {
+
+        auto c = bedhint.shape.circ;
+        auto cc = to_lnCircle(c);
+
+        // Create the arranger for the box shaped bed
+        AutoArranger<lnCircle> arrange(cc, min_obj_distance);
+
+        if(!preshapes.front().empty()) { // If there is something on the plate
+            arrange.preload(preshapes);
+            try_first_to_center(
+                [&arrange](const Item& itm) {return arrange.is_colliding(itm);},
+                [&arrange](Item& itm) { arrange.preload({{itm}}); }
+            );
+        }
+
+        // Arrange and return the items with their respective indices within the
+        // input sequence.
+        result = arrange(shapes.begin(), shapes.end());
+        break;
+    }
+    case BedShapeType::IRREGULAR:
+    case BedShapeType::WHO_KNOWS: {
+        using P = libnest2d::PolygonImpl;
+
+        auto ctour = Slic3rMultiPoint_to_ClipperPath(bed);
+        P irrbed = sl::create<PolygonImpl>(std::move(ctour));
+
+        AutoArranger<P> arrange(irrbed, min_obj_distance);
+
+        if(!preshapes.front().empty()) { // If there is something on the plate
+            arrange.preload(preshapes);
+            try_first_to_center(
+                [&arrange](const Item& itm) {return arrange.is_colliding(itm);},
+                [&arrange](Item& itm) { arrange.preload({{itm}}); }
+            );
+        }
+
+        // Arrange and return the items with their respective indices within the
+        // input sequence.
+        result = arrange(shapes.begin(), shapes.end());
+        break;
+    }
+    };
+
+    // Now we go through the result which will contain the fixed and the moving
+    // polygons as well. We will have to search for our item.
+
+    const auto STRIDE_PADDING = 1.2;
+    Coord stride = Coord(STRIDE_PADDING*binbb.width()*SCALING_FACTOR);
+    Coord batch_offset = 0;
+
+    for(auto& group : result) {
+        for(auto& r : group) if(r.first < shapes.size()) {
+            Item& resultitem = r.second;
+            unsigned idx = r.first;
+            auto offset = resultitem.translation();
+            Radians rot = resultitem.rotation();
+            ModelInstance *minst = shapes_ptr[idx];
+            Vec3d foffset(offset.X*SCALING_FACTOR + batch_offset,
+                          offset.Y*SCALING_FACTOR,
+                          minst->get_offset()(Z));
+
+            // write the transformation data into the model instance
+            minst->set_rotation(Z, rot);
+            minst->set_offset(foffset);
+        }
+        batch_offset += stride;
+    }
 }
+
+}
+
+
 }
diff --git a/src/libslic3r/ModelArrange.hpp b/src/libslic3r/ModelArrange.hpp
index d62e0df30..d76769081 100644
--- a/src/libslic3r/ModelArrange.hpp
+++ b/src/libslic3r/ModelArrange.hpp
@@ -73,7 +73,13 @@ bool arrange(Model &model, coord_t min_obj_distance,
              std::function<void(unsigned)> progressind,
              std::function<bool(void)> stopcondition);
 
-}
+/// This will find a suitable position for a new object instance and leave the
+/// old items untouched.
+void find_new_position(const Model& model,
+                       ModelInstancePtrs instances_to_add,
+                       coord_t min_obj_distance,
+                       const Slic3r::Polyline& bed);
 
-}
+}   // arr
+}   // Slic3r
 #endif // MODELARRANGE_HPP
diff --git a/src/slic3r/GUI/Plater.cpp b/src/slic3r/GUI/Plater.cpp
index 683e354f5..cb08ea3f3 100644
--- a/src/slic3r/GUI/Plater.cpp
+++ b/src/slic3r/GUI/Plater.cpp
@@ -1473,24 +1473,19 @@ std::vector<size_t> Plater::priv::load_model_objects(const ModelObjectPtrs &mode
     const BoundingBoxf bed_shape = bed_shape_bb();
     const Vec3d bed_size = Slic3r::to_3d(bed_shape.size().cast<double>(), 1.0) - 2.0 * Vec3d::Ones();
 
-    bool need_arrange = false;
     bool scaled_down = false;
     std::vector<size_t> obj_idxs;
     unsigned int obj_count = model.objects.size();
 
+    ModelInstancePtrs new_instances;
     for (ModelObject *model_object : model_objects) {
         auto *object = model.add_object(*model_object);
         std::string object_name = object->name.empty() ? fs::path(object->input_file).filename().string() : object->name;
         obj_idxs.push_back(obj_count++);
 
         if (model_object->instances.empty()) {
-            // if object has no defined position(s) we need to rearrange everything after loading
-            need_arrange = true;
-
-            // add a default instance and center object around origin
-            object->center_around_origin();  // also aligns object to Z = 0
-            ModelInstance* instance = object->add_instance();
-            instance->set_offset(Slic3r::to_3d(bed_shape.center().cast<double>(), -object->origin_translation(2)));
+            object->center_around_origin();
+            new_instances.emplace_back(object->add_instance());
         }
 
         const Vec3d size = object->bounding_box().size();
@@ -1518,6 +1513,17 @@ std::vector<size_t> Plater::priv::load_model_objects(const ModelObjectPtrs &mode
         // print.add_model_object(object);
     }
 
+    // FIXME distance should be a config value /////////////////////////////////
+    auto min_obj_distance = static_cast<coord_t>(6/SCALING_FACTOR);
+    const auto *bed_shape_opt = config->opt<ConfigOptionPoints>("bed_shape");
+    assert(bed_shape_opt);
+    auto& bedpoints = bed_shape_opt->values;
+    Polyline bed; bed.points.reserve(bedpoints.size());
+    for(auto& v : bedpoints) bed.append(Point::new_scale(v(0), v(1)));
+
+    arr::find_new_position(model, new_instances, min_obj_distance, bed);
+    // /////////////////////////////////////////////////////////////////////////
+
     if (scaled_down) {
         GUI::show_info(q,
             _(L("Your object appears to be too large, so it was automatically scaled down to fit your print bed.")),