Incorporating performance optimizations from libnest2d

xs/src/libnest2d/CMakeLists.txt

@@ -35,6 +35,7 @@ set(LIBNEST2D_SRCFILES
     libnest2d/geometry_traits.hpp
     libnest2d/geometries_io.hpp
     libnest2d/common.hpp
+    libnest2d/optimizer.hpp
     libnest2d/placers/placer_boilerplate.hpp
     libnest2d/placers/bottomleftplacer.hpp
     libnest2d/placers/nfpplacer.hpp
@@ -43,6 +44,10 @@ set(LIBNEST2D_SRCFILES
     libnest2d/selections/filler.hpp
     libnest2d/selections/firstfit.hpp
     libnest2d/selections/djd_heuristic.hpp
+    libnest2d/optimizers/simplex.hpp
+    libnest2d/optimizers/subplex.hpp
+    libnest2d/optimizers/genetic.hpp
+    libnest2d/optimizers/nlopt_boilerplate.hpp
     )
 
 if((NOT LIBNEST2D_GEOMETRIES_TARGET) OR (LIBNEST2D_GEOMETRIES_TARGET STREQUAL ""))
@@ -70,15 +75,22 @@ endif()
 
 message(STATUS "clipper lib is: ${LIBNEST2D_GEOMETRIES_TARGET}")
 
+find_package(NLopt 1.4)
+if(NOT NLopt_FOUND)
+    message(STATUS "NLopt not found so downloading and automatic build is performed...")
+    include(DownloadNLopt)
+endif()
+find_package(Threads REQUIRED)
+
 add_library(libnest2d_static STATIC ${LIBNEST2D_SRCFILES} )
-target_link_libraries(libnest2d_static PRIVATE ${LIBNEST2D_GEOMETRIES_TARGET})
+target_link_libraries(libnest2d_static PRIVATE ${LIBNEST2D_GEOMETRIES_TARGET} ${NLopt_LIBS})
-target_include_directories(libnest2d_static PUBLIC ${CMAKE_SOURCE_DIR})
+target_include_directories(libnest2d_static PUBLIC ${CMAKE_SOURCE_DIR} ${NLopt_INCLUDE_DIR})
 set_target_properties(libnest2d_static PROPERTIES PREFIX "")
 
 if(LIBNEST2D_BUILD_SHARED_LIB)
     add_library(libnest2d SHARED ${LIBNEST2D_SRCFILES} )
-    target_link_libraries(libnest2d PRIVATE ${LIBNEST2D_GEOMETRIES_TARGET})
+    target_link_libraries(libnest2d PRIVATE ${LIBNEST2D_GEOMETRIES_TARGET} ${NLopt_LIBS})
-    target_include_directories(libnest2d PUBLIC ${CMAKE_SOURCE_DIR})
+    target_include_directories(libnest2d PUBLIC ${CMAKE_SOURCE_DIR}  ${NLopt_INCLUDE_DIR})
     set_target_properties(libnest2d PROPERTIES PREFIX "")
 endif()
 
@@ -98,6 +110,6 @@ if(LIBNEST2D_BUILD_EXAMPLES)
                            tests/svgtools.hpp
                            tests/printer_parts.cpp
                            tests/printer_parts.h)
-    target_link_libraries(example libnest2d_static)
+    target_link_libraries(example libnest2d_static Threads::Threads)
     target_include_directories(example PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
 endif()

xs/src/libnest2d/README.md

@@ -1,29 +1,26 @@
 # Introduction
 
 Libnest2D is a library and framework for the 2D bin packaging problem. 
-Inspired from the [SVGNest](svgnest.com) Javascript library the project is is 
+Inspired from the [SVGNest](svgnest.com) Javascript library the project is  
 built from scratch in C++11. The library is written with a policy that it should
 be usable out of the box with a very simple interface but has to be customizable
-to the very core as well. This has led to a design where the algorithms are 
+to the very core as well. The algorithms are defined in a header only fashion 
-defined in a header only fashion with template only geometry types. These 
+with templated geometry types. These geometries can have custom or already 
-geometries can have custom or already existing implementation to avoid copying 
+existing implementation to avoid copying or having unnecessary dependencies.
-or having unnecessary dependencies.
 
-A default backend is provided if a user just wants to use the library out of the
+A default backend is provided if the user of the library just wants to use it 
-box without implementing the interface of these geometry types. The default 
+out of the box without additional integration. The default backend is reasonably 
-backend is built on top of boost geometry and the 
+fast and robust, being built on top of boost geometry and the 
-[polyclipping](http://www.angusj.com/delphi/clipper.php) library and implies the
+[polyclipping](http://www.angusj.com/delphi/clipper.php) library. Usage of 
-dependency on these packages as well as the compilation of the backend (although
+this default backend implies the dependency on these packages as well as the 
-I may find a solution in the future to make the backend header only as well).
+compilation of the backend itself (The default backend is not yet header only).
 
-This software is currently under heavy construction and lacks a throughout 
+This software is currently under construction and lacks a throughout 
-documentation and some essential algorithms as well. At this point a fairly 
+documentation and some essential algorithms as well. At this stage it works well
-untested version of the DJD selection heuristic is working with a bottom-left 
+for rectangles and convex closed polygons without considering holes and 
-placing strategy which may produce usable arrangements in most cases. 
+concavities.
 
-The no-fit polygon based placement strategy will be implemented in the very near
+Holes and non-convex polygons will be usable in the near future as well.
-future which should produce high quality results for convex and non convex 
-polygons with holes as well.  
 
 # References
 - [SVGNest](https://github.com/Jack000/SVGnest)

xs/src/libnest2d/cmake_modules/DownloadNLopt.cmake

@@ -0,0 +1,31 @@
+include(DownloadProject)
+
+if (CMAKE_VERSION VERSION_LESS 3.2)
+    set(UPDATE_DISCONNECTED_IF_AVAILABLE "")
+else()
+    set(UPDATE_DISCONNECTED_IF_AVAILABLE "UPDATE_DISCONNECTED 1")
+endif()
+
+# set(NLopt_DIR ${CMAKE_BINARY_DIR}/nlopt)
+include(DownloadProject)
+download_project(   PROJ                nlopt
+                    GIT_REPOSITORY      https://github.com/stevengj/nlopt.git
+                    GIT_TAG             1fcbcbf2fe8e34234e016cc43a6c41d3e8453e1f #master #nlopt-2.4.2
+                    # CMAKE_CACHE_ARGS    -DBUILD_SHARED_LIBS:BOOL=OFF -DCMAKE_BUILD_TYPE:STRING=${CMAKE_BUILD_TYPE} -DCMAKE_INSTALL_PREFIX=${NLopt_DIR}
+                    ${UPDATE_DISCONNECTED_IF_AVAILABLE}
+)
+
+set(SHARED_LIBS_STATE BUILD_SHARED_LIBS)
+set(BUILD_SHARED_LIBS OFF CACHE BOOL "" FORCE)
+set(NLOPT_PYTHON OFF CACHE BOOL "" FORCE)
+set(NLOPT_OCTAVE OFF CACHE BOOL "" FORCE)
+set(NLOPT_MATLAB OFF CACHE BOOL "" FORCE)
+set(NLOPT_GUILE OFF CACHE BOOL "" FORCE)
+set(NLOPT_SWIG OFF CACHE BOOL "" FORCE)
+set(NLOPT_LINK_PYTHON OFF CACHE BOOL "" FORCE)
+
+add_subdirectory(${nlopt_SOURCE_DIR} ${nlopt_BINARY_DIR})
+
+set(NLopt_LIBS nlopt)
+set(NLopt_INCLUDE_DIR ${nlopt_BINARY_DIR})
+set(SHARED_LIBS_STATE ${SHARED_STATE})

xs/src/libnest2d/cmake_modules/FindNLopt.cmake

@@ -0,0 +1,125 @@
+#///////////////////////////////////////////////////////////////////////////
+#//-------------------------------------------------------------------------
+#//
+#// Description:
+#//      cmake module for finding NLopt installation
+#//      NLopt installation location is defined by environment variable $NLOPT
+#//
+#//      following variables are defined:
+#//      NLopt_DIR              - NLopt installation directory
+#//      NLopt_INCLUDE_DIR      - NLopt header directory
+#//      NLopt_LIBRARY_DIR      - NLopt library directory
+#//      NLopt_LIBS             - NLopt library files
+#//
+#//      Example usage:
+#//          find_package(NLopt 1.4 REQUIRED)
+#//
+#//
+#//-------------------------------------------------------------------------
+
+
+set(NLopt_FOUND        FALSE)
+set(NLopt_ERROR_REASON "")
+set(NLopt_DEFINITIONS  "")
+set(NLopt_LIBS)
+
+
+set(NLopt_DIR $ENV{NLOPT})
+if(NOT NLopt_DIR)
+
+	set(NLopt_FOUND TRUE)
+
+	set(_NLopt_LIB_NAMES "nlopt")
+	find_library(NLopt_LIBS
+		NAMES ${_NLopt_LIB_NAMES})
+	if(NOT NLopt_LIBS)
+		set(NLopt_FOUND FALSE)
+		set(NLopt_ERROR_REASON "${NLopt_ERROR_REASON} Cannot find NLopt library '${_NLopt_LIB_NAMES}'.")
+	else()
+		get_filename_component(NLopt_DIR ${NLopt_LIBS} PATH)
+	endif()
+	unset(_NLopt_LIB_NAMES)
+
+	set(_NLopt_HEADER_FILE_NAME "nlopt.hpp")
+	find_file(_NLopt_HEADER_FILE
+		NAMES ${_NLopt_HEADER_FILE_NAME})
+	if(NOT _NLopt_HEADER_FILE)
+		set(NLopt_FOUND FALSE)
+		set(NLopt_ERROR_REASON "${NLopt_ERROR_REASON} Cannot find NLopt header file '${_NLopt_HEADER_FILE_NAME}'.")
+	endif()
+	unset(_NLopt_HEADER_FILE_NAME)
+	unset(_NLopt_HEADER_FILE)
+
+	if(NOT NLopt_FOUND)
+		set(NLopt_ERROR_REASON "${NLopt_ERROR_REASON} NLopt not found in system directories (and environment variable NLOPT is not set).")
+	else()
+	get_filename_component(NLopt_INCLUDE_DIR ${_NLopt_HEADER_FILE} DIRECTORY )
+	endif()
+
+
+
+else()
+
+	set(NLopt_FOUND TRUE)
+
+	set(NLopt_INCLUDE_DIR "${NLopt_DIR}/include")
+	if(NOT EXISTS "${NLopt_INCLUDE_DIR}")
+		set(NLopt_FOUND FALSE)
+		set(NLopt_ERROR_REASON "${NLopt_ERROR_REASON} Directory '${NLopt_INCLUDE_DIR}' does not exist.")
+	endif()
+
+	set(NLopt_LIBRARY_DIR "${NLopt_DIR}/lib")
+	if(NOT EXISTS "${NLopt_LIBRARY_DIR}")
+		set(NLopt_FOUND FALSE)
+		set(NLopt_ERROR_REASON "${NLopt_ERROR_REASON} Directory '${NLopt_LIBRARY_DIR}' does not exist.")
+	endif()
+
+	set(_NLopt_LIB_NAMES "nlopt_cxx")
+	find_library(NLopt_LIBS
+		NAMES ${_NLopt_LIB_NAMES}
+		PATHS ${NLopt_LIBRARY_DIR}
+		NO_DEFAULT_PATH)
+	if(NOT NLopt_LIBS)
+		set(NLopt_FOUND FALSE)
+		set(NLopt_ERROR_REASON "${NLopt_ERROR_REASON} Cannot find NLopt library '${_NLopt_LIB_NAMES}' in '${NLopt_LIBRARY_DIR}'.")
+	endif()
+	unset(_NLopt_LIB_NAMES)
+
+	set(_NLopt_HEADER_FILE_NAME "nlopt.hpp")
+	find_file(_NLopt_HEADER_FILE
+		NAMES ${_NLopt_HEADER_FILE_NAME}
+		PATHS ${NLopt_INCLUDE_DIR}
+		NO_DEFAULT_PATH)
+	if(NOT _NLopt_HEADER_FILE)
+		set(NLopt_FOUND FALSE)
+		set(NLopt_ERROR_REASON "${NLopt_ERROR_REASON} Cannot find NLopt header file '${_NLopt_HEADER_FILE_NAME}' in '${NLopt_INCLUDE_DIR}'.")
+	endif()
+	unset(_NLopt_HEADER_FILE_NAME)
+	unset(_NLopt_HEADER_FILE)
+
+endif()
+
+
+# make variables changeable
+mark_as_advanced(
+	NLopt_INCLUDE_DIR
+	NLopt_LIBRARY_DIR
+	NLopt_LIBS
+	NLopt_DEFINITIONS
+	)
+
+
+# report result
+if(NLopt_FOUND)
+	message(STATUS "Found NLopt in '${NLopt_DIR}'.")
+	message(STATUS "Using NLopt include directory '${NLopt_INCLUDE_DIR}'.")
+	message(STATUS "Using NLopt library '${NLopt_LIBS}'.")
+else()
+	if(NLopt_FIND_REQUIRED)
+		message(FATAL_ERROR "Unable to find requested NLopt installation:${NLopt_ERROR_REASON}")
+	else()
+		if(NOT NLopt_FIND_QUIETLY)
+			message(STATUS "NLopt was not found:${NLopt_ERROR_REASON}")
+		endif()
+	endif()
+endif()

xs/src/libnest2d/libnest2d/boost_alg.hpp

@@ -5,6 +5,9 @@
     #include <sstream>
 #endif
 
+#ifdef __clang__
+#undef _MSC_EXTENSIONS
+#endif
 #include <boost/geometry.hpp>
 
 // this should be removed to not confuse the compiler
@@ -152,7 +155,7 @@ template<> struct indexed_access<bp2d::Segment, 0, 0> {
         return bp2d::getX(seg.first());
     }
     static inline void set(bp2d::Segment &seg, bp2d::Coord const& coord) {
-        bp2d::setX(seg.first(), coord);
+        auto p = seg.first(); bp2d::setX(p, coord); seg.first(p);
     }
 };
 
@@ -161,7 +164,7 @@ template<> struct indexed_access<bp2d::Segment, 0, 1> {
         return bp2d::getY(seg.first());
     }
     static inline void set(bp2d::Segment &seg, bp2d::Coord const& coord) {
-        bp2d::setY(seg.first(), coord);
+        auto p = seg.first(); bp2d::setY(p, coord); seg.first(p);
     }
 };
 
@@ -170,7 +173,7 @@ template<> struct indexed_access<bp2d::Segment, 1, 0> {
         return bp2d::getX(seg.second());
     }
     static inline void set(bp2d::Segment &seg, bp2d::Coord const& coord) {
-        bp2d::setX(seg.second(), coord);
+        auto p = seg.second();  bp2d::setX(p, coord); seg.second(p);
     }
 };
 
@@ -179,7 +182,7 @@ template<> struct indexed_access<bp2d::Segment, 1, 1> {
         return bp2d::getY(seg.second());
     }
     static inline void set(bp2d::Segment &seg, bp2d::Coord const& coord) {
-        bp2d::setY(seg.second(), coord);
+        auto p = seg.second(); bp2d::setY(p, coord); seg.second(p);
     }
 };
 
@@ -325,20 +328,23 @@ inline bool ShapeLike::intersects(const PathImpl& sh1,
 // Tell libnest2d how to make string out of a ClipperPolygon object
 template<>
 inline bool ShapeLike::intersects(const PolygonImpl& sh1,
-                                  const PolygonImpl& sh2) {
+                                  const PolygonImpl& sh2)
+{
     return boost::geometry::intersects(sh1, sh2);
 }
 
 // Tell libnest2d how to make string out of a ClipperPolygon object
 template<>
 inline bool ShapeLike::intersects(const bp2d::Segment& s1,
-                                  const bp2d::Segment& s2) {
+                                  const bp2d::Segment& s2)
+{
     return boost::geometry::intersects(s1, s2);
 }
 
 #ifndef DISABLE_BOOST_AREA
 template<>
-inline double ShapeLike::area(const PolygonImpl& shape) {
+inline double ShapeLike::area(const PolygonImpl& shape)
+{
     return boost::geometry::area(shape);
 }
 #endif
@@ -364,16 +370,25 @@ inline bool ShapeLike::touches( const PolygonImpl& sh1,
     return boost::geometry::touches(sh1, sh2);
 }
 
+template<>
+inline bool ShapeLike::touches( const PointImpl& point,
+                                const PolygonImpl& shape)
+{
+    return boost::geometry::touches(point, shape);
+}
+
 #ifndef DISABLE_BOOST_BOUNDING_BOX
 template<>
-inline bp2d::Box ShapeLike::boundingBox(const PolygonImpl& sh) {
+inline bp2d::Box ShapeLike::boundingBox(const PolygonImpl& sh)
+{
     bp2d::Box b;
     boost::geometry::envelope(sh, b);
     return b;
 }
 
 template<>
-inline bp2d::Box ShapeLike::boundingBox(const bp2d::Shapes& shapes) {
+inline bp2d::Box ShapeLike::boundingBox<PolygonImpl>(const bp2d::Shapes& shapes)
+{
     bp2d::Box b;
     boost::geometry::envelope(shapes, b);
     return b;
@@ -398,17 +413,19 @@ inline PolygonImpl ShapeLike::convexHull(const bp2d::Shapes& shapes)
 }
 #endif
 
+#ifndef DISABLE_BOOST_ROTATE
 template<>
 inline void ShapeLike::rotate(PolygonImpl& sh, const Radians& rads)
 {
     namespace trans = boost::geometry::strategy::transform;
 
     PolygonImpl cpy = sh;
-
     trans::rotate_transformer<boost::geometry::radian, Radians, 2, 2>
             rotate(rads);
+
     boost::geometry::transform(cpy, sh, rotate);
 }
+#endif
 
 #ifndef DISABLE_BOOST_TRANSLATE
 template<>

xs/src/libnest2d/libnest2d/clipper_backend/clipper_backend.cpp

@@ -1,15 +1,35 @@
 #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());

xs/src/libnest2d/libnest2d/clipper_backend/clipper_backend.hpp

@@ -103,7 +103,7 @@ inline TCoord<PointImpl>& PointLike::y(PointImpl& p)
 }
 
 template<>
-inline void ShapeLike::reserve(PolygonImpl& sh, unsigned long vertex_capacity)
+inline void ShapeLike::reserve(PolygonImpl& sh, size_t vertex_capacity)
 {
     return sh.Contour.reserve(vertex_capacity);
 }
@@ -164,6 +164,95 @@ inline void ShapeLike::offset(PolygonImpl& sh, TCoord<PointImpl> distance) {
 
 }
 
+// TODO : make it convex hull right and faster than boost
+
+//#define DISABLE_BOOST_CONVEX_HULL
+//inline TCoord<PointImpl> cross( const PointImpl &O,
+//                                const PointImpl &A,
+//                                const PointImpl &B)
+//{
+//    return (A.X - O.X) * (B.Y - O.Y) - (A.Y - O.Y) * (B.X - O.X);
+//}
+
+//template<>
+//inline PolygonImpl ShapeLike::convexHull(const PolygonImpl& sh)
+//{
+//    auto& P = sh.Contour;
+//    PolygonImpl ret;
+
+//    size_t n = P.size(), k = 0;
+//    if (n <= 3) return ret;
+
+//    auto& H = ret.Contour;
+//    H.resize(2*n);
+//    std::vector<unsigned> indices(P.size(), 0);
+//    std::iota(indices.begin(), indices.end(), 0);
+
+//    // Sort points lexicographically
+//    std::sort(indices.begin(), indices.end(), [&P](unsigned i1, unsigned i2){
+//        auto& p1 = P[i1], &p2 = P[i2];
+//        return p1.X < p2.X || (p1.X == p2.X && p1.Y < p2.Y);
+//    });
+
+//    // Build lower hull
+//    for (size_t i = 0; i < n; ++i) {
+//        while (k >= 2 && cross(H[k-2], H[k-1], P[i]) <= 0) k--;
+//        H[k++] = P[i];
+//    }
+
+//    // Build upper hull
+//    for (size_t i = n-1, t = k+1; i > 0; --i) {
+//        while (k >= t && cross(H[k-2], H[k-1], P[i-1]) <= 0) k--;
+//        H[k++] = P[i-1];
+//    }
+
+//    H.resize(k-1);
+//    return ret;
+//}
+
+//template<>
+//inline PolygonImpl ShapeLike::convexHull(
+//        const ShapeLike::Shapes<PolygonImpl>& shapes)
+//{
+//    PathImpl P;
+//    PolygonImpl ret;
+
+//    size_t n = 0, k = 0;
+//    for(auto& sh : shapes) { n += sh.Contour.size(); }
+
+//    P.reserve(n);
+//    for(auto& sh : shapes) { P.insert(P.end(),
+//                                      sh.Contour.begin(), sh.Contour.end()); }
+
+//    if (n <= 3) { ret.Contour = P; return ret; }
+
+//    auto& H = ret.Contour;
+//    H.resize(2*n);
+//    std::vector<unsigned> indices(P.size(), 0);
+//    std::iota(indices.begin(), indices.end(), 0);
+
+//    // Sort points lexicographically
+//    std::sort(indices.begin(), indices.end(), [&P](unsigned i1, unsigned i2){
+//        auto& p1 = P[i1], &p2 = P[i2];
+//        return p1.X < p2.X || (p1.X == p2.X && p1.Y < p2.Y);
+//    });
+
+//    // Build lower hull
+//    for (size_t i = 0; i < n; ++i) {
+//        while (k >= 2 && cross(H[k-2], H[k-1], P[i]) <= 0) k--;
+//        H[k++] = P[i];
+//    }
+
+//    // Build upper hull
+//    for (size_t i = n-1, t = k+1; i > 0; --i) {
+//        while (k >= t && cross(H[k-2], H[k-1], P[i-1]) <= 0) k--;
+//        H[k++] = P[i-1];
+//    }
+
+//    H.resize(k-1);
+//    return ret;
+//}
+
 template<>
 inline PolygonImpl& Nfp::minkowskiAdd(PolygonImpl& sh,
                                             const PolygonImpl& other)
@@ -263,8 +352,30 @@ inline void ShapeLike::translate(PolygonImpl& sh, const PointImpl& offs)
     for(auto& hole : sh.Childs) for(auto& p : hole->Contour) { p += offs; }
 }
 
-#define DISABLE_BOOST_NFP_MERGE
+#define DISABLE_BOOST_ROTATE
+template<>
+inline void ShapeLike::rotate(PolygonImpl& sh, const Radians& rads)
+{
+    using Coord = TCoord<PointImpl>;
 
+    auto cosa = rads.cos();//std::cos(rads);
+    auto sina = rads.sin(); //std::sin(rads);
+
+    for(auto& p : sh.Contour) {
+        p = {
+                static_cast<Coord>(p.X * cosa - p.Y * sina),
+                static_cast<Coord>(p.X * sina + p.Y * cosa)
+             };
+    }
+    for(auto& hole : sh.Childs) for(auto& p : hole->Contour) {
+        p = {
+            static_cast<Coord>(p.X * cosa - p.Y * sina),
+            static_cast<Coord>(p.X * sina + p.Y * cosa)
+             };
+    }
+}
+
+#define DISABLE_BOOST_NFP_MERGE
 template<> inline
 Nfp::Shapes<PolygonImpl> Nfp::merge(const Nfp::Shapes<PolygonImpl>& shapes,
                                     const PolygonImpl& sh)

xs/src/libnest2d/libnest2d/common.hpp

@@ -1,6 +1,15 @@
 #ifndef LIBNEST2D_CONFIG_HPP
 #define LIBNEST2D_CONFIG_HPP
 
+#ifndef NDEBUG
+#include <iostream>
+#endif
+
+#include <stdexcept>
+#include <string>
+#include <cmath>
+#include <type_traits>
+
 #if defined(_MSC_VER) &&  _MSC_VER <= 1800 || __cplusplus < 201103L
     #define BP2D_NOEXCEPT
     #define BP2D_CONSTEXPR
@@ -9,12 +18,50 @@
     #define BP2D_CONSTEXPR constexpr
 #endif
 
-#include <stdexcept>
+
-#include <string>
+/*
-#include <cmath>
+ * Debugging output dout and derr definition
+ */
+//#ifndef NDEBUG
+//#   define dout std::cout
+//#   define derr std::cerr
+//#else
+//#   define dout 0 && std::cout
+//#   define derr 0 && std::cerr
+//#endif
 
 namespace libnest2d {
 
+struct DOut {
+#ifndef NDEBUG
+    std::ostream& out = std::cout;
+#endif
+};
+
+struct DErr {
+#ifndef NDEBUG
+    std::ostream& out = std::cerr;
+#endif
+};
+
+template<class T>
+inline DOut&& operator<<( DOut&& out, T&& d) {
+#ifndef NDEBUG
+    out.out << d;
+#endif
+    return std::move(out);
+}
+
+template<class T>
+inline DErr&& operator<<( DErr&& out, T&& d) {
+#ifndef NDEBUG
+    out.out << d;
+#endif
+    return std::move(out);
+}
+inline DOut dout() { return DOut(); }
+inline DErr derr() { return DErr(); }
+
 template< class T >
 struct remove_cvref {
     using type = typename std::remove_cv<
@@ -24,9 +71,58 @@ struct remove_cvref {
 template< class T >
 using remove_cvref_t = typename remove_cvref<T>::type;
 
+template< class T >
+using remove_ref_t = typename std::remove_reference<T>::type;
+
 template<bool B, class T>
 using enable_if_t = typename std::enable_if<B, T>::type;
 
+template<class F, class...Args>
+struct invoke_result {
+    using type = typename std::result_of<F(Args...)>::type;
+};
+
+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.
@@ -35,12 +131,14 @@ using enable_if_t = typename std::enable_if<B, T>::type;
  */
 template<class T> struct always_false { enum { value = false }; };
 
-const auto BP2D_CONSTEXPR Pi = 3.141592653589793238463; // 2*std::acos(0);
+const double BP2D_CONSTEXPR Pi = 3.141592653589793238463; // 2*std::acos(0);
+const double BP2D_CONSTEXPR Pi_2 = 2*Pi;
 
 /**
  * @brief Only for the Radian and Degrees classes to behave as doubles.
  */
 class Double {
+protected:
   double val_;
 public:
   Double(): val_(double{}) { }
@@ -56,12 +154,29 @@ class Degrees;
  * @brief Data type representing radians. It supports conversion to degrees.
  */
 class Radians: public Double {
+    mutable double sin_ = std::nan(""), cos_ = std::nan("");
 public:
     Radians(double rads = Double() ): Double(rads) {}
     inline Radians(const Degrees& degs);
 
     inline operator Degrees();
     inline double toDegrees();
+
+    inline double sin() const {
+        if(std::isnan(sin_)) {
+            cos_ = std::cos(val_);
+            sin_ = std::sin(val_);
+        }
+        return sin_;
+    }
+
+    inline double cos() const {
+        if(std::isnan(cos_)) {
+            cos_ = std::cos(val_);
+            sin_ = std::sin(val_);
+        }
+        return cos_;
+    }
 };
 
 /**

xs/src/libnest2d/libnest2d/geometries_nfp.hpp

@@ -75,12 +75,12 @@ static RawShape noFitPolygon(const RawShape& sh, const RawShape& other) {
         RawShape rsh;   // Final nfp placeholder
         std::vector<Edge> edgelist;
 
-        size_t cap = ShapeLike::contourVertexCount(sh) +
+        auto cap = ShapeLike::contourVertexCount(sh) +
                 ShapeLike::contourVertexCount(other);
 
         // Reserve the needed memory
         edgelist.reserve(cap);
-        ShapeLike::reserve(rsh, cap);
+        ShapeLike::reserve(rsh, static_cast<unsigned long>(cap));
 
         { // place all edges from sh into edgelist
             auto first = ShapeLike::cbegin(sh);
@@ -208,9 +208,10 @@ static inline bool _vsort(const TPoint<RawShape>& v1,
                           const TPoint<RawShape>& v2)
 {
     using Coord = TCoord<TPoint<RawShape>>;
-    auto diff = getY(v1) - getY(v2);
+    Coord &&x1 = getX(v1), &&x2 = getX(v2), &&y1 = getY(v1), &&y2 = getY(v2);
+    auto diff = y1 - y2;
     if(std::abs(diff) <= std::numeric_limits<Coord>::epsilon())
-        return getX(v1) < getX(v2);
+        return x1 < x2;
 
     return diff < 0;
 }

xs/src/libnest2d/libnest2d/geometry_traits.hpp

@@ -20,17 +20,17 @@ template<class GeomType>
 using TCoord = typename CoordType<remove_cvref_t<GeomType>>::Type;
 
 /// Getting the type of point structure used by a shape.
-template<class Shape> struct PointType { using Type = void; };
+template<class Shape> struct PointType { /*using Type = void;*/ };
 
 /// TPoint<ShapeClass> as shorthand for `typename PointType<ShapeClass>::Type`.
 template<class Shape>
 using TPoint = typename PointType<remove_cvref_t<Shape>>::Type;
 
 /// Getting the VertexIterator type of a shape class.
-template<class Shape> struct VertexIteratorType { using Type = void; };
+template<class Shape> struct VertexIteratorType { /*using Type = void;*/ };
 
 /// Getting the const vertex iterator for a shape class.
-template<class Shape> struct VertexConstIteratorType { using Type = void; };
+template<class Shape> struct VertexConstIteratorType {/* using Type = void;*/ };
 
 /**
  * TVertexIterator<Shape> as shorthand for
@@ -86,23 +86,47 @@ public:
     inline RawPoint center() const BP2D_NOEXCEPT;
 };
 
+/**
+ * \brief An abstraction of a directed line segment with two points.
+ */
 template<class RawPoint>
 class _Segment: PointPair<RawPoint> {
     using PointPair<RawPoint>::p1;
     using PointPair<RawPoint>::p2;
+    mutable Radians angletox_ = std::nan("");
 public:
 
     inline _Segment() {}
+
     inline _Segment(const RawPoint& p, const RawPoint& pp):
         PointPair<RawPoint>({p, pp}) {}
 
+    /**
+     * @brief Get the first point.
+     * @return Returns the starting point.
+     */
     inline const RawPoint& first() const BP2D_NOEXCEPT { return p1; }
+
+    /**
+     * @brief The end point.
+     * @return Returns the end point of the segment.
+     */
     inline const RawPoint& second() const BP2D_NOEXCEPT { return p2; }
 
-    inline RawPoint& first() BP2D_NOEXCEPT { return p1; }
+    inline void first(const RawPoint& p) BP2D_NOEXCEPT
-    inline RawPoint& second() BP2D_NOEXCEPT { return p2; }
+    {
+        angletox_ = std::nan(""); p1 = p;
+    }
 
+    inline void second(const RawPoint& p) BP2D_NOEXCEPT {
+        angletox_ = std::nan(""); p2 = p;
+    }
+
+    /// Returns the angle measured to the X (horizontal) axis.
     inline Radians angleToXaxis() const;
+
+    /// The length of the segment in the measure of the coordinate system.
+    inline double length();
 };
 
 // This struct serves as a namespace. The only difference is that is can be
@@ -204,13 +228,14 @@ struct PointLike {
 };
 
 template<class RawPoint>
-TCoord<RawPoint> _Box<RawPoint>::width() const BP2D_NOEXCEPT {
+TCoord<RawPoint> _Box<RawPoint>::width() const BP2D_NOEXCEPT
+{
     return PointLike::x(maxCorner()) - PointLike::x(minCorner());
 }
 
-
 template<class RawPoint>
-TCoord<RawPoint> _Box<RawPoint>::height() const BP2D_NOEXCEPT {
+TCoord<RawPoint> _Box<RawPoint>::height() const BP2D_NOEXCEPT
+{
     return PointLike::y(maxCorner()) - PointLike::y(minCorner());
 }
 
@@ -221,33 +246,37 @@ template<class RawPoint>
 TCoord<RawPoint> getY(const RawPoint& p) { return PointLike::y<RawPoint>(p); }
 
 template<class RawPoint>
-void setX(RawPoint& p, const TCoord<RawPoint>& val) {
+void setX(RawPoint& p, const TCoord<RawPoint>& val)
+{
     PointLike::x<RawPoint>(p) = val;
 }
 
 template<class RawPoint>
-void setY(RawPoint& p, const TCoord<RawPoint>& val) {
+void setY(RawPoint& p, const TCoord<RawPoint>& val)
+{
     PointLike::y<RawPoint>(p) = val;
 }
 
 template<class RawPoint>
 inline Radians _Segment<RawPoint>::angleToXaxis() const
 {
-    static const double Pi_2 = 2*Pi;
+    if(std::isnan(angletox_)) {
         TCoord<RawPoint> dx = getX(second()) - getX(first());
         TCoord<RawPoint> dy = getY(second()) - getY(first());
 
         double a = std::atan2(dy, dx);
-//    if(dx == 0 && dy >= 0) return Pi/2;
-//    if(dx == 0 && dy < 0) return 3*Pi/2;
-//    if(dy == 0 && dx >= 0) return 0;
-//    if(dy == 0 && dx < 0) return Pi;
-
-//    double ddx = static_cast<double>(dx);
         auto s = std::signbit(a);
-//    double a = std::atan(ddx/dy);
+
         if(s) a += Pi_2;
-    return a;
+        angletox_ = a;
+    }
+    return angletox_;
+}
+
+template<class RawPoint>
+inline double _Segment<RawPoint>::length()
+{
+    return PointLike::distance(first(), second());
 }
 
 template<class RawPoint>
@@ -319,7 +348,7 @@ struct ShapeLike {
 
     // Optional, does nothing by default
     template<class RawShape>
-    static void reserve(RawShape& /*sh*/,  unsigned long /*vertex_capacity*/) {}
+    static void reserve(RawShape& /*sh*/,  size_t /*vertex_capacity*/) {}
 
     template<class RawShape, class...Args>
     static void addVertex(RawShape& sh, Args...args)
@@ -415,6 +444,15 @@ struct ShapeLike {
         return false;
     }
 
+    template<class RawShape>
+    static bool touches( const TPoint<RawShape>& /*point*/,
+                         const RawShape& /*shape*/)
+    {
+        static_assert(always_false<RawShape>::value,
+                      "ShapeLike::touches(point, shape) unimplemented!");
+        return false;
+    }
+
     template<class RawShape>
     static _Box<TPoint<RawShape>> boundingBox(const RawShape& /*sh*/)
     {

xs/src/libnest2d/libnest2d/libnest2d.hpp

@@ -9,6 +9,9 @@
 #include <functional>
 
 #include "geometry_traits.hpp"
+//#include "optimizers/subplex.hpp"
+//#include "optimizers/simplex.hpp"
+#include "optimizers/genetic.hpp"
 
 namespace libnest2d {
 
@@ -48,6 +51,7 @@ class _Item {
     mutable bool area_cache_valid_ = false;
     mutable RawShape offset_cache_;
     mutable bool offset_cache_valid_ = false;
+
 public:
 
     /// The type of the shape which was handed over as the template argument.
@@ -188,7 +192,7 @@ public:
     }
 
     /// The number of the outer ring vertices.
-    inline unsigned long vertexCount() const {
+    inline size_t vertexCount() const {
         return ShapeLike::contourVertexCount(sh_);
     }
 
@@ -495,6 +499,8 @@ public:
     /// Clear the packed items so a new session can be started.
     inline void clearItems() { impl_.clearItems(); }
 
+    inline double filledArea() const { return impl_.filledArea(); }
+
 #ifndef NDEBUG
     inline auto getDebugItems() -> decltype(impl_.debug_items_)&
     {
@@ -629,7 +635,7 @@ template<class PlacementStrategy, class SelectionStrategy >
 class Arranger {
     using TSel = SelectionStrategyLike<SelectionStrategy>;
     TSel selector_;
-
+    bool use_min_bb_rotation_ = false;
 public:
     using Item = typename PlacementStrategy::Item;
     using ItemRef = std::reference_wrapper<Item>;
@@ -713,9 +719,9 @@ public:
     }
 
     /// Set a progress indicatior function object for the selector.
-    inline void progressIndicator(ProgressFunction func)
+    inline Arranger& progressIndicator(ProgressFunction func)
     {
-        selector_.progressIndicator(func);
+        selector_.progressIndicator(func); return *this;
     }
 
     inline PackGroup lastResult() {
@@ -727,6 +733,10 @@ public:
         return ret;
     }
 
+    inline Arranger& useMinimumBoundigBoxRotation(bool s = true) {
+        use_min_bb_rotation_ = s; return *this;
+    }
+
 private:
 
     template<class TIterator,
@@ -814,16 +824,42 @@ private:
         return pg;
     }
 
+    Radians findBestRotation(Item& item) {
+        opt::StopCriteria stopcr;
+        stopcr.stoplimit = 0.01;
+        stopcr.max_iterations = 10000;
+        stopcr.type = opt::StopLimitType::RELATIVE;
+        opt::TOptimizer<opt::Method::G_GENETIC> solver(stopcr);
+
+        auto orig_rot = item.rotation();
+
+        auto result = solver.optimize_min([&item, &orig_rot](Radians rot){
+            item.rotation(orig_rot + rot);
+            auto bb = item.boundingBox();
+            return std::sqrt(bb.height()*bb.width());
+        }, opt::initvals(Radians(0)), opt::bound<Radians>(-Pi/2, Pi/2));
+
+        item.rotation(orig_rot);
+
+        return std::get<0>(result.optimum);
+    }
+
     template<class TIter> inline void __arrange(TIter from, TIter to)
     {
         if(min_obj_distance_ > 0) std::for_each(from, to, [this](Item& item) {
             item.addOffset(static_cast<Unit>(std::ceil(min_obj_distance_/2.0)));
         });
 
+        if(use_min_bb_rotation_)
+            std::for_each(from, to, [this](Item& item){
+                Radians rot = findBestRotation(item);
+                item.rotate(rot);
+            });
+
         selector_.template packItems<PlacementStrategy>(
                     from, to, bin_, pconfig_);
 
-        if(min_obj_distance_ > 0) std::for_each(from, to, [this](Item& item) {
+        if(min_obj_distance_ > 0) std::for_each(from, to, [](Item& item) {
             item.removeOffset();
         });
 

xs/src/libnest2d/libnest2d/optimizer.hpp

@@ -0,0 +1,419 @@
+#ifndef OPTIMIZER_HPP
+#define OPTIMIZER_HPP
+
+#include <tuple>
+#include <functional>
+#include <limits>
+#include "common.hpp"
+
+namespace libnest2d { namespace opt {
+
+using std::forward;
+using std::tuple;
+using std::get;
+using std::tuple_element;
+
+/// A Type trait for upper and lower limit of a numeric type.
+template<class T, class B = void >
+struct limits {
+    inline static T min() { return std::numeric_limits<T>::min(); }
+    inline static T max() { return std::numeric_limits<T>::max(); }
+};
+
+template<class T>
+struct limits<T, enable_if_t<std::numeric_limits<T>::has_infinity, void>> {
+    inline static T min() { return -std::numeric_limits<T>::infinity(); }
+    inline static T max() { return std::numeric_limits<T>::infinity(); }
+};
+
+/// An interval of possible input values for optimization
+template<class T>
+class Bound {
+    T min_;
+    T max_;
+public:
+    Bound(const T& min = limits<T>::min(),
+          const T& max = limits<T>::max()): min_(min), max_(max) {}
+    inline const T min() const BP2D_NOEXCEPT { return min_; }
+    inline const T max() const BP2D_NOEXCEPT { return max_; }
+};
+
+/**
+ * Helper function to make a Bound object with its type deduced automatically.
+ */
+template<class T>
+inline Bound<T> bound(const T& min, const T& max) { return Bound<T>(min, max); }
+
+/**
+ * This is the type of an input tuple for the object function. It holds the
+ * values and their type in each dimension.
+ */
+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)...);
+}
+
+};
+
+/**
+ * @brief Specific optimization methods for which a default optimizer
+ * implementation can be instantiated.
+ */
+enum class Method {
+    L_SIMPLEX,
+    L_SUBPLEX,
+    G_GENETIC,
+    //...
+};
+
+/**
+ * @brief Info about result of an optimization. These codes are exactly the same
+ * as the nlopt codes for convinience.
+ */
+enum ResultCodes {
+    FAILURE = -1, /* generic failure code */
+    INVALID_ARGS = -2,
+    OUT_OF_MEMORY = -3,
+    ROUNDOFF_LIMITED = -4,
+    FORCED_STOP = -5,
+    SUCCESS = 1, /* generic success code */
+    STOPVAL_REACHED = 2,
+    FTOL_REACHED = 3,
+    XTOL_REACHED = 4,
+    MAXEVAL_REACHED = 5,
+    MAXTIME_REACHED = 6
+};
+
+/**
+ * \brief A type to hold the complete result of the optimization.
+ */
+template<class...Args>
+struct Result {
+    ResultCodes resultcode;
+    std::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;
+
+    /// The error value that is interpredted depending on the type property.
+    double stoplimit = 0.0001;
+
+    unsigned max_iterations = 0;
+};
+
+/**
+ * \brief The Optimizer base class with CRTP pattern.
+ */
+template<class Subclass>
+class Optimizer {
+protected:
+    enum class OptDir{
+        MIN,
+        MAX
+    } dir_;
+
+    StopCriteria stopcr_;
+
+public:
+
+    inline explicit Optimizer(const StopCriteria& scr = {}): stopcr_(scr) {}
+
+    /**
+     * \brief Optimize for minimum value of the provided objectfunction.
+     * \param objectfunction The function that will be searched for the minimum
+     * return value.
+     * \param initvals A tuple with the initial values for the search
+     * \param bounds A parameter pack with the bounds for each dimension.
+     * \return Returns a Result<Args...> structure.
+     * An example call would be:
+     *     auto result = opt.optimize_min(
+     *           [](std::tuple<double> x) // object function
+     *           {
+     *               return std::pow(std::get<0>(x), 2);
+     *           },
+     *           std::make_tuple(-0.5),  // initial value
+     *           {-1.0, 1.0}           // search space bounds
+     *       );
+     */
+    template<class Func, class...Args>
+    inline Result<Args...> optimize_min(Func&& objectfunction,
+                                        Input<Args...> initvals,
+                                        Bound<Args>... bounds)
+    {
+        dir_ = OptDir::MIN;
+        return static_cast<Subclass*>(this)->template optimize<Func, Args...>(
+                    forward<Func>(objectfunction), initvals, bounds... );
+    }
+
+    template<class Func, class...Args>
+    inline Result<Args...> optimize_min(Func&& objectfunction,
+                                        Input<Args...> initvals)
+    {
+        dir_ = OptDir::MIN;
+        return static_cast<Subclass*>(this)->template optimize<Func, Args...>(
+                    objectfunction, initvals, Bound<Args>()... );
+    }
+
+    template<class...Args, class Func>
+    inline Result<Args...> optimize_min(Func&& objectfunction)
+    {
+        dir_ = OptDir::MIN;
+        return static_cast<Subclass*>(this)->template optimize<Func, Args...>(
+                    objectfunction,
+                    Input<Args...>(),
+                    Bound<Args>()... );
+    }
+
+    /// Same as optimize_min but optimizes for maximum function value.
+    template<class Func, class...Args>
+    inline Result<Args...> optimize_max(Func&& objectfunction,
+                                        Input<Args...> initvals,
+                                        Bound<Args>... bounds)
+    {
+        dir_ = OptDir::MAX;
+        return static_cast<Subclass*>(this)->template optimize<Func, Args...>(
+                    objectfunction, initvals, bounds... );
+    }
+
+    template<class Func, class...Args>
+    inline Result<Args...> optimize_max(Func&& objectfunction,
+                                        Input<Args...> initvals)
+    {
+        dir_ = OptDir::MAX;
+        return static_cast<Subclass*>(this)->template optimize<Func, Args...>(
+                    objectfunction, initvals, Bound<Args>()... );
+    }
+
+    template<class...Args, class Func>
+    inline Result<Args...> optimize_max(Func&& objectfunction)
+    {
+        dir_ = OptDir::MAX;
+        return static_cast<Subclass*>(this)->template optimize<Func, Args...>(
+                    objectfunction,
+                    Input<Args...>(),
+                    Bound<Args>()... );
+    }
+
+};
+
+// Just to be able to instantiate an unimplemented method and generate compile
+// error.
+template<class T = void>
+class DummyOptimizer : public Optimizer<DummyOptimizer<T>> {
+    friend class Optimizer<DummyOptimizer<T>>;
+
+public:
+    DummyOptimizer() {
+        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)
+    {
+        return Result<Args...>();
+    }
+};
+
+// Specializing this struct will tell what kind of optimizer to generate for
+// a given method
+template<Method m> struct OptimizerSubclass { using Type = DummyOptimizer<>; };
+
+/// Optimizer type based on the method provided in parameter m.
+template<Method m> using TOptimizer = typename OptimizerSubclass<m>::Type;
+
+/// Global optimizer with an explicitly specified local method.
+template<Method m>
+inline TOptimizer<m> GlobalOptimizer(Method, const StopCriteria& scr = {})
+{ // Need to be specialized in order to do anything useful.
+    return TOptimizer<m>(scr);
+}
+
+}
+}
+
+#endif // OPTIMIZER_HPP

xs/src/libnest2d/libnest2d/optimizers/genetic.hpp

@@ -0,0 +1,31 @@
+#ifndef GENETIC_HPP
+#define GENETIC_HPP
+
+#include "nlopt_boilerplate.hpp"
+
+namespace libnest2d { namespace opt {
+
+class GeneticOptimizer: public NloptOptimizer {
+public:
+    inline explicit GeneticOptimizer(const StopCriteria& scr = {}):
+        NloptOptimizer(method2nloptAlg(Method::G_GENETIC), scr) {}
+
+    inline GeneticOptimizer& localMethod(Method m) {
+        localmethod_ = m;
+        return *this;
+    }
+};
+
+template<>
+struct OptimizerSubclass<Method::G_GENETIC> { using Type = GeneticOptimizer; };
+
+template<> TOptimizer<Method::G_GENETIC> GlobalOptimizer<Method::G_GENETIC>(
+        Method localm, const StopCriteria& scr )
+{
+    return GeneticOptimizer (scr).localMethod(localm);
+}
+
+}
+}
+
+#endif // GENETIC_HPP

xs/src/libnest2d/libnest2d/optimizers/nlopt_boilerplate.hpp

@@ -0,0 +1,176 @@
+#ifndef NLOPT_BOILERPLATE_HPP
+#define NLOPT_BOILERPLATE_HPP
+
+#include <nlopt.hpp>
+#include <libnest2d/optimizer.hpp>
+#include <cassert>
+
+#include <utility>
+
+namespace libnest2d { namespace opt {
+
+nlopt::algorithm method2nloptAlg(Method m) {
+
+    switch(m) {
+    case Method::L_SIMPLEX: return nlopt::LN_NELDERMEAD;
+    case Method::L_SUBPLEX: return nlopt::LN_SBPLX;
+    case Method::G_GENETIC: return nlopt::GN_ESCH;
+    default: assert(false); throw(m);
+    }
+}
+
+/**
+ * Optimizer based on NLopt.
+ *
+ * All the optimized types have to be convertible to double.
+ */
+class NloptOptimizer: public Optimizer<NloptOptimizer> {
+protected:
+    nlopt::opt opt_;
+    std::vector<double> lower_bounds_;
+    std::vector<double> upper_bounds_;
+    std::vector<double> initvals_;
+    nlopt::algorithm alg_;
+    Method localmethod_;
+
+    using Base = Optimizer<NloptOptimizer>;
+
+    friend Base;
+
+    // ********************************************************************** */
+
+    // TODO: CHANGE FOR LAMBDAS WHEN WE WILL MOVE TO C++14
+
+    struct BoundsFunc {
+        NloptOptimizer& self;
+        inline explicit BoundsFunc(NloptOptimizer& o): self(o) {}
+
+        template<class T> void operator()(int N, T& bounds)
+        {
+            self.lower_bounds_[N] = bounds.min();
+            self.upper_bounds_[N] = bounds.max();
+        }
+    };
+
+    struct InitValFunc {
+        NloptOptimizer& self;
+        inline explicit InitValFunc(NloptOptimizer& o): self(o) {}
+
+        template<class T> void operator()(int N, T& initval)
+        {
+            self.initvals_[N] = initval;
+        }
+    };
+
+    struct ResultCopyFunc {
+        NloptOptimizer& self;
+        inline explicit ResultCopyFunc(NloptOptimizer& o): self(o) {}
+
+        template<class T> void operator()(int N, T& resultval)
+        {
+            resultval = self.initvals_[N];
+        }
+    };
+
+    struct FunvalCopyFunc {
+        using D = const std::vector<double>;
+        D& params;
+        inline explicit FunvalCopyFunc(D& p): params(p) {}
+
+        template<class T> void operator()(int N, T& resultval)
+        {
+            resultval = params[N];
+        }
+    };
+
+    /* ********************************************************************** */
+
+    template<class Fn, class...Args>
+    static double optfunc(const std::vector<double>& params,
+                          std::vector<double>& grad,
+                          void *data)
+    {
+        auto fnptr = static_cast<remove_ref_t<Fn>*>(data);
+        auto funval = std::tuple<Args...>();
+
+        // copy the obtained objectfunction arguments to the funval tuple.
+        metaloop::apply(FunvalCopyFunc(params), funval);
+
+        auto ret = metaloop::callFunWithTuple(*fnptr, funval,
+                                              index_sequence_for<Args...>());
+
+        return ret;
+    }
+
+    template<class Func, class...Args>
+    Result<Args...> optimize(Func&& func,
+                             std::tuple<Args...> initvals,
+                             Bound<Args>... args)
+    {
+        lower_bounds_.resize(sizeof...(Args));
+        upper_bounds_.resize(sizeof...(Args));
+        initvals_.resize(sizeof...(Args));
+
+        opt_ = nlopt::opt(alg_, sizeof...(Args) );
+
+        // Copy the bounds which is obtained as a parameter pack in args into
+        // lower_bounds_ and upper_bounds_
+        metaloop::apply(BoundsFunc(*this), args...);
+
+        opt_.set_lower_bounds(lower_bounds_);
+        opt_.set_upper_bounds(upper_bounds_);
+
+        nlopt::opt localopt;
+        switch(opt_.get_algorithm()) {
+        case nlopt::GN_MLSL:
+        case nlopt::GN_MLSL_LDS:
+            localopt = nlopt::opt(method2nloptAlg(localmethod_),
+                                  sizeof...(Args));
+            localopt.set_lower_bounds(lower_bounds_);
+            localopt.set_upper_bounds(upper_bounds_);
+            opt_.set_local_optimizer(localopt);
+        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;
+        }
+
+        if(this->stopcr_.max_iterations > 0)
+            opt_.set_maxeval(this->stopcr_.max_iterations );
+
+        // Take care of the initial values, copy them to initvals_
+        metaloop::apply(InitValFunc(*this), initvals);
+
+        switch(dir_) {
+        case OptDir::MIN:
+            opt_.set_min_objective(optfunc<Func, Args...>, &func); break;
+        case OptDir::MAX:
+            opt_.set_max_objective(optfunc<Func, Args...>, &func); break;
+        }
+
+        Result<Args...> result;
+
+        auto rescode = opt_.optimize(initvals_, result.score);
+        result.resultcode = static_cast<ResultCodes>(rescode);
+
+        metaloop::apply(ResultCopyFunc(*this), result.optimum);
+
+        return result;
+    }
+
+public:
+    inline explicit NloptOptimizer(nlopt::algorithm alg,
+                                   StopCriteria stopcr = {}):
+        Base(stopcr), alg_(alg), localmethod_(Method::L_SIMPLEX) {}
+
+};
+
+}
+}
+
+
+#endif // NLOPT_BOILERPLATE_HPP

xs/src/libnest2d/libnest2d/optimizers/simplex.hpp

@@ -0,0 +1,20 @@
+#ifndef SIMPLEX_HPP
+#define SIMPLEX_HPP
+
+#include "nlopt_boilerplate.hpp"
+
+namespace libnest2d { namespace opt {
+
+class SimplexOptimizer: public NloptOptimizer {
+public:
+    inline explicit SimplexOptimizer(const StopCriteria& scr = {}):
+        NloptOptimizer(method2nloptAlg(Method::L_SIMPLEX), scr) {}
+};
+
+template<>
+struct OptimizerSubclass<Method::L_SIMPLEX> { using Type = SimplexOptimizer; };
+
+}
+}
+
+#endif // SIMPLEX_HPP

xs/src/libnest2d/libnest2d/optimizers/subplex.hpp

@@ -0,0 +1,20 @@
+#ifndef SUBPLEX_HPP
+#define SUBPLEX_HPP
+
+#include "nlopt_boilerplate.hpp"
+
+namespace libnest2d { namespace opt {
+
+class SubplexOptimizer: public NloptOptimizer {
+public:
+    inline explicit SubplexOptimizer(const StopCriteria& scr = {}):
+        NloptOptimizer(method2nloptAlg(Method::L_SUBPLEX), scr) {}
+};
+
+template<>
+struct OptimizerSubclass<Method::L_SUBPLEX> { using Type = SubplexOptimizer; };
+
+}
+}
+
+#endif // SUBPLEX_HPP

xs/src/libnest2d/libnest2d/placers/bottomleftplacer.hpp

@@ -348,8 +348,9 @@ protected:
         };*/
 
         auto reverseAddOthers = [&rsh, finish, start, &item](){
-            for(size_t i = finish-1; i > start; i--)
+            for(auto i = finish-1; i > start; i--)
-                ShapeLike::addVertex(rsh, item.vertex(i));
+                ShapeLike::addVertex(rsh, item.vertex(
+                                         static_cast<unsigned long>(i)));
         };
 
         // Final polygon construction...

xs/src/libnest2d/libnest2d/placers/nfpplacer.hpp

@@ -1,8 +1,13 @@
 #ifndef NOFITPOLY_HPP
 #define NOFITPOLY_HPP
 
+#ifndef NDEBUG
+#include <iostream>
+#endif
 #include "placer_boilerplate.hpp"
 #include "../geometries_nfp.hpp"
+#include <libnest2d/optimizers/subplex.hpp>
+//#include <libnest2d/optimizers/genetic.hpp>
 
 namespace libnest2d { namespace strategies {
 
@@ -17,10 +22,183 @@ struct NfpPConfig {
         TOP_RIGHT,
     };
 
-    bool allow_rotations = false;
+    /// Which angles to try out for better results
+    std::vector<Radians> rotations;
+
+    /// Where to align the resulting packed pile
     Alignment alignment;
+
+    NfpPConfig(): rotations({0.0, Pi/2.0, Pi, 3*Pi/2}),
+        alignment(Alignment::CENTER) {}
 };
 
+// A class for getting a point on the circumference of the polygon (in log time)
+template<class RawShape> class EdgeCache {
+    using Vertex = TPoint<RawShape>;
+    using Coord = TCoord<Vertex>;
+    using Edge = _Segment<Vertex>;
+
+    enum Corners {
+        BOTTOM,
+        LEFT,
+        RIGHT,
+        TOP,
+        NUM_CORNERS
+    };
+
+    mutable std::vector<double> corners_;
+
+    std::vector<Edge> emap_;
+    std::vector<double> distances_;
+    double full_distance_ = 0;
+
+    void createCache(const RawShape& sh) {
+        auto first = ShapeLike::cbegin(sh);
+        auto next = first + 1;
+        auto endit = ShapeLike::cend(sh);
+
+        distances_.reserve(ShapeLike::contourVertexCount(sh));
+
+        while(next != endit) {
+            emap_.emplace_back(*(first++), *(next++));
+            full_distance_ += emap_.back().length();
+            distances_.push_back(full_distance_);
+        }
+    }
+
+    void fetchCorners() const {
+        if(!corners_.empty()) return;
+
+        corners_ = std::vector<double>(NUM_CORNERS, 0.0);
+
+        std::vector<unsigned> idx_ud(emap_.size(), 0);
+        std::vector<unsigned> idx_lr(emap_.size(), 0);
+
+        std::iota(idx_ud.begin(), idx_ud.end(), 0);
+        std::iota(idx_lr.begin(), idx_lr.end(), 0);
+
+        std::sort(idx_ud.begin(), idx_ud.end(),
+                  [this](unsigned idx1, unsigned idx2)
+        {
+            const Vertex& v1 = emap_[idx1].first();
+            const Vertex& v2 = emap_[idx2].first();
+            auto diff = getY(v1) - getY(v2);
+            if(std::abs(diff) <= std::numeric_limits<Coord>::epsilon())
+              return getX(v1) < getX(v2);
+
+            return diff < 0;
+        });
+
+        std::sort(idx_lr.begin(), idx_lr.end(),
+                  [this](unsigned idx1, unsigned idx2)
+        {
+            const Vertex& v1 = emap_[idx1].first();
+            const Vertex& v2 = emap_[idx2].first();
+
+            auto diff = getX(v1) - getX(v2);
+            if(std::abs(diff) <= std::numeric_limits<Coord>::epsilon())
+                return getY(v1) < getY(v2);
+
+            return diff < 0;
+        });
+
+        corners_[BOTTOM] = distances_[idx_ud.front()]/full_distance_;
+        corners_[TOP] = distances_[idx_ud.back()]/full_distance_;
+        corners_[LEFT] = distances_[idx_lr.front()]/full_distance_;
+        corners_[RIGHT] = distances_[idx_lr.back()]/full_distance_;
+    }
+
+public:
+
+    using iterator = std::vector<double>::iterator;
+    using const_iterator = std::vector<double>::const_iterator;
+
+    inline EdgeCache() = default;
+
+    inline EdgeCache(const _Item<RawShape>& item)
+    {
+        createCache(item.transformedShape());
+    }
+
+    inline EdgeCache(const RawShape& sh)
+    {
+        createCache(sh);
+    }
+
+    /**
+     * @brief Get a point on the circumference of a polygon.
+     * @param distance A relative distance from the starting point to the end.
+     * Can be from 0.0 to 1.0 where 0.0 is the starting point and 1.0 is the
+     * closing point (which should be eqvivalent with the starting point with
+     * closed polygons).
+     * @return Returns the coordinates of the point lying on the polygon
+     * circumference.
+     */
+    inline Vertex coords(double distance) {
+        assert(distance >= .0 && distance <= 1.0);
+
+        // distance is from 0.0 to 1.0, we scale it up to the full length of
+        // the circumference
+        double d = distance*full_distance_;
+
+        // Magic: we find the right edge in log time
+        auto it = std::lower_bound(distances_.begin(), distances_.end(), d);
+        auto idx = it - distances_.begin();     // get the index of the edge
+        auto edge = emap_[idx];   // extrac the edge
+
+        // Get the remaining distance on the target edge
+        auto ed = d - (idx > 0 ? *std::prev(it) : 0 );
+        auto angle = edge.angleToXaxis();
+        Vertex ret = edge.first();
+
+        // Get the point on the edge which lies in ed distance from the start
+        ret += { static_cast<Coord>(std::round(ed*std::cos(angle))),
+                 static_cast<Coord>(std::round(ed*std::sin(angle))) };
+
+        return ret;
+    }
+
+    inline double circumference() const BP2D_NOEXCEPT { return full_distance_; }
+
+    inline double corner(Corners c) const BP2D_NOEXCEPT {
+        assert(c < NUM_CORNERS);
+        fetchCorners();
+        return corners_[c];
+    }
+
+    inline const std::vector<double>& corners() const BP2D_NOEXCEPT {
+        fetchCorners();
+        return corners_;
+    }
+
+};
+
+// Nfp for a bunch of polygons. If the polygons are convex, the nfp calculated
+// for trsh can be the union of nfp-s calculated with each polygon
+template<class RawShape, class Container>
+Nfp::Shapes<RawShape> nfp(const Container& polygons, const RawShape& trsh )
+{
+    using Item = _Item<RawShape>;
+
+    Nfp::Shapes<RawShape> nfps;
+
+    for(Item& sh : polygons) {
+        auto subnfp = Nfp::noFitPolygon(sh.transformedShape(),
+                                        trsh);
+        #ifndef NDEBUG
+            auto vv = ShapeLike::isValid(sh.transformedShape());
+            assert(vv.first);
+
+            auto vnfp = ShapeLike::isValid(subnfp);
+            assert(vnfp.first);
+        #endif
+
+        nfps = Nfp::merge(nfps, subnfp);
+    }
+
+    return nfps;
+}
+
 template<class RawShape>
 class _NofitPolyPlacer: public PlacerBoilerplate<_NofitPolyPlacer<RawShape>,
         RawShape, _Box<TPoint<RawShape>>, NfpPConfig<RawShape>> {
@@ -32,9 +210,30 @@ class _NofitPolyPlacer: public PlacerBoilerplate<_NofitPolyPlacer<RawShape>,
 
     using Box = _Box<TPoint<RawShape>>;
 
+    const double norm_;
+    const double penality_;
+
+    bool static wouldFit(const RawShape& chull, const RawShape& bin) {
+        auto bbch = ShapeLike::boundingBox<RawShape>(chull);
+        auto bbin = ShapeLike::boundingBox<RawShape>(bin);
+        auto d = bbin.minCorner() - bbch.minCorner();
+        auto chullcpy = chull;
+        ShapeLike::translate(chullcpy, d);
+        return ShapeLike::isInside<RawShape>(chullcpy, bbin);
+    }
+
+    bool static wouldFit(const RawShape& chull, const Box& bin)
+    {
+        auto bbch = ShapeLike::boundingBox<RawShape>(chull);
+        return bbch.width() <= bin.width() && bbch.height() <= bin.height();
+    }
+
 public:
 
-    inline explicit _NofitPolyPlacer(const BinType& bin): Base(bin) {}
+    inline explicit _NofitPolyPlacer(const BinType& bin):
+        Base(bin),
+        norm_(std::sqrt(ShapeLike::area<RawShape>(bin))),
+        penality_(1e6*norm_) {}
 
     PackResult trypack(Item& item) {
 
@@ -47,88 +246,148 @@ public:
             can_pack = item.isInside(bin_);
         } else {
 
-            // place the new item outside of the print bed to make sure it is
+            double global_score = penality_;
-            // disjuct from the current merged pile
+
+            auto initial_tr = item.translation();
+            auto initial_rot = item.rotation();
+            Vertex final_tr = {0, 0};
+            Radians final_rot = initial_rot;
+            Nfp::Shapes<RawShape> nfps;
+
+            for(auto rot : config_.rotations) {
+
+                item.translation(initial_tr);
+                item.rotation(initial_rot + rot);
+
+                // place the new item outside of the print bed to make sure
+                // it is disjuct from the current merged pile
                 placeOutsideOfBin(item);
 
                 auto trsh = item.transformedShape();
-            Nfp::Shapes<RawShape> nfps;
-
-#ifndef NDEBUG
-#ifdef DEBUG_EXPORT_NFP
-            Base::debug_items_.clear();
-#endif
-            auto v = ShapeLike::isValid(trsh);
-            assert(v.first);
-#endif
-            for(Item& sh : items_) {
-                auto subnfp = Nfp::noFitPolygon(sh.transformedShape(),
-                                                trsh);
-#ifndef NDEBUG
-#ifdef DEBUG_EXPORT_NFP
-                Base::debug_items_.emplace_back(subnfp);
-#endif
-                auto vv = ShapeLike::isValid(sh.transformedShape());
-                assert(vv.first);
-
-                auto vnfp = ShapeLike::isValid(subnfp);
-                assert(vnfp.first);
-#endif
-                nfps = Nfp::merge(nfps, subnfp);
-            }
-
-            double min_area = std::numeric_limits<double>::max();
-            Vertex tr = {0, 0};
 
+                nfps = nfp(items_, trsh);
                 auto iv = Nfp::referenceVertex(trsh);
 
-            // place item on each the edge of this nfp
+                auto startpos = item.translation();
-            for(auto& nfp : nfps)
-            ShapeLike::foreachContourVertex(nfp, [&]
-                                            (Vertex& v)
-            {
-                Coord dx = getX(v) - getX(iv);
-                Coord dy = getY(v) - getY(iv);
 
-                Item placeditem(trsh);
+                std::vector<EdgeCache<RawShape>> ecache;
-                placeditem.translate(Vertex(dx, dy));
+                ecache.reserve(nfps.size());
 
-                if( placeditem.isInside(bin_) ) {
+                for(auto& nfp : nfps ) ecache.emplace_back(nfp);
-                    Nfp::Shapes<RawShape> m;
-                    m.reserve(items_.size());
 
-                    for(Item& pi : items_)
+                auto getNfpPoint = [&ecache](double relpos) {
-                        m.emplace_back(pi.transformedShape());
+                    auto relpfloor = std::floor(relpos);
+                    auto nfp_idx = static_cast<unsigned>(relpfloor);
+                    if(nfp_idx >= ecache.size()) nfp_idx--;
+                    auto p = relpos - relpfloor;
+                    return ecache[nfp_idx].coords(p);
+                };
 
-                    m.emplace_back(placeditem.transformedShape());
+                Nfp::Shapes<RawShape> pile;
-
+                pile.reserve(items_.size()+1);
-//                    auto b = ShapeLike::boundingBox(m);
+                double pile_area = 0;
-
+                for(Item& mitem : items_) {
-//                    auto a = static_cast<double>(std::max(b.height(),
+                    pile.emplace_back(mitem.transformedShape());
-//                                                          b.width()));
+                    pile_area += mitem.area();
-
-                    auto b = ShapeLike::convexHull(m);
-                    auto a = ShapeLike::area(b);
-
-                    if(a < min_area) {
-                        can_pack = true;
-                        min_area = a;
-                        tr = {dx, dy};
                 }
+
+                // Our object function for placement
+                auto objfunc = [&] (double relpos)
+                {
+                    Vertex v = getNfpPoint(relpos);
+                    auto d = v - iv;
+                    d += startpos;
+                    item.translation(d);
+
+                    pile.emplace_back(item.transformedShape());
+
+                    double occupied_area = pile_area + item.area();
+                    auto ch = ShapeLike::convexHull(pile);
+
+                    pile.pop_back();
+
+                    // The pack ratio -- how much is the convex hull occupied
+                    double pack_rate = occupied_area/ShapeLike::area(ch);
+
+                    // ratio of waste
+                    double waste = 1.0 - pack_rate;
+
+                    // Score is the square root of waste. This will extend the
+                    // range of good (lower) values and shring the range of bad
+                    // (larger) values.
+                    auto score = std::sqrt(waste);
+
+                    if(!wouldFit(ch, bin_)) score = 2*penality_ - score;
+
+                    return score;
+                };
+
+                opt::StopCriteria stopcr;
+                stopcr.max_iterations = 1000;
+                stopcr.stoplimit = 0.01;
+                stopcr.type = opt::StopLimitType::RELATIVE;
+                opt::TOptimizer<opt::Method::L_SUBPLEX> solver(stopcr);
+
+                double optimum = 0;
+                double best_score = penality_;
+
+                // double max_bound = 1.0*nfps.size();
+                // Genetic should look like this:
+                /*auto result = solver.optimize_min(objfunc,
+                                opt::initvals<double>(0.0),
+                                opt::bound(0.0, max_bound)
+                                );
+
+                if(result.score < penality_) {
+                    best_score = result.score;
+                    optimum = std::get<0>(result.optimum);
+                }*/
+
+                // Local optimization with the four polygon corners as
+                // starting points
+                for(unsigned ch = 0; ch < ecache.size(); ch++) {
+                    auto& cache = ecache[ch];
+
+                    std::for_each(cache.corners().begin(),
+                                  cache.corners().end(),
+                                  [ch, &solver, &objfunc,
+                                  &best_score, &optimum]
+                                  (double pos)
+                    {
+                        try {
+                            auto result = solver.optimize_min(objfunc,
+                                            opt::initvals<double>(ch+pos),
+                                            opt::bound<double>(ch, 1.0 + ch)
+                                            );
+
+                            if(result.score < best_score) {
+                                best_score = result.score;
+                                optimum = std::get<0>(result.optimum);
+                            }
+                        } catch(std::exception&
+                        #ifndef NDEBUG
+                                e
+                        #endif
+                                ) {
+                        #ifndef NDEBUG
+                            std::cerr << "ERROR " << e.what() << std::endl;
+                        #endif
                         }
                     });
-
-#ifndef NDEBUG
-            for(auto&nfp : nfps) {
-                auto val = ShapeLike::isValid(nfp);
-                if(!val.first) std::cout << val.second << std::endl;
-#ifdef DEBUG_EXPORT_NFP
-                Base::debug_items_.emplace_back(nfp);
-#endif
                 }
-#endif
 
-            item.translate(tr);
+                if( best_score < global_score ) {
+                    auto d = getNfpPoint(optimum) - iv;
+                    d += startpos;
+                    final_tr = d;
+                    final_rot = initial_rot + rot;
+                    can_pack = true;
+                    global_score = best_score;
+                }
+            }
+
+            item.translation(final_tr);
+            item.rotation(final_rot);
         }
 
         if(can_pack) {
@@ -138,10 +397,13 @@ public:
         return ret;
     }
 
-private:
+    ~_NofitPolyPlacer() {
+        Nfp::Shapes<RawShape> m;
+        m.reserve(items_.size());
+
+        for(Item& item : items_) m.emplace_back(item.transformedShape());
+        auto&& bb = ShapeLike::boundingBox<RawShape>(m);
 
-    void setInitialPosition(Item& item) {
-        Box&& bb = item.boundingBox();
         Vertex ci, cb;
 
         switch(config_.alignment) {
@@ -173,11 +435,23 @@ private:
         }
 
         auto d = cb - ci;
+        for(Item& item : items_) item.translate(d);
+    }
+
+private:
+
+    void setInitialPosition(Item& item) {
+        Box&& bb = item.boundingBox();
+
+        Vertex ci = bb.minCorner();
+        Vertex cb = bin_.minCorner();
+
+        auto&& d = cb - ci;
         item.translate(d);
     }
 
     void placeOutsideOfBin(Item& item) {
-        auto bb = item.boundingBox();
+        auto&& bb = item.boundingBox();
         Box binbb = ShapeLike::boundingBox<RawShape>(bin_);
 
         Vertex v = { getX(bb.maxCorner()), getY(bb.minCorner()) };

xs/src/libnest2d/libnest2d/placers/placer_boilerplate.hpp

@@ -12,6 +12,8 @@ template<class Subclass, class RawShape, class TBin,
          class Store = std::vector<std::reference_wrapper<_Item<RawShape>>>
          >
 class PlacerBoilerplate {
+    mutable bool farea_valid_ = false;
+    mutable double farea_ = 0.0;
 public:
     using Item = _Item<RawShape>;
     using Vertex = TPoint<RawShape>;
@@ -39,7 +41,10 @@ public:
 
     using ItemGroup = const Container&;
 
-    inline PlacerBoilerplate(const BinType& bin): bin_(bin) {}
+    inline PlacerBoilerplate(const BinType& bin, unsigned cap = 50): bin_(bin)
+    {
+        items_.reserve(cap);
+    }
 
     inline const BinType& bin() const BP2D_NOEXCEPT { return bin_; }
 
@@ -53,7 +58,10 @@ public:
 
     bool pack(Item& item) {
         auto&& r = static_cast<Subclass*>(this)->trypack(item);
-        if(r) items_.push_back(*(r.item_ptr_));
+        if(r) {
+            items_.push_back(*(r.item_ptr_));
+            farea_valid_ = false;
+        }
         return r;
     }
 
@@ -62,14 +70,38 @@ public:
             r.item_ptr_->translation(r.move_);
             r.item_ptr_->rotation(r.rot_);
             items_.push_back(*(r.item_ptr_));
+            farea_valid_ = false;
         }
     }
 
-    void unpackLast() { items_.pop_back(); }
+    void unpackLast() {
+        items_.pop_back();
+        farea_valid_ = false;
+    }
 
     inline ItemGroup getItems() const { return items_; }
 
-    inline void clearItems() { items_.clear(); }
+    inline void clearItems() {
+        items_.clear();
+        farea_valid_ = false;
+#ifndef NDEBUG
+        debug_items_.clear();
+#endif
+    }
+
+    inline double filledArea() const {
+        if(farea_valid_) return farea_;
+        else {
+            farea_ = .0;
+            std::for_each(items_.begin(), items_.end(),
+                          [this] (Item& item) {
+                farea_ += item.area();
+            });
+            farea_valid_ = true;
+        }
+
+        return farea_;
+    }
 
 #ifndef NDEBUG
     std::vector<Item> debug_items_;

xs/src/libnest2d/libnest2d/selections/djd_heuristic.hpp

@@ -2,6 +2,10 @@
 #define DJD_HEURISTIC_HPP
 
 #include <list>
+#include <future>
+#include <atomic>
+#include <functional>
+
 #include "selection_boilerplate.hpp"
 
 namespace libnest2d { namespace strategies {
@@ -13,6 +17,20 @@ namespace libnest2d { namespace strategies {
 template<class RawShape>
 class _DJDHeuristic: public SelectionBoilerplate<RawShape> {
     using Base = SelectionBoilerplate<RawShape>;
+
+    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); }
+    };
+
 public:
     using typename Base::Item;
     using typename Base::ItemRef;
@@ -21,27 +39,54 @@ public:
      * @brief The Config for DJD heuristic.
      */
     struct Config {
-        /// Max number of bins.
-        unsigned max_bins = 0;
 
         /**
          * If true, the algorithm will try to place pair and driplets in all
          * possible order.
          */
         bool try_reverse_order = true;
+
+        /**
+         * The initial fill proportion of the bin area that will be filled before
+         * trying items one by one, or pairs or triplets.
+         *
+         * The initial fill proportion suggested by
+         * [López-Camacho]\
+         * (http://www.cs.stir.ac.uk/~goc/papers/EffectiveHueristic2DAOR2013.pdf)
+         * is one third of the area of bin.
+         */
+        double initial_fill_proportion = 1.0/3.0;
+
+        /**
+         * @brief How much is the acceptable waste incremented at each iteration
+         */
+        double waste_increment = 0.1;
+
+        /**
+         * @brief Allow parallel jobs for filling multiple bins.
+         *
+         * This will decrease the soution quality but can greatly boost up
+         * performance for large number of items.
+         */
+        bool allow_parallel = true;
+
+        /**
+         * @brief Always use parallel processing if the items don't fit into
+         * one bin.
+         */
+        bool force_parallel = false;
     };
 
 private:
     using Base::packed_bins_;
     using ItemGroup = typename Base::ItemGroup;
 
-    using Container = ItemGroup;//typename std::vector<Item>;
+    using Container = ItemGroup;
     Container store_;
     Config config_;
 
-    // The initial fill proportion of the bin area that will be filled before
+    static const unsigned MAX_ITEMS_SEQUENTIALLY = 30;
-    // trying items one by one, or pairs or triplets.
+    static const unsigned MAX_VERTICES_SEQUENTIALLY = MAX_ITEMS_SEQUENTIALLY*20;
-    static const double INITIAL_FILL_PROPORTION;
 
 public:
 
@@ -61,7 +106,9 @@ public:
         using ItemList = std::list<ItemRef>;
 
         const double bin_area = ShapeLike::area<RawShape>(bin);
-        const double w = bin_area * 0.1;
+        const double w = bin_area * config_.waste_increment;
+
+        const double INITIAL_FILL_PROPORTION = config_.initial_fill_proportion;
         const double INITIAL_FILL_AREA = bin_area*INITIAL_FILL_PROPORTION;
 
         store_.clear();
@@ -74,24 +121,22 @@ public:
             return i1.area() > i2.area();
         });
 
-        ItemList not_packed(store_.begin(), store_.end());
+        size_t glob_vertex_count = 0;
+        std::for_each(store_.begin(), store_.end(),
+                      [&glob_vertex_count](const Item& item) {
+             glob_vertex_count += item.vertexCount();
+        });
 
         std::vector<Placer> placers;
 
-        double free_area = 0;
-        double filled_area = 0;
-        double waste = 0;
         bool try_reverse = config_.try_reverse_order;
 
         // Will use a subroutine to add a new bin
-        auto addBin = [this, &placers, &free_area,
+        auto addBin = [this, &placers, &bin, &pconfig]()
-                       &filled_area, &bin, &pconfig]()
         {
             placers.emplace_back(bin);
             packed_bins_.emplace_back();
             placers.back().configure(pconfig);
-            free_area = ShapeLike::area<RawShape>(bin);
-            filled_area = 0;
         };
 
         // Types for pairs and triplets
@@ -136,8 +181,11 @@ public:
         };
 
         auto tryOneByOne = // Subroutine to try adding items one by one.
-                [&not_packed,  &bin_area, &free_area, &filled_area]
+                [&bin_area]
-                (Placer& placer, double waste)
+                (Placer& placer, ItemList& not_packed,
+                double waste,
+                double& free_area,
+                double& filled_area)
         {
             double item_area = 0;
             bool ret = false;
@@ -160,9 +208,12 @@ public:
         };
 
         auto tryGroupsOfTwo = // Try adding groups of two items into the bin.
-                [&not_packed, &bin_area, &free_area, &filled_area, &check_pair,
+                [&bin_area, &check_pair,
                  try_reverse]
-                (Placer& placer, double waste)
+                (Placer& placer, ItemList& not_packed,
+                double waste,
+                double& free_area,
+                double& filled_area)
         {
             double item_area = 0, largest_area = 0, smallest_area = 0;
             double second_largest = 0, second_smallest = 0;
@@ -259,12 +310,15 @@ public:
         };
 
         auto tryGroupsOfThree = // Try adding groups of three items.
-                [&not_packed, &bin_area, &free_area, &filled_area,
+                [&bin_area,
                  &check_pair, &check_triplet, try_reverse]
-                (Placer& placer, double waste)
+                (Placer& placer, ItemList& not_packed,
+                double waste,
+                double& free_area,
+                double& filled_area)
         {
-
+            auto np_size = not_packed.size();
-            if(not_packed.size() < 3) return false;
+            if(np_size < 3) return false;
 
             auto it = not_packed.begin();           // from
             const auto endit = not_packed.end();    // to
@@ -275,6 +329,10 @@ public:
             std::vector<TPair> wrong_pairs;
             std::vector<TTriplet> wrong_triplets;
 
+            auto cap = np_size*np_size / 2 ;
+            wrong_pairs.reserve(cap);
+            wrong_triplets.reserve(cap);
+
             // Will be true if a succesfull pack can be made.
             bool ret = false;
 
@@ -445,33 +503,68 @@ public:
             return ret;
         };
 
-        addBin();
-
         // Safety test: try to pack each item into an empty bin. If it fails
         // then it should be removed from the not_packed list
-        { auto it = not_packed.begin();
+        { auto it = store_.begin();
-            while (it != not_packed.end()) {
+            while (it != store_.end()) {
                 Placer p(bin);
                 if(!p.pack(*it)) {
                     auto itmp = it++;
-                    not_packed.erase(itmp);
+                    store_.erase(itmp);
                 } else it++;
             }
         }
 
-        auto makeProgress = [this, &not_packed](Placer& placer) {
+        int acounter = int(store_.size());
-            packed_bins_.back() = placer.getItems();
+        std::atomic_flag flg = ATOMIC_FLAG_INIT;
+        SpinLock slock(flg);
+
+        auto makeProgress = [this, &acounter, &slock]
+                (Placer& placer, size_t idx, int packednum)
+        {
+
+            packed_bins_[idx] = placer.getItems();
 #ifndef NDEBUG
-            packed_bins_.back().insert(packed_bins_.back().end(),
+            packed_bins_[idx].insert(packed_bins_[idx].end(),
                                        placer.getDebugItems().begin(),
                                        placer.getDebugItems().end());
 #endif
-            this->progress_(not_packed.size());
+            // TODO here should be a spinlock
+            slock.lock();
+            acounter -= packednum;
+            this->progress_(acounter);
+            slock.unlock();
         };
 
-        while(!not_packed.empty()) {
+        double items_area = 0;
+        for(Item& item : store_) items_area += item.area();
 
-            auto& placer = placers.back();
+        // Number of bins that will definitely be needed
+        auto bincount_guess = unsigned(std::ceil(items_area / bin_area));
+
+        // Do parallel if feasible
+        bool do_parallel = config_.allow_parallel && bincount_guess > 1 &&
+                ((glob_vertex_count >  MAX_VERTICES_SEQUENTIALLY ||
+                 store_.size() > MAX_ITEMS_SEQUENTIALLY) ||
+                config_.force_parallel);
+
+        if(do_parallel) dout() << "Parallel execution..." << "\n";
+
+        // The DJD heuristic algorithm itself:
+        auto packjob = [INITIAL_FILL_AREA, bin_area, w,
+                        &tryOneByOne,
+                        &tryGroupsOfTwo,
+                        &tryGroupsOfThree,
+                        &makeProgress]
+                        (Placer& placer, ItemList& not_packed, size_t idx)
+        {
+            bool can_pack = true;
+
+            double filled_area = placer.filledArea();
+            double free_area = bin_area - filled_area;
+            double waste = .0;
+
+            while(!not_packed.empty() && can_pack) {
 
                 {// Fill the bin up to INITIAL_FILL_PROPORTION of its capacity
                     auto it = not_packed.begin();
@@ -484,48 +577,97 @@ public:
                             free_area = bin_area - filled_area;
                             auto itmp = it++;
                             not_packed.erase(itmp);
-                        makeProgress(placer);
+                            makeProgress(placer, idx, 1);
                         } else it++;
                     }
                 }
 
                 // try pieses one by one
-            while(tryOneByOne(placer, waste)) {
+                while(tryOneByOne(placer, not_packed, waste, free_area,
+                                  filled_area)) {
                     waste = 0;
-                makeProgress(placer);
+                    makeProgress(placer, idx, 1);
                 }
 
                 // try groups of 2 pieses
-            while(tryGroupsOfTwo(placer, waste)) {
+                while(tryGroupsOfTwo(placer, not_packed, waste, free_area,
+                                     filled_area)) {
                     waste = 0;
-                makeProgress(placer);
+                    makeProgress(placer, idx, 2);
                 }
 
                 // try groups of 3 pieses
-            while(tryGroupsOfThree(placer, waste)) {
+                while(tryGroupsOfThree(placer, not_packed, waste, free_area,
+                                       filled_area)) {
                     waste = 0;
-                makeProgress(placer);
+                    makeProgress(placer, idx, 3);
                 }
 
                 if(waste < free_area) waste += w;
-            else if(!not_packed.empty()) addBin();
+                else if(!not_packed.empty()) can_pack = false;
             }
 
-//        std::for_each(placers.begin(), placers.end(),
+            return can_pack;
-//                      [this](Placer& placer){
-//            packed_bins_.push_back(placer.getItems());
-//        });
-    }
         };
 
-/*
+        size_t idx = 0;
- * The initial fill proportion suggested by
+        ItemList remaining;
- * [López-Camacho]\
+
- * (http://www.cs.stir.ac.uk/~goc/papers/EffectiveHueristic2DAOR2013.pdf)
+        if(do_parallel) {
- * is one third of the area of bin.
+            std::vector<ItemList> not_packeds(bincount_guess);
- */
+
-template<class RawShape>
+            // Preallocating the bins
-const double _DJDHeuristic<RawShape>::INITIAL_FILL_PROPORTION = 1.0/3.0;
+            for(unsigned b = 0; b < bincount_guess; b++) {
+                addBin();
+                ItemList& not_packed = not_packeds[b];
+                for(unsigned idx = b; idx < store_.size(); idx+=bincount_guess) {
+                    not_packed.push_back(store_[idx]);
+                }
+            }
+
+            // The parallel job
+            auto job = [&placers, &not_packeds, &packjob](unsigned idx) {
+                Placer& placer = placers[idx];
+                ItemList& not_packed = not_packeds[idx];
+                return packjob(placer, not_packed, idx);
+            };
+
+            // We will create jobs for each bin
+            std::vector<std::future<bool>> rets(bincount_guess);
+
+            for(unsigned b = 0; b < bincount_guess; b++) { // launch the jobs
+                rets[b] = std::async(std::launch::async, job, b);
+            }
+
+            for(unsigned fi = 0; fi < rets.size(); ++fi) {
+                rets[fi].wait();
+
+                // Collect remaining items while waiting for the running jobs
+                remaining.merge( not_packeds[fi], [](Item& i1, Item& i2) {
+                    return i1.area() > i2.area();
+                });
+
+            }
+
+            idx = placers.size();
+
+            // Try to put the remaining items into one of the packed bins
+            if(remaining.size() <= placers.size())
+            for(size_t j = 0; j < idx && !remaining.empty(); j++) {
+                packjob(placers[j], remaining, j);
+            }
+
+        } else {
+            remaining = ItemList(store_.begin(), store_.end());
+        }
+
+        while(!remaining.empty()) {
+            addBin();
+            packjob(placers[idx], remaining, idx); idx++;
+        }
+
+    }
+};
 
 }
 }

xs/src/libnest2d/libnest2d/selections/filler.hpp

@@ -33,7 +33,17 @@ public:
 
         store_.clear();
         store_.reserve(last-first);
-        packed_bins_.clear();
+        packed_bins_.emplace_back();
+
+        auto makeProgress = [this](PlacementStrategyLike<TPlacer>& placer) {
+            packed_bins_.back() = placer.getItems();
+#ifndef NDEBUG
+            packed_bins_.back().insert(packed_bins_.back().end(),
+                                       placer.getDebugItems().begin(),
+                                       placer.getDebugItems().end());
+#endif
+            this->progress_(packed_bins_.back().size());
+        };
 
         std::copy(first, last, std::back_inserter(store_));
 
@@ -50,18 +60,22 @@ public:
         PlacementStrategyLike<TPlacer> placer(bin);
         placer.configure(pconfig);
 
-        bool was_packed = false;
+        auto it = store_.begin();
-        for(auto& item : store_ ) {
+        while(it != store_.end()) {
-            if(!placer.pack(item))  {
+            if(!placer.pack(*it))  {
-                packed_bins_.push_back(placer.getItems());
+                if(packed_bins_.back().empty()) ++it;
+                makeProgress(placer);
                 placer.clearItems();
-                was_packed = placer.pack(item);
+                packed_bins_.emplace_back();
-            } else was_packed = true;
+            } else {
+                makeProgress(placer);
+                ++it;
+            }
         }
 
-        if(was_packed) {
+//        if(was_packed) {
-            packed_bins_.push_back(placer.getItems());
+//            packed_bins_.push_back(placer.getItems());
-        }
+//        }
     }
 };
 

xs/src/libnest2d/tests/main.cpp

@@ -1,6 +1,8 @@
 #include <iostream>
-#include <fstream>
 #include <string>
+#include <fstream>
+
+//#define DEBUG_EXPORT_NFP
 
 #include <libnest2d.h>
 #include <libnest2d/geometries_io.hpp>
@@ -8,6 +10,8 @@
 #include "printer_parts.h"
 #include "benchmark.h"
 #include "svgtools.hpp"
+//#include <libnest2d/optimizer.hpp>
+//#include <libnest2d/optimizers/simplex.hpp>
 
 using namespace libnest2d;
 using ItemGroup = std::vector<std::reference_wrapper<Item>>;
@@ -36,51 +40,122 @@ std::vector<Item>& stegoParts() {
 void arrangeRectangles() {
     using namespace libnest2d;
 
-    auto input = stegoParts();
-
     const int SCALE = 1000000;
+    std::vector<Rectangle> rects = {
+        {80*SCALE, 80*SCALE},
+        {60*SCALE, 90*SCALE},
+        {70*SCALE, 30*SCALE},
+        {80*SCALE, 60*SCALE},
+        {60*SCALE, 60*SCALE},
+        {60*SCALE, 40*SCALE},
+        {40*SCALE, 40*SCALE},
+        {10*SCALE, 10*SCALE},
+        {10*SCALE, 10*SCALE},
+        {10*SCALE, 10*SCALE},
+        {10*SCALE, 10*SCALE},
+        {10*SCALE, 10*SCALE},
+        {5*SCALE, 5*SCALE},
+        {5*SCALE, 5*SCALE},
+        {5*SCALE, 5*SCALE},
+        {5*SCALE, 5*SCALE},
+        {5*SCALE, 5*SCALE},
+        {5*SCALE, 5*SCALE},
+        {5*SCALE, 5*SCALE},
+        {20*SCALE, 20*SCALE}
+       };
 
-    Box bin(210*SCALE, 250*SCALE);
+//    std::vector<Rectangle> rects = {
+//        {20*SCALE, 10*SCALE},
+//        {20*SCALE, 10*SCALE},
+//        {20*SCALE, 20*SCALE},
+//    };
 
-    Coord min_obj_distance = 0; //6*SCALE;
+//    std::vector<Item> input {
+//        {{0, 0}, {0, 20*SCALE}, {10*SCALE, 0}, {0, 0}}
+//    };
 
-    NfpPlacer::Config pconf;
+    std::vector<Item> input;
-    pconf.alignment = NfpPlacer::Config::Alignment::TOP_LEFT;
+    input.insert(input.end(), prusaParts().begin(), prusaParts().end());
-    Arranger<NfpPlacer, DJDHeuristic> arrange(bin, min_obj_distance, pconf);
+//    input.insert(input.end(), stegoParts().begin(), stegoParts().end());
+//    input.insert(input.end(), rects.begin(), rects.end());
 
-//    arrange.progressIndicator([&arrange, &bin](unsigned r){
+    Box bin(250*SCALE, 210*SCALE);
+
+    Coord min_obj_distance = 6*SCALE;
+
+    using Packer = Arranger<NfpPlacer, DJDHeuristic>;
+
+    Packer::PlacementConfig pconf;
+    pconf.alignment = NfpPlacer::Config::Alignment::CENTER;
+//    pconf.rotations = {0.0, Pi/2.0, Pi, 3*Pi/2};
+    Packer::SelectionConfig sconf;
+    sconf.allow_parallel = true;
+    sconf.force_parallel = false;
+    sconf.try_reverse_order = false;
+    Packer arrange(bin, min_obj_distance, pconf, sconf);
+
+    arrange.progressIndicator([&](unsigned r){
 //        svg::SVGWriter::Config conf;
 //        conf.mm_in_coord_units = SCALE;
 //        svg::SVGWriter svgw(conf);
 //        svgw.setSize(bin);
 //        svgw.writePackGroup(arrange.lastResult());
-//        svgw.save("out");
+//        svgw.save("debout");
-//        std::cout << "Remaining items: " << r << std::endl;
+        std::cout << "Remaining items: " << r << std::endl;
-//    });
+    }).useMinimumBoundigBoxRotation();
 
     Benchmark bench;
 
     bench.start();
-    auto result = arrange(input.begin(),
+    auto result = arrange.arrange(input.begin(),
                           input.end());
 
     bench.stop();
 
-    std::cout << bench.getElapsedSec() << std::endl;
+    std::vector<double> eff;
+    eff.reserve(result.size());
+
+    auto bin_area = double(bin.height()*bin.width());
+    for(auto& r : result) {
+        double a = 0;
+        std::for_each(r.begin(), r.end(), [&a] (Item& e ){ a += e.area(); });
+        eff.emplace_back(a/bin_area);
+    };
+
+    std::cout << bench.getElapsedSec() << " bin count: " << result.size()
+              << std::endl;
+
+    std::cout << "Bin efficiency: (";
+    for(double e : eff) std::cout << e*100.0 << "% ";
+    std::cout << ") Average: "
+              << std::accumulate(eff.begin(), eff.end(), 0.0)*100.0/result.size()
+              << " %" << std::endl;
+
+    std::cout << "Bin usage: (";
+    unsigned total = 0;
+    for(auto& r : result) { std::cout << r.size() << " "; total += r.size(); }
+    std::cout << ") Total: " << total << std::endl;
 
     for(auto& it : input) {
         auto ret = ShapeLike::isValid(it.transformedShape());
         std::cout << ret.second << std::endl;
     }
 
+    if(total != input.size()) std::cout << "ERROR " << "could not pack "
+                                        << input.size() - total << " elements!"
+                                        << std::endl;
+
     svg::SVGWriter::Config conf;
     conf.mm_in_coord_units = SCALE;
     svg::SVGWriter svgw(conf);
     svgw.setSize(bin);
     svgw.writePackGroup(result);
+//    std::for_each(input.begin(), input.end(), [&svgw](Item& item){ svgw.writeItem(item);});
     svgw.save("out");
 }
 
+
+
 int main(void /*int argc, char **argv*/) {
     arrangeRectangles();
 //    findDegenerateCase();

xs/src/libnest2d/tests/svgtools.hpp

@@ -50,7 +50,9 @@ public:
         if(conf_.origo_location == BOTTOMLEFT)
         for(unsigned i = 0; i < tsh.vertexCount(); i++) {
             auto v = tsh.vertex(i);
-            setY(v, -getY(v) + conf_.height*conf_.mm_in_coord_units);
+            auto d = static_cast<Coord>(
+                        std::round(conf_.height*conf_.mm_in_coord_units) );
+            setY(v, -getY(v) + d);
             tsh.setVertex(i, v);
         }
         currentLayer() += ShapeLike::serialize<Formats::SVG>(tsh.rawShape(),
@@ -78,8 +80,8 @@ public:
     }
 
     void save(const std::string& filepath) {
-        unsigned lyrc = svg_layers_.size() > 1? 1 : 0;
+        size_t lyrc = svg_layers_.size() > 1? 1 : 0;
-        unsigned last = svg_layers_.size() > 1? svg_layers_.size() : 0;
+        size_t last = svg_layers_.size() > 1? svg_layers_.size() : 0;
 
         for(auto& lyr : svg_layers_) {
             std::fstream out(filepath + (lyrc > 0? std::to_string(lyrc) : "") +

xs/src/libnest2d/tests/test.cpp

@@ -253,7 +253,7 @@ TEST(GeometryAlgorithms, LeftAndDownPolygon)
     ASSERT_TRUE(ShapeLike::isValid(leftp.rawShape()).first);
     ASSERT_EQ(leftp.vertexCount(), leftControl.vertexCount());
 
-    for(size_t i = 0; i < leftControl.vertexCount(); i++) {
+    for(unsigned long i = 0; i < leftControl.vertexCount(); i++) {
         ASSERT_EQ(getX(leftp.vertex(i)), getX(leftControl.vertex(i)));
         ASSERT_EQ(getY(leftp.vertex(i)), getY(leftControl.vertex(i)));
     }
@@ -263,7 +263,7 @@ TEST(GeometryAlgorithms, LeftAndDownPolygon)
     ASSERT_TRUE(ShapeLike::isValid(downp.rawShape()).first);
     ASSERT_EQ(downp.vertexCount(), downControl.vertexCount());
 
-    for(size_t i = 0; i < downControl.vertexCount(); i++) {
+    for(unsigned long i = 0; i < downControl.vertexCount(); i++) {
         ASSERT_EQ(getX(downp.vertex(i)), getX(downControl.vertex(i)));
         ASSERT_EQ(getY(downp.vertex(i)), getY(downControl.vertex(i)));
     }
@@ -696,6 +696,27 @@ TEST(GeometryAlgorithms, nfpConvexConvex) {
     }
 }
 
+TEST(GeometryAlgorithms, pointOnPolygonContour) {
+    using namespace libnest2d;
+
+    Rectangle input(10, 10);
+
+    strategies::EdgeCache<PolygonImpl> ecache(input);
+
+    auto first = *input.begin();
+    ASSERT_TRUE(getX(first) == getX(ecache.coords(0)));
+    ASSERT_TRUE(getY(first) == getY(ecache.coords(0)));
+
+    auto last = *std::prev(input.end());
+    ASSERT_TRUE(getX(last) == getX(ecache.coords(1.0)));
+    ASSERT_TRUE(getY(last) == getY(ecache.coords(1.0)));
+
+    for(int i = 0; i <= 100; i++) {
+        auto v = ecache.coords(i*(0.01));
+        ASSERT_TRUE(ShapeLike::touches(v, input.transformedShape()));
+    }
+}
+
 int main(int argc, char **argv) {
   ::testing::InitGoogleTest(&argc, argv);
   return RUN_ALL_TESTS();

xs/src/libslic3r/Model.cpp

@@ -331,14 +331,20 @@ ShapeData2D projectModelFromTop(const Slic3r::Model &model) {
             for(auto objinst : objptr->instances) {
                 if(objinst) {
                     Slic3r::TriangleMesh tmpmesh = rmesh;
-                    objinst->transform_mesh(&tmpmesh);
+//                    objinst->transform_mesh(&tmpmesh);
                     ClipperLib::PolyNode pn;
                     auto p = tmpmesh.convex_hull();
                     p.make_clockwise();
                     p.append(p.first_point());
                     pn.Contour = Slic3rMultiPoint_to_ClipperPath( p );
 
-                    ret.emplace_back(objinst, Item(std::move(pn)));
+                    Item item(std::move(pn));
+                    item.rotation(objinst->rotation);
+                    item.translation( {
+                        ClipperLib::cInt(objinst->offset.x/SCALING_FACTOR),
+                        ClipperLib::cInt(objinst->offset.y/SCALING_FACTOR)
+                    });
+                    ret.emplace_back(objinst, item);
                 }
             }
         }
@@ -407,7 +413,7 @@ bool arrange(Model &model, coordf_t dist, const Slic3r::BoundingBoxf* bb,
 //    std::cout << "}" << std::endl;
 //    return true;
 
-    double area = 0;
+    bool hasbin = bb != nullptr && bb->defined;
     double area_max = 0;
     Item *biggest = nullptr;
 
@@ -416,9 +422,10 @@ bool arrange(Model &model, coordf_t dist, const Slic3r::BoundingBoxf* bb,
     std::vector<std::reference_wrapper<Item>> shapes;
     shapes.reserve(shapemap.size());
     std::for_each(shapemap.begin(), shapemap.end(),
-                  [&shapes, &area, min_obj_distance, &area_max, &biggest]
+                  [&shapes, min_obj_distance, &area_max, &biggest,hasbin]
                   (ShapeData2D::value_type& it)
     {
+        if(!hasbin) {
             Item& item = it.second;
             item.addOffset(min_obj_distance);
             auto b = ShapeLike::boundingBox(item.transformedShape());
@@ -427,13 +434,13 @@ bool arrange(Model &model, coordf_t dist, const Slic3r::BoundingBoxf* bb,
                 area_max = static_cast<double>(a);
                 biggest = &item;
             }
-        area += b.width()*b.height();
+        }
         shapes.push_back(std::ref(it.second));
     });
 
     Box bin;
 
-    if(bb != nullptr && bb->defined) {
+    if(hasbin) {
         // Scale up the bounding box to clipper scale.
         BoundingBoxf bbb = *bb;
         bbb.scale(1.0/SCALING_FACTOR);
@@ -456,8 +463,13 @@ bool arrange(Model &model, coordf_t dist, const Slic3r::BoundingBoxf* bb,
     using Arranger = Arranger<NfpPlacer, DJDHeuristic>;
 
     Arranger::PlacementConfig pcfg;
-    pcfg.alignment = Arranger::PlacementConfig::Alignment::BOTTOM_LEFT;
+    Arranger::SelectionConfig scfg;
-    Arranger arranger(bin, min_obj_distance, pcfg);
+
+    scfg.try_reverse_order = false;
+    scfg.force_parallel = true;
+    pcfg.alignment = Arranger::PlacementConfig::Alignment::CENTER;
+    Arranger arranger(bin, min_obj_distance, pcfg, scfg);
+    arranger.useMinimumBoundigBoxRotation();
 
     std::cout << "Arranging model..." << std::endl;
     bench.start();
@@ -483,18 +495,15 @@ bool arrange(Model &model, coordf_t dist, const Slic3r::BoundingBoxf* bb,
 
             // Get the tranformation data from the item object and scale it
             // appropriately
-            Radians rot = item.rotation();
             auto off = item.translation();
-            Pointf foff(off.X*SCALING_FACTOR + batch_offset,
+            Radians rot = item.rotation();
+            Pointf foff(off.X*SCALING_FACTOR,
                         off.Y*SCALING_FACTOR);
 
             // write the tranformation data into the model instance
-            inst_ptr->rotation += rot;
+            inst_ptr->rotation = rot;
-            inst_ptr->offset += foff;
+            inst_ptr->offset = foff;
-
+            inst_ptr->offset_z = -batch_offset;
-            // Debug
-            /*std::cout << "item " << idx << ": \n" << "\toffset_x: "
-             * << foff.x << "\n\toffset_y: " << foff.y << std::endl;*/
         }
     };
 
@@ -503,14 +512,23 @@ bool arrange(Model &model, coordf_t dist, const Slic3r::BoundingBoxf* bb,
     if(first_bin_only) {
         applyResult(result.front(), 0);
     } else {
+
+        const auto STRIDE_PADDING = 1.2;
+        const auto MIN_STRIDE = 100;
+
+        auto h = STRIDE_PADDING * model.bounding_box().size().z;
+        h = h < MIN_STRIDE ? MIN_STRIDE : h;
+        Coord stride = static_cast<Coord>(h);
+
         Coord batch_offset = 0;
+
         for(auto& group : result) {
             applyResult(group, batch_offset);
 
             // Only the first pack group can be placed onto the print bed. The
             // other objects which could not fit will be placed next to the
             // print bed
-            batch_offset += static_cast<Coord>(2*bin.width()*SCALING_FACTOR);
+            batch_offset += stride;
         }
     }
     bench.stop();
@@ -1236,7 +1254,7 @@ void ModelInstance::transform_mesh(TriangleMesh* mesh, bool dont_translate) cons
     mesh->rotate_z(this->rotation);                 // rotate around mesh origin
     mesh->scale(this->scaling_factor);              // scale around mesh origin
     if (!dont_translate)
-        mesh->translate(this->offset.x, this->offset.y, 0);
+        mesh->translate(this->offset.x, this->offset.y, this->offset_z);
 }
 
 BoundingBoxf3 ModelInstance::transform_mesh_bounding_box(const TriangleMesh* mesh, bool dont_translate) const

xs/src/libslic3r/Model.hpp

@@ -201,8 +201,9 @@ public:
     double rotation;            // Rotation around the Z axis, in radians around mesh center point
     double scaling_factor;
     Pointf offset;              // in unscaled coordinates
+    double offset_z = 0;
     
-    ModelObject* get_object() const { return this->object; };
+    ModelObject* get_object() const { return this->object; }
 
     // To be called on an external mesh
     void transform_mesh(TriangleMesh* mesh, bool dont_translate = false) const;

xs/src/libslic3r/Print.cpp

@@ -444,7 +444,7 @@ bool Print::apply_config(DynamicPrintConfig config)
                         const ModelVolume &volume = *object->model_object()->volumes[volume_id];
                         if (this_region_config_set) {
                             // If the new config for this volume differs from the other
-                            // volume configs currently associated to this region, it means
+                            // volume configs currently associated to this region, it means
                             // the region subdivision does not make sense anymore.
                             if (! this_region_config.equals(this->_region_config_from_model_volume(volume))) {
                                 rearrange_regions = true;