diff --git a/xs/src/libnest2d/CMakeLists.txt b/xs/src/libnest2d/CMakeLists.txt
index bfdb551fc..835e8311d 100644
--- a/xs/src/libnest2d/CMakeLists.txt
+++ b/xs/src/libnest2d/CMakeLists.txt
@@ -2,8 +2,6 @@ cmake_minimum_required(VERSION 2.8)
 
 project(Libnest2D)
 
-enable_testing()
-
 if(CMAKE_COMPILER_IS_GNUCC OR CMAKE_COMPILER_IS_GNUCXX)
     # Update if necessary
     set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wno-long-long ")
@@ -32,6 +30,7 @@ set(LIBNEST2D_SRCFILES
     ${CMAKE_CURRENT_SOURCE_DIR}/libnest2d/geometry_traits.hpp
     ${CMAKE_CURRENT_SOURCE_DIR}/libnest2d/common.hpp
     ${CMAKE_CURRENT_SOURCE_DIR}/libnest2d/optimizer.hpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/libnest2d/metaloop.hpp
     ${CMAKE_CURRENT_SOURCE_DIR}/libnest2d/placers/placer_boilerplate.hpp
     ${CMAKE_CURRENT_SOURCE_DIR}/libnest2d/placers/bottomleftplacer.hpp
     ${CMAKE_CURRENT_SOURCE_DIR}/libnest2d/placers/nfpplacer.hpp
@@ -60,8 +59,7 @@ if(LIBNEST2D_GEOMETRIES_BACKEND STREQUAL "clipper")
     include_directories(BEFORE ${CLIPPER_INCLUDE_DIRS})
     include_directories(${Boost_INCLUDE_DIRS})
 
-    list(APPEND LIBNEST2D_SRCFILES ${CMAKE_CURRENT_SOURCE_DIR}/libnest2d/clipper_backend/clipper_backend.cpp
-                                   ${CMAKE_CURRENT_SOURCE_DIR}/libnest2d/clipper_backend/clipper_backend.hpp
+    list(APPEND LIBNEST2D_SRCFILES ${CMAKE_CURRENT_SOURCE_DIR}/libnest2d/clipper_backend/clipper_backend.hpp
                                    ${CMAKE_CURRENT_SOURCE_DIR}/libnest2d/boost_alg.hpp)
     list(APPEND LIBNEST2D_LIBRARIES ${CLIPPER_LIBRARIES})
     list(APPEND LIBNEST2D_HEADERS ${CLIPPER_INCLUDE_DIRS}
@@ -81,22 +79,12 @@ if(LIBNEST2D_OPTIMIZER_BACKEND STREQUAL "nlopt")
                                     ${CMAKE_CURRENT_SOURCE_DIR}/libnest2d/optimizers/subplex.hpp
                                     ${CMAKE_CURRENT_SOURCE_DIR}/libnest2d/optimizers/genetic.hpp
                                     ${CMAKE_CURRENT_SOURCE_DIR}/libnest2d/optimizers/nlopt_boilerplate.hpp)
-    list(APPEND LIBNEST2D_LIBRARIES ${NLopt_LIBS}
-#                                    Threads::Threads
-    )
+    list(APPEND LIBNEST2D_LIBRARIES ${NLopt_LIBS})
     list(APPEND LIBNEST2D_HEADERS ${NLopt_INCLUDE_DIR})
 endif()
 
-# Currently we are outsourcing the non-convex NFP implementation from
-# libnfporb and it needs libgmp to work
-#find_package(GMP)
-#if(GMP_FOUND)
-#    list(APPEND LIBNEST2D_LIBRARIES ${GMP_LIBRARIES})
-#    list(APPEND LIBNEST2D_HEADERS ${GMP_INCLUDE_DIR})
-#    add_definitions(-DLIBNFP_USE_RATIONAL)
-#endif()
-
 if(LIBNEST2D_UNITTESTS)
+    enable_testing()
     add_subdirectory(tests)
 endif()
 
diff --git a/xs/src/libnest2d/cmake_modules/FindGMP.cmake b/xs/src/libnest2d/cmake_modules/FindGMP.cmake
deleted file mode 100644
index db173bc90..000000000
--- a/xs/src/libnest2d/cmake_modules/FindGMP.cmake
+++ /dev/null
@@ -1,35 +0,0 @@
-# Try to find the GMP libraries:
-# GMP_FOUND - System has GMP lib
-# GMP_INCLUDE_DIR - The GMP include directory
-# GMP_LIBRARIES - Libraries needed to use GMP
-
-if (GMP_INCLUDE_DIR AND GMP_LIBRARIES)
-	# Force search at every time, in case configuration changes
-	unset(GMP_INCLUDE_DIR CACHE)
-	unset(GMP_LIBRARIES CACHE)
-endif (GMP_INCLUDE_DIR AND GMP_LIBRARIES)
-
-find_path(GMP_INCLUDE_DIR NAMES gmp.h)
-
-if(WIN32)
-	find_library(GMP_LIBRARIES NAMES libgmp.a gmp gmp.lib mpir mpir.lib)
-else(WIN32)
-	if(STBIN)
-		message(STATUS "STBIN: ${STBIN}")
-		find_library(GMP_LIBRARIES NAMES libgmp.a gmp)
-	else(STBIN)
-		find_library(GMP_LIBRARIES NAMES libgmp.so gmp)
-	endif(STBIN)
-endif(WIN32)
-
-if(GMP_INCLUDE_DIR AND GMP_LIBRARIES)
-   set(GMP_FOUND TRUE)
-endif(GMP_INCLUDE_DIR AND GMP_LIBRARIES)
-
-if(GMP_FOUND)
-	message(STATUS "Configured GMP: ${GMP_LIBRARIES}")
-else(GMP_FOUND)
-	message(STATUS "Could NOT find GMP")
-endif(GMP_FOUND)
-
-mark_as_advanced(GMP_INCLUDE_DIR GMP_LIBRARIES)
\ No newline at end of file
diff --git a/xs/src/libnest2d/examples/main.cpp b/xs/src/libnest2d/examples/main.cpp
index e5a47161e..a97618578 100644
--- a/xs/src/libnest2d/examples/main.cpp
+++ b/xs/src/libnest2d/examples/main.cpp
@@ -535,17 +535,18 @@ void arrangeRectangles() {
     proba[0].rotate(Pi/3);
     proba[1].rotate(Pi-Pi/3);
 
+//    std::vector<Item> input(25, Rectangle(70*SCALE, 10*SCALE));
     std::vector<Item> input;
     input.insert(input.end(), prusaParts().begin(), prusaParts().end());
 //    input.insert(input.end(), prusaExParts().begin(), prusaExParts().end());
-    input.insert(input.end(), stegoParts().begin(), stegoParts().end());
+//    input.insert(input.end(), stegoParts().begin(), stegoParts().end());
 //    input.insert(input.end(), rects.begin(), rects.end());
-    input.insert(input.end(), proba.begin(), proba.end());
+//    input.insert(input.end(), proba.begin(), proba.end());
 //    input.insert(input.end(), crasher.begin(), crasher.end());
 
     Box bin(250*SCALE, 210*SCALE);
 
-    Coord min_obj_distance = 6*SCALE;
+    auto min_obj_distance = static_cast<Coord>(0*SCALE);
 
     using Placer = NfpPlacer;
     using Packer = Arranger<Placer, FirstFitSelection>;
@@ -554,21 +555,45 @@ void arrangeRectangles() {
 
     Packer::PlacementConfig pconf;
     pconf.alignment = Placer::Config::Alignment::CENTER;
-    pconf.starting_point = Placer::Config::Alignment::CENTER;
+    pconf.starting_point = Placer::Config::Alignment::BOTTOM_LEFT;
     pconf.rotations = {0.0/*, Pi/2.0, Pi, 3*Pi/2*/};
-    pconf.object_function = [&bin](Placer::Pile pile, double area,
-                               double norm, double penality) {
+
+    double norm_2 = std::nan("");
+    pconf.object_function = [&bin, &norm_2](Placer::Pile pile, const Item& item,
+            double /*area*/, double norm, double penality) {
+
+        using pl = PointLike;
 
         auto bb = ShapeLike::boundingBox(pile);
+        auto ibb = item.boundingBox();
+        auto minc = ibb.minCorner();
+        auto maxc = ibb.maxCorner();
 
-        auto& sh = pile.back();
-        auto rv = Nfp::referenceVertex(sh);
-        auto c = bin.center();
-        auto d = PointLike::distance(rv, c);
-        double score = double(d)/norm;
+        if(std::isnan(norm_2)) norm_2 = pow(norm, 2);
+
+        // We get the distance of the reference point from the center of the
+        // heat bed
+        auto cc = bb.center();
+        auto top_left = PointImpl{getX(minc), getY(maxc)};
+        auto bottom_right = PointImpl{getX(maxc), getY(minc)};
+
+        auto a = pl::distance(ibb.maxCorner(), cc);
+        auto b = pl::distance(ibb.minCorner(), cc);
+        auto c = pl::distance(ibb.center(), cc);
+        auto d = pl::distance(top_left, cc);
+        auto e = pl::distance(bottom_right, cc);
+
+        auto area = bb.width() * bb.height() / norm_2;
+
+        auto min_dist = std::min({a, b, c, d, e}) / norm;
+
+        // The score will be the normalized distance which will be minimized,
+        // effectively creating a circle shaped pile of items
+        double score = 0.8*min_dist  + 0.2*area;
 
         // If it does not fit into the print bed we will beat it
-        // with a large penality
+        // with a large penality. If we would not do this, there would be only
+        // one big pile that doesn't care whether it fits onto the print bed.
         if(!NfpPlacer::wouldFit(bb, bin)) score = 2*penality - score;
 
         return score;
@@ -577,7 +602,7 @@ void arrangeRectangles() {
     Packer::SelectionConfig sconf;
 //    sconf.allow_parallel = false;
 //    sconf.force_parallel = false;
-//    sconf.try_triplets = false;
+//    sconf.try_triplets = true;
 //    sconf.try_reverse_order = true;
 //    sconf.waste_increment = 0.005;
 
@@ -630,7 +655,7 @@ void arrangeRectangles() {
               << " %" << std::endl;
 
     std::cout << "Bin usage: (";
-    unsigned total = 0;
+    size_t total = 0;
     for(auto& r : result) { std::cout << r.size() << " "; total += r.size(); }
     std::cout << ") Total: " << total << std::endl;
 
@@ -643,9 +668,11 @@ void arrangeRectangles() {
                                         << input.size() - total << " elements!"
                                         << std::endl;
 
-    svg::SVGWriter::Config conf;
+    using SVGWriter = svg::SVGWriter<PolygonImpl>;
+
+    SVGWriter::Config conf;
     conf.mm_in_coord_units = SCALE;
-    svg::SVGWriter svgw(conf);
+    SVGWriter svgw(conf);
     svgw.setSize(bin);
     svgw.writePackGroup(result);
 //    std::for_each(input.begin(), input.end(), [&svgw](Item& item){ svgw.writeItem(item);});
diff --git a/xs/src/libnest2d/libnest2d.h b/xs/src/libnest2d/libnest2d.h
index 1e0a98f6a..e0ad05c41 100644
--- a/xs/src/libnest2d/libnest2d.h
+++ b/xs/src/libnest2d/libnest2d.h
@@ -6,7 +6,7 @@
 #include <libnest2d/clipper_backend/clipper_backend.hpp>
 
 // We include the stock optimizers for local and global optimization
-#include <libnest2d/optimizers/simplex.hpp>     // Local subplex for NfpPlacer
+#include <libnest2d/optimizers/simplex.hpp>     // Local simplex for NfpPlacer
 #include <libnest2d/optimizers/genetic.hpp>     // Genetic for min. bounding box
 
 #include <libnest2d/libnest2d.hpp>
diff --git a/xs/src/libnest2d/libnest2d/boost_alg.hpp b/xs/src/libnest2d/libnest2d/boost_alg.hpp
index a50b397d3..422616d20 100644
--- a/xs/src/libnest2d/libnest2d/boost_alg.hpp
+++ b/xs/src/libnest2d/libnest2d/boost_alg.hpp
@@ -8,8 +8,16 @@
 #ifdef __clang__
 #undef _MSC_EXTENSIONS
 #endif
-#include <boost/geometry.hpp>
 
+#ifdef _MSC_VER
+#pragma warning(push)
+#pragma warning(disable: 4244)
+#pragma warning(disable: 4267)
+#endif
+#include <boost/geometry.hpp>
+#ifdef _MSC_VER
+#pragma warning(pop)
+#endif
 // this should be removed to not confuse the compiler
 // #include <libnest2d.h>
 
@@ -461,15 +469,6 @@ inline bp2d::Shapes Nfp::merge(const bp2d::Shapes& shapes,
 }
 #endif
 
-//#ifndef DISABLE_BOOST_MINKOWSKI_ADD
-//template<>
-//inline PolygonImpl& Nfp::minkowskiAdd(PolygonImpl& sh,
-//                                      const PolygonImpl& /*other*/)
-//{
-//    return sh;
-//}
-//#endif
-
 #ifndef DISABLE_BOOST_SERIALIZE
 template<> inline std::string ShapeLike::serialize<libnest2d::Formats::SVG>(
         const PolygonImpl& sh, double scale)
diff --git a/xs/src/libnest2d/libnest2d/clipper_backend/clipper_backend.cpp b/xs/src/libnest2d/libnest2d/clipper_backend/clipper_backend.cpp
deleted file mode 100644
index 830d235a3..000000000
--- a/xs/src/libnest2d/libnest2d/clipper_backend/clipper_backend.cpp
+++ /dev/null
@@ -1,58 +0,0 @@
-//#include "clipper_backend.hpp"
-//#include <atomic>
-
-//namespace libnest2d {
-
-//namespace  {
-
-//class SpinLock {
-//    std::atomic_flag& lck_;
-//public:
-
-//    inline SpinLock(std::atomic_flag& flg): lck_(flg) {}
-
-//    inline void lock() {
-//        while(lck_.test_and_set(std::memory_order_acquire)) {}
-//    }
-
-//    inline void unlock() { lck_.clear(std::memory_order_release); }
-//};
-
-//class HoleCache {
-//    friend struct libnest2d::ShapeLike;
-
-//    std::unordered_map< const PolygonImpl*, ClipperLib::Paths> map;
-
-//    ClipperLib::Paths& _getHoles(const PolygonImpl* p) {
-//        static std::atomic_flag flg = ATOMIC_FLAG_INIT;
-//        SpinLock lock(flg);
-
-//        lock.lock();
-//        ClipperLib::Paths& paths = map[p];
-//        lock.unlock();
-
-//        if(paths.size() != p->Childs.size()) {
-//            paths.reserve(p->Childs.size());
-
-//            for(auto np : p->Childs) {
-//                paths.emplace_back(np->Contour);
-//            }
-//        }
-
-//        return paths;
-//    }
-
-//    ClipperLib::Paths& getHoles(PolygonImpl& p) {
-//        return _getHoles(&p);
-//    }
-
-//    const ClipperLib::Paths& getHoles(const PolygonImpl& p) {
-//        return _getHoles(&p);
-//    }
-//};
-//}
-
-//HoleCache holeCache;
-
-//}
-
diff --git a/xs/src/libnest2d/libnest2d/clipper_backend/clipper_backend.hpp b/xs/src/libnest2d/libnest2d/clipper_backend/clipper_backend.hpp
index 8cc27573a..15ceb1576 100644
--- a/xs/src/libnest2d/libnest2d/clipper_backend/clipper_backend.hpp
+++ b/xs/src/libnest2d/libnest2d/clipper_backend/clipper_backend.hpp
@@ -21,7 +21,7 @@ struct PolygonImpl {
     PathImpl Contour;
     HoleStore Holes;
 
-    inline PolygonImpl() {}
+    inline PolygonImpl() = default;
 
     inline explicit PolygonImpl(const PathImpl& cont): Contour(cont) {}
     inline explicit PolygonImpl(const HoleStore& holes):
@@ -66,6 +66,19 @@ inline PointImpl operator-(const PointImpl& p1, const PointImpl& p2) {
     ret -= p2;
     return ret;
 }
+
+inline PointImpl& operator *=(PointImpl& p, const PointImpl& pa ) {
+    p.X *= pa.X;
+    p.Y *= pa.Y;
+    return p;
+}
+
+inline PointImpl operator*(const PointImpl& p1, const PointImpl& p2) {
+    PointImpl ret = p1;
+    ret *= p2;
+    return ret;
+}
+
 }
 
 namespace libnest2d {
@@ -135,7 +148,7 @@ inline void ShapeLike::reserve(PolygonImpl& sh, size_t vertex_capacity)
 
 namespace _smartarea {
 template<Orientation o>
-inline double area(const PolygonImpl& sh) {
+inline double area(const PolygonImpl& /*sh*/) {
     return std::nan("");
 }
 
@@ -220,22 +233,6 @@ inline void ShapeLike::offset(PolygonImpl& sh, TCoord<PointImpl> distance) {
     }
 }
 
-//template<> // TODO make it support holes if this method will ever be needed.
-//inline PolygonImpl Nfp::minkowskiDiff(const PolygonImpl& sh,
-//                                      const PolygonImpl& other)
-//{
-//    #define DISABLE_BOOST_MINKOWSKI_ADD
-
-//    ClipperLib::Paths solution;
-
-//    ClipperLib::MinkowskiDiff(sh.Contour, other.Contour, solution);
-
-//    PolygonImpl ret;
-//    ret.Contour = solution.front();
-
-//    return sh;
-//}
-
 // Tell libnest2d how to make string out of a ClipperPolygon object
 template<> inline std::string ShapeLike::toString(const PolygonImpl& sh) {
     std::stringstream ss;
@@ -406,35 +403,12 @@ inline void ShapeLike::rotate(PolygonImpl& sh, const Radians& rads)
 }
 
 #define DISABLE_BOOST_NFP_MERGE
-template<> inline Nfp::Shapes<PolygonImpl>
-Nfp::merge(const Nfp::Shapes<PolygonImpl>& shapes, const PolygonImpl& sh)
-{
+inline Nfp::Shapes<PolygonImpl> _merge(ClipperLib::Clipper& clipper) {
     Nfp::Shapes<PolygonImpl> retv;
 
-    ClipperLib::Clipper clipper(ClipperLib::ioReverseSolution);
-
-    bool closed = true;
-    bool valid = false;
-
-    valid = clipper.AddPath(sh.Contour, ClipperLib::ptSubject, closed);
-
-    for(auto& hole : sh.Holes) {
-        valid &= clipper.AddPath(hole, ClipperLib::ptSubject, closed);
-    }
-
-    for(auto& path : shapes) {
-        valid &= clipper.AddPath(path.Contour, ClipperLib::ptSubject, closed);
-
-        for(auto& hole : path.Holes) {
-            valid &= clipper.AddPath(hole, ClipperLib::ptSubject, closed);
-        }
-    }
-
-    if(!valid) throw GeometryException(GeomErr::MERGE);
-
     ClipperLib::PolyTree result;
-    clipper.Execute(ClipperLib::ctUnion, result, ClipperLib::pftNonZero);
-    retv.reserve(result.Total());
+    clipper.Execute(ClipperLib::ctUnion, result, ClipperLib::pftNegative);
+    retv.reserve(static_cast<size_t>(result.Total()));
 
     std::function<void(ClipperLib::PolyNode*, PolygonImpl&)> processHole;
 
@@ -445,7 +419,8 @@ Nfp::merge(const Nfp::Shapes<PolygonImpl>& shapes, const PolygonImpl& sh)
         retv.push_back(poly);
     };
 
-    processHole = [&processPoly](ClipperLib::PolyNode *pptr, PolygonImpl& poly) {
+    processHole = [&processPoly](ClipperLib::PolyNode *pptr, PolygonImpl& poly)
+    {
         poly.Holes.push_back(pptr->Contour);
         poly.Holes.back().push_back(poly.Holes.back().front());
         for(auto c : pptr->Childs) processPoly(c);
@@ -463,6 +438,27 @@ Nfp::merge(const Nfp::Shapes<PolygonImpl>& shapes, const PolygonImpl& sh)
     return retv;
 }
 
+template<> inline Nfp::Shapes<PolygonImpl>
+Nfp::merge(const Nfp::Shapes<PolygonImpl>& shapes)
+{
+    ClipperLib::Clipper clipper(ClipperLib::ioReverseSolution);
+
+    bool closed = true;
+    bool valid = true;
+
+    for(auto& path : shapes) {
+        valid &= clipper.AddPath(path.Contour, ClipperLib::ptSubject, closed);
+
+        for(auto& hole : path.Holes) {
+            valid &= clipper.AddPath(hole, ClipperLib::ptSubject, closed);
+        }
+    }
+
+    if(!valid) throw GeometryException(GeomErr::MERGE);
+
+    return _merge(clipper);
+}
+
 }
 
 //#define DISABLE_BOOST_SERIALIZE
diff --git a/xs/src/libnest2d/libnest2d/common.hpp b/xs/src/libnest2d/libnest2d/common.hpp
index 18f313712..6867f76f3 100644
--- a/xs/src/libnest2d/libnest2d/common.hpp
+++ b/xs/src/libnest2d/libnest2d/common.hpp
@@ -13,6 +13,7 @@
 #if defined(_MSC_VER) &&  _MSC_VER <= 1800 || __cplusplus < 201103L
     #define BP2D_NOEXCEPT
     #define BP2D_CONSTEXPR
+    #define BP2D_COMPILER_MSVC12
 #elif __cplusplus >= 201103L
     #define BP2D_NOEXCEPT noexcept
     #define BP2D_CONSTEXPR constexpr
@@ -84,44 +85,6 @@ struct invoke_result {
 template<class F, class...Args>
 using invoke_result_t = typename invoke_result<F, Args...>::type;
 
-/* ************************************************************************** */
-/* C++14 std::index_sequence implementation:                                  */
-/* ************************************************************************** */
-
-/**
- * \brief C++11 conformant implementation of the index_sequence type from C++14
- */
-template<size_t...Ints> struct index_sequence {
-    using value_type = size_t;
-    BP2D_CONSTEXPR value_type size() const { return sizeof...(Ints); }
-};
-
-// A Help structure to generate the integer list
-template<size_t...Nseq> struct genSeq;
-
-// Recursive template to generate the list
-template<size_t I, size_t...Nseq> struct genSeq<I, Nseq...> {
-    // Type will contain a genSeq with Nseq appended by one element
-    using Type = typename genSeq< I - 1, I - 1, Nseq...>::Type;
-};
-
-// Terminating recursion
-template <size_t ... Nseq> struct genSeq<0, Nseq...> {
-    // If I is zero, Type will contain index_sequence with the fuly generated
-    // integer list.
-    using Type = index_sequence<Nseq...>;
-};
-
-/// Helper alias to make an index sequence from 0 to N
-template<size_t N> using make_index_sequence = typename genSeq<N>::Type;
-
-/// Helper alias to make an index sequence for a parameter pack
-template<class...Args>
-using index_sequence_for = make_index_sequence<sizeof...(Args)>;
-
-
-/* ************************************************************************** */
-
 /**
  * A useful little tool for triggering static_assert error messages e.g. when
  * a mandatory template specialization (implementation) is missing.
@@ -229,7 +192,7 @@ public:
 
     GeomErr errcode() const { return errcode_; }
 
-    virtual const char * what() const BP2D_NOEXCEPT override {
+    const char * what() const BP2D_NOEXCEPT override {
         return errorstr(errcode_).c_str();
     }
 };
diff --git a/xs/src/libnest2d/libnest2d/geometry_traits.hpp b/xs/src/libnest2d/libnest2d/geometry_traits.hpp
index dbd609201..568c0a766 100644
--- a/xs/src/libnest2d/libnest2d/geometry_traits.hpp
+++ b/xs/src/libnest2d/libnest2d/geometry_traits.hpp
@@ -68,7 +68,7 @@ class _Box: PointPair<RawPoint> {
     using PointPair<RawPoint>::p2;
 public:
 
-    inline _Box() {}
+    inline _Box() = default;
     inline _Box(const RawPoint& p, const RawPoint& pp):
         PointPair<RawPoint>({p, pp}) {}
 
@@ -97,7 +97,7 @@ class _Segment: PointPair<RawPoint> {
     mutable Radians angletox_ = std::nan("");
 public:
 
-    inline _Segment() {}
+    inline _Segment() = default;
 
     inline _Segment(const RawPoint& p, const RawPoint& pp):
         PointPair<RawPoint>({p, pp}) {}
@@ -188,7 +188,7 @@ struct PointLike {
 
         if( (y < y1 && y < y2) || (y > y1 && y > y2) )
             return {0, false};
-        else if ((y == y1 && y == y2) && (x > x1 && x > x2))
+        if ((y == y1 && y == y2) && (x > x1 && x > x2))
             ret = std::min( x-x1, x -x2);
         else if( (y == y1 && y == y2) && (x < x1 && x < x2))
             ret = -std::min(x1 - x, x2 - x);
@@ -214,7 +214,7 @@ struct PointLike {
 
         if( (x < x1 && x < x2) || (x > x1 && x > x2) )
             return {0, false};
-        else if ((x == x1 && x == x2) && (y > y1 && y > y2))
+        if ((x == x1 && x == x2) && (y > y1 && y > y2))
             ret = std::min( y-y1, y -y2);
         else if( (x == x1 && x == x2) && (y < y1 && y < y2))
             ret = -std::min(y1 - y, y2 - y);
@@ -329,7 +329,7 @@ enum class Formats {
 };
 
 // This struct serves as a namespace. The only difference is that it can be
-// used in friend declarations.
+// used in friend declarations and can be aliased at class scope.
 struct ShapeLike {
 
     template<class RawShape>
@@ -361,6 +361,51 @@ struct ShapeLike {
         return create<RawShape>(contour, {});
     }
 
+    template<class RawShape>
+    static THolesContainer<RawShape>& holes(RawShape& /*sh*/)
+    {
+        static THolesContainer<RawShape> empty;
+        return empty;
+    }
+
+    template<class RawShape>
+    static const THolesContainer<RawShape>& holes(const RawShape& /*sh*/)
+    {
+        static THolesContainer<RawShape> empty;
+        return empty;
+    }
+
+    template<class RawShape>
+    static TContour<RawShape>& getHole(RawShape& sh, unsigned long idx)
+    {
+        return holes(sh)[idx];
+    }
+
+    template<class RawShape>
+    static const TContour<RawShape>& getHole(const RawShape& sh,
+                                              unsigned long idx)
+    {
+        return holes(sh)[idx];
+    }
+
+    template<class RawShape>
+    static size_t holeCount(const RawShape& sh)
+    {
+        return holes(sh).size();
+    }
+
+    template<class RawShape>
+    static TContour<RawShape>& getContour(RawShape& sh)
+    {
+        return sh;
+    }
+
+    template<class RawShape>
+    static const TContour<RawShape>& getContour(const RawShape& sh)
+    {
+        return sh;
+    }
+
     // Optional, does nothing by default
     template<class RawShape>
     static void reserve(RawShape& /*sh*/,  size_t /*vertex_capacity*/) {}
@@ -402,7 +447,7 @@ struct ShapeLike {
     }
 
     template<Formats, class RawShape>
-    static std::string serialize(const RawShape& /*sh*/, double scale=1)
+    static std::string serialize(const RawShape& /*sh*/, double /*scale*/=1)
     {
         static_assert(always_false<RawShape>::value,
                       "ShapeLike::serialize() unimplemented!");
@@ -498,51 +543,6 @@ struct ShapeLike {
         return RawShape();
     }
 
-    template<class RawShape>
-    static THolesContainer<RawShape>& holes(RawShape& /*sh*/)
-    {
-        static THolesContainer<RawShape> empty;
-        return empty;
-    }
-
-    template<class RawShape>
-    static const THolesContainer<RawShape>& holes(const RawShape& /*sh*/)
-    {
-        static THolesContainer<RawShape> empty;
-        return empty;
-    }
-
-    template<class RawShape>
-    static TContour<RawShape>& getHole(RawShape& sh, unsigned long idx)
-    {
-        return holes(sh)[idx];
-    }
-
-    template<class RawShape>
-    static const TContour<RawShape>& getHole(const RawShape& sh,
-                                              unsigned long idx)
-    {
-        return holes(sh)[idx];
-    }
-
-    template<class RawShape>
-    static size_t holeCount(const RawShape& sh)
-    {
-        return holes(sh).size();
-    }
-
-    template<class RawShape>
-    static TContour<RawShape>& getContour(RawShape& sh)
-    {
-        return sh;
-    }
-
-    template<class RawShape>
-    static const TContour<RawShape>& getContour(const RawShape& sh)
-    {
-        return sh;
-    }
-
     template<class RawShape>
     static void rotate(RawShape& /*sh*/, const Radians& /*rads*/)
     {
@@ -621,14 +621,12 @@ struct ShapeLike {
     }
 
     template<class RawShape>
-    static double area(const Shapes<RawShape>& shapes)
+    static inline double area(const Shapes<RawShape>& shapes)
     {
-        double ret = 0;
-        std::accumulate(shapes.first(), shapes.end(),
-                        [](const RawShape& a, const RawShape& b) {
-            return area(a) + area(b);
+        return std::accumulate(shapes.begin(), shapes.end(), 0.0,
+                        [](double a, const RawShape& b) {
+            return a += area(b);
         });
-        return ret;
     }
 
     template<class RawShape> // Potential O(1) implementation may exist
diff --git a/xs/src/libnest2d/libnest2d/geometry_traits_nfp.hpp b/xs/src/libnest2d/libnest2d/geometry_traits_nfp.hpp
index 581b6bed0..90cf21be5 100644
--- a/xs/src/libnest2d/libnest2d/geometry_traits_nfp.hpp
+++ b/xs/src/libnest2d/libnest2d/geometry_traits_nfp.hpp
@@ -3,7 +3,9 @@
 
 #include "geometry_traits.hpp"
 #include <algorithm>
+#include <functional>
 #include <vector>
+#include <iterator>
 
 namespace libnest2d {
 
@@ -23,64 +25,22 @@ struct Nfp {
 template<class RawShape>
 using Shapes = typename ShapeLike::Shapes<RawShape>;
 
-/// Minkowski addition (not used yet)
+/**
+ * Merge a bunch of polygons with the specified additional polygon.
+ *
+ * \tparam RawShape the Polygon data type.
+ * \param shc The pile of polygons that will be unified with sh.
+ * \param sh A single polygon to unify with shc.
+ *
+ * \return A set of polygons that is the union of the input polygons. Note that
+ * mostly it will be a set containing only one big polygon but if the input
+ * polygons are disjuct than the resulting set will contain more polygons.
+ */
 template<class RawShape>
-static RawShape minkowskiDiff(const RawShape& sh, const RawShape& cother)
+static Shapes<RawShape> merge(const Shapes<RawShape>& /*shc*/)
 {
-    using Vertex = TPoint<RawShape>;
-    //using Coord = TCoord<Vertex>;
-    using Edge = _Segment<Vertex>;
-    using sl = ShapeLike;
-    using std::signbit;
-
-    // Copy the orbiter (controur only), we will have to work on it
-    RawShape orbiter = sl::create(sl::getContour(cother));
-
-    // Make the orbiter reverse oriented
-    for(auto &v : sl::getContour(orbiter)) v = -v;
-
-    // An egde with additional data for marking it
-    struct MarkedEdge { Edge e; Radians turn_angle; bool is_turning_point; };
-
-    // Container for marked edges
-    using EdgeList = std::vector<MarkedEdge>;
-
-    EdgeList A, B;
-
-    auto fillEdgeList = [](EdgeList& L, const RawShape& poly) {
-        L.reserve(sl::contourVertexCount(poly));
-
-        auto it = sl::cbegin(poly);
-        auto nextit = std::next(it);
-
-        L.emplace_back({Edge(*it, *nextit), 0, false});
-        it++; nextit++;
-
-        while(nextit != sl::cend(poly)) {
-            Edge e(*it, *nextit);
-            auto& L_prev = L.back();
-            auto phi = L_prev.e.angleToXaxis();
-            auto phi_prev = e.angleToXaxis();
-            auto turn_angle = phi-phi_prev;
-            if(turn_angle > Pi) turn_angle -= 2*Pi;
-            L.emplace_back({
-                              e,
-                              turn_angle,
-                              signbit(turn_angle) != signbit(L_prev.turn_angle)
-                           });
-            it++; nextit++;
-        }
-
-        L.front().turn_angle = L.front().e.angleToXaxis() -
-                               L.back().e.angleToXaxis();
-
-        if(L.front().turn_angle > Pi) L.front().turn_angle -= 2*Pi;
-    };
-
-    fillEdgeList(A, sh);
-    fillEdgeList(B, orbiter);
-
-    return sh;
+    static_assert(always_false<RawShape>::value,
+                  "Nfp::merge(shapes, shape) unimplemented!");
 }
 
 /**
@@ -95,10 +55,12 @@ static RawShape minkowskiDiff(const RawShape& sh, const RawShape& cother)
  * polygons are disjuct than the resulting set will contain more polygons.
  */
 template<class RawShape>
-static Shapes<RawShape> merge(const Shapes<RawShape>& shc, const RawShape& sh)
+static Shapes<RawShape> merge(const Shapes<RawShape>& shc,
+                              const RawShape& sh)
 {
-    static_assert(always_false<RawShape>::value,
-                  "Nfp::merge(shapes, shape) unimplemented!");
+    auto m = merge(shc);
+    m.push_back(sh);
+    return merge(m);
 }
 
 /**
@@ -139,16 +101,20 @@ template<class RawShape>
 static TPoint<RawShape> rightmostUpVertex(const RawShape& sh)
 {
 
-    // find min x and min y vertex
+    // find max x and max y vertex
     auto it = std::max_element(ShapeLike::cbegin(sh), ShapeLike::cend(sh),
                                _vsort<RawShape>);
 
     return *it;
 }
 
+template<class RawShape>
+using NfpResult = std::pair<RawShape, TPoint<RawShape>>;
+
 /// Helper function to get the NFP
 template<NfpLevel nfptype, class RawShape>
-static RawShape noFitPolygon(const RawShape& sh, const RawShape& other)
+static NfpResult<RawShape> noFitPolygon(const RawShape& sh,
+                                        const RawShape& other)
 {
     NfpImpl<RawShape, nfptype> nfp;
     return nfp(sh, other);
@@ -167,44 +133,46 @@ static RawShape noFitPolygon(const RawShape& sh, const RawShape& other)
  * \tparam RawShape the Polygon data type.
  * \param sh The stationary polygon
  * \param cother The orbiting polygon
- * \return Returns the NFP of the two input polygons which have to be strictly
- * convex. The resulting NFP is proven to be convex as well in this case.
+ * \return Returns a pair of the NFP and its reference vertex of the two input
+ * polygons which have to be strictly convex. The resulting NFP is proven to be
+ * convex as well in this case.
  *
  */
 template<class RawShape>
-static RawShape nfpConvexOnly(const RawShape& sh, const RawShape& cother)
+static NfpResult<RawShape> nfpConvexOnly(const RawShape& sh,
+                                         const RawShape& other)
 {
     using Vertex = TPoint<RawShape>; using Edge = _Segment<Vertex>;
-
-    RawShape other = cother;
-
-    // Make the other polygon counter-clockwise
-    std::reverse(ShapeLike::begin(other), ShapeLike::end(other));
+    using sl = ShapeLike;
 
     RawShape rsh;   // Final nfp placeholder
+    Vertex top_nfp;
     std::vector<Edge> edgelist;
 
-    auto cap = ShapeLike::contourVertexCount(sh) +
-            ShapeLike::contourVertexCount(other);
+    auto cap = sl::contourVertexCount(sh) + sl::contourVertexCount(other);
 
     // Reserve the needed memory
     edgelist.reserve(cap);
-    ShapeLike::reserve(rsh, static_cast<unsigned long>(cap));
+    sl::reserve(rsh, static_cast<unsigned long>(cap));
 
     { // place all edges from sh into edgelist
-        auto first = ShapeLike::cbegin(sh);
-        auto next = first + 1;
-        auto endit = ShapeLike::cend(sh);
+        auto first = sl::cbegin(sh);
+        auto next = std::next(first);
 
-        while(next != endit) edgelist.emplace_back(*(first++), *(next++));
+        while(next != sl::cend(sh)) {
+            edgelist.emplace_back(*(first), *(next));
+            ++first; ++next;
+        }
     }
 
     { // place all edges from other into edgelist
-        auto first = ShapeLike::cbegin(other);
-        auto next = first + 1;
-        auto endit = ShapeLike::cend(other);
+        auto first = sl::cbegin(other);
+        auto next = std::next(first);
 
-        while(next != endit) edgelist.emplace_back(*(first++), *(next++));
+        while(next != sl::cend(other)) {
+            edgelist.emplace_back(*(next), *(first));
+            ++first; ++next;
+        }
     }
 
     // Sort the edges by angle to X axis.
@@ -215,10 +183,16 @@ static RawShape nfpConvexOnly(const RawShape& sh, const RawShape& cother)
     });
 
     // Add the two vertices from the first edge into the final polygon.
-    ShapeLike::addVertex(rsh, edgelist.front().first());
-    ShapeLike::addVertex(rsh, edgelist.front().second());
+    sl::addVertex(rsh, edgelist.front().first());
+    sl::addVertex(rsh, edgelist.front().second());
 
-    auto tmp = std::next(ShapeLike::begin(rsh));
+    // Sorting function for the nfp reference vertex search
+    auto& cmp = _vsort<RawShape>;
+
+    // the reference (rightmost top) vertex so far
+    top_nfp = *std::max_element(sl::cbegin(rsh), sl::cend(rsh), cmp );
+
+    auto tmp = std::next(sl::begin(rsh));
 
     // Construct final nfp by placing each edge to the end of the previous
     for(auto eit = std::next(edgelist.begin());
@@ -226,56 +200,325 @@ static RawShape nfpConvexOnly(const RawShape& sh, const RawShape& cother)
         ++eit)
     {
         auto d = *tmp - eit->first();
-        auto p = eit->second() + d;
+        Vertex p = eit->second() + d;
 
-        ShapeLike::addVertex(rsh, p);
+        sl::addVertex(rsh, p);
+
+        // Set the new reference vertex
+        if(cmp(top_nfp, p)) top_nfp = p;
 
         tmp = std::next(tmp);
     }
 
-    // Now we have an nfp somewhere in the dark. We need to get it
-    // to the right position around the stationary shape.
-    // This is done by choosing the leftmost lowest vertex of the
-    // orbiting polygon to be touched with the rightmost upper
-    // vertex of the stationary polygon. In this configuration, the
-    // reference vertex of the orbiting polygon (which can be dragged around
-    // the nfp) will be its rightmost upper vertex that coincides with the
-    // rightmost upper vertex of the nfp. No proof provided other than Jonas
-    // Lindmark's reasoning about the reference vertex of nfp in his thesis
-    // ("No fit polygon problem" - section 2.1.9)
+    return {rsh, top_nfp};
+}
 
-    // TODO: dont do this here. Cache the rmu and lmd in Item and get translate
-    // the nfp after this call
+template<class RawShape>
+static NfpResult<RawShape> nfpSimpleSimple(const RawShape& cstationary,
+                                           const RawShape& cother)
+{
 
-    auto csh = sh;  // Copy sh, we will sort the verices in the copy
-    auto& cmp = _vsort<RawShape>;
-    std::sort(ShapeLike::begin(csh), ShapeLike::end(csh), cmp);
-    std::sort(ShapeLike::begin(other), ShapeLike::end(other), cmp);
+    // Algorithms are from the original algorithm proposed in paper:
+    // https://eprints.soton.ac.uk/36850/1/CORMSIS-05-05.pdf
 
-    // leftmost lower vertex of the stationary polygon
-    auto& touch_sh = *(std::prev(ShapeLike::end(csh)));
-    // rightmost upper vertex of the orbiting polygon
-    auto& touch_other = *(ShapeLike::begin(other));
+    // /////////////////////////////////////////////////////////////////////////
+    // Algorithm 1: Obtaining the minkowski sum
+    // /////////////////////////////////////////////////////////////////////////
 
-    // Calculate the difference and move the orbiter to the touch position.
-    auto dtouch = touch_sh - touch_other;
-    auto top_other = *(std::prev(ShapeLike::end(other))) + dtouch;
+    // 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.
 
-    // Get the righmost upper vertex of the nfp and move it to the RMU of
-    // the orbiter because they should coincide.
-    auto&& top_nfp = rightmostUpVertex(rsh);
-    auto dnfp = top_other - top_nfp;
-    std::for_each(ShapeLike::begin(rsh), ShapeLike::end(rsh),
-                  [&dnfp](Vertex& v) { v+= dnfp; } );
+    using Result = NfpResult<RawShape>;
+    using Vertex = TPoint<RawShape>;
+    using Coord = TCoord<Vertex>;
+    using Edge = _Segment<Vertex>;
+    using sl = ShapeLike;
+    using std::signbit;
+    using std::sort;
+    using std::vector;
+    using std::ref;
+    using std::reference_wrapper;
 
-    return rsh;
+    // 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::getContour(cstationary);
+    std::reverse(stcont.begin(), stcont.end());
+    RawShape stationary;
+    sl::getContour(stationary) = stcont;
+
+    // Reverse the orbiter contour to counter clockwise
+    auto orbcont = sl::getContour(cother);
+
+    std::reverse(orbcont.begin(), orbcont.end());
+
+    // Copy the orbiter (contour only), we will have to work on it
+    RawShape orbiter;
+    sl::getContour(orbiter) = orbcont;
+
+    // Step 1: Make the orbiter reverse oriented
+    for(auto &v : sl::getContour(orbiter)) v = -v;
+
+    // An egde 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) {}
+    };
+
+    // 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& poly, int dir) {
+        L.reserve(sl::contourVertexCount(poly));
+
+        auto it = sl::cbegin(poly);
+        auto nextit = std::next(it);
+
+        double turn_angle = 0;
+        bool is_turn_point = false;
+
+        while(nextit != sl::cend(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 TwoPi = 2.0*Pi;
+            if(phi > Pi) phi -= TwoPi;
+            if(phi_prev > Pi) phi_prev -= TwoPi;
+            auto turn_angle = phi-phi_prev;
+            if(turn_angle > Pi) turn_angle -= TwoPi;
+            return phi-phi_prev;
+        };
+
+        if(dir > 0) {
+            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);
+                enext->is_turning_point =
+                        signbit(enext->turn_angle) != signbit(eit->turn_angle);
+                ++eit; ++enext;
+            }
+
+            L.front().is_turning_point = signbit(L.front().turn_angle) !=
+                                         signbit(L.back().turn_angle);
+        } else {
+            std::cout << L.size() << std::endl;
+
+            auto eit = L.rbegin();
+            auto enext = std::next(eit);
+
+            eit->turn_angle = getTurnAngle(L.back().e, L.front().e);
+
+            while(enext != L.rend()) {
+                enext->turn_angle = getTurnAngle(enext->e, eit->e);
+                enext->is_turning_point =
+                        signbit(enext->turn_angle) != signbit(eit->turn_angle);
+                std::cout << enext->is_turning_point << " " << enext->turn_angle << std::endl;
+
+                ++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);
+
+    // 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) );
+    });
+
+    struct EdgeGroup { typename EdgeRefList::const_iterator first, last; };
+
+    auto mink = [sortfn] // the Mink(Q, R, direction) sub-procedure
+            (const EdgeGroup& Q, const EdgeGroup& 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.last - Q.first) + (R.last - R.first));
+
+        merged.insert(merged.end(), Q.first, Q.last);
+        merged.insert(merged.end(), R.first, R.last);
+        sort(merged.begin(), merged.end(), sortfn);
+
+        // Step 2 "set i = 1, k = 1, direction = 1, s1 = q1"
+        // we dont 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.first->container.get();
+        const auto& Qcont = Q.first->container.get();
+
+        // Set the intial direction
+        Coord dir = positive? 1 : -1;
+
+        // roughly i = 1 (so q = Q.first) and s1 = q1 so S[0] = q;
+        auto q = Q.first;
+        S.push_back(*q++);
+
+        // Roughly step 3
+        while(q != Q.last) {
+            auto it = merged.begin();
+            while(it != merged.end() && !(it->eq(*(Q.first))) ) {
+                if(it->isFrom(Rcont)) {
+                    auto s = *it;
+                    s.dir = dir;
+                    S.push_back(s);
+                }
+                if(it->eq(*q)) {
+                    S.push_back(*q);
+                    if(it->isTurningPoint()) dir = -dir;
+                    if(q != Q.first) it += dir;
+                }
+                else it += dir;
+            }
+            ++q; // "Set i = i + 1"
+        }
+
+        // Step 4:
+
+        // "Let starting edge r1 be in position si in sequence"
+        // whaaat? I guess this means the following:
+        S[0] = *R.first;
+        auto it = S.begin();
+
+        // "Set j = 1, next = 2, direction = 1, seq1 = si"
+        // we dont use j, seq is expanded dynamically.
+        dir = 1; auto next = std::next(R.first);
+
+        // Step 5:
+        // "If all si edges have been allocated to seqj" should mean that
+        // we loop until seq has equal size with S
+        while(seq.size() < S.size()) {
+            ++it; if(it == S.end()) it = S.begin();
+
+            if(it->isFrom(Qcont)) {
+                seq.push_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; }
+            }
+
+            if(it->eq(*next) && dir == next->dir) { // "If si = direction.rnext"
+                // "j = j + 1, seqj = si, next = next + direction"
+                seq.push_back(*it); next += dir;
+            }
+        }
+
+        return seq;
+    };
+
+    EdgeGroup R{ Bref.begin(), Bref.begin() }, Q{ Aref.begin(), Aref.end() };
+    auto it = Bref.begin();
+    bool orientation = true;
+    EdgeRefList seqlist;
+    seqlist.reserve(3*(Aref.size() + Bref.size()));
+
+    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
+    // /////////////////////////////////////////////////////////////////////////
+
+
+
+    return Result(stationary, Vertex());
 }
 
 // Specializable NFP implementation class. Specialize it if you have a faster
 // or better NFP implementation
 template<class RawShape, NfpLevel nfptype>
 struct NfpImpl {
-    RawShape operator()(const RawShape& sh, const RawShape& other) {
+    NfpResult<RawShape> operator()(const RawShape& sh, const RawShape& other)
+    {
         static_assert(nfptype == NfpLevel::CONVEX_ONLY,
                       "Nfp::noFitPolygon() unimplemented!");
 
diff --git a/xs/src/libnest2d/libnest2d/libnest2d.hpp b/xs/src/libnest2d/libnest2d/libnest2d.hpp
index 96316c344..37b5fea95 100644
--- a/xs/src/libnest2d/libnest2d/libnest2d.hpp
+++ b/xs/src/libnest2d/libnest2d/libnest2d.hpp
@@ -9,6 +9,7 @@
 #include <functional>
 
 #include "geometry_traits.hpp"
+#include "optimizer.hpp"
 
 namespace libnest2d {
 
@@ -27,6 +28,7 @@ class _Item {
     using Coord = TCoord<TPoint<RawShape>>;
     using Vertex = TPoint<RawShape>;
     using Box = _Box<Vertex>;
+    using sl = ShapeLike;
 
     // The original shape that gets encapsulated.
     RawShape sh_;
@@ -56,6 +58,13 @@ class _Item {
     };
 
     mutable Convexity convexity_ = Convexity::UNCHECKED;
+    mutable TVertexConstIterator<RawShape> rmt_;    // rightmost top vertex
+    mutable TVertexConstIterator<RawShape> lmb_;    // leftmost bottom vertex
+    mutable bool rmt_valid_ = false, lmb_valid_ = false;
+    mutable struct BBCache {
+        Box bb; bool valid; Vertex tr;
+        BBCache(): valid(false), tr(0, 0) {}
+    } bb_cache_;
 
 public:
 
@@ -104,15 +113,15 @@ public:
      * @param il The initializer list of vertices.
      */
     inline _Item(const std::initializer_list< Vertex >& il):
-        sh_(ShapeLike::create<RawShape>(il)) {}
+        sh_(sl::create<RawShape>(il)) {}
 
     inline _Item(const TContour<RawShape>& contour,
                  const THolesContainer<RawShape>& holes = {}):
-        sh_(ShapeLike::create<RawShape>(contour, holes)) {}
+        sh_(sl::create<RawShape>(contour, holes)) {}
 
     inline _Item(TContour<RawShape>&& contour,
                  THolesContainer<RawShape>&& holes):
-        sh_(ShapeLike::create<RawShape>(std::move(contour),
+        sh_(sl::create<RawShape>(std::move(contour),
                                         std::move(holes))) {}
 
     /**
@@ -122,31 +131,31 @@ public:
      */
     inline std::string toString() const
     {
-        return ShapeLike::toString(sh_);
+        return sl::toString(sh_);
     }
 
     /// Iterator tho the first contour vertex in the polygon.
     inline Iterator begin() const
     {
-        return ShapeLike::cbegin(sh_);
+        return sl::cbegin(sh_);
     }
 
     /// Alias to begin()
     inline Iterator cbegin() const
     {
-        return ShapeLike::cbegin(sh_);
+        return sl::cbegin(sh_);
     }
 
     /// Iterator to the last contour vertex.
     inline Iterator end() const
     {
-        return ShapeLike::cend(sh_);
+        return sl::cend(sh_);
     }
 
     /// Alias to end()
     inline Iterator cend() const
     {
-        return ShapeLike::cend(sh_);
+        return sl::cend(sh_);
     }
 
     /**
@@ -161,7 +170,7 @@ public:
      */
     inline Vertex vertex(unsigned long idx) const
     {
-        return ShapeLike::vertex(sh_, idx);
+        return sl::vertex(sh_, idx);
     }
 
     /**
@@ -176,7 +185,7 @@ public:
     inline void setVertex(unsigned long idx, const Vertex& v )
     {
         invalidateCache();
-        ShapeLike::vertex(sh_, idx) = v;
+        sl::vertex(sh_, idx) = v;
     }
 
     /**
@@ -191,7 +200,7 @@ public:
         double ret ;
         if(area_cache_valid_) ret = area_cache_;
         else {
-            ret = ShapeLike::area(offsettedShape());
+            ret = sl::area(offsettedShape());
             area_cache_ = ret;
             area_cache_valid_ = true;
         }
@@ -203,7 +212,7 @@ public:
 
         switch(convexity_) {
         case Convexity::UNCHECKED:
-            ret = ShapeLike::isConvex<RawShape>(ShapeLike::getContour(transformedShape()));
+            ret = sl::isConvex<RawShape>(sl::getContour(transformedShape()));
             convexity_ = ret? Convexity::TRUE : Convexity::FALSE;
             break;
         case Convexity::TRUE: ret = true; break;
@@ -213,7 +222,7 @@ public:
         return ret;
     }
 
-    inline bool isHoleConvex(unsigned holeidx) const {
+    inline bool isHoleConvex(unsigned /*holeidx*/) const {
         return false;
     }
 
@@ -223,11 +232,11 @@ public:
 
     /// The number of the outer ring vertices.
     inline size_t vertexCount() const {
-        return ShapeLike::contourVertexCount(sh_);
+        return sl::contourVertexCount(sh_);
     }
 
     inline size_t holeCount() const {
-        return ShapeLike::holeCount(sh_);
+        return sl::holeCount(sh_);
     }
 
     /**
@@ -235,36 +244,33 @@ public:
      * @param p
      * @return
      */
-    inline bool isPointInside(const Vertex& p)
+    inline bool isPointInside(const Vertex& p) const
     {
-        return ShapeLike::isInside(p, sh_);
+        return sl::isInside(p, transformedShape());
     }
 
     inline bool isInside(const _Item& sh) const
     {
-        return ShapeLike::isInside(transformedShape(), sh.transformedShape());
+        return sl::isInside(transformedShape(), sh.transformedShape());
     }
 
-    inline bool isInside(const _Box<TPoint<RawShape>>& box);
+    inline bool isInside(const _Box<TPoint<RawShape>>& box) const;
 
     inline void translate(const Vertex& d) BP2D_NOEXCEPT
     {
-        translation_ += d; has_translation_ = true;
-        tr_cache_valid_ = false;
+        translation(translation() + d);
     }
 
     inline void rotate(const Radians& rads) BP2D_NOEXCEPT
     {
-        rotation_ += rads;
-        has_rotation_ = true;
-        tr_cache_valid_ = false;
+        rotation(rotation() + rads);
     }
 
     inline void addOffset(Coord distance) BP2D_NOEXCEPT
     {
         offset_distance_ = distance;
         has_offset_ = true;
-        offset_cache_valid_ = false;
+        invalidateCache();
     }
 
     inline void removeOffset() BP2D_NOEXCEPT {
@@ -286,6 +292,8 @@ public:
     {
         if(rotation_ != rot) {
             rotation_ = rot; has_rotation_ = true; tr_cache_valid_ = false;
+            rmt_valid_ = false; lmb_valid_ = false;
+            bb_cache_.valid = false;
         }
     }
 
@@ -293,6 +301,7 @@ public:
     {
         if(translation_ != tr) {
             translation_ = tr; has_translation_ = true; tr_cache_valid_ = false;
+            bb_cache_.valid = false;
         }
     }
 
@@ -301,9 +310,10 @@ public:
         if(tr_cache_valid_) return tr_cache_;
 
         RawShape cpy = offsettedShape();
-        if(has_rotation_) ShapeLike::rotate(cpy, rotation_);
-        if(has_translation_) ShapeLike::translate(cpy, translation_);
+        if(has_rotation_) sl::rotate(cpy, rotation_);
+        if(has_translation_) sl::translate(cpy, translation_);
         tr_cache_ = cpy; tr_cache_valid_ = true;
+        rmt_valid_ = false; lmb_valid_ = false;
 
         return tr_cache_;
     }
@@ -321,23 +331,53 @@ public:
     inline void resetTransformation() BP2D_NOEXCEPT
     {
         has_translation_ = false; has_rotation_ = false; has_offset_ = false;
+        invalidateCache();
     }
 
     inline Box boundingBox() const {
-        return ShapeLike::boundingBox(transformedShape());
+        if(!bb_cache_.valid) {
+            bb_cache_.bb = sl::boundingBox(transformedShape());
+            bb_cache_.tr = {0, 0};
+            bb_cache_.valid = true;
+        }
+
+        auto &bb = bb_cache_.bb; auto &tr = bb_cache_.tr;
+        return {bb.minCorner() + tr, bb.maxCorner() + tr};
+    }
+
+    inline Vertex referenceVertex() const {
+        return rightmostTopVertex();
+    }
+
+    inline Vertex rightmostTopVertex() const {
+        if(!rmt_valid_ || !tr_cache_valid_) {  // find max x and max y vertex
+            auto& tsh = transformedShape();
+            rmt_ = std::max_element(sl::cbegin(tsh), sl::cend(tsh), vsort);
+            rmt_valid_ = true;
+        }
+        return *rmt_;
+    }
+
+    inline Vertex leftmostBottomVertex() const {
+        if(!lmb_valid_ || !tr_cache_valid_) {  // find min x and min y vertex
+            auto& tsh = transformedShape();
+            lmb_ = std::min_element(sl::cbegin(tsh), sl::cend(tsh), vsort);
+            lmb_valid_ = true;
+        }
+        return *lmb_;
     }
 
     //Static methods:
 
     inline static bool intersects(const _Item& sh1, const _Item& sh2)
     {
-        return ShapeLike::intersects(sh1.transformedShape(),
+        return sl::intersects(sh1.transformedShape(),
                                      sh2.transformedShape());
     }
 
     inline static bool touches(const _Item& sh1, const _Item& sh2)
     {
-        return ShapeLike::touches(sh1.transformedShape(),
+        return sl::touches(sh1.transformedShape(),
                                   sh2.transformedShape());
     }
 
@@ -346,12 +386,11 @@ private:
     inline const RawShape& offsettedShape() const {
         if(has_offset_ ) {
             if(offset_cache_valid_) return offset_cache_;
-            else {
-                offset_cache_ = sh_;
-                ShapeLike::offset(offset_cache_, offset_distance_);
-                offset_cache_valid_ = true;
-                return offset_cache_;
-            }
+
+            offset_cache_ = sh_;
+            sl::offset(offset_cache_, offset_distance_);
+            offset_cache_valid_ = true;
+            return offset_cache_;
         }
         return sh_;
     }
@@ -359,10 +398,23 @@ private:
     inline void invalidateCache() const BP2D_NOEXCEPT
     {
         tr_cache_valid_ = false;
+        lmb_valid_ = false; rmt_valid_ = false;
         area_cache_valid_ = false;
         offset_cache_valid_ = false;
+        bb_cache_.valid = false;
         convexity_ = Convexity::UNCHECKED;
     }
+
+    static inline bool vsort(const Vertex& v1, const Vertex& v2)
+    {
+        Coord &&x1 = getX(v1), &&x2 = getX(v2);
+        Coord &&y1 = getY(v1), &&y2 = getY(v2);
+        auto diff = y1 - y2;
+        if(std::abs(diff) <= std::numeric_limits<Coord>::epsilon())
+            return x1 < x2;
+
+        return diff < 0;
+    }
 };
 
 /**
@@ -370,7 +422,6 @@ private:
  */
 template<class RawShape>
 class _Rectangle: public _Item<RawShape> {
-    RawShape sh_;
     using _Item<RawShape>::vertex;
     using TO = Orientation;
 public:
@@ -415,7 +466,7 @@ public:
 };
 
 template<class RawShape>
-inline bool _Item<RawShape>::isInside(const _Box<TPoint<RawShape>>& box) {
+inline bool _Item<RawShape>::isInside(const _Box<TPoint<RawShape>>& box) const {
     _Rectangle<RawShape> rect(box.width(), box.height());
     return _Item<RawShape>::isInside(rect);
 }
@@ -874,9 +925,8 @@ private:
 
     Radians findBestRotation(Item& item) {
         opt::StopCriteria stopcr;
-        stopcr.stoplimit = 0.01;
+        stopcr.absolute_score_difference = 0.01;
         stopcr.max_iterations = 10000;
-        stopcr.type = opt::StopLimitType::RELATIVE;
         opt::TOptimizer<opt::Method::G_GENETIC> solver(stopcr);
 
         auto orig_rot = item.rotation();
@@ -910,7 +960,6 @@ private:
         if(min_obj_distance_ > 0) std::for_each(from, to, [](Item& item) {
             item.removeOffset();
         });
-
     }
 };
 
diff --git a/xs/src/libnest2d/libnest2d/metaloop.hpp b/xs/src/libnest2d/libnest2d/metaloop.hpp
new file mode 100644
index 000000000..18755525c
--- /dev/null
+++ b/xs/src/libnest2d/libnest2d/metaloop.hpp
@@ -0,0 +1,227 @@
+#ifndef METALOOP_HPP
+#define METALOOP_HPP
+
+#include "common.hpp"
+#include <tuple>
+#include <functional>
+
+namespace libnest2d {
+
+/* ************************************************************************** */
+/* C++14 std::index_sequence implementation:                                  */
+/* ************************************************************************** */
+
+/**
+ * \brief C++11 conformant implementation of the index_sequence type from C++14
+ */
+template<size_t...Ints> struct index_sequence {
+    using value_type = size_t;
+    BP2D_CONSTEXPR value_type size() const { return sizeof...(Ints); }
+};
+
+// A Help structure to generate the integer list
+template<size_t...Nseq> struct genSeq;
+
+// Recursive template to generate the list
+template<size_t I, size_t...Nseq> struct genSeq<I, Nseq...> {
+    // Type will contain a genSeq with Nseq appended by one element
+    using Type = typename genSeq< I - 1, I - 1, Nseq...>::Type;
+};
+
+// Terminating recursion
+template <size_t ... Nseq> struct genSeq<0, Nseq...> {
+    // If I is zero, Type will contain index_sequence with the fuly generated
+    // integer list.
+    using Type = index_sequence<Nseq...>;
+};
+
+/// Helper alias to make an index sequence from 0 to N
+template<size_t N> using make_index_sequence = typename genSeq<N>::Type;
+
+/// Helper alias to make an index sequence for a parameter pack
+template<class...Args>
+using index_sequence_for = make_index_sequence<sizeof...(Args)>;
+
+
+/* ************************************************************************** */
+
+namespace opt {
+
+using std::forward;
+using std::tuple;
+using std::get;
+using std::tuple_element;
+
+/**
+ * @brief Helper class to be able to loop over a parameter pack's elements.
+ */
+class metaloop {
+
+// The implementation is based on partial struct template specializations.
+// Basically we need a template type that is callable and takes an integer
+// non-type template parameter which can be used to implement recursive calls.
+//
+// C++11 will not allow the usage of a plain template function that is why we
+// use struct with overloaded call operator. At the same time C++11 prohibits
+// partial template specialization with a non type parameter such as int. We
+// need to wrap that in a type (see metaloop::Int).
+
+/*
+ * A helper alias to create integer values wrapped as a type. It is nessecary
+ * because a non type template parameter (such as int) would be prohibited in
+ * a partial specialization. Also for the same reason we have to use a class
+ * _Metaloop instead of a simple function as a functor. A function cannot be
+ * partially specialized in a way that is neccesary for this trick.
+ */
+template<int N> using Int = std::integral_constant<int, N>;
+
+/*
+ * Helper class to implement in-place functors.
+ *
+ * We want to be able to use inline functors like a lambda to keep the code
+ * as clear as possible.
+ */
+template<int N, class Fn> class MapFn {
+    Fn&& fn_;
+public:
+
+    // It takes the real functor that can be specified in-place but only
+    // with C++14 because the second parameter's type will depend on the
+    // type of the parameter pack element that is processed. In C++14 we can
+    // specify this second parameter type as auto in the lamda parameter list.
+    inline MapFn(Fn&& fn): fn_(forward<Fn>(fn)) {}
+
+    template<class T> void operator ()(T&& pack_element) {
+        // We provide the index as the first parameter and the pack (or tuple)
+        // element as the second parameter to the functor.
+        fn_(N, forward<T>(pack_element));
+    }
+};
+
+/*
+ * Implementation of the template loop trick.
+ * We create a mechanism for looping over a parameter pack in compile time.
+ * \tparam Idx is the loop index which will be decremented at each recursion.
+ * \tparam Args The parameter pack that will be processed.
+ *
+ */
+template <typename Idx, class...Args>
+class _MetaLoop {};
+
+// Implementation for the first element of Args...
+template <class...Args>
+class _MetaLoop<Int<0>, Args...> {
+public:
+
+    const static BP2D_CONSTEXPR int N = 0;
+    const static BP2D_CONSTEXPR int ARGNUM = sizeof...(Args)-1;
+
+    template<class Tup, class Fn>
+    void run( Tup&& valtup, Fn&& fn) {
+        MapFn<ARGNUM-N, Fn> {forward<Fn>(fn)} (get<ARGNUM-N>(valtup));
+    }
+};
+
+// Implementation for the N-th element of Args...
+template <int N, class...Args>
+class _MetaLoop<Int<N>, Args...> {
+public:
+
+    const static BP2D_CONSTEXPR int ARGNUM = sizeof...(Args)-1;
+
+    template<class Tup, class Fn>
+    void run(Tup&& valtup, Fn&& fn) {
+        MapFn<ARGNUM-N, Fn> {forward<Fn>(fn)} (std::get<ARGNUM-N>(valtup));
+
+        // Recursive call to process the next element of Args
+        _MetaLoop<Int<N-1>, Args...> ().run(forward<Tup>(valtup),
+                                            forward<Fn>(fn));
+    }
+};
+
+/*
+ * Instantiation: We must instantiate the template with the last index because
+ * the generalized version calls the decremented instantiations recursively.
+ * Once the instantiation with the first index is called, the terminating
+ * version of run is called which does not call itself anymore.
+ *
+ * If you are utterly annoyed, at least you have learned a super crazy
+ * functional metaprogramming pattern.
+ */
+template<class...Args>
+using MetaLoop = _MetaLoop<Int<sizeof...(Args)-1>, Args...>;
+
+public:
+
+/**
+ * \brief The final usable function template.
+ *
+ * This is similar to what varags was on C but in compile time C++11.
+ * You can call:
+ * apply(<the mapping function>, <arbitrary number of arguments of any type>);
+ * For example:
+ *
+ *      struct mapfunc {
+ *          template<class T> void operator()(int N, T&& element) {
+ *              std::cout << "The value of the parameter "<< N <<": "
+ *                        << element << std::endl;
+ *          }
+ *      };
+ *
+ *      apply(mapfunc(), 'a', 10, 151.545);
+ *
+ * C++14:
+ *      apply([](int N, auto&& element){
+ *          std::cout << "The value of the parameter "<< N <<": "
+ *                        << element << std::endl;
+ *      }, 'a', 10, 151.545);
+ *
+ * This yields the output:
+ * The value of the parameter 0: a
+ * The value of the parameter 1: 10
+ * The value of the parameter 2: 151.545
+ *
+ * As an addition, the function can be called with a tuple as the second
+ * parameter holding the arguments instead of a parameter pack.
+ *
+ */
+template<class...Args, class Fn>
+inline static void apply(Fn&& fn, Args&&...args) {
+    MetaLoop<Args...>().run(tuple<Args&&...>(forward<Args>(args)...),
+                            forward<Fn>(fn));
+}
+
+/// The version of apply with a tuple rvalue reference.
+template<class...Args, class Fn>
+inline static void apply(Fn&& fn, tuple<Args...>&& tup) {
+    MetaLoop<Args...>().run(std::move(tup), forward<Fn>(fn));
+}
+
+/// The version of apply with a tuple lvalue reference.
+template<class...Args, class Fn>
+inline static void apply(Fn&& fn, tuple<Args...>& tup) {
+    MetaLoop<Args...>().run(tup, forward<Fn>(fn));
+}
+
+/// The version of apply with a tuple const reference.
+template<class...Args, class Fn>
+inline static void apply(Fn&& fn, const tuple<Args...>& tup) {
+    MetaLoop<Args...>().run(tup, forward<Fn>(fn));
+}
+
+/**
+ * Call a function with its arguments encapsualted in a tuple.
+ */
+template<class Fn, class Tup, std::size_t...Is>
+inline static auto
+callFunWithTuple(Fn&& fn, Tup&& tup, index_sequence<Is...>) ->
+    decltype(fn(std::get<Is>(tup)...))
+{
+    return fn(std::get<Is>(tup)...);
+}
+
+};
+}
+}
+
+#endif // METALOOP_HPP
diff --git a/xs/src/libnest2d/libnest2d/optimizer.hpp b/xs/src/libnest2d/libnest2d/optimizer.hpp
index 52e67f7d5..c8ed2e378 100644
--- a/xs/src/libnest2d/libnest2d/optimizer.hpp
+++ b/xs/src/libnest2d/libnest2d/optimizer.hpp
@@ -10,8 +10,7 @@ namespace libnest2d { namespace opt {
 
 using std::forward;
 using std::tuple;
-using std::get;
-using std::tuple_element;
+using std::make_tuple;
 
 /// A Type trait for upper and lower limit of a numeric type.
 template<class T, class B = void >
@@ -51,176 +50,7 @@ inline Bound<T> bound(const T& min, const T& max) { return Bound<T>(min, max); }
 template<class...Args> using Input = tuple<Args...>;
 
 template<class...Args>
-inline tuple<Args...> initvals(Args...args) { return std::make_tuple(args...); }
-
-/**
- * @brief Helper class to be able to loop over a parameter pack's elements.
- */
-class metaloop {
-// The implementation is based on partial struct template specializations.
-// Basically we need a template type that is callable and takes an integer
-// non-type template parameter which can be used to implement recursive calls.
-//
-// C++11 will not allow the usage of a plain template function that is why we
-// use struct with overloaded call operator. At the same time C++11 prohibits
-// partial template specialization with a non type parameter such as int. We
-// need to wrap that in a type (see metaloop::Int).
-
-/*
- * A helper alias to create integer values wrapped as a type. It is nessecary
- * because a non type template parameter (such as int) would be prohibited in
- * a partial specialization. Also for the same reason we have to use a class
- * _Metaloop instead of a simple function as a functor. A function cannot be
- * partially specialized in a way that is neccesary for this trick.
- */
-template<int N> using Int = std::integral_constant<int, N>;
-
-/*
- * Helper class to implement in-place functors.
- *
- * We want to be able to use inline functors like a lambda to keep the code
- * as clear as possible.
- */
-template<int N, class Fn> class MapFn {
-    Fn&& fn_;
-public:
-
-    // It takes the real functor that can be specified in-place but only
-    // with C++14 because the second parameter's type will depend on the
-    // type of the parameter pack element that is processed. In C++14 we can
-    // specify this second parameter type as auto in the lamda parameter list.
-    inline MapFn(Fn&& fn): fn_(forward<Fn>(fn)) {}
-
-    template<class T> void operator ()(T&& pack_element) {
-        // We provide the index as the first parameter and the pack (or tuple)
-        // element as the second parameter to the functor.
-        fn_(N, forward<T>(pack_element));
-    }
-};
-
-/*
- * Implementation of the template loop trick.
- * We create a mechanism for looping over a parameter pack in compile time.
- * \tparam Idx is the loop index which will be decremented at each recursion.
- * \tparam Args The parameter pack that will be processed.
- *
- */
-template <typename Idx, class...Args>
-class _MetaLoop {};
-
-// Implementation for the first element of Args...
-template <class...Args>
-class _MetaLoop<Int<0>, Args...> {
-public:
-
-    const static BP2D_CONSTEXPR int N = 0;
-    const static BP2D_CONSTEXPR int ARGNUM = sizeof...(Args)-1;
-
-    template<class Tup, class Fn>
-    void run( Tup&& valtup, Fn&& fn) {
-        MapFn<ARGNUM-N, Fn> {forward<Fn>(fn)} (get<ARGNUM-N>(valtup));
-    }
-};
-
-// Implementation for the N-th element of Args...
-template <int N, class...Args>
-class _MetaLoop<Int<N>, Args...> {
-public:
-
-    const static BP2D_CONSTEXPR int ARGNUM = sizeof...(Args)-1;
-
-    template<class Tup, class Fn>
-    void run(Tup&& valtup, Fn&& fn) {
-        MapFn<ARGNUM-N, Fn> {forward<Fn>(fn)} (std::get<ARGNUM-N>(valtup));
-
-        // Recursive call to process the next element of Args
-        _MetaLoop<Int<N-1>, Args...> ().run(forward<Tup>(valtup),
-                                            forward<Fn>(fn));
-    }
-};
-
-/*
- * Instantiation: We must instantiate the template with the last index because
- * the generalized version calls the decremented instantiations recursively.
- * Once the instantiation with the first index is called, the terminating
- * version of run is called which does not call itself anymore.
- *
- * If you are utterly annoyed, at least you have learned a super crazy
- * functional metaprogramming pattern.
- */
-template<class...Args>
-using MetaLoop = _MetaLoop<Int<sizeof...(Args)-1>, Args...>;
-
-public:
-
-/**
- * \brief The final usable function template.
- *
- * This is similar to what varags was on C but in compile time C++11.
- * You can call:
- * apply(<the mapping function>, <arbitrary number of arguments of any type>);
- * For example:
- *
- *      struct mapfunc {
- *          template<class T> void operator()(int N, T&& element) {
- *              std::cout << "The value of the parameter "<< N <<": "
- *                        << element << std::endl;
- *          }
- *      };
- *
- *      apply(mapfunc(), 'a', 10, 151.545);
- *
- * C++14:
- *      apply([](int N, auto&& element){
- *          std::cout << "The value of the parameter "<< N <<": "
- *                        << element << std::endl;
- *      }, 'a', 10, 151.545);
- *
- * This yields the output:
- * The value of the parameter 0: a
- * The value of the parameter 1: 10
- * The value of the parameter 2: 151.545
- *
- * As an addition, the function can be called with a tuple as the second
- * parameter holding the arguments instead of a parameter pack.
- *
- */
-template<class...Args, class Fn>
-inline static void apply(Fn&& fn, Args&&...args) {
-    MetaLoop<Args...>().run(tuple<Args&&...>(forward<Args>(args)...),
-                            forward<Fn>(fn));
-}
-
-/// The version of apply with a tuple rvalue reference.
-template<class...Args, class Fn>
-inline static void apply(Fn&& fn, tuple<Args...>&& tup) {
-    MetaLoop<Args...>().run(std::move(tup), forward<Fn>(fn));
-}
-
-/// The version of apply with a tuple lvalue reference.
-template<class...Args, class Fn>
-inline static void apply(Fn&& fn, tuple<Args...>& tup) {
-    MetaLoop<Args...>().run(tup, forward<Fn>(fn));
-}
-
-/// The version of apply with a tuple const reference.
-template<class...Args, class Fn>
-inline static void apply(Fn&& fn, const tuple<Args...>& tup) {
-    MetaLoop<Args...>().run(tup, forward<Fn>(fn));
-}
-
-/**
- * Call a function with its arguments encapsualted in a tuple.
- */
-template<class Fn, class Tup, std::size_t...Is>
-inline static auto
-callFunWithTuple(Fn&& fn, Tup&& tup, index_sequence<Is...>) ->
-    decltype(fn(std::get<Is>(tup)...))
-{
-    return fn(std::get<Is>(tup)...);
-}
-
-};
+inline tuple<Args...> initvals(Args...args) { return make_tuple(args...); }
 
 /**
  * @brief Specific optimization methods for which a default optimizer
@@ -257,29 +87,20 @@ enum ResultCodes {
 template<class...Args>
 struct Result {
     ResultCodes resultcode;
-    std::tuple<Args...> optimum;
+    tuple<Args...> optimum;
     double score;
 };
 
-/**
- * @brief The stop limit can be specified as the absolute error or as the
- * relative error, just like in nlopt.
- */
-enum class StopLimitType {
-    ABSOLUTE,
-    RELATIVE
-};
-
 /**
  * @brief A type for specifying the stop criteria.
  */
 struct StopCriteria {
 
-    /// Relative or absolute termination error
-    StopLimitType type = StopLimitType::RELATIVE;
+    /// If the absolute value difference between two scores.
+    double absolute_score_difference = std::nan("");
 
-    /// The error value that is interpredted depending on the type property.
-    double stoplimit = 0.0001;
+    /// If the relative value difference between two scores.
+    double relative_score_difference = std::nan("");
 
     unsigned max_iterations = 0;
 };
@@ -310,11 +131,11 @@ public:
      * \return Returns a Result<Args...> structure.
      * An example call would be:
      *     auto result = opt.optimize_min(
-     *           [](std::tuple<double> x) // object function
+     *           [](tuple<double> x) // object function
      *           {
      *               return std::pow(std::get<0>(x), 2);
      *           },
-     *           std::make_tuple(-0.5),  // initial value
+     *           make_tuple(-0.5),  // initial value
      *           {-1.0, 1.0}           // search space bounds
      *       );
      */
@@ -390,10 +211,14 @@ public:
         static_assert(always_false<T>::value, "Optimizer unimplemented!");
     }
 
+    DummyOptimizer(const StopCriteria&) {
+        static_assert(always_false<T>::value, "Optimizer unimplemented!");
+    }
+
     template<class Func, class...Args>
-    Result<Args...> optimize(Func&& func,
-                             std::tuple<Args...> initvals,
-                             Bound<Args>... args)
+    Result<Args...> optimize(Func&& /*func*/,
+                             tuple<Args...> /*initvals*/,
+                             Bound<Args>...  /*args*/)
     {
         return Result<Args...>();
     }
diff --git a/xs/src/libnest2d/libnest2d/optimizers/nlopt_boilerplate.hpp b/xs/src/libnest2d/libnest2d/optimizers/nlopt_boilerplate.hpp
index 798bf9622..737ca6e3c 100644
--- a/xs/src/libnest2d/libnest2d/optimizers/nlopt_boilerplate.hpp
+++ b/xs/src/libnest2d/libnest2d/optimizers/nlopt_boilerplate.hpp
@@ -1,15 +1,25 @@
 #ifndef NLOPT_BOILERPLATE_HPP
 #define NLOPT_BOILERPLATE_HPP
 
+#ifdef _MSC_VER
+#pragma warning(push)
+#pragma warning(disable: 4244)
+#pragma warning(disable: 4267)
+#endif
 #include <nlopt.hpp>
+#ifdef _MSC_VER
+#pragma warning(pop)
+#endif
+
 #include <libnest2d/optimizer.hpp>
 #include <cassert>
+#include "libnest2d/metaloop.hpp"
 
 #include <utility>
 
 namespace libnest2d { namespace opt {
 
-nlopt::algorithm method2nloptAlg(Method m) {
+inline nlopt::algorithm method2nloptAlg(Method m) {
 
     switch(m) {
     case Method::L_SIMPLEX: return nlopt::LN_NELDERMEAD;
@@ -87,7 +97,7 @@ protected:
 
     template<class Fn, class...Args>
     static double optfunc(const std::vector<double>& params,
-                          std::vector<double>& grad,
+                          std::vector<double>& /*grad*/,
                           void *data)
     {
         auto fnptr = static_cast<remove_ref_t<Fn>*>(data);
@@ -132,12 +142,10 @@ protected:
         default: ;
         }
 
-        switch(this->stopcr_.type) {
-        case StopLimitType::ABSOLUTE:
-            opt_.set_ftol_abs(stopcr_.stoplimit); break;
-        case StopLimitType::RELATIVE:
-            opt_.set_ftol_rel(stopcr_.stoplimit); break;
-        }
+        auto abs_diff = stopcr_.absolute_score_difference;
+        auto rel_diff = stopcr_.relative_score_difference;
+        if(!std::isnan(abs_diff)) opt_.set_ftol_abs(abs_diff);
+        if(!std::isnan(rel_diff)) opt_.set_ftol_rel(rel_diff);
 
         if(this->stopcr_.max_iterations > 0)
             opt_.set_maxeval(this->stopcr_.max_iterations );
diff --git a/xs/src/libnest2d/libnest2d/placers/nfpplacer.hpp b/xs/src/libnest2d/libnest2d/placers/nfpplacer.hpp
index 89848eb53..61d923b87 100644
--- a/xs/src/libnest2d/libnest2d/placers/nfpplacer.hpp
+++ b/xs/src/libnest2d/libnest2d/placers/nfpplacer.hpp
@@ -6,6 +6,10 @@
 #endif
 #include "placer_boilerplate.hpp"
 #include "../geometry_traits_nfp.hpp"
+#include "libnest2d/optimizer.hpp"
+#include <cassert>
+
+#include "tools/svgtools.hpp"
 
 namespace libnest2d { namespace strategies {
 
@@ -20,15 +24,62 @@ struct NfpPConfig {
         TOP_RIGHT,
     };
 
-    /// Which angles to try out for better results
+    /// Which angles to try out for better results.
     std::vector<Radians> rotations;
 
-    /// Where to align the resulting packed pile
+    /// Where to align the resulting packed pile.
     Alignment alignment;
 
+    /// Where to start putting objects in the bin.
     Alignment starting_point;
 
-    std::function<double(const Nfp::Shapes<RawShape>&, double, double, double)>
+    /**
+     * @brief A function object representing the fitting function in the
+     * placement optimization process. (Optional)
+     *
+     * This is the most versatile tool to configure the placer. The fitting
+     * function is evaluated many times when a new item is being placed into the
+     * bin. The output should be a rated score of the new item's position.
+     *
+     * This is not a mandatory option as there is a default fitting function
+     * that will optimize for the best pack efficiency. With a custom fitting
+     * function you can e.g. influence the shape of the arranged pile.
+     *
+     * \param shapes The first parameter is a container with all the placed
+     * polygons including the current candidate. You can calculate a bounding
+     * box or convex hull on this pile of polygons.
+     *
+     * \param item The second parameter is the candidate item. Note that
+     * calling transformedShape() on this second argument returns an identical
+     * shape as calling shapes.back(). These would not be the same objects only
+     * identical shapes! Using the second parameter is a lot faster due to
+     * caching some properties of the polygon (area, etc...)
+     *
+     * \param occupied_area The third parameter is the sum of areas of the
+     * items in the first parameter so you don't have to iterate through them
+     * if you only need their area.
+     *
+     * \param norm A norming factor for physical dimensions. E.g. if your score
+     * is the distance between the item and the bin center, you should divide
+     * that distance with the norming factor. If the score is an area than
+     * divide it with the square of the norming factor. Imagine it as a unit of
+     * distance.
+     *
+     * \param penality The fifth parameter is the amount of minimum penality if
+     * the arranged pile would't fit into the bin. You can use the wouldFit()
+     * function to check this. Note that the pile can be outside the bin's
+     * boundaries while the placement algorithm is running. Your job is only to
+     * check if the pile could be translated into a position in the bin where
+     * all the items would be inside. For a box shaped bin you can use the
+     * pile's bounding box to check whether it's width and height is small
+     * enough. If the pile would not fit, you have to make sure that the
+     * resulting score will be higher then the penality value. A good solution
+     * would be to set score = 2*penality-score in case the pile wouldn't fit
+     * into the bin.
+     *
+     */
+    std::function<double(Nfp::Shapes<RawShape>&, const _Item<RawShape>&,
+                         double, double, double)>
     object_function;
 
     /**
@@ -38,11 +89,30 @@ struct NfpPConfig {
      */
     float accuracy = 1.0;
 
+    /**
+     * @brief If you want to see items inside other item's holes, you have to
+     * turn this switch on.
+     *
+     * This will only work if a suitable nfp implementation is provided.
+     * The library has no such implementation right now.
+     */
+    bool explore_holes = false;
+
     NfpPConfig(): rotations({0.0, Pi/2.0, Pi, 3*Pi/2}),
         alignment(Alignment::CENTER), starting_point(Alignment::CENTER) {}
 };
 
-// A class for getting a point on the circumference of the polygon (in log time)
+/**
+ * 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.
+ */
 template<class RawShape> class EdgeCache {
     using Vertex = TPoint<RawShape>;
     using Coord = TCoord<Vertex>;
@@ -93,18 +163,22 @@ template<class RawShape> class EdgeCache {
     void fetchCorners() const {
         if(!contour_.corners.empty()) return;
 
-        // TODO Accuracy
-        contour_.corners = contour_.distances;
-        for(auto& d : contour_.corners) d /= contour_.full_distance;
+        contour_.corners.reserve(contour_.distances.size() / 3 + 1);
+        for(size_t i = 0; i < contour_.distances.size() - 1; i += 3) {
+            contour_.corners.emplace_back(
+                    contour_.distances.at(i) / contour_.full_distance);
+        }
     }
 
     void fetchHoleCorners(unsigned hidx) const {
         auto& hc = holes_[hidx];
         if(!hc.corners.empty()) return;
 
-        // TODO Accuracy
-        hc.corners = hc.distances;
-        for(auto& d : hc.corners) d /= hc.full_distance;
+        hc.corners.reserve(hc.distances.size() / 3 + 1);
+        for(size_t i = 0; i < hc.distances.size() - 1; i += 3) {
+            hc.corners.emplace_back(
+                    hc.distances.at(i) / hc.full_distance);
+        }
     }
 
     inline Vertex coords(const ContourCache& cache, double distance) const {
@@ -176,24 +250,64 @@ public:
         return holes_[hidx].full_distance;
     }
 
+    /// Get the normalized distance values for each vertex
     inline const std::vector<double>& corners() const BP2D_NOEXCEPT {
         fetchCorners();
         return contour_.corners;
     }
 
+    /// corners for a specific hole
     inline const std::vector<double>&
     corners(unsigned holeidx) const BP2D_NOEXCEPT {
         fetchHoleCorners(holeidx);
         return holes_[holeidx].corners;
     }
 
-    inline unsigned holeCount() const BP2D_NOEXCEPT { return holes_.size(); }
+    /// The number of holes in the abstracted polygon
+    inline size_t holeCount() const BP2D_NOEXCEPT { return holes_.size(); }
 
 };
 
 template<NfpLevel lvl>
 struct Lvl { static const NfpLevel value = lvl; };
 
+template<class RawShape>
+inline void correctNfpPosition(Nfp::NfpResult<RawShape>& nfp,
+                               const _Item<RawShape>& stationary,
+                               const _Item<RawShape>& orbiter)
+{
+    // The provided nfp is somewhere in the dark. We need to get it
+    // to the right position around the stationary shape.
+    // This is done by choosing the leftmost lowest vertex of the
+    // orbiting polygon to be touched with the rightmost upper
+    // vertex of the stationary polygon. In this configuration, the
+    // reference vertex of the orbiting polygon (which can be dragged around
+    // the nfp) will be its rightmost upper vertex that coincides with the
+    // rightmost upper vertex of the nfp. No proof provided other than Jonas
+    // Lindmark's reasoning about the reference vertex of nfp in his thesis
+    // ("No fit polygon problem" - section 2.1.9)
+
+    auto touch_sh = stationary.rightmostTopVertex();
+    auto touch_other = orbiter.leftmostBottomVertex();
+    auto dtouch = touch_sh - touch_other;
+    auto top_other = orbiter.rightmostTopVertex() + dtouch;
+    auto dnfp = top_other - nfp.second; // nfp.second is the nfp reference point
+    ShapeLike::translate(nfp.first, dnfp);
+}
+
+template<class RawShape>
+inline void correctNfpPosition(Nfp::NfpResult<RawShape>& nfp,
+                               const RawShape& stationary,
+                               const _Item<RawShape>& orbiter)
+{
+    auto touch_sh = Nfp::rightmostUpVertex(stationary);
+    auto touch_other = orbiter.leftmostBottomVertex();
+    auto dtouch = touch_sh - touch_other;
+    auto top_other = orbiter.rightmostTopVertex() + dtouch;
+    auto dnfp = top_other - nfp.second;
+    ShapeLike::translate(nfp.first, dnfp);
+}
+
 template<class RawShape, class Container>
 Nfp::Shapes<RawShape> nfp( const Container& polygons,
                            const _Item<RawShape>& trsh,
@@ -203,18 +317,35 @@ Nfp::Shapes<RawShape> nfp( const Container& polygons,
 
     Nfp::Shapes<RawShape> nfps;
 
+    //int pi = 0;
     for(Item& sh : polygons) {
-        auto subnfp = Nfp::noFitPolygon<NfpLevel::CONVEX_ONLY>(
-                    sh.transformedShape(), trsh.transformedShape());
+        auto subnfp_r = Nfp::noFitPolygon<NfpLevel::CONVEX_ONLY>(
+                            sh.transformedShape(), trsh.transformedShape());
         #ifndef NDEBUG
             auto vv = ShapeLike::isValid(sh.transformedShape());
             assert(vv.first);
 
-            auto vnfp = ShapeLike::isValid(subnfp);
+            auto vnfp = ShapeLike::isValid(subnfp_r.first);
             assert(vnfp.first);
         #endif
 
-        nfps = Nfp::merge(nfps, subnfp);
+        correctNfpPosition(subnfp_r, sh, trsh);
+
+        nfps = Nfp::merge(nfps, subnfp_r.first);
+
+//        double SCALE = 1000000;
+//        using SVGWriter = svg::SVGWriter<RawShape>;
+//        SVGWriter::Config conf;
+//        conf.mm_in_coord_units = SCALE;
+//        SVGWriter svgw(conf);
+//        Box bin(250*SCALE, 210*SCALE);
+//        svgw.setSize(bin);
+//        for(int i = 0; i <= pi; i++) svgw.writeItem(polygons[i]);
+//        svgw.writeItem(trsh);
+////        svgw.writeItem(Item(subnfp_r.first));
+//        for(auto& n : nfps) svgw.writeItem(Item(n));
+//        svgw.save("nfpout");
+//        pi++;
     }
 
     return nfps;
@@ -227,42 +358,65 @@ Nfp::Shapes<RawShape> nfp( const Container& polygons,
 {
     using Item = _Item<RawShape>;
 
-    Nfp::Shapes<RawShape> nfps, stationary;
+    Nfp::Shapes<RawShape> nfps;
+
+    auto& orb = trsh.transformedShape();
+    bool orbconvex = trsh.isContourConvex();
 
     for(Item& sh : polygons) {
-        stationary = Nfp::merge(stationary, sh.transformedShape());
-    }
+        Nfp::NfpResult<RawShape> subnfp;
+        auto& stat = sh.transformedShape();
 
-    std::cout << "pile size: " << stationary.size() << std::endl;
-    for(RawShape& sh : stationary) {
+        if(sh.isContourConvex() && orbconvex)
+            subnfp = Nfp::noFitPolygon<NfpLevel::CONVEX_ONLY>(stat, orb);
+        else if(orbconvex)
+            subnfp = Nfp::noFitPolygon<NfpLevel::ONE_CONVEX>(stat, orb);
+        else
+            subnfp = Nfp::noFitPolygon<Level::value>(stat, orb);
 
-        RawShape subnfp;
-//        if(sh.isContourConvex() && trsh.isContourConvex()) {
-//            subnfp = Nfp::noFitPolygon<NfpLevel::CONVEX_ONLY>(
-//                        sh.transformedShape(), trsh.transformedShape());
-//        } else {
-            subnfp = Nfp::noFitPolygon<Level::value>( sh/*.transformedShape()*/,
-                                                      trsh.transformedShape());
-//        }
+        correctNfpPosition(subnfp, sh, trsh);
 
-//        #ifndef NDEBUG
-//            auto vv = ShapeLike::isValid(sh.transformedShape());
-//            assert(vv.first);
-
-//            auto vnfp = ShapeLike::isValid(subnfp);
-//            assert(vnfp.first);
-//        #endif
-
-//            auto vnfp = ShapeLike::isValid(subnfp);
-//            if(!vnfp.first) {
-//                std::cout << vnfp.second << std::endl;
-//                std::cout << ShapeLike::toString(subnfp) << std::endl;
-//            }
-
-        nfps = Nfp::merge(nfps, subnfp);
+        nfps = Nfp::merge(nfps, subnfp.first);
     }
 
     return nfps;
+
+
+//    using Item = _Item<RawShape>;
+//    using sl = ShapeLike;
+
+//    Nfp::Shapes<RawShape> nfps, stationary;
+
+//    for(Item& sh : polygons) {
+//        stationary = Nfp::merge(stationary, sh.transformedShape());
+//    }
+
+//    for(RawShape& sh : stationary) {
+
+////        auto vv = sl::isValid(sh);
+////        std::cout << vv.second << std::endl;
+
+
+//        Nfp::NfpResult<RawShape> subnfp;
+//        bool shconvex = sl::isConvex<RawShape>(sl::getContour(sh));
+//        if(shconvex && trsh.isContourConvex()) {
+//            subnfp = Nfp::noFitPolygon<NfpLevel::CONVEX_ONLY>(
+//                        sh, trsh.transformedShape());
+//        } else if(trsh.isContourConvex()) {
+//            subnfp = Nfp::noFitPolygon<NfpLevel::ONE_CONVEX>(
+//                        sh, trsh.transformedShape());
+//        }
+//        else {
+//            subnfp = Nfp::noFitPolygon<Level::value>( sh,
+//                                                      trsh.transformedShape());
+//        }
+
+//        correctNfpPosition(subnfp, sh, trsh);
+
+//        nfps = Nfp::merge(nfps, subnfp.first);
+//    }
+
+//    return nfps;
 }
 
 template<class RawShape>
@@ -283,13 +437,21 @@ class _NofitPolyPlacer: public PlacerBoilerplate<_NofitPolyPlacer<RawShape>,
 
 public:
 
-    using Pile = const Nfp::Shapes<RawShape>&;
+    using Pile = Nfp::Shapes<RawShape>;
 
     inline explicit _NofitPolyPlacer(const BinType& bin):
         Base(bin),
         norm_(std::sqrt(ShapeLike::area<RawShape>(bin))),
         penality_(1e6*norm_) {}
 
+    _NofitPolyPlacer(const _NofitPolyPlacer&) = default;
+    _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;
+#endif
+
     bool static inline wouldFit(const RawShape& chull, const RawShape& bin) {
         auto bbch = ShapeLike::boundingBox<RawShape>(chull);
         auto bbin = ShapeLike::boundingBox<RawShape>(bin);
@@ -363,7 +525,7 @@ public:
                 auto getNfpPoint = [&ecache](const Optimum& opt)
                 {
                     return opt.hidx < 0? ecache[opt.nfpidx].coords(opt.relpos) :
-                            ecache[opt.nfpidx].coords(opt.nfpidx, opt.relpos);
+                            ecache[opt.nfpidx].coords(opt.hidx, opt.relpos);
                 };
 
                 Nfp::Shapes<RawShape> pile;
@@ -378,8 +540,9 @@ public:
                 // customizable by the library client
                 auto _objfunc = config_.object_function?
                             config_.object_function :
-                [this](const Nfp::Shapes<RawShape>& pile, double occupied_area,
-                            double /*norm*/, double penality)
+                [this](Nfp::Shapes<RawShape>& pile, Item,
+                            double occupied_area, double /*norm*/,
+                            double penality)
                 {
                     auto ch = ShapeLike::convexHull(pile);
 
@@ -406,22 +569,22 @@ public:
                     d += startpos;
                     item.translation(d);
 
-                    pile.emplace_back(item.transformedShape());
+//                    pile.emplace_back(item.transformedShape());
 
                     double occupied_area = pile_area + item.area();
 
-                    double score = _objfunc(pile, occupied_area,
+                    double score = _objfunc(pile, item, occupied_area,
                                             norm_, penality_);
 
-                    pile.pop_back();
+//                    pile.pop_back();
 
                     return score;
                 };
 
                 opt::StopCriteria stopcr;
                 stopcr.max_iterations = 1000;
-                stopcr.stoplimit = 0.001;
-                stopcr.type = opt::StopLimitType::RELATIVE;
+                stopcr.absolute_score_difference = 1e-20*norm_;
+//                stopcr.relative_score_difference = 1e-20;
                 opt::TOptimizer<opt::Method::L_SIMPLEX> solver(stopcr);
 
                 Optimum optimum(0, 0);
@@ -458,6 +621,14 @@ public:
                         } catch(std::exception& e) {
                             derr() << "ERROR: " << e.what() << "\n";
                         }
+
+//                        auto sc = contour_ofn(pos);
+//                        if(sc < best_score) {
+//                            best_score = sc;
+//                            optimum.relpos = pos;
+//                            optimum.nfpidx = ch;
+//                            optimum.hidx = -1;
+//                        }
                     });
 
                     for(unsigned hidx = 0; hidx < cache.holeCount(); ++hidx) {
@@ -490,6 +661,13 @@ public:
                             } catch(std::exception& e) {
                                 derr() << "ERROR: " << e.what() << "\n";
                             }
+//                            auto sc = hole_ofn(pos);
+//                            if(sc < best_score) {
+//                                best_score = sc;
+//                                optimum.relpos = pos;
+//                                optimum.nfpidx = ch;
+//                                optimum.hidx = hidx;
+//                            }
                         });
                     }
                 }
diff --git a/xs/src/libnest2d/libnest2d/selections/djd_heuristic.hpp b/xs/src/libnest2d/libnest2d/selections/djd_heuristic.hpp
index 1d233cf35..17ac1167d 100644
--- a/xs/src/libnest2d/libnest2d/selections/djd_heuristic.hpp
+++ b/xs/src/libnest2d/libnest2d/selections/djd_heuristic.hpp
@@ -256,14 +256,14 @@ public:
 
             if(not_packed.size() < 2)
                 return false; // No group of two items
-            else {
-                double largest_area = not_packed.front().get().area();
-                auto itmp = not_packed.begin(); itmp++;
-                double second_largest = itmp->get().area();
-                if( free_area - second_largest - largest_area > waste)
-                    return false; // If even the largest two items do not fill
-                    // the bin to the desired waste than we can end here.
-            }
+
+            double largest_area = not_packed.front().get().area();
+            auto itmp = not_packed.begin(); itmp++;
+            double second_largest = itmp->get().area();
+            if( free_area - second_largest - largest_area > waste)
+                return false; // If even the largest two items do not fill
+                // the bin to the desired waste than we can end here.
+
 
             bool ret = false;
             auto it = not_packed.begin();
@@ -481,7 +481,7 @@ public:
                             {
                                 std::array<bool, 3> packed = {false};
 
-                                for(auto id : idx) packed[id] =
+                                for(auto id : idx) packed.at(id) =
                                         placer.pack(candidates[id]);
 
                                 bool check =
@@ -537,8 +537,7 @@ public:
             while (it != store_.end()) {
                 Placer p(bin);
                 if(!p.pack(*it)) {
-                    auto itmp = it++;
-                    store_.erase(itmp);
+                    it = store_.erase(it);
                 } else it++;
             }
         }
@@ -605,8 +604,7 @@ public:
                         if(placer.pack(*it)) {
                             filled_area += it->get().area();
                             free_area = bin_area - filled_area;
-                            auto itmp = it++;
-                            not_packed.erase(itmp);
+                            it = not_packed.erase(it);
                             makeProgress(placer, idx, 1);
                         } else it++;
                     }
diff --git a/xs/src/libnest2d/libnest2d/selections/firstfit.hpp b/xs/src/libnest2d/libnest2d/selections/firstfit.hpp
index b6e80520c..f34961f80 100644
--- a/xs/src/libnest2d/libnest2d/selections/firstfit.hpp
+++ b/xs/src/libnest2d/libnest2d/selections/firstfit.hpp
@@ -52,7 +52,7 @@ public:
         auto total = last-first;
         auto makeProgress = [this, &total](Placer& placer, size_t idx) {
             packed_bins_[idx] = placer.getItems();
-            this->progress_(--total);
+            this->progress_(static_cast<unsigned>(--total));
         };
 
         // Safety test: try to pack each item into an empty bin. If it fails
@@ -61,8 +61,7 @@ public:
             while (it != store_.end()) {
                 Placer p(bin);
                 if(!p.pack(*it)) {
-                    auto itmp = it++;
-                    store_.erase(itmp);
+                    it = store_.erase(it);
                 } else it++;
             }
         }
diff --git a/xs/src/libnest2d/tests/test.cpp b/xs/src/libnest2d/tests/test.cpp
index b37274f84..39315ff1a 100644
--- a/xs/src/libnest2d/tests/test.cpp
+++ b/xs/src/libnest2d/tests/test.cpp
@@ -682,7 +682,9 @@ void testNfp(const std::vector<ItemPair>& testdata) {
         auto&& nfp = Nfp::noFitPolygon<lvl>(stationary.rawShape(),
                                             orbiter.transformedShape());
 
-        auto v = ShapeLike::isValid(nfp);
+        strategies::correctNfpPosition(nfp, stationary, orbiter);
+
+        auto v = ShapeLike::isValid(nfp.first);
 
         if(!v.first) {
             std::cout << v.second << std::endl;
@@ -690,7 +692,7 @@ void testNfp(const std::vector<ItemPair>& testdata) {
 
         ASSERT_TRUE(v.first);
 
-        Item infp(nfp);
+        Item infp(nfp.first);
 
         int i = 0;
         auto rorbiter = orbiter.transformedShape();
@@ -742,6 +744,15 @@ TEST(GeometryAlgorithms, nfpConvexConvex) {
 //    testNfp<NfpLevel::BOTH_CONCAVE, 1000>(nfp_concave_testdata);
 //}
 
+TEST(GeometryAlgorithms, nfpConcaveConcave) {
+    using namespace libnest2d;
+
+//    Rectangle r1(10, 10);
+//    Rectangle r2(20, 20);
+//    auto result = Nfp::nfpSimpleSimple(r1.transformedShape(),
+//                                       r2.transformedShape());
+}
+
 TEST(GeometryAlgorithms, pointOnPolygonContour) {
     using namespace libnest2d;
 
diff --git a/xs/src/libnest2d/tools/libnfpglue.cpp b/xs/src/libnest2d/tools/libnfpglue.cpp
index 4cbdb5442..18656fd40 100644
--- a/xs/src/libnest2d/tools/libnfpglue.cpp
+++ b/xs/src/libnest2d/tools/libnfpglue.cpp
@@ -49,18 +49,18 @@ libnfporb::point_t scale(const libnfporb::point_t& p, long double factor) {
     long double px = p.x_.val();
     long double py = p.y_.val();
 #endif
-    return libnfporb::point_t(px*factor, py*factor);
+    return {px*factor, py*factor};
 }
 
 }
 
-PolygonImpl _nfp(const PolygonImpl &sh, const PolygonImpl &cother)
+NfpR _nfp(const PolygonImpl &sh, const PolygonImpl &cother)
 {
     using Vertex = PointImpl;
 
-    PolygonImpl ret;
+    NfpR ret;
 
-//    try {
+    try {
         libnfporb::polygon_t pstat, porb;
 
         boost::geometry::convert(sh, pstat);
@@ -85,7 +85,7 @@ PolygonImpl _nfp(const PolygonImpl &sh, const PolygonImpl &cother)
         // this can throw
         auto nfp = libnfporb::generateNFP(pstat, porb, true);
 
-        auto &ct = ShapeLike::getContour(ret);
+        auto &ct = ShapeLike::getContour(ret.first);
         ct.reserve(nfp.front().size()+1);
         for(auto v : nfp.front()) {
             v = scale(v, refactor);
@@ -94,10 +94,10 @@ PolygonImpl _nfp(const PolygonImpl &sh, const PolygonImpl &cother)
         ct.push_back(ct.front());
         std::reverse(ct.begin(), ct.end());
 
-        auto &rholes = ShapeLike::holes(ret);
+        auto &rholes = ShapeLike::holes(ret.first);
         for(size_t hidx = 1; hidx < nfp.size(); ++hidx) {
             if(nfp[hidx].size() >= 3) {
-                rholes.push_back({});
+                rholes.emplace_back();
                 auto& h = rholes.back();
                 h.reserve(nfp[hidx].size()+1);
 
@@ -110,73 +110,48 @@ PolygonImpl _nfp(const PolygonImpl &sh, const PolygonImpl &cother)
             }
         }
 
-        auto& cmp = vsort;
-        std::sort(pstat.outer().begin(), pstat.outer().end(), cmp);
-        std::sort(porb.outer().begin(), porb.outer().end(), cmp);
+        ret.second = Nfp::referenceVertex(ret.first);
 
-        // leftmost lower vertex of the stationary polygon
-        auto& touch_sh = scale(pstat.outer().back(), refactor);
-        // rightmost upper vertex of the orbiting polygon
-        auto& touch_other = scale(porb.outer().front(), refactor);
-
-        // Calculate the difference and move the orbiter to the touch position.
-        auto dtouch = touch_sh - touch_other;
-        auto _top_other = scale(porb.outer().back(), refactor) + dtouch;
-
-        Vertex top_other(getX(_top_other), getY(_top_other));
-
-        // Get the righmost upper vertex of the nfp and move it to the RMU of
-        // the orbiter because they should coincide.
-        auto&& top_nfp = Nfp::rightmostUpVertex(ret);
-        auto dnfp = top_other - top_nfp;
-
-        std::for_each(ShapeLike::begin(ret), ShapeLike::end(ret),
-                      [&dnfp](Vertex& v) { v+= dnfp; } );
-
-        for(auto& h : ShapeLike::holes(ret))
-            std::for_each( h.begin(), h.end(),
-                           [&dnfp](Vertex& v) { v += dnfp; } );
-
-//    } catch(std::exception& e) {
-//        std::cout << "Error: " << e.what() << "\nTrying with convex hull..." << std::endl;
+    } catch(std::exception& e) {
+        std::cout << "Error: " << e.what() << "\nTrying with convex hull..." << std::endl;
 //        auto ch_stat = ShapeLike::convexHull(sh);
 //        auto ch_orb = ShapeLike::convexHull(cother);
-//        ret = Nfp::nfpConvexOnly(ch_stat, ch_orb);
-//    }
+        ret = Nfp::nfpConvexOnly(sh, cother);
+    }
 
     return ret;
 }
 
-PolygonImpl Nfp::NfpImpl<PolygonImpl, NfpLevel::CONVEX_ONLY>::operator()(
+NfpR Nfp::NfpImpl<PolygonImpl, NfpLevel::CONVEX_ONLY>::operator()(
         const PolygonImpl &sh, const ClipperLib::PolygonImpl &cother)
 {
     return _nfp(sh, cother);//nfpConvexOnly(sh, cother);
 }
 
-PolygonImpl Nfp::NfpImpl<PolygonImpl, NfpLevel::ONE_CONVEX>::operator()(
+NfpR Nfp::NfpImpl<PolygonImpl, NfpLevel::ONE_CONVEX>::operator()(
         const PolygonImpl &sh, const ClipperLib::PolygonImpl &cother)
 {
     return _nfp(sh, cother);
 }
 
-PolygonImpl Nfp::NfpImpl<PolygonImpl, NfpLevel::BOTH_CONCAVE>::operator()(
+NfpR Nfp::NfpImpl<PolygonImpl, NfpLevel::BOTH_CONCAVE>::operator()(
         const PolygonImpl &sh, const ClipperLib::PolygonImpl &cother)
 {
     return _nfp(sh, cother);
 }
 
-PolygonImpl
-Nfp::NfpImpl<PolygonImpl, NfpLevel::ONE_CONVEX_WITH_HOLES>::operator()(
-        const PolygonImpl &sh, const ClipperLib::PolygonImpl &cother)
-{
-    return _nfp(sh, cother);
-}
+//PolygonImpl
+//Nfp::NfpImpl<PolygonImpl, NfpLevel::ONE_CONVEX_WITH_HOLES>::operator()(
+//        const PolygonImpl &sh, const ClipperLib::PolygonImpl &cother)
+//{
+//    return _nfp(sh, cother);
+//}
 
-PolygonImpl
-Nfp::NfpImpl<PolygonImpl, NfpLevel::BOTH_CONCAVE_WITH_HOLES>::operator()(
-        const PolygonImpl &sh, const ClipperLib::PolygonImpl &cother)
-{
-    return _nfp(sh, cother);
-}
+//PolygonImpl
+//Nfp::NfpImpl<PolygonImpl, NfpLevel::BOTH_CONCAVE_WITH_HOLES>::operator()(
+//        const PolygonImpl &sh, const ClipperLib::PolygonImpl &cother)
+//{
+//    return _nfp(sh, cother);
+//}
 
 }
diff --git a/xs/src/libnest2d/tools/libnfpglue.hpp b/xs/src/libnest2d/tools/libnfpglue.hpp
index 87b0e0833..75f639445 100644
--- a/xs/src/libnest2d/tools/libnfpglue.hpp
+++ b/xs/src/libnest2d/tools/libnfpglue.hpp
@@ -5,37 +5,39 @@
 
 namespace libnest2d {
 
-PolygonImpl _nfp(const PolygonImpl& sh, const PolygonImpl& cother);
+using NfpR = Nfp::NfpResult<PolygonImpl>;
+
+NfpR _nfp(const PolygonImpl& sh, const PolygonImpl& cother);
 
 template<>
 struct Nfp::NfpImpl<PolygonImpl, NfpLevel::CONVEX_ONLY> {
-    PolygonImpl operator()(const PolygonImpl& sh, const PolygonImpl& cother);
+    NfpR operator()(const PolygonImpl& sh, const PolygonImpl& cother);
 };
 
 template<>
 struct Nfp::NfpImpl<PolygonImpl, NfpLevel::ONE_CONVEX> {
-    PolygonImpl operator()(const PolygonImpl& sh, const PolygonImpl& cother);
+    NfpR operator()(const PolygonImpl& sh, const PolygonImpl& cother);
 };
 
 template<>
 struct Nfp::NfpImpl<PolygonImpl, NfpLevel::BOTH_CONCAVE> {
-    PolygonImpl operator()(const PolygonImpl& sh, const PolygonImpl& cother);
+    NfpR operator()(const PolygonImpl& sh, const PolygonImpl& cother);
 };
 
-template<>
-struct Nfp::NfpImpl<PolygonImpl, NfpLevel::ONE_CONVEX_WITH_HOLES> {
-    PolygonImpl operator()(const PolygonImpl& sh, const PolygonImpl& cother);
-};
+//template<>
+//struct Nfp::NfpImpl<PolygonImpl, NfpLevel::ONE_CONVEX_WITH_HOLES> {
+//    NfpResult operator()(const PolygonImpl& sh, const PolygonImpl& cother);
+//};
 
-template<>
-struct Nfp::NfpImpl<PolygonImpl, NfpLevel::BOTH_CONCAVE_WITH_HOLES> {
-    PolygonImpl operator()(const PolygonImpl& sh, const PolygonImpl& cother);
-};
+//template<>
+//struct Nfp::NfpImpl<PolygonImpl, NfpLevel::BOTH_CONCAVE_WITH_HOLES> {
+//    NfpResult operator()(const PolygonImpl& sh, const PolygonImpl& cother);
+//};
 
 template<> struct Nfp::MaxNfpLevel<PolygonImpl> {
     static const BP2D_CONSTEXPR NfpLevel value =
 //            NfpLevel::CONVEX_ONLY;
-            NfpLevel::BOTH_CONCAVE_WITH_HOLES;
+            NfpLevel::BOTH_CONCAVE;
 };
 
 }
diff --git a/xs/src/libnest2d/tools/svgtools.hpp b/xs/src/libnest2d/tools/svgtools.hpp
index 273ecabac..3a83caa70 100644
--- a/xs/src/libnest2d/tools/svgtools.hpp
+++ b/xs/src/libnest2d/tools/svgtools.hpp
@@ -5,11 +5,17 @@
 #include <fstream>
 #include <string>
 
-#include <libnest2d.h>
+#include <libnest2d/libnest2d.hpp>
 
 namespace libnest2d { namespace svg {
 
+template<class RawShape>
 class SVGWriter {
+    using Item = _Item<RawShape>;
+    using Coord = TCoord<TPoint<RawShape>>;
+    using Box = _Box<TPoint<RawShape>>;
+    using PackGroup = _PackGroup<RawShape>;
+
 public:
 
     enum OrigoLocation {
diff --git a/xs/src/libslic3r/Model.cpp b/xs/src/libslic3r/Model.cpp
index 2f73cad20..6e5cadd0d 100644
--- a/xs/src/libslic3r/Model.cpp
+++ b/xs/src/libslic3r/Model.cpp
@@ -538,32 +538,98 @@ bool arrange(Model &model, coordf_t dist, const Slic3r::BoundingBoxf* bb,
     // effect, the arrange procedure is a lot faster (we do not need to
     // calculate the convex hulls)
     pcfg.object_function = [bin, hasbin](
-            NfpPlacer::Pile pile,   // The currently arranged pile
+            NfpPlacer::Pile& pile,   // The currently arranged pile
+            Item item,
             double /*area*/,        // Sum area of items (not needed)
             double norm,            // A norming factor for physical dimensions
             double penality)        // Min penality in case of bad arrangement
     {
-        auto bb = ShapeLike::boundingBox(pile);
+        using pl = PointLike;
 
-        // We get the current item that's being evaluated.
-        auto& sh = pile.back();
+        static const double BIG_ITEM_TRESHOLD = 0.2;
+        static const double GRAVITY_RATIO = 0.5;
+        static const double DENSITY_RATIO = 1.0 - GRAVITY_RATIO;
 
-        // We retrieve the reference point of this item
-        auto rv = Nfp::referenceVertex(sh);
+        // We will treat big items (compared to the print bed) differently
+        NfpPlacer::Pile bigs;
+        bigs.reserve(pile.size());
+        for(auto& p : pile) {
+            auto pbb = ShapeLike::boundingBox(p);
+            auto na = std::sqrt(pbb.width()*pbb.height())/norm;
+            if(na > BIG_ITEM_TRESHOLD) bigs.emplace_back(p);
+        }
 
-        // We get the distance of the reference point from the center of the
-        // heat bed
-        auto c = bin.center();
-        auto d = PointLike::distance(rv, c);
+        // Candidate item bounding box
+        auto ibb = item.boundingBox();
 
-        // The score will be the normalized distance which will be minimized,
-        // effectively creating a circle shaped pile of items
-        double score = double(d)/norm;
+        // Calculate the full bounding box of the pile with the candidate item
+        pile.emplace_back(item.transformedShape());
+        auto fullbb = ShapeLike::boundingBox(pile);
+        pile.pop_back();
+
+        // The bounding box of the big items (they will accumulate in the center
+        // of the pile
+        auto bigbb = bigs.empty()? fullbb : ShapeLike::boundingBox(bigs);
+
+        // The size indicator of the candidate item. This is not the area,
+        // but almost...
+        auto itemnormarea = std::sqrt(ibb.width()*ibb.height())/norm;
+
+        // Will hold the resulting score
+        double score = 0;
+
+        if(itemnormarea > BIG_ITEM_TRESHOLD) {
+            // This branch is for the bigger items..
+            // Here we will use the closest point of the item bounding box to
+            // the already arranged pile. So not the bb center nor the a choosen
+            // corner but whichever is the closest to the center. This will
+            // prevent unwanted strange arrangements.
+
+            auto minc = ibb.minCorner(); // bottom left corner
+            auto maxc = ibb.maxCorner(); // top right corner
+
+            // top left and bottom right corners
+            auto top_left = PointImpl{getX(minc), getY(maxc)};
+            auto bottom_right = PointImpl{getX(maxc), getY(minc)};
+
+            auto cc = fullbb.center(); // The gravity center
+
+            // Now the distnce of the gravity center will be calculated to the
+            // five anchor points and the smallest will be chosen.
+            std::array<double, 5> dists;
+            dists[0] = pl::distance(minc, cc);
+            dists[1] = pl::distance(maxc, cc);
+            dists[2] = pl::distance(ibb.center(), cc);
+            dists[3] = pl::distance(top_left, cc);
+            dists[4] = pl::distance(bottom_right, cc);
+
+            auto dist = *(std::min_element(dists.begin(), dists.end())) / norm;
+
+            // Density is the pack density: how big is the arranged pile
+            auto density = std::sqrt(fullbb.width()*fullbb.height()) / norm;
+
+            // The score is a weighted sum of the distance from pile center
+            // and the pile size
+            score = GRAVITY_RATIO * dist + DENSITY_RATIO * density;
+
+        } else if(itemnormarea < BIG_ITEM_TRESHOLD && bigs.empty()) {
+            // If there are no big items, only small, we should consider the
+            // density here as well to not get silly results
+            auto bindist = pl::distance(ibb.center(), bin.center()) / norm;
+            auto density = std::sqrt(fullbb.width()*fullbb.height()) / norm;
+            score = GRAVITY_RATIO * bindist + DENSITY_RATIO * density;
+        } else {
+            // 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 = pl::distance(ibb.center(), bigbb.center()) / norm;
+        }
 
         // If it does not fit into the print bed we will beat it
         // with a large penality. If we would not do this, there would be only
         // one big pile that doesn't care whether it fits onto the print bed.
-        if(hasbin && !NfpPlacer::wouldFit(bb, bin)) score = 2*penality - score;
+        if(!NfpPlacer::wouldFit(fullbb, bin)) score = 2*penality - score;
 
         return score;
     };