Fixing nesting crash in debug mode.

Also commented out unnecessary bloat from AppController

src/libnest2d/include/libnest2d.h

@@ -0,0 +1,175 @@
+#ifndef LIBNEST2D_H
+#define LIBNEST2D_H
+
+// The type of backend should be set conditionally by the cmake configuriation
+// for now we set it statically to clipper backend
+#ifdef LIBNEST2D_BACKEND_CLIPPER
+#include <libnest2d/backends/clipper/geometries.hpp>
+#endif
+
+#ifdef LIBNEST2D_OPTIMIZER_NLOPT
+// We include the stock optimizers for local and global optimization
+#include <libnest2d/optimizers/nlopt/subplex.hpp>     // Local subplex for NfpPlacer
+#include <libnest2d/optimizers/nlopt/genetic.hpp>     // Genetic for min. bounding box
+#endif
+
+#include <libnest2d/libnest2d.hpp>
+#include <libnest2d/placers/bottomleftplacer.hpp>
+#include <libnest2d/placers/nfpplacer.hpp>
+#include <libnest2d/selections/firstfit.hpp>
+#include <libnest2d/selections/filler.hpp>
+#include <libnest2d/selections/djd_heuristic.hpp>
+
+namespace libnest2d {
+
+using Point = PointImpl;
+using Coord = TCoord<PointImpl>;
+using Box = _Box<PointImpl>;
+using Segment = _Segment<PointImpl>;
+using Circle = _Circle<PointImpl>;
+
+using Item = _Item<PolygonImpl>;
+using Rectangle = _Rectangle<PolygonImpl>;
+
+using PackGroup = _PackGroup<PolygonImpl>;
+using IndexedPackGroup = _IndexedPackGroup<PolygonImpl>;
+
+using FillerSelection = selections::_FillerSelection<PolygonImpl>;
+using FirstFitSelection = selections::_FirstFitSelection<PolygonImpl>;
+using DJDHeuristic  = selections::_DJDHeuristic<PolygonImpl>;
+
+template<class Bin> // Generic placer for arbitrary bin types
+using _NfpPlacer = placers::_NofitPolyPlacer<PolygonImpl, Bin>;
+
+// NfpPlacer is with Box bin
+using NfpPlacer = _NfpPlacer<Box>;
+
+// This supports only box shaped bins
+using BottomLeftPlacer = placers::_BottomLeftPlacer<PolygonImpl>;
+
+template<class Placer = NfpPlacer,
+         class Selector = FirstFitSelection,
+         class Iterator = std::vector<Item>::iterator>
+PackGroup nest(Iterator from, Iterator to,
+               const typename Placer::BinType& bin,
+               Coord dist = 0,
+               const typename Placer::Config& pconf = {},
+               const typename Selector::Config& sconf = {})
+{
+    Nester<Placer, Selector> nester(bin, dist, pconf, sconf);
+    return nester.execute(from, to);
+}
+
+template<class Placer = NfpPlacer,
+         class Selector = FirstFitSelection,
+         class Container = std::vector<Item>>
+PackGroup nest(Container&& cont,
+               const typename Placer::BinType& bin,
+               Coord dist = 0,
+               const typename Placer::Config& pconf = {},
+               const typename Selector::Config& sconf = {})
+{
+    return nest<Placer, Selector>(cont.begin(), cont.end(),
+                                  bin, dist, pconf, sconf);
+}
+
+template<class Placer = NfpPlacer,
+         class Selector = FirstFitSelection,
+         class Iterator = std::vector<Item>::iterator>
+PackGroup nest(Iterator from, Iterator to,
+               const typename Placer::BinType& bin,
+               ProgressFunction prg,
+               StopCondition scond = []() { return false; },
+               Coord dist = 0,
+               const typename Placer::Config& pconf = {},
+               const typename Selector::Config& sconf = {})
+{
+    Nester<Placer, Selector> nester(bin, dist, pconf, sconf);
+    if(prg) nester.progressIndicator(prg);
+    if(scond) nester.stopCondition(scond);
+    return nester.execute(from, to);
+}
+
+template<class Placer = NfpPlacer,
+         class Selector = FirstFitSelection,
+         class Container = std::vector<Item>>
+PackGroup nest(Container&& cont,
+               const typename Placer::BinType& bin,
+               ProgressFunction prg,
+               StopCondition scond = []() { return false; },
+               Coord dist = 0,
+               const typename Placer::Config& pconf = {},
+               const typename Selector::Config& sconf = {})
+{
+    return nest<Placer, Selector>(cont.begin(), cont.end(),
+                                  bin, prg, scond, dist, pconf, sconf);
+}
+
+#ifdef LIBNEST2D_STATIC
+extern template
+PackGroup nest<NfpPlacer, FirstFitSelection, std::vector<Item>&>(
+    std::vector<Item>& cont,
+    const Box& bin,
+    Coord dist,
+    const NfpPlacer::Config& pcfg,
+    const FirstFitSelection::Config& scfg
+);
+
+extern template
+PackGroup nest<NfpPlacer, FirstFitSelection, std::vector<Item>&>(
+    std::vector<Item>& cont,
+    const Box& bin,
+    ProgressFunction prg,
+    StopCondition scond,
+    Coord dist,
+    const NfpPlacer::Config& pcfg,
+    const FirstFitSelection::Config& scfg
+);
+
+extern template
+PackGroup nest<NfpPlacer, FirstFitSelection, std::vector<Item>>(
+    std::vector<Item>&& cont,
+    const Box& bin,
+    Coord dist,
+    const NfpPlacer::Config& pcfg,
+    const FirstFitSelection::Config& scfg
+);
+
+extern template
+PackGroup nest<NfpPlacer, FirstFitSelection, std::vector<Item>>(
+    std::vector<Item>&& cont,
+    const Box& bin,
+    ProgressFunction prg,
+    StopCondition scond,
+    Coord dist,
+    const NfpPlacer::Config& pcfg,
+    const FirstFitSelection::Config& scfg
+);
+
+extern template
+PackGroup nest<NfpPlacer, FirstFitSelection, std::vector<Item>::iterator>(
+    std::vector<Item>::iterator from,
+    std::vector<Item>::iterator to,
+    const Box& bin,
+    Coord dist,
+    const NfpPlacer::Config& pcfg,
+    const FirstFitSelection::Config& scfg
+);
+
+extern template
+PackGroup nest<NfpPlacer, FirstFitSelection, std::vector<Item>::iterator>(
+    std::vector<Item>::iterator from,
+    std::vector<Item>::iterator to,
+    const Box& bin,
+    ProgressFunction prg,
+    StopCondition scond,
+    Coord dist,
+    const NfpPlacer::Config& pcfg,
+    const FirstFitSelection::Config& scfg
+);
+
+#endif
+
+}
+
+#endif // LIBNEST2D_H

src/libnest2d/include/libnest2d/backends/clipper/CMakeLists.txt

@@ -0,0 +1,73 @@
+if(NOT TARGET clipper) # If there is a clipper target in the parent project we are good to go.
+
+    find_package(Clipper 6.1)
+
+    if(NOT CLIPPER_FOUND)
+        find_package(Subversion QUIET)
+        if(Subversion_FOUND)        
+
+            set(URL_CLIPPER "https://svn.code.sf.net/p/polyclipping/code/trunk/cpp" 
+            CACHE STRING "Clipper source code repository location.")
+
+            message(STATUS "Clipper not found so it will be downloaded.")
+            # Silently download and build the library in the build dir
+
+            if (CMAKE_VERSION VERSION_LESS 3.2)
+                set(UPDATE_DISCONNECTED_IF_AVAILABLE "")
+            else()
+                set(UPDATE_DISCONNECTED_IF_AVAILABLE "UPDATE_DISCONNECTED 1")
+            endif()
+
+            include(DownloadProject)
+            download_project(   PROJ                clipper_library
+                                SVN_REPOSITORY      ${URL_CLIPPER}
+                                SVN_REVISION        -r540
+                                #SOURCE_SUBDIR       cpp
+                                INSTALL_COMMAND     ""
+                                CONFIGURE_COMMAND   "" # Not working, I will just add the source files
+                                ${UPDATE_DISCONNECTED_IF_AVAILABLE}
+            )
+
+            # This is not working and I dont have time to fix it
+            # add_subdirectory(${clipper_library_SOURCE_DIR}/cpp
+            #                  ${clipper_library_BINARY_DIR}
+            # )
+
+            add_library(ClipperBackend STATIC
+                ${clipper_library_SOURCE_DIR}/clipper.cpp
+                ${clipper_library_SOURCE_DIR}/clipper.hpp)
+
+            target_include_directories(ClipperBackend INTERFACE 
+                ${clipper_library_SOURCE_DIR})
+        else()
+            message(FATAL_ERROR "Can't find clipper library and no SVN client found to download.
+                You can download the clipper sources and define a clipper target in your project, that will work for libnest2d.")
+        endif()
+    else()
+        add_library(ClipperBackend INTERFACE)
+        target_link_libraries(ClipperBackend INTERFACE Clipper::Clipper)
+    endif()
+else()
+    # set(CLIPPER_INCLUDE_DIRS "" PARENT_SCOPE)
+    # set(CLIPPER_LIBRARIES clipper PARENT_SCOPE)
+    add_library(ClipperBackend INTERFACE)
+    target_link_libraries(ClipperBackend INTERFACE clipper)
+endif()
+
+# Clipper backend is not enough on its own, it still needs some functions
+# from Boost geometry
+if(NOT Boost_INCLUDE_DIRS_FOUND)
+    find_package(Boost 1.58 REQUIRED)
+    # TODO automatic download of boost geometry headers
+endif()
+
+target_include_directories(ClipperBackend INTERFACE ${Boost_INCLUDE_DIRS} )
+target_sources(ClipperBackend INTERFACE
+    ${CMAKE_CURRENT_SOURCE_DIR}/geometries.hpp
+    ${SRC_DIR}/libnest2d/utils/boost_alg.hpp )
+
+target_compile_definitions(ClipperBackend INTERFACE LIBNEST2D_BACKEND_CLIPPER)
+
+# And finally plug the ClipperBackend into libnest2d
+target_link_libraries(libnest2d INTERFACE ClipperBackend)
+

src/libnest2d/include/libnest2d/backends/clipper/geometries.hpp

@@ -0,0 +1,460 @@
+#ifndef CLIPPER_BACKEND_HPP
+#define CLIPPER_BACKEND_HPP
+
+#include <sstream>
+#include <unordered_map>
+#include <cassert>
+#include <vector>
+#include <iostream>
+
+#include <libnest2d/geometry_traits.hpp>
+#include <libnest2d/geometry_traits_nfp.hpp>
+
+#include <clipper.hpp>
+
+namespace ClipperLib {
+using PointImpl = IntPoint;
+using PathImpl = Path;
+using HoleStore = std::vector<PathImpl>;
+
+struct PolygonImpl {
+    PathImpl Contour;
+    HoleStore Holes;
+
+    inline PolygonImpl() = default;
+
+    inline explicit PolygonImpl(const PathImpl& cont): Contour(cont) {}
+    inline explicit PolygonImpl(const HoleStore& holes):
+        Holes(holes) {}
+    inline PolygonImpl(const Path& cont, const HoleStore& holes):
+        Contour(cont), Holes(holes) {}
+
+    inline explicit PolygonImpl(PathImpl&& cont): Contour(std::move(cont)) {}
+    inline explicit PolygonImpl(HoleStore&& holes): Holes(std::move(holes)) {}
+    inline PolygonImpl(Path&& cont, HoleStore&& holes):
+        Contour(std::move(cont)), Holes(std::move(holes)) {}
+};
+
+inline PointImpl& operator +=(PointImpl& p, const PointImpl& pa ) {
+    // This could be done with SIMD
+    p.X += pa.X;
+    p.Y += pa.Y;
+    return p;
+}
+
+inline PointImpl operator+(const PointImpl& p1, const PointImpl& p2) {
+    PointImpl ret = p1;
+    ret += p2;
+    return ret;
+}
+
+inline PointImpl& operator -=(PointImpl& p, const PointImpl& pa ) {
+    p.X -= pa.X;
+    p.Y -= pa.Y;
+    return p;
+}
+
+inline PointImpl operator -(PointImpl& p ) {
+    PointImpl ret = p;
+    ret.X = -ret.X;
+    ret.Y = -ret.Y;
+    return ret;
+}
+
+inline PointImpl operator-(const PointImpl& p1, const PointImpl& p2) {
+    PointImpl ret = p1;
+    ret -= p2;
+    return ret;
+}
+
+inline PointImpl& operator *=(PointImpl& p, const PointImpl& pa ) {
+    p.X *= pa.X;
+    p.Y *= pa.Y;
+    return p;
+}
+
+inline PointImpl operator*(const PointImpl& p1, const PointImpl& p2) {
+    PointImpl ret = p1;
+    ret *= p2;
+    return ret;
+}
+
+}
+
+namespace libnest2d {
+
+// Aliases for convinience
+using ClipperLib::PointImpl;
+using ClipperLib::PathImpl;
+using ClipperLib::PolygonImpl;
+using ClipperLib::HoleStore;
+
+// Type of coordinate units used by Clipper
+template<> struct CoordType<PointImpl> {
+    using Type = ClipperLib::cInt;
+};
+
+// Type of point used by Clipper
+template<> struct PointType<PolygonImpl> {
+    using Type = PointImpl;
+};
+
+template<> struct PointType<PathImpl> {
+    using Type = PointImpl;
+};
+
+template<> struct PointType<PointImpl> {
+    using Type = PointImpl;
+};
+
+template<> struct CountourType<PolygonImpl> {
+    using Type = PathImpl;
+};
+
+template<> struct ShapeTag<PolygonImpl> { using Type = PolygonTag; };
+template<> struct ShapeTag<PathImpl> { using Type = PathTag; };
+
+template<> struct ShapeTag<TMultiShape<PolygonImpl>> {
+    using Type = MultiPolygonTag;
+};
+
+template<> struct PointType<TMultiShape<PolygonImpl>> {
+    using Type = PointImpl;
+};
+
+template<> struct HolesContainer<PolygonImpl> {
+    using Type = ClipperLib::Paths;
+};
+
+namespace pointlike {
+
+// Tell libnest2d how to extract the X coord from a ClipperPoint object
+template<> inline TCoord<PointImpl> x(const PointImpl& p)
+{
+    return p.X;
+}
+
+// Tell libnest2d how to extract the Y coord from a ClipperPoint object
+template<> inline TCoord<PointImpl> y(const PointImpl& p)
+{
+    return p.Y;
+}
+
+// Tell libnest2d how to extract the X coord from a ClipperPoint object
+template<> inline TCoord<PointImpl>& x(PointImpl& p)
+{
+    return p.X;
+}
+
+// Tell libnest2d how to extract the Y coord from a ClipperPoint object
+template<> inline TCoord<PointImpl>& y(PointImpl& p)
+{
+    return p.Y;
+}
+
+}
+
+#define DISABLE_BOOST_AREA
+
+namespace _smartarea {
+template<Orientation o>
+inline double area(const PolygonImpl& /*sh*/) {
+    return std::nan("");
+}
+
+template<>
+inline double area<Orientation::CLOCKWISE>(const PolygonImpl& sh) {
+    double a = 0;
+
+    std::for_each(sh.Holes.begin(), sh.Holes.end(), [&a](const PathImpl& h)
+    {
+        a -= ClipperLib::Area(h);
+    });
+
+    return -ClipperLib::Area(sh.Contour) + a;
+}
+
+template<>
+inline double area<Orientation::COUNTER_CLOCKWISE>(const PolygonImpl& sh) {
+    double a = 0;
+
+    std::for_each(sh.Holes.begin(), sh.Holes.end(), [&a](const PathImpl& h)
+    {
+        a += ClipperLib::Area(h);
+    });
+
+    return ClipperLib::Area(sh.Contour) + a;
+}
+
+}
+
+namespace shapelike {
+
+// Tell libnest2d how to make string out of a ClipperPolygon object
+template<> inline double area(const PolygonImpl& sh, const PolygonTag&)
+{
+    return _smartarea::area<OrientationType<PolygonImpl>::Value>(sh);
+}
+
+template<> inline void offset(PolygonImpl& sh, TCoord<PointImpl> distance)
+{
+    #define DISABLE_BOOST_OFFSET
+
+    using ClipperLib::ClipperOffset;
+    using ClipperLib::jtMiter;
+    using ClipperLib::etClosedPolygon;
+    using ClipperLib::Paths;
+
+    // If the input is not at least a triangle, we can not do this algorithm
+    if(sh.Contour.size() <= 3 ||
+       std::any_of(sh.Holes.begin(), sh.Holes.end(),
+                   [](const PathImpl& p) { return p.size() <= 3; })
+       ) throw GeometryException(GeomErr::OFFSET);
+
+    ClipperOffset offs;
+    Paths result;
+    offs.AddPath(sh.Contour, jtMiter, etClosedPolygon);
+    offs.AddPaths(sh.Holes, jtMiter, etClosedPolygon);
+    offs.Execute(result, static_cast<double>(distance));
+
+    // Offsetting reverts the orientation and also removes the last vertex
+    // so boost will not have a closed polygon.
+
+    bool found_the_contour = false;
+    for(auto& r : result) {
+        if(ClipperLib::Orientation(r)) {
+            // We don't like if the offsetting generates more than one contour
+            // but throwing would be an overkill. Instead, we should warn the
+            // caller about the inability to create correct geometries
+            if(!found_the_contour) {
+                sh.Contour = r;
+                ClipperLib::ReversePath(sh.Contour);
+                sh.Contour.push_back(sh.Contour.front());
+                found_the_contour = true;
+            } else {
+                dout() << "Warning: offsetting result is invalid!";
+                /* TODO warning */
+            }
+        } else {
+            // TODO If there are multiple contours we can't be sure which hole
+            // belongs to the first contour. (But in this case the situation is
+            // bad enough to let it go...)
+            sh.Holes.push_back(r);
+            ClipperLib::ReversePath(sh.Holes.back());
+            sh.Holes.back().push_back(sh.Holes.back().front());
+        }
+    }
+}
+
+// Tell libnest2d how to make string out of a ClipperPolygon object
+template<> inline std::string toString(const PolygonImpl& sh)
+{
+    std::stringstream ss;
+
+    ss << "Contour {\n";
+    for(auto p : sh.Contour) {
+        ss << "\t" << p.X << " " << p.Y << "\n";
+    }
+    ss << "}\n";
+
+    for(auto& h : sh.Holes) {
+        ss << "Holes {\n";
+        for(auto p : h)  {
+            ss << "\t{\n";
+            ss << "\t\t" << p.X << " " << p.Y << "\n";
+            ss << "\t}\n";
+        }
+        ss << "}\n";
+    }
+
+    return ss.str();
+}
+
+template<>
+inline PolygonImpl create(const PathImpl& path, const HoleStore& holes)
+{
+    PolygonImpl p;
+    p.Contour = path;
+
+    // Expecting that the coordinate system Y axis is positive in upwards
+    // direction
+    if(ClipperLib::Orientation(p.Contour)) {
+        // Not clockwise then reverse the b*tch
+        ClipperLib::ReversePath(p.Contour);
+    }
+
+    p.Holes = holes;
+    for(auto& h : p.Holes) {
+        if(!ClipperLib::Orientation(h)) {
+            ClipperLib::ReversePath(h);
+        }
+    }
+
+    return p;
+}
+
+template<> inline PolygonImpl create( PathImpl&& path, HoleStore&& holes) {
+    PolygonImpl p;
+    p.Contour.swap(path);
+
+    // Expecting that the coordinate system Y axis is positive in upwards
+    // direction
+    if(ClipperLib::Orientation(p.Contour)) {
+        // Not clockwise then reverse the b*tch
+        ClipperLib::ReversePath(p.Contour);
+    }
+
+    p.Holes.swap(holes);
+
+    for(auto& h : p.Holes) {
+        if(!ClipperLib::Orientation(h)) {
+            ClipperLib::ReversePath(h);
+        }
+    }
+
+    return p;
+}
+
+template<>
+inline const THolesContainer<PolygonImpl>& holes(const PolygonImpl& sh)
+{
+    return sh.Holes;
+}
+
+template<> inline THolesContainer<PolygonImpl>& holes(PolygonImpl& sh)
+{
+    return sh.Holes;
+}
+
+template<>
+inline TContour<PolygonImpl>& hole(PolygonImpl& sh, unsigned long idx)
+{
+    return sh.Holes[idx];
+}
+
+template<>
+inline const TContour<PolygonImpl>& hole(const PolygonImpl& sh,
+                                            unsigned long idx)
+{
+    return sh.Holes[idx];
+}
+
+template<> inline size_t holeCount(const PolygonImpl& sh)
+{
+    return sh.Holes.size();
+}
+
+template<> inline PathImpl& contour(PolygonImpl& sh)
+{
+    return sh.Contour;
+}
+
+template<>
+inline const PathImpl& contour(const PolygonImpl& sh)
+{
+    return sh.Contour;
+}
+
+#define DISABLE_BOOST_TRANSLATE
+template<>
+inline void translate(PolygonImpl& sh, const PointImpl& offs)
+{
+    for(auto& p : sh.Contour) { p += offs; }
+    for(auto& hole : sh.Holes) for(auto& p : hole) { p += offs; }
+}
+
+#define DISABLE_BOOST_ROTATE
+template<>
+inline void rotate(PolygonImpl& sh, const Radians& rads)
+{
+    using Coord = TCoord<PointImpl>;
+
+    auto cosa = rads.cos();
+    auto sina = rads.sin();
+
+    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.Holes) for(auto& p : hole) {
+        p = {
+                static_cast<Coord>(p.X * cosa - p.Y * sina),
+                static_cast<Coord>(p.X * sina + p.Y * cosa)
+            };
+    }
+}
+
+} // namespace shapelike
+
+#define DISABLE_BOOST_NFP_MERGE
+inline std::vector<PolygonImpl> _merge(ClipperLib::Clipper& clipper) {
+    shapelike::Shapes<PolygonImpl> retv;
+
+    ClipperLib::PolyTree result;
+    clipper.Execute(ClipperLib::ctUnion, result, ClipperLib::pftNegative);
+    retv.reserve(static_cast<size_t>(result.Total()));
+
+    std::function<void(ClipperLib::PolyNode*, PolygonImpl&)> processHole;
+
+    auto processPoly = [&retv, &processHole](ClipperLib::PolyNode *pptr) {
+        PolygonImpl poly(pptr->Contour);
+        poly.Contour.push_back(poly.Contour.front());
+        for(auto h : pptr->Childs) { processHole(h, poly); }
+        retv.push_back(poly);
+    };
+
+    processHole = [&processPoly](ClipperLib::PolyNode *pptr, PolygonImpl& poly)
+    {
+        poly.Holes.push_back(pptr->Contour);
+        poly.Holes.back().push_back(poly.Holes.back().front());
+        for(auto c : pptr->Childs) processPoly(c);
+    };
+
+    auto traverse = [&processPoly] (ClipperLib::PolyNode *node)
+    {
+        for(auto ch : node->Childs) {
+            processPoly(ch);
+        }
+    };
+
+    traverse(&result);
+
+    return retv;
+}
+
+namespace nfp {
+
+template<> inline std::vector<PolygonImpl>
+merge(const std::vector<PolygonImpl>& shapes)
+{
+    ClipperLib::Clipper clipper(ClipperLib::ioReverseSolution);
+
+    bool closed = true;
+    bool valid = true;
+
+    for(auto& path : shapes) {
+        valid &= clipper.AddPath(path.Contour, ClipperLib::ptSubject, closed);
+
+        for(auto& hole : path.Holes) {
+            valid &= clipper.AddPath(hole, ClipperLib::ptSubject, closed);
+        }
+    }
+
+    if(!valid) throw GeometryException(GeomErr::MERGE);
+
+    return _merge(clipper);
+}
+
+}
+
+}
+
+//#define DISABLE_BOOST_SERIALIZE
+//#define DISABLE_BOOST_UNSERIALIZE
+
+// All other operators and algorithms are implemented with boost
+#include <libnest2d/utils/boost_alg.hpp>
+
+#endif // CLIPPER_BACKEND_HPP

src/libnest2d/include/libnest2d/common.hpp

@@ -0,0 +1,202 @@
+#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
+    #define BP2D_COMPILER_MSVC12
+#elif __cplusplus >= 201103L
+    #define BP2D_NOEXCEPT noexcept
+    #define BP2D_CONSTEXPR constexpr
+#endif
+
+/*
+ * 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<
+        typename std::remove_reference<T>::type>::type;
+};
+
+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;
+
+/**
+ * A useful little tool for triggering static_assert error messages e.g. when
+ * a mandatory template specialization (implementation) is missing.
+ *
+ * \tparam T A template argument that may come from and outer template method.
+ */
+template<class T> struct always_false { enum { value = false }; };
+
+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{}) { }
+  Double(double d) : val_(d) { }
+
+  operator double() const BP2D_NOEXCEPT { return val_; }
+  operator double&() BP2D_NOEXCEPT { return val_; }
+};
+
+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_;
+    }
+};
+
+/**
+ * @brief Data type representing degrees. It supports conversion to radians.
+ */
+class Degrees: public Double {
+public:
+    Degrees(double deg = Double()): Double(deg) {}
+    Degrees(const Radians& rads): Double( rads * 180/Pi ) {}
+    inline double toRadians() { return Radians(*this);}
+};
+
+inline bool operator==(const Degrees& deg, const Radians& rads) {
+    Degrees deg2 = rads;
+    auto diff = std::abs(deg - deg2);
+    return diff < 0.0001;
+}
+
+inline bool operator==(const Radians& rads, const Degrees& deg) {
+    return deg == rads;
+}
+
+inline Radians::operator Degrees() { return *this * 180/Pi; }
+
+inline Radians::Radians(const Degrees &degs): Double( degs * Pi/180) {}
+
+inline double Radians::toDegrees() { return operator Degrees(); }
+
+enum class GeomErr : std::size_t {
+    OFFSET,
+    MERGE,
+    NFP
+};
+
+const std::string ERROR_STR[] = {
+    "Offsetting could not be done! An invalid geometry may have been added.",
+    "Error while merging geometries!",
+    "No fit polygon cannot be calculated."
+};
+
+class GeometryException: public std::exception {
+
+    virtual const std::string& errorstr(GeomErr errcode) const BP2D_NOEXCEPT {
+        return ERROR_STR[static_cast<std::size_t>(errcode)];
+    }
+
+    GeomErr errcode_;
+public:
+
+    GeometryException(GeomErr code): errcode_(code) {}
+
+    GeomErr errcode() const { return errcode_; }
+
+    const char * what() const BP2D_NOEXCEPT override {
+        return errorstr(errcode_).c_str();
+    }
+};
+
+
+}
+#endif // LIBNEST2D_CONFIG_HPP

src/libnest2d/include/libnest2d/geometry_traits.hpp

@@ -0,0 +1,965 @@
+#ifndef GEOMETRY_TRAITS_HPP
+#define GEOMETRY_TRAITS_HPP
+
+#include <string>
+#include <type_traits>
+#include <algorithm>
+#include <array>
+#include <vector>
+#include <numeric>
+#include <limits>
+#include <iterator>
+#include <cmath>
+
+#include "common.hpp"
+
+namespace libnest2d {
+
+/// Getting the coordinate data type for a geometry class.
+template<class GeomClass> struct CoordType { using Type = long; };
+
+/// TCoord<GeomType> as shorthand for typename `CoordType<GeomType>::Type`.
+template<class GeomType>
+using TCoord = typename CoordType<remove_cvref_t<GeomType>>::Type;
+
+
+/// Getting the type of point structure used by a shape.
+template<class Sh> struct PointType { using Type = typename Sh::PointType; };
+
+/// TPoint<ShapeClass> as shorthand for `typename PointType<ShapeClass>::Type`.
+template<class Shape>
+using TPoint = typename PointType<remove_cvref_t<Shape>>::Type;
+
+
+template<class RawShape> struct CountourType { using Type = RawShape; };
+
+template<class RawShape>
+using TContour = typename CountourType<remove_cvref_t<RawShape>>::Type;
+
+
+template<class RawShape>
+struct HolesContainer { using Type = std::vector<TContour<RawShape>>;  };
+
+template<class RawShape>
+using THolesContainer = typename HolesContainer<remove_cvref_t<RawShape>>::Type;
+
+
+template<class RawShape>
+struct LastPointIsFirst { static const bool Value = true; };
+
+enum class Orientation {
+    CLOCKWISE,
+    COUNTER_CLOCKWISE
+};
+
+template<class RawShape>
+struct OrientationType {
+
+    // Default Polygon orientation that the library expects
+    static const Orientation Value = Orientation::CLOCKWISE;
+};
+
+/**
+ * \brief A point pair base class for other point pairs (segment, box, ...).
+ * \tparam RawPoint The actual point type to use.
+ */
+template<class RawPoint>
+struct PointPair {
+    RawPoint p1;
+    RawPoint p2;
+};
+
+struct PolygonTag {};
+struct PathTag {};
+struct MultiPolygonTag {};
+struct BoxTag {};
+struct CircleTag {};
+
+template<class Shape> struct ShapeTag { using Type = typename Shape::Tag; };
+template<class S> using Tag = typename ShapeTag<remove_cvref_t<S>>::Type;
+
+template<class S> struct MultiShape { using Type = std::vector<S>; };
+template<class S>
+using TMultiShape =typename MultiShape<remove_cvref_t<S>>::Type;
+
+/**
+ * \brief An abstraction of a box;
+ */
+template<class RawPoint>
+class _Box: PointPair<RawPoint> {
+    using PointPair<RawPoint>::p1;
+    using PointPair<RawPoint>::p2;
+public:
+
+    using Tag = BoxTag;
+    using PointType = RawPoint;
+
+    inline _Box() = default;
+    inline _Box(const RawPoint& p, const RawPoint& pp):
+        PointPair<RawPoint>({p, pp}) {}
+
+    inline _Box(TCoord<RawPoint> width, TCoord<RawPoint> height):
+        _Box(RawPoint{0, 0}, RawPoint{width, height}) {}
+
+    inline const RawPoint& minCorner() const BP2D_NOEXCEPT { return p1; }
+    inline const RawPoint& maxCorner() const BP2D_NOEXCEPT { return p2; }
+
+    inline RawPoint& minCorner() BP2D_NOEXCEPT { return p1; }
+    inline RawPoint& maxCorner() BP2D_NOEXCEPT { return p2; }
+
+    inline TCoord<RawPoint> width() const BP2D_NOEXCEPT;
+    inline TCoord<RawPoint> height() const BP2D_NOEXCEPT;
+
+    inline RawPoint center() const BP2D_NOEXCEPT;
+
+    inline double area() const BP2D_NOEXCEPT {
+        return double(width()*height());
+    }
+};
+
+template<class RawPoint>
+class _Circle {
+    RawPoint center_;
+    double radius_ = 0;
+public:
+
+    using Tag = CircleTag;
+    using PointType = RawPoint;
+
+    _Circle() = default;
+
+    _Circle(const RawPoint& center, double r): center_(center), radius_(r) {}
+
+    inline const RawPoint& center() const BP2D_NOEXCEPT { return center_; }
+    inline const void center(const RawPoint& c) { center_ = c; }
+
+    inline double radius() const BP2D_NOEXCEPT { return radius_; }
+    inline void radius(double r) { radius_ = r; }
+
+    inline double area() const BP2D_NOEXCEPT {
+        return 2.0*Pi*radius_*radius_;
+    }
+};
+
+/**
+ * \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:
+
+    using PointType = RawPoint;
+
+    inline _Segment() = default;
+
+    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 void first(const RawPoint& p) BP2D_NOEXCEPT
+    {
+        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 almost as a namespace. The only difference is that is can
+// used in friend declarations.
+namespace pointlike {
+
+template<class RawPoint>
+inline TCoord<RawPoint> x(const RawPoint& p)
+{
+    return p.x();
+}
+
+template<class RawPoint>
+inline TCoord<RawPoint> y(const RawPoint& p)
+{
+    return p.y();
+}
+
+template<class RawPoint>
+inline TCoord<RawPoint>& x(RawPoint& p)
+{
+    return p.x();
+}
+
+template<class RawPoint>
+inline TCoord<RawPoint>& y(RawPoint& p)
+{
+    return p.y();
+}
+
+template<class RawPoint>
+inline double distance(const RawPoint& /*p1*/, const RawPoint& /*p2*/)
+{
+    static_assert(always_false<RawPoint>::value,
+                  "PointLike::distance(point, point) unimplemented!");
+    return 0;
+}
+
+template<class RawPoint>
+inline double distance(const RawPoint& /*p1*/,
+                       const _Segment<RawPoint>& /*s*/)
+{
+    static_assert(always_false<RawPoint>::value,
+                  "PointLike::distance(point, segment) unimplemented!");
+    return 0;
+}
+
+template<class RawPoint>
+inline std::pair<TCoord<RawPoint>, bool> horizontalDistance(
+        const RawPoint& p, const _Segment<RawPoint>& s)
+{
+    using Unit = TCoord<RawPoint>;
+    auto x = pointlike::x(p), y = pointlike::y(p);
+    auto x1 = pointlike::x(s.first()), y1 = pointlike::y(s.first());
+    auto x2 = pointlike::x(s.second()), y2 = pointlike::y(s.second());
+
+    TCoord<RawPoint> ret;
+
+    if( (y < y1 && y < y2) || (y > y1 && y > y2) )
+        return {0, false};
+    if ((y == y1 && y == y2) && (x > x1 && x > x2))
+        ret = std::min( x-x1, x -x2);
+    else if( (y == y1 && y == y2) && (x < x1 && x < x2))
+        ret = -std::min(x1 - x, x2 - x);
+    else if(std::abs(y - y1) <= std::numeric_limits<Unit>::epsilon() &&
+            std::abs(y - y2) <= std::numeric_limits<Unit>::epsilon())
+        ret = 0;
+    else
+        ret = x - x1 + (x1 - x2)*(y1 - y)/(y1 - y2);
+
+    return {ret, true};
+}
+
+template<class RawPoint>
+inline std::pair<TCoord<RawPoint>, bool> verticalDistance(
+        const RawPoint& p, const _Segment<RawPoint>& s)
+{
+    using Unit = TCoord<RawPoint>;
+    auto x = pointlike::x(p), y = pointlike::y(p);
+    auto x1 = pointlike::x(s.first()), y1 = pointlike::y(s.first());
+    auto x2 = pointlike::x(s.second()), y2 = pointlike::y(s.second());
+
+    TCoord<RawPoint> ret;
+
+    if( (x < x1 && x < x2) || (x > x1 && x > x2) )
+        return {0, false};
+    if ((x == x1 && x == x2) && (y > y1 && y > y2))
+        ret = std::min( y-y1, y -y2);
+    else if( (x == x1 && x == x2) && (y < y1 && y < y2))
+        ret = -std::min(y1 - y, y2 - y);
+    else if(std::abs(x - x1) <= std::numeric_limits<Unit>::epsilon() &&
+            std::abs(x - x2) <= std::numeric_limits<Unit>::epsilon())
+        ret = 0;
+    else
+        ret = y - y1 + (y1 - y2)*(x1 - x)/(x1 - x2);
+
+    return {ret, true};
+}
+}
+
+template<class RawPoint>
+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
+{
+    return pointlike::y(maxCorner()) - pointlike::y(minCorner());
+}
+
+template<class RawPoint>
+TCoord<RawPoint> getX(const RawPoint& p) { return pointlike::x<RawPoint>(p); }
+
+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)
+{
+    pointlike::x<RawPoint>(p) = val;
+}
+
+template<class RawPoint>
+void setY(RawPoint& p, const TCoord<RawPoint>& val)
+{
+    pointlike::y<RawPoint>(p) = val;
+}
+
+template<class RawPoint>
+inline Radians _Segment<RawPoint>::angleToXaxis() const
+{
+    if(std::isnan(static_cast<double>(angletox_))) {
+        TCoord<RawPoint> dx = getX(second()) - getX(first());
+        TCoord<RawPoint> dy = getY(second()) - getY(first());
+
+        double a = std::atan2(dy, dx);
+        auto s = std::signbit(a);
+
+        if(s) a += Pi_2;
+        angletox_ = a;
+    }
+    return angletox_;
+}
+
+template<class RawPoint>
+inline double _Segment<RawPoint>::length()
+{
+    return pointlike::distance(first(), second());
+}
+
+template<class RawPoint>
+inline RawPoint _Box<RawPoint>::center() const BP2D_NOEXCEPT {
+    auto& minc = minCorner();
+    auto& maxc = maxCorner();
+
+    using Coord = TCoord<RawPoint>;
+
+    RawPoint ret =  { // No rounding here, we dont know if these are int coords
+        static_cast<Coord>( (getX(minc) + getX(maxc))/2.0 ),
+        static_cast<Coord>( (getY(minc) + getY(maxc))/2.0 )
+    };
+
+    return ret;
+}
+
+enum class Formats {
+    WKT,
+    SVG
+};
+
+// This struct serves as a namespace. The only difference is that it can be
+// used in friend declarations and can be aliased at class scope.
+namespace shapelike {
+
+template<class RawShape>
+using Shapes = TMultiShape<RawShape>;
+
+template<class RawShape>
+inline RawShape create(const TContour<RawShape>& contour,
+                       const THolesContainer<RawShape>& holes)
+{
+    return RawShape(contour, holes);
+}
+
+template<class RawShape>
+inline RawShape create(TContour<RawShape>&& contour,
+                       THolesContainer<RawShape>&& holes)
+{
+    return RawShape(contour, holes);
+}
+
+template<class RawShape>
+inline RawShape create(const TContour<RawShape>& contour)
+{
+    return create<RawShape>(contour, {});
+}
+
+template<class RawShape>
+inline RawShape create(TContour<RawShape>&& contour)
+{
+    return create<RawShape>(contour, {});
+}
+
+template<class RawShape>
+inline THolesContainer<RawShape>& holes(RawShape& /*sh*/)
+{
+    static THolesContainer<RawShape> empty;
+    return empty;
+}
+
+template<class RawShape>
+inline const THolesContainer<RawShape>& holes(const RawShape& /*sh*/)
+{
+    static THolesContainer<RawShape> empty;
+    return empty;
+}
+
+template<class RawShape>
+inline TContour<RawShape>& hole(RawShape& sh, unsigned long idx)
+{
+    return holes(sh)[idx];
+}
+
+template<class RawShape>
+inline const TContour<RawShape>& hole(const RawShape& sh, unsigned long idx)
+{
+    return holes(sh)[idx];
+}
+
+template<class RawShape>
+inline size_t holeCount(const RawShape& sh)
+{
+    return holes(sh).size();
+}
+
+template<class RawShape>
+inline TContour<RawShape>& contour(RawShape& sh)
+{
+    static_assert(always_false<RawShape>::value,
+                  "shapelike::contour() unimplemented!");
+    return sh;
+}
+
+template<class RawShape>
+inline const TContour<RawShape>& contour(const RawShape& sh)
+{
+    static_assert(always_false<RawShape>::value,
+                  "shapelike::contour() unimplemented!");
+    return sh;
+}
+
+// Optional, does nothing by default
+template<class RawPath>
+inline void reserve(RawPath& p, size_t vertex_capacity, const PathTag&)
+{
+    p.reserve(vertex_capacity);
+}
+
+template<class RawShape, class...Args>
+inline void addVertex(RawShape& sh, const PathTag&, Args...args)
+{
+    return sh.emplace_back(std::forward<Args>(args)...);
+}
+
+template<class RawShape, class Fn>
+inline void foreachVertex(RawShape& sh, Fn fn, const PathTag&) {
+    std::for_each(sh.begin(), sh.end(), fn);
+}
+
+template<class RawShape>
+inline typename RawShape::iterator begin(RawShape& sh, const PathTag&)
+{
+    return sh.begin();
+}
+
+template<class RawShape>
+inline typename RawShape::iterator end(RawShape& sh, const PathTag&)
+{
+    return sh.end();
+}
+
+template<class RawShape>
+inline typename RawShape::const_iterator
+cbegin(const RawShape& sh, const PathTag&)
+{
+    return sh.cbegin();
+}
+
+template<class RawShape>
+inline typename RawShape::const_iterator
+cend(const RawShape& sh, const PathTag&)
+{
+    return sh.cend();
+}
+
+template<class RawShape>
+inline std::string toString(const RawShape& /*sh*/)
+{
+    return "";
+}
+
+template<Formats, class RawShape>
+inline std::string serialize(const RawShape& /*sh*/, double /*scale*/=1)
+{
+    static_assert(always_false<RawShape>::value,
+                  "shapelike::serialize() unimplemented!");
+    return "";
+}
+
+template<Formats, class RawShape>
+inline void unserialize(RawShape& /*sh*/, const std::string& /*str*/)
+{
+    static_assert(always_false<RawShape>::value,
+                  "shapelike::unserialize() unimplemented!");
+}
+
+template<class RawShape>
+inline double area(const RawShape& /*sh*/, const PolygonTag&)
+{
+    static_assert(always_false<RawShape>::value,
+                  "shapelike::area() unimplemented!");
+    return 0;
+}
+
+template<class RawShape>
+inline bool intersects(const RawShape& /*sh*/, const RawShape& /*sh*/)
+{
+    static_assert(always_false<RawShape>::value,
+                  "shapelike::intersects() unimplemented!");
+    return false;
+}
+
+template<class RawShape>
+inline bool isInside(const TPoint<RawShape>& /*point*/,
+                     const RawShape& /*shape*/)
+{
+    static_assert(always_false<RawShape>::value,
+                  "shapelike::isInside(point, shape) unimplemented!");
+    return false;
+}
+
+template<class RawShape>
+inline bool isInside(const RawShape& /*shape*/,
+                     const RawShape& /*shape*/)
+{
+    static_assert(always_false<RawShape>::value,
+                  "shapelike::isInside(shape, shape) unimplemented!");
+    return false;
+}
+
+template<class RawShape>
+inline bool touches( const RawShape& /*shape*/,
+                     const RawShape& /*shape*/)
+{
+    static_assert(always_false<RawShape>::value,
+                  "shapelike::touches(shape, shape) unimplemented!");
+    return false;
+}
+
+template<class RawShape>
+inline 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>
+inline _Box<TPoint<RawShape>> boundingBox(const RawShape& /*sh*/,
+                                          const PolygonTag&)
+{
+    static_assert(always_false<RawShape>::value,
+                  "shapelike::boundingBox(shape) unimplemented!");
+}
+
+template<class RawShapes>
+inline _Box<TPoint<typename RawShapes::value_type>>
+boundingBox(const RawShapes& /*sh*/, const MultiPolygonTag&)
+{
+    static_assert(always_false<RawShapes>::value,
+                  "shapelike::boundingBox(shapes) unimplemented!");
+}
+
+template<class RawShape>
+inline RawShape convexHull(const RawShape& /*sh*/, const PolygonTag&)
+{
+    static_assert(always_false<RawShape>::value,
+                  "shapelike::convexHull(shape) unimplemented!");
+    return RawShape();
+}
+
+template<class RawShapes>
+inline typename RawShapes::value_type
+convexHull(const RawShapes& /*sh*/, const MultiPolygonTag&)
+{
+    static_assert(always_false<RawShapes>::value,
+                  "shapelike::convexHull(shapes) unimplemented!");
+    return typename RawShapes::value_type();
+}
+
+template<class RawShape>
+inline void rotate(RawShape& /*sh*/, const Radians& /*rads*/)
+{
+    static_assert(always_false<RawShape>::value,
+                  "shapelike::rotate() unimplemented!");
+}
+
+template<class RawShape, class RawPoint>
+inline void translate(RawShape& /*sh*/, const RawPoint& /*offs*/)
+{
+    static_assert(always_false<RawShape>::value,
+                  "shapelike::translate() unimplemented!");
+}
+
+template<class RawShape>
+inline void offset(RawShape& /*sh*/, TCoord<TPoint<RawShape>> /*distance*/)
+{
+    dout() << "The current geometry backend does not support offsetting!\n";
+}
+
+template<class RawShape>
+inline std::pair<bool, std::string> isValid(const RawShape& /*sh*/)
+{
+    return {false, "shapelike::isValid() unimplemented!"};
+}
+
+template<class RawPath> inline bool isConvex(const RawPath& sh, const PathTag&)
+{
+    using Vertex = TPoint<RawPath>;
+    auto first = begin(sh);
+    auto middle = std::next(first);
+    auto last = std::next(middle);
+    using CVrRef = const Vertex&;
+
+    auto zcrossproduct = [](CVrRef k, CVrRef k1, CVrRef k2) {
+        auto dx1 = getX(k1) - getX(k);
+        auto dy1 = getY(k1) - getY(k);
+        auto dx2 = getX(k2) - getX(k1);
+        auto dy2 = getY(k2) - getY(k1);
+        return dx1*dy2 - dy1*dx2;
+    };
+
+    auto firstprod = zcrossproduct( *(std::prev(std::prev(end(sh)))),
+                                    *first,
+                                    *middle );
+
+    bool ret = true;
+    bool frsign = firstprod > 0;
+    while(last != end(sh)) {
+        auto &k = *first, &k1 = *middle, &k2 = *last;
+        auto zc = zcrossproduct(k, k1, k2);
+        ret &= frsign == (zc > 0);
+        ++first; ++middle; ++last;
+    }
+
+    return ret;
+}
+
+// *****************************************************************************
+// No need to implement these
+// *****************************************************************************
+
+template<class RawShape>
+inline typename TContour<RawShape>::iterator
+begin(RawShape& sh, const PolygonTag& t)
+{
+    return begin(contour(sh), PathTag());
+}
+
+template<class RawShape> // Tag dispatcher
+inline auto begin(RawShape& sh) -> decltype(begin(sh, Tag<RawShape>()))
+{
+    return begin(sh, Tag<RawShape>());
+}
+
+template<class RawShape>
+inline typename TContour<RawShape>::const_iterator
+cbegin(const RawShape& sh, const PolygonTag&)
+{
+    return cbegin(contour(sh), PathTag());
+}
+
+template<class RawShape> // Tag dispatcher
+inline auto cbegin(const RawShape& sh) -> decltype(cbegin(sh, Tag<RawShape>()))
+{
+    return cbegin(sh, Tag<RawShape>());
+}
+
+template<class RawShape>
+inline typename TContour<RawShape>::iterator
+end(RawShape& sh, const PolygonTag&)
+{
+    return end(contour(sh), PathTag());
+}
+
+template<class RawShape> // Tag dispatcher
+inline auto end(RawShape& sh) -> decltype(begin(sh, Tag<RawShape>()))
+{
+    return end(sh, Tag<RawShape>());
+}
+
+template<class RawShape>
+inline typename TContour<RawShape>::const_iterator
+cend(const RawShape& sh, const PolygonTag&)
+{
+    return cend(contour(sh), PathTag());
+}
+
+template<class RawShape> // Tag dispatcher
+inline auto cend(const RawShape& sh) -> decltype(cend(sh, Tag<RawShape>()))
+{
+    return cend(sh, Tag<RawShape>());
+}
+
+template<class It> std::reverse_iterator<It> _backward(It iter) {
+    return std::reverse_iterator<It>(iter);
+}
+
+template<class P> auto rbegin(P& p) -> decltype(_backward(end(p)))
+{
+    return _backward(end(p));
+}
+
+template<class P> auto rcbegin(const P& p) -> decltype(_backward(end(p)))
+{
+    return _backward(end(p));
+}
+
+template<class P> auto rend(P& p) -> decltype(_backward(begin(p)))
+{
+    return _backward(begin(p));
+}
+
+template<class P> auto rcend(const P& p) -> decltype(_backward(cbegin(p)))
+{
+    return _backward(cbegin(p));
+}
+
+template<class P> TPoint<P> front(const P& p) { return *shapelike::cbegin(p); }
+template<class P> TPoint<P> back (const P& p) {
+    return *backward(shapelike::cend(p));
+}
+
+// Optional, does nothing by default
+template<class RawShape>
+inline void reserve(RawShape& sh, size_t vertex_capacity, const PolygonTag&)
+{
+    reserve(contour(sh), vertex_capacity, PathTag());
+}
+
+template<class T> // Tag dispatcher
+inline void reserve(T& sh, size_t vertex_capacity) {
+    reserve(sh, vertex_capacity, Tag<T>());
+}
+
+template<class RawShape, class...Args>
+inline void addVertex(RawShape& sh, const PolygonTag&, Args...args)
+{
+    return addVertex(contour(sh), PathTag(), std::forward<Args>(args)...);
+}
+
+template<class RawShape, class...Args> // Tag dispatcher
+inline void addVertex(RawShape& sh, Args...args)
+{
+    return addVertex(sh, Tag<RawShape>(), std::forward<Args>(args)...);
+}
+
+template<class Box>
+inline Box boundingBox(const Box& box, const BoxTag& )
+{
+    return box;
+}
+
+template<class Circle>
+inline _Box<typename Circle::PointType> boundingBox(
+        const Circle& circ, const CircleTag&)
+{
+    using Point = typename Circle::PointType;
+    using Coord = TCoord<Point>;
+    Point pmin = {
+        static_cast<Coord>(getX(circ.center()) - circ.radius()),
+        static_cast<Coord>(getY(circ.center()) - circ.radius()) };
+
+    Point pmax = {
+        static_cast<Coord>(getX(circ.center()) + circ.radius()),
+        static_cast<Coord>(getY(circ.center()) + circ.radius()) };
+
+    return {pmin, pmax};
+}
+
+template<class S> // Dispatch function
+inline _Box<TPoint<S>> boundingBox(const S& sh)
+{
+    return boundingBox(sh, Tag<S>() );
+}
+
+template<class Box>
+inline double area(const Box& box, const BoxTag& )
+{
+    return box.area();
+}
+
+template<class Circle>
+inline double area(const Circle& circ, const CircleTag& )
+{
+    return circ.area();
+}
+
+template<class RawShape> // Dispatching function
+inline double area(const RawShape& sh)
+{
+    return area(sh, Tag<RawShape>());
+}
+
+template<class RawShapes>
+inline double area(const RawShapes& shapes, const MultiPolygonTag&)
+{
+    using RawShape = typename RawShapes::value_type;
+    return std::accumulate(shapes.begin(), shapes.end(), 0.0,
+                    [](double a, const RawShape& b) {
+        return a += area(b);
+    });
+}
+
+template<class RawShape>
+inline auto convexHull(const RawShape& sh)
+    -> decltype(convexHull(sh, Tag<RawShape>())) // TODO: C++14 could deduce
+{
+    return convexHull(sh, Tag<RawShape>());
+}
+
+template<class RawShape>
+inline bool isInside(const TPoint<RawShape>& point,
+                     const _Circle<TPoint<RawShape>>& circ)
+{
+    return pointlike::distance(point, circ.center()) < circ.radius();
+}
+
+template<class RawShape>
+inline bool isInside(const TPoint<RawShape>& point,
+                     const _Box<TPoint<RawShape>>& box)
+{
+    auto px = getX(point);
+    auto py = getY(point);
+    auto minx = getX(box.minCorner());
+    auto miny = getY(box.minCorner());
+    auto maxx = getX(box.maxCorner());
+    auto maxy = getY(box.maxCorner());
+
+    return px > minx && px < maxx && py > miny && py < maxy;
+}
+
+template<class RawShape>
+inline bool isInside(const RawShape& sh,
+                     const _Circle<TPoint<RawShape>>& circ)
+{
+    return std::all_of(cbegin(sh), cend(sh),
+                       [&circ](const TPoint<RawShape>& p){
+        return isInside<RawShape>(p, circ);
+    });
+}
+
+template<class RawShape>
+inline bool isInside(const _Box<TPoint<RawShape>>& box,
+                     const _Circle<TPoint<RawShape>>& circ)
+{
+    return isInside<RawShape>(box.minCorner(), circ) &&
+            isInside<RawShape>(box.maxCorner(), circ);
+}
+
+template<class RawShape>
+inline bool isInside(const _Box<TPoint<RawShape>>& ibb,
+                     const _Box<TPoint<RawShape>>& box)
+{
+    auto iminX = getX(ibb.minCorner());
+    auto imaxX = getX(ibb.maxCorner());
+    auto iminY = getY(ibb.minCorner());
+    auto imaxY = getY(ibb.maxCorner());
+
+    auto minX = getX(box.minCorner());
+    auto maxX = getX(box.maxCorner());
+    auto minY = getY(box.minCorner());
+    auto maxY = getY(box.maxCorner());
+
+    return iminX > minX && imaxX < maxX && iminY > minY && imaxY < maxY;
+}
+
+template<class RawShape> // Potential O(1) implementation may exist
+inline TPoint<RawShape>& vertex(RawShape& sh, unsigned long idx,
+                                const PolygonTag&)
+{
+    return *(shapelike::begin(contour(sh)) + idx);
+}
+
+template<class RawShape> // Potential O(1) implementation may exist
+inline TPoint<RawShape>& vertex(RawShape& sh, unsigned long idx,
+                                const PathTag&)
+{
+    return *(shapelike::begin(sh) + idx);
+}
+
+template<class RawShape> // Potential O(1) implementation may exist
+inline TPoint<RawShape>& vertex(RawShape& sh, unsigned long idx)
+{
+    return vertex(sh, idx, Tag<RawShape>());
+}
+
+template<class RawShape> // Potential O(1) implementation may exist
+inline const TPoint<RawShape>& vertex(const RawShape& sh,
+                                      unsigned long idx,
+                                      const PolygonTag&)
+{
+    return *(shapelike::cbegin(contour(sh)) + idx);
+}
+
+template<class RawShape> // Potential O(1) implementation may exist
+inline const TPoint<RawShape>& vertex(const RawShape& sh,
+                                      unsigned long idx,
+                                      const PathTag&)
+{
+    return *(shapelike::cbegin(sh) + idx);
+}
+
+
+template<class RawShape> // Potential O(1) implementation may exist
+inline const TPoint<RawShape>& vertex(const RawShape& sh,
+                                      unsigned long idx)
+{
+    return vertex(sh, idx, Tag<RawShape>());
+}
+
+template<class RawShape>
+inline size_t contourVertexCount(const RawShape& sh)
+{
+    return shapelike::cend(sh) - shapelike::cbegin(sh);
+}
+
+template<class RawShape, class Fn>
+inline void foreachVertex(RawShape& sh, Fn fn, const PolygonTag&) {
+    foreachVertex(contour(sh), fn, PathTag());
+    for(auto& h : holes(sh)) foreachVertex(h, fn, PathTag());
+}
+
+template<class RawShape, class Fn>
+inline void foreachVertex(RawShape& sh, Fn fn) {
+    foreachVertex(sh, fn, Tag<RawShape>());
+}
+
+template<class Poly> inline bool isConvex(const Poly& sh, const PolygonTag&)
+{
+    bool convex = true;
+    convex &= isConvex(contour(sh), PathTag());
+    convex &= holeCount(sh) == 0;
+    return convex;
+}
+
+template<class RawShape> inline bool isConvex(const RawShape& sh) // dispatch
+{
+    return isConvex(sh, Tag<RawShape>());
+}
+
+}
+
+#define DECLARE_MAIN_TYPES(T)        \
+    using Polygon = T;               \
+    using Point   = TPoint<T>;       \
+    using Coord   = TCoord<Point>;   \
+    using Contour = TContour<T>;     \
+    using Box     = _Box<Point>;     \
+    using Circle  = _Circle<Point>;  \
+    using Segment = _Segment<Point>; \
+    using Polygons = TMultiShape<T>
+
+}
+
+#endif // GEOMETRY_TRAITS_HPP

src/libnest2d/include/libnest2d/geometry_traits_nfp.hpp

@@ -0,0 +1,736 @@
+#ifndef GEOMETRIES_NOFITPOLYGON_HPP
+#define GEOMETRIES_NOFITPOLYGON_HPP
+
+#include "geometry_traits.hpp"
+#include <algorithm>
+#include <functional>
+#include <vector>
+#include <iterator>
+
+namespace libnest2d {
+
+namespace __nfp {
+// Do not specialize this...
+template<class RawShape>
+inline bool _vsort(const TPoint<RawShape>& v1, const TPoint<RawShape>& v2)
+{
+    using Coord = TCoord<TPoint<RawShape>>;
+    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 x1 < x2;
+
+    return diff < 0;
+}
+
+template<class EdgeList, class RawShape, class Vertex = TPoint<RawShape>>
+inline void buildPolygon(const EdgeList& edgelist,
+                         RawShape& rpoly,
+                         Vertex& top_nfp)
+{
+    namespace sl = shapelike;
+
+    auto& rsh = sl::contour(rpoly);
+
+    sl::reserve(rsh, 2*edgelist.size());
+
+    // Add the two vertices from the first edge into the final polygon.
+    sl::addVertex(rsh, edgelist.front().first());
+    sl::addVertex(rsh, edgelist.front().second());
+
+    // Sorting function for the nfp reference vertex search
+    auto& cmp = _vsort<RawShape>;
+
+    // the reference (rightmost top) vertex so far
+    top_nfp = *std::max_element(sl::cbegin(rsh), sl::cend(rsh), cmp );
+
+    auto tmp = std::next(sl::begin(rsh));
+
+    // Construct final nfp by placing each edge to the end of the previous
+    for(auto eit = std::next(edgelist.begin());
+        eit != edgelist.end();
+        ++eit)
+    {
+        auto d = *tmp - eit->first();
+        Vertex p = eit->second() + d;
+
+        sl::addVertex(rsh, p);
+
+        // Set the new reference vertex
+        if(cmp(top_nfp, p)) top_nfp = p;
+
+        tmp = std::next(tmp);
+    }
+
+}
+
+template<class Container, class Iterator = typename Container::iterator>
+void advance(Iterator& it, Container& cont, bool direction)
+{
+    int dir = direction ? 1 : -1;
+    if(dir < 0 && it == cont.begin()) it = std::prev(cont.end());
+    else it += dir;
+    if(dir > 0 && it == cont.end()) it = cont.begin();
+}
+
+}
+
+/// A collection of static methods for handling the no fit polygon creation.
+namespace nfp {
+
+const double BP2D_CONSTEXPR TwoPi = 2*Pi;
+
+/// The complexity level of a polygon that an NFP implementation can handle.
+enum class NfpLevel: unsigned {
+    CONVEX_ONLY,
+    ONE_CONVEX,
+    BOTH_CONCAVE,
+    ONE_CONVEX_WITH_HOLES,
+    BOTH_CONCAVE_WITH_HOLES
+};
+
+template<class RawShape>
+using NfpResult = std::pair<RawShape, TPoint<RawShape>>;
+
+template<class RawShape> struct MaxNfpLevel {
+    static const BP2D_CONSTEXPR NfpLevel value = NfpLevel::CONVEX_ONLY;
+};
+
+
+// Shorthand for a pile of polygons
+template<class RawShape>
+using Shapes = TMultiShape<RawShape>;
+
+/**
+ * Merge a bunch of polygons with the specified additional polygon.
+ *
+ * \tparam RawShape the Polygon data type.
+ * \param shc The pile of polygons that will be unified with sh.
+ * \param sh A single polygon to unify with shc.
+ *
+ * \return A set of polygons that is the union of the input polygons. Note that
+ * mostly it will be a set containing only one big polygon but if the input
+ * polygons are disjunct than the resulting set will contain more polygons.
+ */
+template<class RawShapes>
+inline RawShapes merge(const RawShapes& /*shc*/)
+{
+    static_assert(always_false<RawShapes>::value,
+                  "Nfp::merge(shapes, shape) unimplemented!");
+}
+
+/**
+ * Merge a bunch of polygons with the specified additional polygon.
+ *
+ * \tparam RawShape the Polygon data type.
+ * \param shc The pile of polygons that will be unified with sh.
+ * \param sh A single polygon to unify with shc.
+ *
+ * \return A set of polygons that is the union of the input polygons. Note that
+ * mostly it will be a set containing only one big polygon but if the input
+ * polygons are disjunct than the resulting set will contain more polygons.
+ */
+template<class RawShape>
+inline TMultiShape<RawShape> merge(const TMultiShape<RawShape>& shc,
+                                   const RawShape& sh)
+{
+    auto m = nfp::merge(shc);
+    m.emplace_back(sh);
+    return nfp::merge(m);
+}
+
+/**
+ * Get the vertex of the polygon that is at the lowest values (bottom) in the Y
+ * axis and if there are more than one vertices on the same Y coordinate than
+ * the result will be the leftmost (with the highest X coordinate).
+ */
+template<class RawShape>
+inline TPoint<RawShape> leftmostDownVertex(const RawShape& sh)
+{
+
+    // find min x and min y vertex
+    auto it = std::min_element(shapelike::cbegin(sh), shapelike::cend(sh),
+                               __nfp::_vsort<RawShape>);
+
+    return it == shapelike::cend(sh) ? TPoint<RawShape>() : *it;;
+}
+
+/**
+ * Get the vertex of the polygon that is at the highest values (top) in the Y
+ * axis and if there are more than one vertices on the same Y coordinate than
+ * the result will be the rightmost (with the lowest X coordinate).
+ */
+template<class RawShape>
+TPoint<RawShape> rightmostUpVertex(const RawShape& sh)
+{
+
+    // find max x and max y vertex
+    auto it = std::max_element(shapelike::cbegin(sh), shapelike::cend(sh),
+                               __nfp::_vsort<RawShape>);
+
+    return it == shapelike::cend(sh) ? TPoint<RawShape>() : *it;
+}
+
+/**
+ * A method to get a vertex from a polygon that always maintains a relative
+ * position to the coordinate system: It is always the rightmost top vertex.
+ *
+ * This way it does not matter in what order the vertices are stored, the
+ * reference will be always the same for the same polygon.
+ */
+template<class RawShape>
+inline TPoint<RawShape> referenceVertex(const RawShape& sh)
+{
+    return rightmostUpVertex(sh);
+}
+
+/**
+ * The "trivial" Cuninghame-Green implementation of NFP for convex polygons.
+ *
+ * You can use this even if you provide implementations for the more complex
+ * cases (Through specializing the the NfpImpl struct). Currently, no other
+ * cases are covered in the library.
+ *
+ * Complexity should be no more than nlogn (std::sort) in the number of edges
+ * of the input polygons.
+ *
+ * \tparam RawShape the Polygon data type.
+ * \param sh The stationary polygon
+ * \param cother The orbiting polygon
+ * \return Returns a pair of the NFP and its reference vertex of the two input
+ * polygons which have to be strictly convex. The resulting NFP is proven to be
+ * convex as well in this case.
+ *
+ */
+template<class RawShape>
+inline NfpResult<RawShape> nfpConvexOnly(const RawShape& sh,
+                                         const RawShape& other)
+{
+    using Vertex = TPoint<RawShape>; using Edge = _Segment<Vertex>;
+    namespace sl = shapelike;
+
+    RawShape rsh;   // Final nfp placeholder
+    Vertex top_nfp;
+    std::vector<Edge> edgelist;
+
+    auto cap = sl::contourVertexCount(sh) + sl::contourVertexCount(other);
+
+    // Reserve the needed memory
+    edgelist.reserve(cap);
+    sl::reserve(rsh, static_cast<unsigned long>(cap));
+
+    { // place all edges from sh into edgelist
+        auto first = sl::cbegin(sh);
+        auto next = std::next(first);
+
+        while(next != sl::cend(sh)) {
+            edgelist.emplace_back(*(first), *(next));
+            ++first; ++next;
+        }
+    }
+
+    { // place all edges from other into edgelist
+        auto first = sl::cbegin(other);
+        auto next = std::next(first);
+
+        while(next != sl::cend(other)) {
+            edgelist.emplace_back(*(next), *(first));
+            ++first; ++next;
+        }
+    }
+
+    // Sort the edges by angle to X axis.
+    std::sort(edgelist.begin(), edgelist.end(),
+              [](const Edge& e1, const Edge& e2)
+    {
+        return e1.angleToXaxis() > e2.angleToXaxis();
+    });
+
+    __nfp::buildPolygon(edgelist, rsh, top_nfp);
+
+    return {rsh, top_nfp};
+}
+
+template<class RawShape>
+NfpResult<RawShape> nfpSimpleSimple(const RawShape& cstationary,
+                                           const RawShape& cother)
+{
+
+    // Algorithms are from the original algorithm proposed in paper:
+    // https://eprints.soton.ac.uk/36850/1/CORMSIS-05-05.pdf
+
+    // /////////////////////////////////////////////////////////////////////////
+    // Algorithm 1: Obtaining the minkowski sum
+    // /////////////////////////////////////////////////////////////////////////
+
+    // I guess this is not a full minkowski sum of the two input polygons by
+    // definition. This yields a subset that is compatible with the next 2
+    // algorithms.
+
+    using Result = NfpResult<RawShape>;
+    using Vertex = TPoint<RawShape>;
+    using Coord = TCoord<Vertex>;
+    using Edge = _Segment<Vertex>;
+    namespace sl = shapelike;
+    using std::signbit;
+    using std::sort;
+    using std::vector;
+    using std::ref;
+    using std::reference_wrapper;
+
+    // TODO The original algorithms expects the stationary polygon in
+    // counter clockwise and the orbiter in clockwise order.
+    // So for preventing any further complication, I will make the input
+    // the way it should be, than make my way around the orientations.
+
+    // Reverse the stationary contour to counter clockwise
+    auto stcont = sl::contour(cstationary);
+    {
+        std::reverse(sl::begin(stcont), sl::end(stcont));
+        stcont.pop_back();
+        auto it = std::min_element(sl::begin(stcont), sl::end(stcont),
+                               [](const Vertex& v1, const Vertex& v2) {
+            return getY(v1) < getY(v2);
+        });
+        std::rotate(sl::begin(stcont), it, sl::end(stcont));
+        sl::addVertex(stcont, sl::front(stcont));
+    }
+    RawShape stationary;
+    sl::contour(stationary) = stcont;
+
+    // Reverse the orbiter contour to counter clockwise
+    auto orbcont = sl::contour(cother);
+    {
+        std::reverse(orbcont.begin(), orbcont.end());
+
+        // Step 1: Make the orbiter reverse oriented
+
+        orbcont.pop_back();
+        auto it = std::min_element(orbcont.begin(), orbcont.end(),
+                              [](const Vertex& v1, const Vertex& v2) {
+            return getY(v1) < getY(v2);
+        });
+
+        std::rotate(orbcont.begin(), it, orbcont.end());
+        orbcont.emplace_back(orbcont.front());
+
+        for(auto &v : orbcont) v = -v;
+
+    }
+
+    // Copy the orbiter (contour only), we will have to work on it
+    RawShape orbiter;
+    sl::contour(orbiter) = orbcont;
+
+    // An edge with additional data for marking it
+    struct MarkedEdge {
+        Edge e; Radians turn_angle = 0; bool is_turning_point = false;
+        MarkedEdge() = default;
+        MarkedEdge(const Edge& ed, Radians ta, bool tp):
+            e(ed), turn_angle(ta), is_turning_point(tp) {}
+
+        // debug
+        std::string label;
+    };
+
+    // Container for marked edges
+    using EdgeList = vector<MarkedEdge>;
+
+    EdgeList A, B;
+
+    // This is how an edge list is created from the polygons
+    auto fillEdgeList = [](EdgeList& L, const RawShape& ppoly, int dir) {
+        auto& poly = sl::contour(ppoly);
+
+        L.reserve(sl::contourVertexCount(poly));
+
+        if(dir > 0) {
+            auto it = poly.begin();
+            auto nextit = std::next(it);
+
+            double turn_angle = 0;
+            bool is_turn_point = false;
+
+            while(nextit != poly.end()) {
+                L.emplace_back(Edge(*it, *nextit), turn_angle, is_turn_point);
+                it++; nextit++;
+            }
+        } else {
+            auto it = sl::rbegin(poly);
+            auto nextit = std::next(it);
+
+            double turn_angle = 0;
+            bool is_turn_point = false;
+
+            while(nextit != sl::rend(poly)) {
+                L.emplace_back(Edge(*it, *nextit), turn_angle, is_turn_point);
+                it++; nextit++;
+            }
+        }
+
+        auto getTurnAngle = [](const Edge& e1, const Edge& e2) {
+            auto phi = e1.angleToXaxis();
+            auto phi_prev = e2.angleToXaxis();
+            auto turn_angle = phi-phi_prev;
+            if(turn_angle > Pi) turn_angle -= TwoPi;
+            if(turn_angle < -Pi) turn_angle += TwoPi;
+            return turn_angle;
+        };
+
+        auto eit = L.begin();
+        auto enext = std::next(eit);
+
+        eit->turn_angle = getTurnAngle(L.front().e, L.back().e);
+
+        while(enext != L.end()) {
+            enext->turn_angle = getTurnAngle( enext->e, eit->e);
+            eit->is_turning_point =
+                    signbit(enext->turn_angle) != signbit(eit->turn_angle);
+            ++eit; ++enext;
+        }
+
+        L.back().is_turning_point = signbit(L.back().turn_angle) !=
+                                    signbit(L.front().turn_angle);
+
+    };
+
+    // Step 2: Fill the edgelists
+    fillEdgeList(A, stationary, 1);
+    fillEdgeList(B, orbiter, 1);
+
+    int i = 1;
+    for(MarkedEdge& me : A) {
+        std::cout << "a" << i << ":\n\t"
+                  << getX(me.e.first()) << " " << getY(me.e.first()) << "\n\t"
+                  << getX(me.e.second()) << " " << getY(me.e.second()) << "\n\t"
+                  << "Turning point: " << (me.is_turning_point ? "yes" : "no")
+                  << std::endl;
+
+        me.label = "a"; me.label += std::to_string(i);
+        i++;
+    }
+
+    i = 1;
+    for(MarkedEdge& me : B) {
+        std::cout << "b" << i << ":\n\t"
+                  << getX(me.e.first()) << " " << getY(me.e.first()) << "\n\t"
+                  << getX(me.e.second()) << " " << getY(me.e.second()) << "\n\t"
+                  << "Turning point: " << (me.is_turning_point ? "yes" : "no")
+                  << std::endl;
+        me.label = "b"; me.label += std::to_string(i);
+        i++;
+    }
+
+    // A reference to a marked edge that also knows its container
+    struct MarkedEdgeRef {
+        reference_wrapper<MarkedEdge> eref;
+        reference_wrapper<vector<MarkedEdgeRef>> container;
+        Coord dir = 1;  // Direction modifier
+
+        inline Radians angleX() const { return eref.get().e.angleToXaxis(); }
+        inline const Edge& edge() const { return eref.get().e; }
+        inline Edge& edge() { return eref.get().e; }
+        inline bool isTurningPoint() const {
+            return eref.get().is_turning_point;
+        }
+        inline bool isFrom(const vector<MarkedEdgeRef>& cont ) {
+            return &(container.get()) == &cont;
+        }
+        inline bool eq(const MarkedEdgeRef& mr) {
+            return &(eref.get()) == &(mr.eref.get());
+        }
+
+        MarkedEdgeRef(reference_wrapper<MarkedEdge> er,
+                      reference_wrapper<vector<MarkedEdgeRef>> ec):
+            eref(er), container(ec), dir(1) {}
+
+        MarkedEdgeRef(reference_wrapper<MarkedEdge> er,
+                      reference_wrapper<vector<MarkedEdgeRef>> ec,
+                      Coord d):
+            eref(er), container(ec), dir(d) {}
+    };
+
+    using EdgeRefList = vector<MarkedEdgeRef>;
+
+    // Comparing two marked edges
+    auto sortfn = [](const MarkedEdgeRef& e1, const MarkedEdgeRef& e2) {
+        return e1.angleX() < e2.angleX();
+    };
+
+    EdgeRefList Aref, Bref;     // We create containers for the references
+    Aref.reserve(A.size()); Bref.reserve(B.size());
+
+    // Fill reference container for the stationary polygon
+    std::for_each(A.begin(), A.end(), [&Aref](MarkedEdge& me) {
+        Aref.emplace_back( ref(me), ref(Aref) );
+    });
+
+    // Fill reference container for the orbiting polygon
+    std::for_each(B.begin(), B.end(), [&Bref](MarkedEdge& me) {
+        Bref.emplace_back( ref(me), ref(Bref) );
+    });
+
+    auto mink = [sortfn] // the Mink(Q, R, direction) sub-procedure
+            (const EdgeRefList& Q, const EdgeRefList& R, bool positive)
+    {
+
+        // Step 1 "merge sort_list(Q) and sort_list(R) to form merge_list(Q,R)"
+        // Sort the containers of edge references and merge them.
+        // Q could be sorted only once and be reused here but we would still
+        // need to merge it with sorted(R).
+
+        EdgeRefList merged;
+        EdgeRefList S, seq;
+        merged.reserve(Q.size() + R.size());
+
+        merged.insert(merged.end(), R.begin(), R.end());
+        std::stable_sort(merged.begin(), merged.end(), sortfn);
+        merged.insert(merged.end(), Q.begin(), Q.end());
+        std::stable_sort(merged.begin(), merged.end(), sortfn);
+
+        // Step 2 "set i = 1, k = 1, direction = 1, s1 = q1"
+        // we don't use i, instead, q is an iterator into Q. k would be an index
+        // into the merged sequence but we use "it" as an iterator for that
+
+        // here we obtain references for the containers for later comparisons
+        const auto& Rcont = R.begin()->container.get();
+        const auto& Qcont = Q.begin()->container.get();
+
+        // Set the initial direction
+        Coord dir = 1;
+
+        // roughly i = 1 (so q = Q.begin()) and s1 = q1 so S[0] = q;
+        if(positive) {
+            auto q = Q.begin();
+            S.emplace_back(*q);
+
+            // Roughly step 3
+
+            std::cout << "merged size: " << merged.size() << std::endl;
+            auto mit = merged.begin();
+            for(bool finish = false; !finish && q != Q.end();) {
+                ++q; // "Set i = i + 1"
+
+                while(!finish && mit != merged.end()) {
+                    if(mit->isFrom(Rcont)) {
+                        auto s = *mit;
+                        s.dir = dir;
+                        S.emplace_back(s);
+                    }
+
+                    if(mit->eq(*q)) {
+                        S.emplace_back(*q);
+                        if(mit->isTurningPoint()) dir = -dir;
+                        if(q == Q.begin()) finish = true;
+                        break;
+                    }
+
+                    mit += dir;
+    //                __nfp::advance(mit, merged, dir > 0);
+                }
+            }
+        } else {
+            auto q = Q.rbegin();
+            S.emplace_back(*q);
+
+            // Roughly step 3
+
+            std::cout << "merged size: " << merged.size() << std::endl;
+            auto mit = merged.begin();
+            for(bool finish = false; !finish && q != Q.rend();) {
+                ++q; // "Set i = i + 1"
+
+                while(!finish && mit != merged.end()) {
+                    if(mit->isFrom(Rcont)) {
+                        auto s = *mit;
+                        s.dir = dir;
+                        S.emplace_back(s);
+                    }
+
+                    if(mit->eq(*q)) {
+                        S.emplace_back(*q);
+                        S.back().dir = -1;
+                        if(mit->isTurningPoint()) dir = -dir;
+                        if(q == Q.rbegin()) finish = true;
+                        break;
+                    }
+
+                    mit += dir;
+            //                __nfp::advance(mit, merged, dir > 0);
+                }
+            }
+        }
+
+
+        // Step 4:
+
+        // "Let starting edge r1 be in position si in sequence"
+        // whaaat? I guess this means the following:
+        auto it = S.begin();
+        while(!it->eq(*R.begin())) ++it;
+
+        // "Set j = 1, next = 2, direction = 1, seq1 = si"
+        // we don't use j, seq is expanded dynamically.
+        dir = 1;
+        auto next = std::next(R.begin()); seq.emplace_back(*it);
+
+        // Step 5:
+        // "If all si edges have been allocated to seqj" should mean that
+        // we loop until seq has equal size with S
+        auto send = it; //it == S.begin() ? it : std::prev(it);
+        while(it != S.end()) {
+            ++it; if(it == S.end()) it = S.begin();
+            if(it == send) break;
+
+            if(it->isFrom(Qcont)) {
+                seq.emplace_back(*it); // "If si is from Q, j = j + 1, seqj = si"
+
+                // "If si is a turning point in Q,
+                // direction = - direction, next = next + direction"
+                if(it->isTurningPoint()) {
+                    dir = -dir;
+                    next += dir;
+//                    __nfp::advance(next, R, dir > 0);
+                }
+            }
+
+            if(it->eq(*next) /*&& dir == next->dir*/) { // "If si = direction.rnext"
+                // "j = j + 1, seqj = si, next = next + direction"
+                seq.emplace_back(*it);
+                next += dir;
+//                __nfp::advance(next, R, dir > 0);
+            }
+        }
+
+        return seq;
+    };
+
+    std::vector<EdgeRefList> seqlist;
+    seqlist.reserve(Bref.size());
+
+    EdgeRefList Bslope = Bref;  // copy Bref, we will make a slope diagram
+
+    // make the slope diagram of B
+    std::sort(Bslope.begin(), Bslope.end(), sortfn);
+
+    auto slopeit = Bslope.begin(); // search for the first turning point
+    while(!slopeit->isTurningPoint() && slopeit != Bslope.end()) slopeit++;
+
+    if(slopeit == Bslope.end()) {
+        // no turning point means convex polygon.
+        seqlist.emplace_back(mink(Aref, Bref, true));
+    } else {
+        int dir = 1;
+
+        auto firstturn = Bref.begin();
+        while(!firstturn->eq(*slopeit)) ++firstturn;
+
+        assert(firstturn != Bref.end());
+
+        EdgeRefList bgroup; bgroup.reserve(Bref.size());
+        bgroup.emplace_back(*slopeit);
+
+        auto b_it = std::next(firstturn);
+        while(b_it != firstturn) {
+            if(b_it == Bref.end()) b_it = Bref.begin();
+
+            while(!slopeit->eq(*b_it)) {
+                __nfp::advance(slopeit, Bslope, dir > 0);
+            }
+
+            if(!slopeit->isTurningPoint()) {
+                bgroup.emplace_back(*slopeit);
+            } else {
+                if(!bgroup.empty()) {
+                    if(dir > 0) bgroup.emplace_back(*slopeit);
+                    for(auto& me : bgroup) {
+                        std::cout << me.eref.get().label << ", ";
+                    }
+                    std::cout << std::endl;
+                    seqlist.emplace_back(mink(Aref, bgroup, dir == 1 ? true : false));
+                    bgroup.clear();
+                    if(dir < 0) bgroup.emplace_back(*slopeit);
+                } else {
+                    bgroup.emplace_back(*slopeit);
+                }
+
+                dir *= -1;
+            }
+            ++b_it;
+        }
+    }
+
+//    while(it != Bref.end()) // This is step 3 and step 4 in one loop
+//        if(it->isTurningPoint()) {
+//            R = {R.last, it++};
+//            auto seq = mink(Q, R, orientation);
+
+//            // TODO step 6 (should be 5 shouldn't it?): linking edges from A
+//            // I don't get this step
+
+//            seqlist.insert(seqlist.end(), seq.begin(), seq.end());
+//            orientation = !orientation;
+//        } else ++it;
+
+//    if(seqlist.empty()) seqlist = mink(Q, {Bref.begin(), Bref.end()}, true);
+
+    // /////////////////////////////////////////////////////////////////////////
+    // Algorithm 2: breaking Minkowski sums into track line trips
+    // /////////////////////////////////////////////////////////////////////////
+
+
+    // /////////////////////////////////////////////////////////////////////////
+    // Algorithm 3: finding the boundary of the NFP from track line trips
+    // /////////////////////////////////////////////////////////////////////////
+
+
+    for(auto& seq : seqlist) {
+        std::cout << "seqlist size: " << seq.size() << std::endl;
+        for(auto& s : seq) {
+            std::cout << (s.dir > 0 ? "" : "-") << s.eref.get().label << ", ";
+        }
+        std::cout << std::endl;
+    }
+
+    auto& seq = seqlist.front();
+    RawShape rsh;
+    Vertex top_nfp;
+    std::vector<Edge> edgelist; edgelist.reserve(seq.size());
+    for(auto& s : seq) {
+        edgelist.emplace_back(s.eref.get().e);
+    }
+
+    __nfp::buildPolygon(edgelist, rsh, top_nfp);
+
+    return Result(rsh, top_nfp);
+}
+
+// Specializable NFP implementation class. Specialize it if you have a faster
+// or better NFP implementation
+template<class RawShape, NfpLevel nfptype>
+struct NfpImpl {
+    NfpResult<RawShape> operator()(const RawShape& sh, const RawShape& other)
+    {
+        static_assert(nfptype == NfpLevel::CONVEX_ONLY,
+                      "Nfp::noFitPolygon() unimplemented!");
+
+        // Libnest2D has a default implementation for convex polygons and will
+        // use it if feasible.
+        return nfpConvexOnly(sh, other);
+    }
+};
+
+/// Helper function to get the NFP
+template<NfpLevel nfptype, class RawShape>
+inline NfpResult<RawShape> noFitPolygon(const RawShape& sh,
+                                        const RawShape& other)
+{
+    NfpImpl<RawShape, nfptype> nfps;
+    return nfps(sh, other);
+}
+
+}
+
+}
+
+#endif // GEOMETRIES_NOFITPOLYGON_HPP

src/libnest2d/include/libnest2d/libnest2d.hpp

@@ -0,0 +1,988 @@
+#ifndef LIBNEST2D_HPP
+#define LIBNEST2D_HPP
+
+#include <memory>
+#include <vector>
+#include <map>
+#include <array>
+#include <algorithm>
+#include <functional>
+
+#include "geometry_traits.hpp"
+
+namespace libnest2d {
+
+namespace sl = shapelike;
+namespace pl = pointlike;
+
+/**
+ * \brief An item to be placed on a bin.
+ *
+ * It holds a copy of the original shape object but supports move construction
+ * from the shape objects if its an rvalue reference. This way we can construct
+ * the items without the cost of copying a potentially large amount of input.
+ *
+ * The results of some calculations are cached for maintaining fast run times.
+ * For this reason, memory demands are much higher but this should pay off.
+ */
+template<class RawShape>
+class _Item {
+    using Coord = TCoord<TPoint<RawShape>>;
+    using Vertex = TPoint<RawShape>;
+    using Box = _Box<Vertex>;
+
+    using VertexConstIterator = typename TContour<RawShape>::const_iterator;
+
+    // The original shape that gets encapsulated.
+    RawShape sh_;
+
+    // Transformation data
+    Vertex translation_;
+    Radians rotation_;
+    Coord offset_distance_;
+
+    // Info about whether the transformations will have to take place
+    // This is needed because if floating point is used, it is hard to say
+    // that a zero angle is not a rotation because of testing for equality.
+    bool has_rotation_ = false, has_translation_ = false, has_offset_ = false;
+
+    // For caching the calculations as they can get pretty expensive.
+    mutable RawShape tr_cache_;
+    mutable bool tr_cache_valid_ = false;
+    mutable double area_cache_ = 0;
+    mutable bool area_cache_valid_ = false;
+    mutable RawShape offset_cache_;
+    mutable bool offset_cache_valid_ = false;
+
+    enum class Convexity: char {
+        UNCHECKED,
+        C_TRUE,
+        C_FALSE
+    };
+
+    mutable Convexity convexity_ = Convexity::UNCHECKED;
+    mutable VertexConstIterator rmt_;    // rightmost top vertex
+    mutable VertexConstIterator lmb_;    // leftmost bottom vertex
+    mutable bool rmt_valid_ = false, lmb_valid_ = false;
+    mutable struct BBCache {
+        Box bb; bool valid;
+        BBCache(): valid(false) {}
+    } bb_cache_;
+
+public:
+
+    /// The type of the shape which was handed over as the template argument.
+    using ShapeType = RawShape;
+
+    /**
+     * \brief Iterator type for the outer vertices.
+     *
+     * Only const iterators can be used. The _Item type is not intended to
+     * modify the carried shapes from the outside. The main purpose of this type
+     * is to cache the calculation results from the various operators it
+     * supports. Giving out a non const iterator would make it impossible to
+     * perform correct cache invalidation.
+     */
+    using Iterator = VertexConstIterator;
+
+    /**
+     * @brief Get the orientation of the polygon.
+     *
+     * The orientation have to be specified as a specialization of the
+     * OrientationType struct which has a Value constant.
+     *
+     * @return The orientation type identifier for the _Item type.
+     */
+    static BP2D_CONSTEXPR Orientation orientation() {
+        return OrientationType<RawShape>::Value;
+    }
+
+    /**
+     * @brief Constructing an _Item form an existing raw shape. The shape will
+     * be copied into the _Item object.
+     * @param sh The original shape object.
+     */
+    explicit inline _Item(const RawShape& sh): sh_(sh) {}
+
+    /**
+     * @brief Construction of an item by moving the content of the raw shape,
+     * assuming that it supports move semantics.
+     * @param sh The original shape object.
+     */
+    explicit inline _Item(RawShape&& sh): sh_(std::move(sh)) {}
+
+    /**
+     * @brief Create an item from an initializer list.
+     * @param il The initializer list of vertices.
+     */
+    inline _Item(const std::initializer_list< Vertex >& il):
+        sh_(sl::create<RawShape>(il)) {}
+
+    inline _Item(const TContour<RawShape>& contour,
+                 const THolesContainer<RawShape>& holes = {}):
+        sh_(sl::create<RawShape>(contour, holes)) {}
+
+    inline _Item(TContour<RawShape>&& contour,
+                 THolesContainer<RawShape>&& holes):
+        sh_(sl::create<RawShape>(std::move(contour),
+                                        std::move(holes))) {}
+
+    /**
+     * @brief Convert the polygon to string representation. The format depends
+     * on the implementation of the polygon.
+     * @return
+     */
+    inline std::string toString() const
+    {
+        return sl::toString(sh_);
+    }
+
+    /// Iterator tho the first contour vertex in the polygon.
+    inline Iterator begin() const
+    {
+        return sl::cbegin(sh_);
+    }
+
+    /// Alias to begin()
+    inline Iterator cbegin() const
+    {
+        return sl::cbegin(sh_);
+    }
+
+    /// Iterator to the last contour vertex.
+    inline Iterator end() const
+    {
+        return sl::cend(sh_);
+    }
+
+    /// Alias to end()
+    inline Iterator cend() const
+    {
+        return sl::cend(sh_);
+    }
+
+    /**
+     * @brief Get a copy of an outer vertex within the carried shape.
+     *
+     * Note that the vertex considered here is taken from the original shape
+     * that this item is constructed from. This means that no transformation is
+     * applied to the shape in this call.
+     *
+     * @param idx The index of the requested vertex.
+     * @return A copy of the requested vertex.
+     */
+    inline Vertex vertex(unsigned long idx) const
+    {
+        return sl::vertex(sh_, idx);
+    }
+
+    /**
+     * @brief Modify a vertex.
+     *
+     * Note that this method will invalidate every cached calculation result
+     * including polygon offset and transformations.
+     *
+     * @param idx The index of the requested vertex.
+     * @param v The new vertex data.
+     */
+    inline void setVertex(unsigned long idx, const Vertex& v )
+    {
+        invalidateCache();
+        sl::vertex(sh_, idx) = v;
+    }
+
+    /**
+     * @brief Calculate the shape area.
+     *
+     * The method returns absolute value and does not reflect polygon
+     * orientation. The result is cached, subsequent calls will have very little
+     * cost.
+     * @return The shape area in floating point double precision.
+     */
+    inline double area() const {
+        double ret ;
+        if(area_cache_valid_) ret = area_cache_;
+        else {
+            ret = sl::area(offsettedShape());
+            area_cache_ = ret;
+            area_cache_valid_ = true;
+        }
+        return ret;
+    }
+
+    inline bool isContourConvex() const {
+        bool ret = false;
+
+        switch(convexity_) {
+        case Convexity::UNCHECKED:
+            ret = sl::isConvex(sl::contour(transformedShape()));
+            convexity_ = ret? Convexity::C_TRUE : Convexity::C_FALSE;
+            break;
+        case Convexity::C_TRUE: ret = true; break;
+        case Convexity::C_FALSE:;
+        }
+
+        return ret;
+    }
+
+    inline bool isHoleConvex(unsigned /*holeidx*/) const {
+        return false;
+    }
+
+    inline bool areHolesConvex() const {
+        return false;
+    }
+
+    /// The number of the outer ring vertices.
+    inline size_t vertexCount() const {
+        return sl::contourVertexCount(sh_);
+    }
+
+    inline size_t holeCount() const {
+        return sl::holeCount(sh_);
+    }
+
+    /**
+     * @brief isPointInside
+     * @param p
+     * @return
+     */
+    inline bool isInside(const Vertex& p) const
+    {
+        return sl::isInside(p, transformedShape());
+    }
+
+    inline bool isInside(const _Item& sh) const
+    {
+        return sl::isInside(transformedShape(), sh.transformedShape());
+    }
+
+    inline bool isInside(const RawShape& sh) const
+    {
+        return sl::isInside(transformedShape(), sh);
+    }
+
+    inline bool isInside(const _Box<TPoint<RawShape>>& box) const;
+    inline bool isInside(const _Circle<TPoint<RawShape>>& box) const;
+
+    inline void translate(const Vertex& d) BP2D_NOEXCEPT
+    {
+        translation(translation() + d);
+    }
+
+    inline void rotate(const Radians& rads) BP2D_NOEXCEPT
+    {
+        rotation(rotation() + rads);
+    }
+
+    inline void addOffset(Coord distance) BP2D_NOEXCEPT
+    {
+        offset_distance_ = distance;
+        has_offset_ = true;
+        invalidateCache();
+    }
+
+    inline void removeOffset() BP2D_NOEXCEPT {
+        has_offset_ = false;
+        invalidateCache();
+    }
+
+    inline Radians rotation() const BP2D_NOEXCEPT
+    {
+        return rotation_;
+    }
+
+    inline TPoint<RawShape> translation() const BP2D_NOEXCEPT
+    {
+        return translation_;
+    }
+
+    inline void rotation(Radians rot) BP2D_NOEXCEPT
+    {
+        if(rotation_ != rot) {
+            rotation_ = rot; has_rotation_ = true; tr_cache_valid_ = false;
+            rmt_valid_ = false; lmb_valid_ = false;
+            bb_cache_.valid = false;
+        }
+    }
+
+    inline void translation(const TPoint<RawShape>& tr) BP2D_NOEXCEPT
+    {
+        if(translation_ != tr) {
+            translation_ = tr; has_translation_ = true; tr_cache_valid_ = false;
+            //bb_cache_.valid = false;
+        }
+    }
+
+    inline const RawShape& transformedShape() const
+    {
+        if(tr_cache_valid_) return tr_cache_;
+
+        RawShape cpy = offsettedShape();
+        if(has_rotation_) sl::rotate(cpy, rotation_);
+        if(has_translation_) sl::translate(cpy, translation_);
+        tr_cache_ = cpy; tr_cache_valid_ = true;
+        rmt_valid_ = false; lmb_valid_ = false;
+
+        return tr_cache_;
+    }
+
+    inline operator RawShape() const
+    {
+        return transformedShape();
+    }
+
+    inline const RawShape& rawShape() const BP2D_NOEXCEPT
+    {
+        return sh_;
+    }
+
+    inline void resetTransformation() BP2D_NOEXCEPT
+    {
+        has_translation_ = false; has_rotation_ = false; has_offset_ = false;
+        invalidateCache();
+    }
+
+    inline Box boundingBox() const {
+        if(!bb_cache_.valid) {
+            if(!has_rotation_)
+                bb_cache_.bb = sl::boundingBox(offsettedShape());
+            else {
+                // TODO make sure this works
+                auto rotsh = offsettedShape();
+                sl::rotate(rotsh, rotation_);
+                bb_cache_.bb = sl::boundingBox(rotsh);
+            }
+            bb_cache_.valid = true;
+        }
+
+        auto &bb = bb_cache_.bb; auto &tr = translation_;
+        return {bb.minCorner() + tr, bb.maxCorner() + tr };
+    }
+
+    inline Vertex referenceVertex() const {
+        return rightmostTopVertex();
+    }
+
+    inline Vertex rightmostTopVertex() const {
+        if(!rmt_valid_ || !tr_cache_valid_) {  // find max x and max y vertex
+            auto& tsh = transformedShape();
+            rmt_ = std::max_element(sl::cbegin(tsh), sl::cend(tsh), vsort);
+            rmt_valid_ = true;
+        }
+        return *rmt_;
+    }
+
+    inline Vertex leftmostBottomVertex() const {
+        if(!lmb_valid_ || !tr_cache_valid_) {  // find min x and min y vertex
+            auto& tsh = transformedShape();
+            lmb_ = std::min_element(sl::cbegin(tsh), sl::cend(tsh), vsort);
+            lmb_valid_ = true;
+        }
+        return *lmb_;
+    }
+
+    //Static methods:
+
+    inline static bool intersects(const _Item& sh1, const _Item& sh2)
+    {
+        return sl::intersects(sh1.transformedShape(),
+                                     sh2.transformedShape());
+    }
+
+    inline static bool touches(const _Item& sh1, const _Item& sh2)
+    {
+        return sl::touches(sh1.transformedShape(),
+                                  sh2.transformedShape());
+    }
+
+private:
+
+    inline const RawShape& offsettedShape() const {
+        if(has_offset_ ) {
+            if(offset_cache_valid_) return offset_cache_;
+
+            offset_cache_ = sh_;
+            sl::offset(offset_cache_, offset_distance_);
+            offset_cache_valid_ = true;
+            return offset_cache_;
+        }
+        return sh_;
+    }
+
+    inline void invalidateCache() const BP2D_NOEXCEPT
+    {
+        tr_cache_valid_ = false;
+        lmb_valid_ = false; rmt_valid_ = false;
+        area_cache_valid_ = false;
+        offset_cache_valid_ = false;
+        bb_cache_.valid = false;
+        convexity_ = Convexity::UNCHECKED;
+    }
+
+    static inline bool vsort(const Vertex& v1, const Vertex& v2)
+    {
+        Coord &&x1 = getX(v1), &&x2 = getX(v2);
+        Coord &&y1 = getY(v1), &&y2 = getY(v2);
+        auto diff = y1 - y2;
+        if(std::abs(diff) <= std::numeric_limits<Coord>::epsilon())
+            return x1 < x2;
+
+        return diff < 0;
+    }
+};
+
+/**
+ * \brief Subclass of _Item for regular rectangle items.
+ */
+template<class RawShape>
+class _Rectangle: public _Item<RawShape> {
+    using _Item<RawShape>::vertex;
+    using TO = Orientation;
+public:
+
+    using Unit = TCoord<TPoint<RawShape>>;
+
+    template<TO o = OrientationType<RawShape>::Value>
+    inline _Rectangle(Unit width, Unit height,
+                      // disable this ctor if o != CLOCKWISE
+                      enable_if_t< o == TO::CLOCKWISE, int> = 0 ):
+        _Item<RawShape>( sl::create<RawShape>( {
+                                                        {0, 0},
+                                                        {0, height},
+                                                        {width, height},
+                                                        {width, 0},
+                                                        {0, 0}
+                                                      } ))
+    {
+    }
+
+    template<TO o = OrientationType<RawShape>::Value>
+    inline _Rectangle(Unit width, Unit height,
+                      // disable this ctor if o != COUNTER_CLOCKWISE
+                      enable_if_t< o == TO::COUNTER_CLOCKWISE, int> = 0 ):
+        _Item<RawShape>( sl::create<RawShape>( {
+                                                        {0, 0},
+                                                        {width, 0},
+                                                        {width, height},
+                                                        {0, height},
+                                                        {0, 0}
+                                                      } ))
+    {
+    }
+
+    inline Unit width() const BP2D_NOEXCEPT {
+        return getX(vertex(2));
+    }
+
+    inline Unit height() const BP2D_NOEXCEPT {
+        return getY(vertex(2));
+    }
+};
+
+template<class RawShape>
+inline bool _Item<RawShape>::isInside(const _Box<TPoint<RawShape>>& box) const {
+    return sl::isInside<RawShape>(boundingBox(), box);
+}
+
+template<class RawShape> inline bool
+_Item<RawShape>::isInside(const _Circle<TPoint<RawShape>>& circ) const {
+    return sl::isInside<RawShape>(transformedShape(), circ);
+}
+
+
+template<class I> using _ItemRef = std::reference_wrapper<I>;
+template<class I> using _ItemGroup = std::vector<_ItemRef<I>>;
+
+template<class Iterator>
+struct ConstItemRange {
+    Iterator from;
+    Iterator to;
+    bool valid = false;
+
+    ConstItemRange() = default;
+    ConstItemRange(Iterator f, Iterator t): from(f), to(t), valid(true) {}
+};
+
+template<class Container>
+inline ConstItemRange<typename Container::const_iterator>
+rem(typename Container::const_iterator it, const Container& cont) {
+    return {std::next(it), cont.end()};
+}
+
+/**
+ * \brief A wrapper interface (trait) class for any placement strategy provider.
+ *
+ * If a client wants to use its own placement algorithm, all it has to do is to
+ * specialize this class template and define all the ten methods it has. It can
+ * use the strategies::PlacerBoilerplace class for creating a new placement
+ * strategy where only the constructor and the trypack method has to be provided
+ * and it will work out of the box.
+ */
+template<class PlacementStrategy>
+class PlacementStrategyLike {
+    PlacementStrategy impl_;
+public:
+
+    /// The item type that the placer works with.
+    using Item = typename PlacementStrategy::Item;
+
+    /// The placer's config type. Should be a simple struct but can be anything.
+    using Config = typename PlacementStrategy::Config;
+
+    /**
+     * \brief The type of the bin that the placer works with.
+     *
+     * Can be a box or an arbitrary shape or just a width or height without a
+     * second dimension if an infinite bin is considered.
+     */
+    using BinType = typename PlacementStrategy::BinType;
+
+    /**
+     * \brief Pack result that can be used to accept or discard it. See trypack
+     * method.
+     */
+    using PackResult = typename PlacementStrategy::PackResult;
+
+    using ItemRef = _ItemRef<Item>;
+    using ItemGroup = _ItemGroup<Item>;
+    using DefaultIterator = typename ItemGroup::const_iterator;
+
+    /**
+     * @brief Constructor taking the bin and an optional configuration.
+     * @param bin The bin object whose type is defined by the placement strategy.
+     * @param config The configuration for the particular placer.
+     */
+    explicit PlacementStrategyLike(const BinType& bin,
+                                   const Config& config = Config()):
+        impl_(bin)
+    {
+        configure(config);
+    }
+
+    /**
+     * @brief Provide a different configuration for the placer.
+     *
+     * Note that it depends on the particular placer implementation how it
+     * reacts to config changes in the middle of a calculation.
+     *
+     * @param config The configuration object defined by the placement strategy.
+     */
+    inline void configure(const Config& config) { impl_.configure(config); }
+
+    /**
+     * Try to pack an item with a result object that contains the packing
+     * information for later accepting it.
+     *
+     * \param item_store A container of items that are intended to be packed
+     * later. Can be used by the placer to switch tactics. When it's knows that
+     * many items will come a greedy strategy may not be the best.
+     * \param from The iterator to the item from which the packing should start,
+     * including the pointed item
+     * \param count How many items should be packed. If the value is 1, than
+     * just the item pointed to by "from" argument should be packed.
+     */
+    template<class Iter = DefaultIterator>
+    inline PackResult trypack(
+            Item& item,
+            const ConstItemRange<Iter>& remaining = ConstItemRange<Iter>())
+    {
+        return impl_.trypack(item, remaining);
+    }
+
+    /**
+     * @brief A method to accept a previously tried item (or items).
+     *
+     * If the pack result is a failure the method should ignore it.
+     * @param r The result of a previous trypack call.
+     */
+    inline void accept(PackResult& r) { impl_.accept(r); }
+
+    /**
+     * @brief pack Try to pack and immediately accept it on success.
+     *
+     * A default implementation would be to call
+     * { auto&& r = trypack(...); accept(r); return r; } but we should let the
+     * implementor of the placement strategy to harvest any optimizations from
+     * the absence of an intermediate step. The above version can still be used
+     * in the implementation.
+     *
+     * @param item The item to pack.
+     * @return Returns true if the item was packed or false if it could not be
+     * packed.
+     */
+    template<class Range = ConstItemRange<DefaultIterator>>
+    inline bool pack(
+            Item& item,
+            const Range& remaining = Range())
+    {
+        return impl_.pack(item, remaining);
+    }
+
+    /// Unpack the last element (remove it from the list of packed items).
+    inline void unpackLast() { impl_.unpackLast(); }
+
+    /// Get the bin object.
+    inline const BinType& bin() const { return impl_.bin(); }
+
+    /// Set a new bin object.
+    inline void bin(const BinType& bin) { impl_.bin(bin); }
+
+    /// Get the packed items.
+    inline ItemGroup getItems() { return impl_.getItems(); }
+
+    /// Clear the packed items so a new session can be started.
+    inline void clearItems() { impl_.clearItems(); }
+
+    inline double filledArea() const { return impl_.filledArea(); }
+
+};
+
+// The progress function will be called with the number of placed items
+using ProgressFunction = std::function<void(unsigned)>;
+using StopCondition = std::function<bool(void)>;
+
+/**
+ * A wrapper interface (trait) class for any selections strategy provider.
+ */
+template<class SelectionStrategy>
+class SelectionStrategyLike {
+    SelectionStrategy impl_;
+public:
+    using Item = typename SelectionStrategy::Item;
+    using Config = typename SelectionStrategy::Config;
+
+    using ItemRef = std::reference_wrapper<Item>;
+    using ItemGroup = std::vector<ItemRef>;
+
+    /**
+     * @brief Provide a different configuration for the selection strategy.
+     *
+     * Note that it depends on the particular placer implementation how it
+     * reacts to config changes in the middle of a calculation.
+     *
+     * @param config The configuration object defined by the selection strategy.
+     */
+    inline void configure(const Config& config) {
+        impl_.configure(config);
+    }
+
+    /**
+     * @brief A function callback which should be called whenever an item or
+     * a group of items where successfully packed.
+     * @param fn A function callback object taking one unsigned integer as the
+     * number of the remaining items to pack.
+     */
+    void progressIndicator(ProgressFunction fn) { impl_.progressIndicator(fn); }
+
+    void stopCondition(StopCondition cond) { impl_.stopCondition(cond); }
+
+    /**
+     * \brief A method to start the calculation on the input sequence.
+     *
+     * \tparam TPlacer The only mandatory template parameter is the type of
+     * placer compatible with the PlacementStrategyLike interface.
+     *
+     * \param first, last The first and last iterator if the input sequence. It
+     * can be only an iterator of a type convertible to Item.
+     * \param bin. The shape of the bin. It has to be supported by the placement
+     * strategy.
+     * \param An optional config object for the placer.
+     */
+    template<class TPlacer, class TIterator,
+             class TBin = typename PlacementStrategyLike<TPlacer>::BinType,
+             class PConfig = typename PlacementStrategyLike<TPlacer>::Config>
+    inline void packItems(
+            TIterator first,
+            TIterator last,
+            TBin&& bin,
+            PConfig&& config = PConfig() )
+    {
+        impl_.template packItems<TPlacer>(first, last,
+                                 std::forward<TBin>(bin),
+                                 std::forward<PConfig>(config));
+    }
+
+    /**
+     * \brief Get the number of bins opened by the selection algorithm.
+     *
+     * Initially it is zero and after the call to packItems it will return
+     * the number of bins opened by the packing procedure.
+     *
+     * \return The number of bins opened.
+     */
+    inline size_t binCount() const { return impl_.binCount(); }
+
+    /**
+     * @brief Get the items for a particular bin.
+     * @param binIndex The index of the requested bin.
+     * @return Returns a list of all items packed into the requested bin.
+     */
+    inline ItemGroup itemsForBin(size_t binIndex) {
+        return impl_.itemsForBin(binIndex);
+    }
+
+    /// Same as itemsForBin but for a const context.
+    inline const ItemGroup itemsForBin(size_t binIndex) const {
+        return impl_.itemsForBin(binIndex);
+    }
+};
+
+
+/**
+ * \brief A list of packed item vectors. Each vector represents a bin.
+ */
+template<class RawShape>
+using _PackGroup = std::vector<
+                        std::vector<
+                            std::reference_wrapper<_Item<RawShape>>
+                        >
+                   >;
+
+/**
+ * \brief A list of packed (index, item) pair vectors. Each vector represents a
+ * bin.
+ *
+ * The index is points to the position of the item in the original input
+ * sequence. This way the caller can use the items as a transformation data
+ * carrier and transform the original objects manually.
+ */
+template<class RawShape>
+using _IndexedPackGroup = std::vector<
+                               std::vector<
+                                   std::pair<
+                                       unsigned,
+                                       std::reference_wrapper<_Item<RawShape>>
+                                   >
+                               >
+                          >;
+
+/**
+ * The Arranger is the front-end class for the libnest2d library. It takes the
+ * input items and outputs the items with the proper transformations to be
+ * inside the provided bin.
+ */
+template<class PlacementStrategy, class SelectionStrategy >
+class Nester {
+    using TSel = SelectionStrategyLike<SelectionStrategy>;
+    TSel selector_;
+public:
+    using Item = typename PlacementStrategy::Item;
+    using ItemRef = std::reference_wrapper<Item>;
+    using TPlacer = PlacementStrategyLike<PlacementStrategy>;
+    using BinType = typename TPlacer::BinType;
+    using PlacementConfig = typename TPlacer::Config;
+    using SelectionConfig = typename TSel::Config;
+
+    using Unit = TCoord<TPoint<typename Item::ShapeType>>;
+
+    using IndexedPackGroup = _IndexedPackGroup<typename Item::ShapeType>;
+    using PackGroup = _PackGroup<typename Item::ShapeType>;
+    using ResultType = PackGroup;
+    using ResultTypeIndexed = IndexedPackGroup;
+
+private:
+    BinType bin_;
+    PlacementConfig pconfig_;
+    Unit min_obj_distance_;
+
+    using SItem =  typename SelectionStrategy::Item;
+    using TPItem = remove_cvref_t<Item>;
+    using TSItem = remove_cvref_t<SItem>;
+
+    std::vector<TPItem> item_cache_;
+
+public:
+
+    /**
+     * \brief Constructor taking the bin as the only mandatory parameter.
+     *
+     * \param bin The bin shape that will be used by the placers. The type
+     * of the bin should be one that is supported by the placer type.
+     */
+    template<class TBinType = BinType,
+             class PConf = PlacementConfig,
+             class SConf = SelectionConfig>
+    Nester( TBinType&& bin,
+              Unit min_obj_distance = 0,
+              const PConf& pconfig = PConf(),
+              const SConf& sconfig = SConf()):
+        bin_(std::forward<TBinType>(bin)),
+        pconfig_(pconfig),
+        min_obj_distance_(min_obj_distance)
+    {
+        static_assert( std::is_same<TPItem, TSItem>::value,
+                       "Incompatible placement and selection strategy!");
+
+        selector_.configure(sconfig);
+    }
+
+    void configure(const PlacementConfig& pconf) { pconfig_ = pconf; }
+    void configure(const SelectionConfig& sconf) { selector_.configure(sconf); }
+    void configure(const PlacementConfig& pconf, const SelectionConfig& sconf) {
+        pconfig_ = pconf;
+        selector_.configure(sconf);
+    }
+    void configure(const SelectionConfig& sconf, const PlacementConfig& pconf) {
+        pconfig_ = pconf;
+        selector_.configure(sconf);
+    }
+
+    /**
+     * \brief Arrange an input sequence and return a PackGroup object with
+     * the packed groups corresponding to the bins.
+     *
+     * The number of groups in the pack group is the number of bins opened by
+     * the selection algorithm.
+     */
+    template<class TIterator>
+    inline PackGroup execute(TIterator from, TIterator to)
+    {
+        return _execute(from, to);
+    }
+
+    /**
+     * A version of the arrange method returning an IndexedPackGroup with
+     * the item indexes into the original input sequence.
+     *
+     * Takes a little longer to collect the indices. Scales linearly with the
+     * input sequence size.
+     */
+    template<class TIterator>
+    inline IndexedPackGroup executeIndexed(TIterator from, TIterator to)
+    {
+        return _executeIndexed(from, to);
+    }
+
+    /// Shorthand to normal arrange method.
+    template<class TIterator>
+    inline PackGroup operator() (TIterator from, TIterator to)
+    {
+        return _execute(from, to);
+    }
+
+    /// Set a progress indicator function object for the selector.
+    inline Nester& progressIndicator(ProgressFunction func)
+    {
+        selector_.progressIndicator(func); return *this;
+    }
+
+    /// Set a predicate to tell when to abort nesting.
+    inline Nester& stopCondition(StopCondition fn) {
+        selector_.stopCondition(fn); return *this;
+    }
+
+    inline PackGroup lastResult() {
+        PackGroup ret;
+        for(size_t i = 0; i < selector_.binCount(); i++) {
+            auto items = selector_.itemsForBin(i);
+            ret.push_back(items);
+        }
+        return ret;
+    }
+
+private:
+
+    template<class TIterator,
+             class IT = remove_cvref_t<typename TIterator::value_type>,
+
+             // This function will be used only if the iterators are pointing to
+             // a type compatible with the libnets2d::_Item template.
+             // This way we can use references to input elements as they will
+             // have to exist for the lifetime of this call.
+             class T = enable_if_t< std::is_convertible<IT, TPItem>::value, IT>
+             >
+    inline PackGroup _execute(TIterator from, TIterator to, bool = false)
+    {
+        __execute(from, to);
+        return lastResult();
+    }
+
+    template<class TIterator,
+             class IT = remove_cvref_t<typename TIterator::value_type>,
+             class T = enable_if_t<!std::is_convertible<IT, TPItem>::value, IT>
+             >
+    inline PackGroup _execute(TIterator from, TIterator to, int = false)
+    {
+        item_cache_ = {from, to};
+
+        __execute(item_cache_.begin(), item_cache_.end());
+        return lastResult();
+    }
+
+    template<class TIterator,
+             class IT = remove_cvref_t<typename TIterator::value_type>,
+
+             // This function will be used only if the iterators are pointing to
+             // a type compatible with the libnest2d::_Item template.
+             // This way we can use references to input elements as they will
+             // have to exist for the lifetime of this call.
+             class T = enable_if_t< std::is_convertible<IT, TPItem>::value, IT>
+             >
+    inline IndexedPackGroup _executeIndexed(TIterator from,
+                                            TIterator to,
+                                            bool = false)
+    {
+        __execute(from, to);
+        return createIndexedPackGroup(from, to, selector_);
+    }
+
+    template<class TIterator,
+             class IT = remove_cvref_t<typename TIterator::value_type>,
+             class T = enable_if_t<!std::is_convertible<IT, TPItem>::value, IT>
+             >
+    inline IndexedPackGroup _executeIndexed(TIterator from,
+                                            TIterator to,
+                                            int = false)
+    {
+        item_cache_ = {from, to};
+        __execute(item_cache_.begin(), item_cache_.end());
+        return createIndexedPackGroup(from, to, selector_);
+    }
+
+    template<class TIterator>
+    static IndexedPackGroup createIndexedPackGroup(TIterator from,
+                                                   TIterator to,
+                                                   TSel& selector)
+    {
+        IndexedPackGroup pg;
+        pg.reserve(selector.binCount());
+
+        for(size_t i = 0; i < selector.binCount(); i++) {
+            auto items = selector.itemsForBin(i);
+            pg.push_back({});
+            pg[i].reserve(items.size());
+
+            for(Item& itemA : items) {
+                auto it = from;
+                unsigned idx = 0;
+                while(it != to) {
+                    Item& itemB = *it;
+                    if(&itemB == &itemA) break;
+                    it++; idx++;
+                }
+                pg[i].emplace_back(idx, itemA);
+            }
+        }
+
+        return pg;
+    }
+
+    template<class TIter> inline void __execute(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)));
+        });
+
+        selector_.template packItems<PlacementStrategy>(
+                    from, to, bin_, pconfig_);
+
+        if(min_obj_distance_ > 0) std::for_each(from, to, [](Item& item) {
+            item.removeOffset();
+        });
+    }
+};
+
+}
+
+#endif // LIBNEST2D_HPP

src/libnest2d/include/libnest2d/optimizer.hpp

@@ -0,0 +1,252 @@
+#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::make_tuple;
+
+/// 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 make_tuple(args...); }
+
+/**
+ * @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;
+    tuple<Args...> optimum;
+    double score;
+};
+
+/**
+ * @brief A type for specifying the stop criteria.
+ */
+struct StopCriteria {
+
+    /// If the absolute value difference between two scores.
+    double absolute_score_difference = std::nan("");
+
+    /// If the relative value difference between two scores.
+    double relative_score_difference = std::nan("");
+
+    /// Stop if this value or better is found.
+    double stop_score = std::nan("");
+
+    /// A predicate that if evaluates to true, the optimization should terminate
+    /// and the best result found prior to termination should be returned.
+    std::function<bool()> stop_condition = [] { return false; };
+
+    /// The max allowed number of iterations.
+    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(
+     *           [](tuple<double> x) // object function
+     *           {
+     *               return std::pow(std::get<0>(x), 2);
+     *           },
+     *           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!");
+    }
+
+    DummyOptimizer(const StopCriteria&) {
+        static_assert(always_false<T>::value, "Optimizer unimplemented!");
+    }
+
+    template<class Func, class...Args>
+    Result<Args...> optimize(Func&& /*func*/,
+                             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

src/libnest2d/include/libnest2d/optimizers/nlopt/CMakeLists.txt

@@ -0,0 +1,61 @@
+find_package(NLopt 1.4)
+
+if(NOT NLopt_FOUND)
+    message(STATUS  "NLopt not found so downloading "
+                    "and automatic build is performed...")
+
+    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(URL_NLOPT "https://github.com/stevengj/nlopt.git" 
+        CACHE STRING "Location of the nlopt git repository")
+
+    # set(NLopt_DIR ${CMAKE_BINARY_DIR}/nlopt)
+    include(DownloadProject)
+    download_project(   PROJ                nlopt
+                        GIT_REPOSITORY      ${URL_NLOPT}
+                        GIT_TAG             v2.5.0
+                        # 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} ${nlopt_BINARY_DIR}/src/api)
+    set(SHARED_LIBS_STATE ${SHARED_STATE})
+
+    add_library(NloptOptimizer INTERFACE)
+    target_link_libraries(NloptOptimizer INTERFACE nlopt)
+    target_include_directories(NloptOptimizer INTERFACE ${NLopt_INCLUDE_DIR})
+
+else()
+    add_library(NloptOptimizer INTERFACE)
+    target_link_libraries(NloptOptimizer INTERFACE Nlopt::Nlopt)
+endif()
+
+target_sources( NloptOptimizer INTERFACE
+${CMAKE_CURRENT_SOURCE_DIR}/simplex.hpp
+${CMAKE_CURRENT_SOURCE_DIR}/subplex.hpp
+${CMAKE_CURRENT_SOURCE_DIR}/genetic.hpp
+${CMAKE_CURRENT_SOURCE_DIR}/nlopt_boilerplate.hpp
+)
+
+target_compile_definitions(NloptOptimizer INTERFACE LIBNEST2D_OPTIMIZER_NLOPT)
+
+# And finally plug the NloptOptimizer into libnest2d
+target_link_libraries(libnest2d INTERFACE NloptOptimizer)

src/libnest2d/include/libnest2d/optimizers/nlopt/genetic.hpp

@@ -0,0 +1,32 @@
+#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<>
+inline TOptimizer<Method::G_GENETIC> GlobalOptimizer<Method::G_GENETIC>(
+        Method localm, const StopCriteria& scr )
+{
+    return GeneticOptimizer (scr).localMethod(localm);
+}
+
+}
+}
+
+#endif // GENETIC_HPP

src/libnest2d/include/libnest2d/optimizers/nlopt/nlopt_boilerplate.hpp

@@ -0,0 +1,200 @@
+#ifndef NLOPT_BOILERPLATE_HPP
+#define NLOPT_BOILERPLATE_HPP
+
+#ifdef _MSC_VER
+#pragma warning(push)
+#pragma warning(disable: 4244)
+#pragma warning(disable: 4267)
+#endif
+#include <nlopt.hpp>
+#ifdef _MSC_VER
+#pragma warning(pop)
+#endif
+
+#include <libnest2d/optimizer.hpp>
+#include <cassert>
+#include <libnest2d/utils/metaloop.hpp>
+
+#include <utility>
+
+namespace libnest2d { namespace opt {
+
+inline 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)
+    {
+        using TData = std::pair<remove_ref_t<Fn>*, NloptOptimizer*>;
+        auto typeddata = static_cast<TData*>(data);
+
+        if(typeddata->second->stopcr_.stop_condition())
+            typeddata->second->opt_.force_stop();
+
+        auto fnptr = typeddata->first;
+        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: ;
+        }
+
+        double abs_diff = stopcr_.absolute_score_difference;
+        double rel_diff = stopcr_.relative_score_difference;
+        double stopval = stopcr_.stop_score;
+        if(!std::isnan(abs_diff)) opt_.set_ftol_abs(abs_diff);
+        if(!std::isnan(rel_diff)) opt_.set_ftol_rel(rel_diff);
+        if(!std::isnan(stopval))  opt_.set_stopval(stopval);
+
+        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);
+
+        std::pair<remove_ref_t<Func>*, NloptOptimizer*> data =
+                std::make_pair(&func, this);
+
+        switch(dir_) {
+        case OptDir::MIN:
+            opt_.set_min_objective(optfunc<Func, Args...>, &data); break;
+        case OptDir::MAX:
+            opt_.set_max_objective(optfunc<Func, Args...>, &data); break;
+        }
+
+        Result<Args...> result;
+        nlopt::result rescode;
+
+        try {
+            rescode = opt_.optimize(initvals_, result.score);
+            result.resultcode = static_cast<ResultCodes>(rescode);
+        } catch( nlopt::forced_stop& ) {
+            result.resultcode = ResultCodes::FORCED_STOP;
+        }
+
+        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

src/libnest2d/include/libnest2d/optimizers/nlopt/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

src/libnest2d/include/libnest2d/optimizers/nlopt/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

src/libnest2d/include/libnest2d/optimizers/optimlib/CMakeLists.txt

@@ -0,0 +1,5 @@
+find_package(Armadillo REQUIRED)
+
+add_library(OptimlibOptimizer INTERFACE)
+target_include_directories(OptimlibOptimizer INTERFACE ${ARMADILLO_INCLUDE_DIRS})
+target_link_libraries(OptimlibOptimizer INTERFACE ${ARMADILLO_LIBRARIES})

src/libnest2d/include/libnest2d/placers/bottomleftplacer.hpp

@@ -0,0 +1,415 @@
+#ifndef BOTTOMLEFT_HPP
+#define BOTTOMLEFT_HPP
+
+#include <limits>
+
+#include "placer_boilerplate.hpp"
+
+namespace libnest2d { namespace placers {
+
+template<class T, class = T> struct Epsilon {};
+
+template<class T>
+struct Epsilon<T, enable_if_t<std::is_integral<T>::value, T> > {
+    static const T Value = 1;
+};
+
+template<class T>
+struct Epsilon<T, enable_if_t<std::is_floating_point<T>::value, T> > {
+    static const T Value = 1e-3;
+};
+
+template<class RawShape>
+struct BLConfig {
+    DECLARE_MAIN_TYPES(RawShape);
+
+    Coord min_obj_distance = 0;
+    Coord epsilon = Epsilon<Coord>::Value;
+    bool allow_rotations = false;
+};
+
+template<class RawShape>
+class _BottomLeftPlacer: public PlacerBoilerplate<
+        _BottomLeftPlacer<RawShape>,
+        RawShape, _Box<TPoint<RawShape>>,
+        BLConfig<RawShape> >
+{
+    using Base = PlacerBoilerplate<_BottomLeftPlacer<RawShape>, RawShape,
+                                    _Box<TPoint<RawShape>>, BLConfig<RawShape>>;
+    DECLARE_PLACER(Base)
+
+public:
+
+    explicit _BottomLeftPlacer(const BinType& bin): Base(bin) {}
+
+    template<class Range = ConstItemRange<typename Base::DefaultIter>>
+    PackResult trypack(Item& item,
+                       const Range& = Range())
+    {
+        auto r = _trypack(item);
+        if(!r && Base::config_.allow_rotations) {
+
+            item.rotate(Degrees(90));
+            r =_trypack(item);
+        }
+        return r;
+    }
+
+    enum class Dir {
+        LEFT,
+        DOWN
+    };
+
+    inline RawShape leftPoly(const Item& item) const {
+        return toWallPoly(item, Dir::LEFT);
+    }
+
+    inline RawShape downPoly(const Item& item) const {
+        return toWallPoly(item, Dir::DOWN);
+    }
+
+    inline Unit availableSpaceLeft(const Item& item) {
+        return availableSpace(item, Dir::LEFT);
+    }
+
+    inline Unit availableSpaceDown(const Item& item) {
+        return availableSpace(item, Dir::DOWN);
+    }
+
+protected:
+
+    PackResult _trypack(Item& item) {
+
+        // Get initial position for item in the top right corner
+        setInitialPosition(item);
+
+        Unit d = availableSpaceDown(item);
+        auto eps = config_.epsilon;
+        bool can_move = d > eps;
+        bool can_be_packed = can_move;
+        bool left = true;
+
+        while(can_move) {
+            if(left) { // write previous down move and go down
+                item.translate({0, -d+eps});
+                d = availableSpaceLeft(item);
+                can_move = d > eps;
+                left = false;
+            } else { // write previous left move and go down
+                item.translate({-d+eps, 0});
+                d = availableSpaceDown(item);
+                can_move = d > eps;
+                left = true;
+            }
+        }
+
+        if(can_be_packed) {
+            Item trsh(item.transformedShape());
+            for(auto& v : trsh) can_be_packed = can_be_packed &&
+                    getX(v) < bin_.width() &&
+                    getY(v) < bin_.height();
+        }
+
+        return can_be_packed? PackResult(item) : PackResult();
+    }
+
+    void setInitialPosition(Item& item) {
+        auto bb = item.boundingBox();
+
+        Vertex v = { getX(bb.maxCorner()), getY(bb.minCorner()) };
+
+
+        Coord dx = getX(bin_.maxCorner()) - getX(v);
+        Coord dy = getY(bin_.maxCorner()) - getY(v);
+
+        item.translate({dx, dy});
+    }
+
+    template<class C = Coord>
+    static enable_if_t<std::is_floating_point<C>::value, bool>
+    isInTheWayOf( const Item& item,
+                  const Item& other,
+                  const RawShape& scanpoly)
+    {
+        auto tsh = other.transformedShape();
+        return ( sl::intersects(tsh, scanpoly) ||
+                 sl::isInside(tsh, scanpoly) ) &&
+               ( !sl::intersects(tsh, item.rawShape()) &&
+                 !sl::isInside(tsh, item.rawShape()) );
+    }
+
+    template<class C = Coord>
+    static enable_if_t<std::is_integral<C>::value, bool>
+    isInTheWayOf( const Item& item,
+                  const Item& other,
+                  const RawShape& scanpoly)
+    {
+        auto tsh = other.transformedShape();
+
+        bool inters_scanpoly = sl::intersects(tsh, scanpoly) &&
+                !sl::touches(tsh, scanpoly);
+        bool inters_item = sl::intersects(tsh, item.rawShape()) &&
+                !sl::touches(tsh, item.rawShape());
+
+        return ( inters_scanpoly ||
+                 sl::isInside(tsh, scanpoly)) &&
+               ( !inters_item &&
+                 !sl::isInside(tsh, item.rawShape())
+                 );
+    }
+
+    ItemGroup itemsInTheWayOf(const Item& item, const Dir dir) {
+        // Get the left or down polygon, that has the same area as the shadow
+        // of input item reflected to the left or downwards
+        auto&& scanpoly = dir == Dir::LEFT? leftPoly(item) :
+                                            downPoly(item);
+
+        ItemGroup ret;    // packed items 'in the way' of item
+        ret.reserve(items_.size());
+
+        // Predicate to find items that are 'in the way' for left (down) move
+        auto predicate = [&scanpoly, &item](const Item& it) {
+            return isInTheWayOf(item, it, scanpoly);
+        };
+
+        // Get the items that are in the way for the left (or down) movement
+        std::copy_if(items_.begin(), items_.end(),
+                     std::back_inserter(ret), predicate);
+
+        return ret;
+    }
+
+    Unit availableSpace(const Item& _item, const Dir dir) {
+
+        Item item (_item.transformedShape());
+
+
+        std::function<Coord(const Vertex&)> getCoord;
+        std::function< std::pair<Coord, bool>(const Segment&, const Vertex&) >
+            availableDistanceSV;
+
+        std::function< std::pair<Coord, bool>(const Vertex&, const Segment&) >
+            availableDistance;
+
+        if(dir == Dir::LEFT) {
+            getCoord = [](const Vertex& v) { return getX(v); };
+            availableDistance = pointlike::horizontalDistance<Vertex>;
+            availableDistanceSV = [](const Segment& s, const Vertex& v) {
+                auto ret = pointlike::horizontalDistance<Vertex>(v, s);
+                if(ret.second) ret.first = -ret.first;
+                return ret;
+            };
+        }
+        else {
+            getCoord = [](const Vertex& v) { return getY(v); };
+            availableDistance = pointlike::verticalDistance<Vertex>;
+            availableDistanceSV = [](const Segment& s, const Vertex& v) {
+                auto ret = pointlike::verticalDistance<Vertex>(v, s);
+                if(ret.second) ret.first = -ret.first;
+                return ret;
+            };
+        }
+
+        auto&& items_in_the_way = itemsInTheWayOf(item, dir);
+
+        // Comparison function for finding min vertex
+        auto cmp = [&getCoord](const Vertex& v1, const Vertex& v2) {
+            return getCoord(v1) < getCoord(v2);
+        };
+
+        // find minimum left or down coordinate of item
+        auto minvertex_it = std::min_element(item.begin(),
+                                                   item.end(),
+                                                   cmp);
+
+        // Get the initial distance in floating point
+        Unit m = getCoord(*minvertex_it);
+
+        // Check available distance for every vertex of item to the objects
+        // in the way for the nearest intersection
+        if(!items_in_the_way.empty()) { // This is crazy, should be optimized...
+            for(Item& pleft : items_in_the_way) {
+                // For all segments in items_to_left
+
+                assert(pleft.vertexCount() > 0);
+
+                auto trpleft_poly = pleft.transformedShape();
+                auto& trpleft = sl::contour(trpleft_poly);
+                auto first = sl::begin(trpleft);
+                auto next = first + 1;
+                auto endit = sl::end(trpleft);
+
+                while(next != endit) {
+                    Segment seg(*(first++), *(next++));
+                    for(auto& v : item) {   // For all vertices in item
+
+                        auto d = availableDistance(v, seg);
+
+                        if(d.second && d.first < m)  m = d.first;
+                    }
+                }
+            }
+
+            auto first = item.begin();
+            auto next = first + 1;
+            auto endit = item.end();
+
+            // For all edges in item:
+            while(next != endit) {
+                Segment seg(*(first++), *(next++));
+
+                // for all shapes in items_to_left
+                for(Item& sh : items_in_the_way) {
+                    assert(sh.vertexCount() > 0);
+
+                    Item tsh(sh.transformedShape());
+                    for(auto& v : tsh) {   // For all vertices in item
+
+                        auto d = availableDistanceSV(seg, v);
+
+                        if(d.second && d.first < m)  m = d.first;
+                    }
+                }
+            }
+        }
+
+        return m;
+    }
+
+    /**
+     * Implementation of the left (and down) polygon as described by
+     * [López-Camacho et al. 2013]\
+     * (http://www.cs.stir.ac.uk/~goc/papers/EffectiveHueristic2DAOR2013.pdf)
+     * see algorithm 8 for details...
+     */
+    RawShape toWallPoly(const Item& _item, const Dir dir) const {
+        // The variable names reflect the case of left polygon calculation.
+        //
+        // We will iterate through the item's vertices and search for the top
+        // and bottom vertices (or right and left if dir==Dir::DOWN).
+        // Save the relevant vertices and their indices into `bottom` and
+        // `top` vectors. In case of left polygon construction these will
+        // contain the top and bottom polygons which have the same vertical
+        // coordinates (in case there is more of them).
+        //
+        // We get the leftmost (or downmost) vertex from the `bottom` and `top`
+        // vectors and construct the final polygon.
+
+        Item item (_item.transformedShape());
+
+        auto getCoord = [dir](const Vertex& v) {
+            return dir == Dir::LEFT? getY(v) : getX(v);
+        };
+
+        Coord max_y = std::numeric_limits<Coord>::min();
+        Coord min_y = std::numeric_limits<Coord>::max();
+
+        using El = std::pair<size_t, std::reference_wrapper<const Vertex>>;
+
+        std::function<bool(const El&, const El&)> cmp;
+
+        if(dir == Dir::LEFT)
+            cmp = [](const El& e1, const El& e2) {
+                return getX(e1.second.get()) < getX(e2.second.get());
+            };
+        else
+            cmp = [](const El& e1, const El& e2) {
+                return getY(e1.second.get()) < getY(e2.second.get());
+            };
+
+        std::vector< El > top;
+        std::vector< El > bottom;
+
+        size_t idx = 0;
+        for(auto& v : item) { // Find the bottom and top vertices and save them
+            auto vref = std::cref(v);
+            auto vy = getCoord(v);
+
+            if( vy > max_y ) {
+                max_y = vy;
+                top.clear();
+                top.emplace_back(idx, vref);
+            }
+            else if(vy == max_y) { top.emplace_back(idx, vref); }
+
+            if(vy < min_y) {
+                min_y = vy;
+                bottom.clear();
+                bottom.emplace_back(idx, vref);
+            }
+            else if(vy == min_y) { bottom.emplace_back(idx, vref); }
+
+            idx++;
+        }
+
+        // Get the top and bottom leftmost vertices, or the right and left
+        // downmost vertices (if dir == Dir::DOWN)
+        auto topleft_it = std::min_element(top.begin(), top.end(), cmp);
+        auto bottomleft_it =
+                std::min_element(bottom.begin(), bottom.end(), cmp);
+
+        auto& topleft_vertex = topleft_it->second.get();
+        auto& bottomleft_vertex = bottomleft_it->second.get();
+
+        // Start and finish positions for the vertices that will be part of the
+        // new polygon
+        auto start = std::min(topleft_it->first, bottomleft_it->first);
+        auto finish = std::max(topleft_it->first, bottomleft_it->first);
+
+        RawShape ret;
+
+        // the return shape
+        auto& rsh = sl::contour(ret);
+
+        // reserve for all vertices plus 2 for the left horizontal wall, 2 for
+        // the additional vertices for maintaning min object distance
+        sl::reserve(rsh, finish-start+4);
+
+        /*auto addOthers = [&rsh, finish, start, &item](){
+            for(size_t i = start+1; i < finish; i++)
+                sl::addVertex(rsh, item.vertex(i));
+        };*/
+
+        auto reverseAddOthers = [&rsh, finish, start, &item](){
+            for(auto i = finish-1; i > start; i--)
+                sl::addVertex(rsh, item.vertex(
+                                         static_cast<unsigned long>(i)));
+        };
+
+        // Final polygon construction...
+
+        static_assert(OrientationType<RawShape>::Value ==
+                      Orientation::CLOCKWISE,
+                      "Counter clockwise toWallPoly() Unimplemented!");
+
+        // Clockwise polygon construction
+
+        sl::addVertex(rsh, topleft_vertex);
+
+        if(dir == Dir::LEFT) reverseAddOthers();
+        else {
+            sl::addVertex(rsh, getX(topleft_vertex), 0);
+            sl::addVertex(rsh, getX(bottomleft_vertex), 0);
+        }
+
+        sl::addVertex(rsh, bottomleft_vertex);
+
+        if(dir == Dir::LEFT) {
+            sl::addVertex(rsh, 0, getY(bottomleft_vertex));
+            sl::addVertex(rsh, 0, getY(topleft_vertex));
+        }
+        else reverseAddOthers();
+
+
+        // Close the polygon
+        sl::addVertex(rsh, topleft_vertex);
+
+        return ret;
+    }
+
+};
+
+}
+}
+
+#endif //BOTTOMLEFT_HPP

src/libnest2d/include/libnest2d/placers/placer_boilerplate.hpp

@@ -0,0 +1,134 @@
+#ifndef PLACER_BOILERPLATE_HPP
+#define PLACER_BOILERPLATE_HPP
+
+#include <libnest2d/libnest2d.hpp>
+
+namespace libnest2d { namespace placers {
+
+struct EmptyConfig {};
+
+template<class Subclass, class RawShape, class TBin, class Cfg = EmptyConfig>
+class PlacerBoilerplate {
+    mutable bool farea_valid_ = false;
+    mutable double farea_ = 0.0;
+public:
+    using Item = _Item<RawShape>;
+    using Vertex = TPoint<RawShape>;
+    using Segment = _Segment<Vertex>;
+    using BinType = TBin;
+    using Coord = TCoord<Vertex>;
+    using Unit = Coord;
+    using Config = Cfg;
+    using ItemGroup = _ItemGroup<Item>;
+    using DefaultIter = typename ItemGroup::const_iterator;
+
+    class PackResult {
+        Item *item_ptr_;
+        Vertex move_;
+        Radians rot_;
+        double overfit_;
+        friend class PlacerBoilerplate;
+        friend Subclass;
+
+        PackResult(Item& item):
+            item_ptr_(&item),
+            move_(item.translation()),
+            rot_(item.rotation()) {}
+
+        PackResult(double overfit = 1.0):
+            item_ptr_(nullptr), overfit_(overfit) {}
+
+    public:
+        operator bool() { return item_ptr_ != nullptr; }
+        double overfit() const { return overfit_; }
+    };
+
+    inline PlacerBoilerplate(const BinType& bin, unsigned cap = 50): bin_(bin)
+    {
+        items_.reserve(cap);
+    }
+
+    inline const BinType& bin() const BP2D_NOEXCEPT { return bin_; }
+
+    template<class TB> inline void bin(TB&& b) {
+        bin_ = std::forward<BinType>(b);
+    }
+
+    inline void configure(const Config& config) BP2D_NOEXCEPT {
+        config_ = config;
+    }
+
+    template<class Range = ConstItemRange<DefaultIter>>
+    bool pack(Item& item,
+              const Range& rem = Range()) {
+        auto&& r = static_cast<Subclass*>(this)->trypack(item, rem);
+        if(r) {
+            items_.push_back(*(r.item_ptr_));
+            farea_valid_ = false;
+        }
+        return r;
+    }
+
+    void accept(PackResult& r) {
+        if(r) {
+            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();
+        farea_valid_ = false;
+    }
+
+    inline const ItemGroup& getItems() const { return items_; }
+
+    inline void clearItems() {
+        items_.clear();
+        farea_valid_ = false;
+    }
+
+    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_;
+    }
+
+protected:
+
+    BinType bin_;
+    ItemGroup items_;
+    Cfg config_;
+};
+
+
+#define DECLARE_PLACER(Base) \
+using Base::bin_;                 \
+using Base::items_;               \
+using Base::config_;              \
+public:                           \
+using typename Base::Item;        \
+using typename Base::ItemGroup;   \
+using typename Base::BinType;     \
+using typename Base::Config;      \
+using typename Base::Vertex;      \
+using typename Base::Segment;     \
+using typename Base::PackResult;  \
+using typename Base::Coord;       \
+using typename Base::Unit;        \
+private:
+
+}
+}
+
+#endif // PLACER_BOILERPLATE_HPP

src/libnest2d/include/libnest2d/selections/djd_heuristic.hpp

@@ -0,0 +1,719 @@
+#ifndef DJD_HEURISTIC_HPP
+#define DJD_HEURISTIC_HPP
+
+#include <list>
+#include <future>
+#include <atomic>
+#include <functional>
+
+#include "selection_boilerplate.hpp"
+
+namespace libnest2d { namespace selections {
+
+/**
+ * Selection heuristic based on [López-Camacho]\
+ * (http://www.cs.stir.ac.uk/~goc/papers/EffectiveHueristic2DAOR2013.pdf)
+ */
+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;
+
+    /**
+     * @brief The Config for DJD heuristic.
+     */
+    struct Config {
+
+        /**
+         * If true, the algorithm will try to place pair and triplets in all
+         * possible order. It will have a hugely negative impact on performance.
+         */
+        bool try_reverse_order = true;
+
+        /**
+         * @brief try_pairs Whether to try pairs of items to pack. It will add
+         * a quadratic component to the complexity.
+         */
+        bool try_pairs = true;
+
+        /**
+         * @brief Whether to try groups of 3 items to pack. This could be very
+         * slow for large number of items (>100) as it adds a cubic component
+         * to the complexity.
+         */
+        bool try_triplets = false;
+
+        /**
+         * 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;
+    Container store_;
+    Config config_;
+
+    static const unsigned MAX_ITEMS_SEQUENTIALLY = 30;
+    static const unsigned MAX_VERTICES_SEQUENTIALLY = MAX_ITEMS_SEQUENTIALLY*20;
+
+public:
+
+    inline void configure(const Config& config) {
+        config_ = config;
+    }
+
+    template<class TPlacer, class TIterator,
+             class TBin = typename PlacementStrategyLike<TPlacer>::BinType,
+             class PConfig = typename PlacementStrategyLike<TPlacer>::Config>
+    void packItems( TIterator first,
+                    TIterator last,
+                    const TBin& bin,
+                    PConfig&& pconfig = PConfig() )
+    {
+        using Placer = PlacementStrategyLike<TPlacer>;
+        using ItemList = std::list<ItemRef>;
+
+        const double bin_area = sl::area(bin);
+        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();
+        store_.reserve(last-first);
+        packed_bins_.clear();
+
+        std::copy(first, last, std::back_inserter(store_));
+
+        std::sort(store_.begin(), store_.end(), [](Item& i1, Item& i2) {
+            return i1.area() > i2.area();
+        });
+
+        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;
+
+        bool try_reverse = config_.try_reverse_order;
+
+        // Will use a subroutine to add a new bin
+        auto addBin = [this, &placers, &bin, &pconfig]()
+        {
+            placers.emplace_back(bin);
+            packed_bins_.emplace_back();
+            placers.back().configure(pconfig);
+        };
+
+        // Types for pairs and triplets
+        using TPair = std::tuple<ItemRef, ItemRef>;
+        using TTriplet = std::tuple<ItemRef, ItemRef, ItemRef>;
+
+
+        // Method for checking a pair whether it was a pack failure.
+        auto check_pair = [](const std::vector<TPair>& wrong_pairs,
+                ItemRef i1, ItemRef i2)
+        {
+            return std::any_of(wrong_pairs.begin(), wrong_pairs.end(),
+                               [&i1, &i2](const TPair& pair)
+            {
+                Item& pi1 = std::get<0>(pair), &pi2 = std::get<1>(pair);
+                Item& ri1 = i1, &ri2 = i2;
+                return (&pi1 == &ri1 && &pi2 == &ri2) ||
+                       (&pi1 == &ri2 && &pi2 == &ri1);
+            });
+        };
+
+        // Method for checking if a triplet was a pack failure
+        auto check_triplet = [](
+                const std::vector<TTriplet>& wrong_triplets,
+                ItemRef i1,
+                ItemRef i2,
+                ItemRef i3)
+        {
+            return std::any_of(wrong_triplets.begin(),
+                               wrong_triplets.end(),
+                               [&i1, &i2, &i3](const TTriplet& tripl)
+            {
+                Item& pi1 = std::get<0>(tripl);
+                Item& pi2 = std::get<1>(tripl);
+                Item& pi3 = std::get<2>(tripl);
+                Item& ri1 = i1, &ri2 = i2, &ri3 = i3;
+                return  (&pi1 == &ri1 && &pi2 == &ri2 && &pi3 == &ri3) ||
+                        (&pi1 == &ri1 && &pi2 == &ri3 && &pi3 == &ri2) ||
+                        (&pi1 == &ri2 && &pi2 == &ri1 && &pi3 == &ri3) ||
+                        (&pi1 == &ri3 && &pi2 == &ri2 && &pi3 == &ri1);
+            });
+        };
+
+        using ItemListIt = typename ItemList::iterator;
+
+        auto largestPiece = [](ItemListIt it, ItemList& not_packed) {
+            return it == not_packed.begin()? std::next(it) : not_packed.begin();
+        };
+
+        auto secondLargestPiece = [&largestPiece](ItemListIt it,
+                ItemList& not_packed) {
+            auto ret = std::next(largestPiece(it, not_packed));
+            return ret == it? std::next(ret) : ret;
+        };
+
+        auto smallestPiece = [](ItemListIt it, ItemList& not_packed) {
+            auto last = std::prev(not_packed.end());
+            return it == last? std::prev(it) : last;
+        };
+
+        auto secondSmallestPiece = [&smallestPiece](ItemListIt it,
+                ItemList& not_packed) {
+            auto ret = std::prev(smallestPiece(it, not_packed));
+            return ret == it? std::prev(ret) : ret;
+        };
+
+        auto tryOneByOne = // Subroutine to try adding items one by one.
+                [&bin_area]
+                (Placer& placer, ItemList& not_packed,
+                double waste,
+                double& free_area,
+                double& filled_area)
+        {
+            double item_area = 0;
+            bool ret = false;
+            auto it = not_packed.begin();
+
+            auto pack = [&placer, &not_packed](ItemListIt it) {
+                return placer.pack(*it, rem(it, not_packed));
+            };
+
+            while(it != not_packed.end() && !ret &&
+                  free_area - (item_area = it->get().area()) <= waste)
+            {
+                if(item_area <= free_area && pack(it) ) {
+                    free_area -= item_area;
+                    filled_area = bin_area - free_area;
+                    ret = true;
+                } else
+                    it++;
+            }
+
+            if(ret) not_packed.erase(it);
+
+            return ret;
+        };
+
+        auto tryGroupsOfTwo = // Try adding groups of two items into the bin.
+                [&bin_area, &check_pair, &largestPiece, &smallestPiece,
+                 try_reverse]
+                (Placer& placer, ItemList& not_packed,
+                double waste,
+                double& free_area,
+                double& filled_area)
+        {
+            double item_area = 0;
+            const auto endit = not_packed.end();
+
+            if(not_packed.size() < 2)
+                return false; // No group of two items
+
+            double largest_area = not_packed.front().get().area();
+            auto itmp = not_packed.begin(); itmp++;
+            double second_largest = itmp->get().area();
+            if( free_area - second_largest - largest_area > waste)
+                return false; // If even the largest two items do not fill
+                // the bin to the desired waste than we can end here.
+
+
+            bool ret = false;
+            auto it = not_packed.begin();
+            auto it2 = it;
+
+            std::vector<TPair> wrong_pairs;
+            using std::placeholders::_1;
+
+            auto trypack = [&placer, &not_packed](ItemListIt it) {
+                return placer.trypack(*it, rem(it, not_packed));
+            };
+
+            while(it != endit && !ret &&
+                  free_area - (item_area = it->get().area()) -
+                  largestPiece(it, not_packed)->get().area() <= waste)
+            {
+                if(item_area + smallestPiece(it, not_packed)->get().area() >
+                        free_area ) { it++; continue; }
+
+                auto pr = trypack(it);
+
+                // First would fit
+                it2 = not_packed.begin();
+                double item2_area = 0;
+                while(it2 != endit && pr && !ret && free_area -
+                      (item2_area = it2->get().area()) - item_area <= waste)
+                {
+                    double area_sum = item_area + item2_area;
+
+                    if(it == it2 || area_sum > free_area ||
+                            check_pair(wrong_pairs, *it, *it2)) {
+                        it2++; continue;
+                    }
+
+                    placer.accept(pr);
+                    auto pr2 = trypack(it2);
+                    if(!pr2) {
+                        placer.unpackLast(); // remove first
+                        if(try_reverse) {
+                            pr2 = trypack(it2);
+                            if(pr2) {
+                                placer.accept(pr2);
+                                auto pr12 = trypack(it);
+                                if(pr12) {
+                                    placer.accept(pr12);
+                                    ret = true;
+                                } else {
+                                    placer.unpackLast();
+                                }
+                            }
+                        }
+                    } else {
+                        placer.accept(pr2); ret = true;
+                    }
+
+                    if(ret)
+                    { // Second fits as well
+                        free_area -= area_sum;
+                        filled_area = bin_area - free_area;
+                    } else {
+                        wrong_pairs.emplace_back(*it, *it2);
+                        it2++;
+                    }
+                }
+
+                if(!ret) it++;
+            }
+
+            if(ret) { not_packed.erase(it); not_packed.erase(it2); }
+
+            return ret;
+        };
+
+        auto tryGroupsOfThree = // Try adding groups of three items.
+                [&bin_area,
+                 &smallestPiece, &largestPiece,
+                 &secondSmallestPiece, &secondLargestPiece,
+                 &check_pair, &check_triplet, try_reverse]
+                (Placer& placer, ItemList& not_packed,
+                double waste,
+                double& free_area,
+                double& filled_area)
+        {
+            auto np_size = not_packed.size();
+            if(np_size < 3) return false;
+
+            auto it = not_packed.begin();           // from
+            const auto endit = not_packed.end();    // to
+            auto it2 = it, it3 = it;
+
+            // Containers for pairs and triplets that were tried before and
+            // do not work.
+            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;
+
+            auto area = [](const ItemListIt& it) {
+                return it->get().area();
+            };
+
+            auto trypack = [&placer, &not_packed](ItemListIt it) {
+                return placer.trypack(*it, rem(it, not_packed));
+            };
+
+            auto pack = [&placer, &not_packed](ItemListIt it) {
+                return placer.pack(*it, rem(it, not_packed));
+            };
+
+            while (it != endit && !ret) { // drill down 1st level
+
+                // We need to determine in each iteration the largest, second
+                // largest, smallest and second smallest item in terms of area.
+
+                Item& largest = *largestPiece(it, not_packed);
+                Item& second_largest = *secondLargestPiece(it, not_packed);
+
+                double area_of_two_largest =
+                        largest.area() + second_largest.area();
+
+                // Check if there is enough free area for the item and the two
+                // largest item
+                if(free_area - area(it) - area_of_two_largest > waste)
+                    break;
+
+                // Determine the area of the two smallest item.
+                Item& smallest = *smallestPiece(it, not_packed);
+                Item& second_smallest = *secondSmallestPiece(it, not_packed);
+
+                // Check if there is enough free area for the item and the two
+                // smallest item.
+                double area_of_two_smallest =
+                        smallest.area() + second_smallest.area();
+
+                if(area(it) + area_of_two_smallest > free_area) {
+                    it++; continue;
+                }
+
+                auto pr = trypack(it);
+
+                // Check for free area and try to pack the 1st item...
+                if(!pr) { it++; continue; }
+
+                it2 = not_packed.begin();
+                double rem2_area = free_area - largest.area();
+                double a2_sum = 0;
+
+                while(it2 != endit && !ret &&
+                      rem2_area - (a2_sum = area(it) + area(it2)) <= waste) {
+                    // Drill down level 2
+
+                    if(a2_sum != area(it) + area(it2)) throw -1;
+
+                    if(it == it2 || check_pair(wrong_pairs, *it, *it2)) {
+                        it2++; continue;
+                    }
+
+                    if(a2_sum + smallest.area() > free_area) {
+                        it2++; continue;
+                    }
+
+                    bool can_pack2 = false;
+
+                    placer.accept(pr);
+                    auto pr2 = trypack(it2);
+                    auto pr12 = pr;
+                    if(!pr2) {
+                        placer.unpackLast(); // remove first
+                        if(try_reverse) {
+                            pr2 = trypack(it2);
+                            if(pr2) {
+                                placer.accept(pr2);
+                                pr12 = trypack(it);
+                                if(pr12) can_pack2 = true;
+                                placer.unpackLast();
+                            }
+                        }
+                    } else {
+                        placer.unpackLast();
+                        can_pack2 = true;
+                    }
+
+                    if(!can_pack2) {
+                        wrong_pairs.emplace_back(*it, *it2);
+                        it2++;
+                        continue;
+                    }
+
+                    // Now we have packed a group of 2 items.
+                    // The 'smallest' variable now could be identical with
+                    // it2 but we don't bother with that
+
+                    it3 = not_packed.begin();
+
+                    double a3_sum = 0;
+
+                    while(it3 != endit && !ret &&
+                          free_area - (a3_sum = a2_sum + area(it3)) <= waste) {
+                        // 3rd level
+
+                        if(it3 == it || it3 == it2 ||
+                                check_triplet(wrong_triplets, *it, *it2, *it3))
+                        { it3++; continue; }
+
+                        if(a3_sum > free_area) { it3++; continue; }
+
+                        placer.accept(pr12); placer.accept(pr2);
+                        bool can_pack3 = pack(it3);
+
+                        if(!can_pack3) {
+                            placer.unpackLast();
+                            placer.unpackLast();
+                        }
+
+                        if(!can_pack3 && try_reverse) {
+
+                            std::array<size_t, 3> indices = {0, 1, 2};
+                            std::array<typename ItemList::iterator, 3>
+                                    candidates = {it, it2, it3};
+
+                            auto tryPack = [&placer, &candidates, &pack](
+                                    const decltype(indices)& idx)
+                            {
+                                std::array<bool, 3> packed = {false};
+
+                                for(auto id : idx) packed.at(id) =
+                                        pack(candidates[id]);
+
+                                bool check =
+                                std::all_of(packed.begin(),
+                                            packed.end(),
+                                            [](bool b) { return b; });
+
+                                if(!check) for(bool b : packed) if(b)
+                                        placer.unpackLast();
+
+                                return check;
+                            };
+
+                            while (!can_pack3 && std::next_permutation(
+                                       indices.begin(),
+                                       indices.end())){
+                                can_pack3 = tryPack(indices);
+                            };
+                        }
+
+                        if(can_pack3) {
+                            // finishit
+                            free_area -= a3_sum;
+                            filled_area = bin_area - free_area;
+                            ret = true;
+                        } else {
+                            wrong_triplets.emplace_back(*it, *it2, *it3);
+                            it3++;
+                        }
+
+                    } // 3rd while
+
+                    if(!ret) it2++;
+
+                } // Second while
+
+                if(!ret) it++;
+
+            } // First while
+
+            if(ret) { // If we eventually succeeded, remove all the packed ones.
+                not_packed.erase(it);
+                not_packed.erase(it2);
+                not_packed.erase(it3);
+            }
+
+            return ret;
+        };
+
+        // 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 = store_.begin();
+            while (it != store_.end() && !this->stopcond_()) {
+                Placer p(bin); p.configure(pconfig);
+                if(!p.pack(*it, rem(it, store_))) {
+                    it = store_.erase(it);
+                } else it++;
+            }
+        }
+
+        int acounter = int(store_.size());
+        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();
+
+            // TODO here should be a spinlock
+            slock.lock();
+            acounter -= packednum;
+            this->progress_(acounter);
+            slock.unlock();
+        };
+
+        double items_area = 0;
+        for(Item& item : store_) items_area += item.area();
+
+        // 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";
+
+        bool do_pairs = config_.try_pairs;
+        bool do_triplets = config_.try_triplets;
+        StopCondition stopcond = this->stopcond_;
+
+        // The DJD heuristic algorithm itself:
+        auto packjob = [INITIAL_FILL_AREA, bin_area, w, do_triplets, do_pairs,
+                        stopcond,
+                        &tryOneByOne,
+                        &tryGroupsOfTwo,
+                        &tryGroupsOfThree,
+                        &makeProgress]
+                        (Placer& placer, ItemList& not_packed, size_t idx)
+        {
+            double filled_area = placer.filledArea();
+            double free_area = bin_area - filled_area;
+            double waste = .0;
+            bool lasttry = false;
+
+            while(!not_packed.empty() && !stopcond()) {
+
+                {// Fill the bin up to INITIAL_FILL_PROPORTION of its capacity
+                    auto it = not_packed.begin();
+
+                    while(it != not_packed.end() && !stopcond() &&
+                          filled_area < INITIAL_FILL_AREA)
+                    {
+                        if(placer.pack(*it, rem(it, not_packed))) {
+                            filled_area += it->get().area();
+                            free_area = bin_area - filled_area;
+                            it = not_packed.erase(it);
+                            makeProgress(placer, idx, 1);
+                        } else it++;
+                    }
+                }
+
+                // try pieces one by one
+                while(tryOneByOne(placer, not_packed, waste, free_area,
+                                  filled_area)) {
+                    waste = 0; lasttry = false;
+                    makeProgress(placer, idx, 1);
+                }
+
+                // try groups of 2 pieces
+                while(do_pairs &&
+                      tryGroupsOfTwo(placer, not_packed, waste, free_area,
+                                     filled_area)) {
+                    waste = 0; lasttry = false;
+                    makeProgress(placer, idx, 2);
+                }
+
+                // try groups of 3 pieces
+                while(do_triplets &&
+                      tryGroupsOfThree(placer, not_packed, waste, free_area,
+                                       filled_area)) {
+                    waste = 0; lasttry = false;
+                    makeProgress(placer, idx, 3);
+                }
+
+                waste += w;
+                if(!lasttry && waste > free_area) lasttry = true;
+                else if(lasttry) break;
+            }
+        };
+
+        size_t idx = 0;
+        ItemList remaining;
+
+        if(do_parallel) {
+            std::vector<ItemList> not_packeds(bincount_guess);
+
+            // Preallocating the bins
+            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<void>> 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++;
+        }
+
+    }
+};
+
+}
+}
+
+#endif // DJD_HEURISTIC_HPP

src/libnest2d/include/libnest2d/selections/filler.hpp

@@ -0,0 +1,80 @@
+#ifndef FILLER_HPP
+#define FILLER_HPP
+
+#include "selection_boilerplate.hpp"
+
+namespace libnest2d { namespace selections {
+
+template<class RawShape>
+class _FillerSelection: public SelectionBoilerplate<RawShape> {
+    using Base = SelectionBoilerplate<RawShape>;
+public:
+    using typename Base::Item;
+    using Config = int; //dummy
+
+private:
+    using Base::packed_bins_;
+    using typename Base::ItemGroup;
+    using Container = ItemGroup;
+    Container store_;
+
+public:
+
+    void configure(const Config& /*config*/) { }
+
+    template<class TPlacer, class TIterator,
+             class TBin = typename PlacementStrategyLike<TPlacer>::BinType,
+             class PConfig = typename PlacementStrategyLike<TPlacer>::Config>
+    void packItems(TIterator first,
+                   TIterator last,
+                   TBin&& bin,
+                   PConfig&& pconfig = PConfig())
+    {
+
+        store_.clear();
+        auto total = last-first;
+        store_.reserve(total);
+        packed_bins_.emplace_back();
+
+        auto makeProgress = [this, &total](
+                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_(--total);
+        };
+
+        std::copy(first, last, std::back_inserter(store_));
+
+        auto sortfunc = [](Item& i1, Item& i2) {
+            return i1.area() > i2.area();
+        };
+
+        std::sort(store_.begin(), store_.end(), sortfunc);
+
+        PlacementStrategyLike<TPlacer> placer(bin);
+        placer.configure(pconfig);
+
+        auto it = store_.begin();
+        while(it != store_.end() && !this->stopcond_()) {
+            if(!placer.pack(*it, {std::next(it), store_.end()}))  {
+                if(packed_bins_.back().empty()) ++it;
+                placer.clearItems();
+                packed_bins_.emplace_back();
+            } else {
+                makeProgress(placer);
+                ++it;
+            }
+        }
+
+    }
+};
+
+}
+}
+
+#endif //BOTTOMLEFT_HPP

src/libnest2d/include/libnest2d/selections/firstfit.hpp

@@ -0,0 +1,99 @@
+#ifndef FIRSTFIT_HPP
+#define FIRSTFIT_HPP
+
+#include "selection_boilerplate.hpp"
+
+namespace libnest2d { namespace selections {
+
+template<class RawShape>
+class _FirstFitSelection: public SelectionBoilerplate<RawShape> {
+    using Base = SelectionBoilerplate<RawShape>;
+public:
+    using typename Base::Item;
+    using Config = int; //dummy
+
+private:
+    using Base::packed_bins_;
+    using typename Base::ItemGroup;
+    using Container = ItemGroup;//typename std::vector<_Item<RawShape>>;
+
+    Container store_;
+
+public:
+
+    void configure(const Config& /*config*/) { }
+
+    template<class TPlacer, class TIterator,
+             class TBin = typename PlacementStrategyLike<TPlacer>::BinType,
+             class PConfig = typename PlacementStrategyLike<TPlacer>::Config>
+    void packItems(TIterator first,
+                   TIterator last,
+                   TBin&& bin,
+                   PConfig&& pconfig = PConfig())
+    {
+
+        using Placer = PlacementStrategyLike<TPlacer>;
+
+        store_.clear();
+        store_.reserve(last-first);
+        packed_bins_.clear();
+
+        std::vector<Placer> placers;
+        placers.reserve(last-first);
+
+        std::copy(first, last, std::back_inserter(store_));
+
+        auto sortfunc = [](Item& i1, Item& i2) {
+            return i1.area() > i2.area();
+        };
+
+        std::sort(store_.begin(), store_.end(), sortfunc);
+
+        auto total = last-first;
+        auto makeProgress = [this, &total](Placer& placer, size_t idx) {
+            packed_bins_[idx] = placer.getItems();
+            this->progress_(static_cast<unsigned>(--total));
+        };
+
+        auto& cancelled = this->stopcond_;
+
+        // Safety test: try to pack each item into an empty bin. If it fails
+        // then it should be removed from the list
+        { auto it = store_.begin();
+            while (it != store_.end() && !cancelled()) {
+                Placer p(bin); p.configure(pconfig);
+                if(!p.pack(*it)) {
+                    it = store_.erase(it);
+                } else it++;
+            }
+        }
+
+
+        auto it = store_.begin();
+
+        while(it != store_.end() && !cancelled()) {
+            bool was_packed = false;
+            size_t j = 0;
+            while(!was_packed && !cancelled()) {
+                for(; j < placers.size() && !was_packed && !cancelled(); j++) {
+                    if((was_packed = placers[j].pack(*it, rem(it, store_) )))
+                            makeProgress(placers[j], j);
+                }
+
+                if(!was_packed) {
+                    placers.emplace_back(bin);
+                    placers.back().configure(pconfig);
+                    packed_bins_.emplace_back();
+                    j = placers.size() - 1;
+                }
+            }
+            ++it;
+        }
+    }
+
+};
+
+}
+}
+
+#endif // FIRSTFIT_HPP

src/libnest2d/include/libnest2d/selections/selection_boilerplate.hpp

@@ -0,0 +1,43 @@
+#ifndef SELECTION_BOILERPLATE_HPP
+#define SELECTION_BOILERPLATE_HPP
+
+#include <atomic>
+#include <libnest2d/libnest2d.hpp>
+
+namespace libnest2d { namespace selections {
+
+template<class RawShape>
+class SelectionBoilerplate {
+public:
+    using Item = _Item<RawShape>;
+    using ItemRef = std::reference_wrapper<Item>;
+    using ItemGroup = std::vector<ItemRef>;
+    using PackGroup = std::vector<ItemGroup>;
+
+    size_t binCount() const { return packed_bins_.size(); }
+
+    ItemGroup itemsForBin(size_t binIndex) {
+        assert(binIndex < packed_bins_.size());
+        return packed_bins_[binIndex];
+    }
+
+    inline const ItemGroup itemsForBin(size_t binIndex) const {
+        assert(binIndex < packed_bins_.size());
+        return packed_bins_[binIndex];
+    }
+
+    inline void progressIndicator(ProgressFunction fn) { progress_ = fn; }
+
+    inline void stopCondition(StopCondition cond) { stopcond_ = cond; }
+
+protected:
+
+    PackGroup packed_bins_;
+    ProgressFunction progress_ = [](unsigned){};
+    StopCondition stopcond_ = [](){ return false; };
+};
+
+}
+}
+
+#endif // SELECTION_BOILERPLATE_HPP

src/libnest2d/include/libnest2d/utils/boost_alg.hpp

@@ -0,0 +1,526 @@
+#ifndef BOOST_ALG_HPP
+#define BOOST_ALG_HPP
+
+#ifndef DISABLE_BOOST_SERIALIZE
+    #include <sstream>
+#endif
+
+#ifdef __clang__
+#undef _MSC_EXTENSIONS
+#endif
+
+#ifdef _MSC_VER
+#pragma warning(push)
+#pragma warning(disable: 4244)
+#pragma warning(disable: 4267)
+#endif
+#include <boost/geometry.hpp>
+#ifdef _MSC_VER
+#pragma warning(pop)
+#endif
+// this should be removed to not confuse the compiler
+// #include <libnest2d.h>
+
+namespace bp2d {
+
+using libnest2d::TCoord;
+using libnest2d::PointImpl;
+using Coord = TCoord<PointImpl>;
+using libnest2d::PolygonImpl;
+using libnest2d::PathImpl;
+using libnest2d::Orientation;
+using libnest2d::OrientationType;
+using libnest2d::getX;
+using libnest2d::getY;
+using libnest2d::setX;
+using libnest2d::setY;
+using Box = libnest2d::_Box<PointImpl>;
+using Segment = libnest2d::_Segment<PointImpl>;
+using Shapes = libnest2d::nfp::Shapes<PolygonImpl>;
+
+}
+
+/**
+ * We have to make all the libnest2d geometry types available to boost. The real
+ * models of the geometries remain the same if a conforming model for libnest2d
+ * was defined by the library client. Boost is used only as an optional
+ * implementer of some algorithms that can be implemented by the model itself
+ * if a faster alternative exists.
+ *
+ * However, boost has its own type traits and we have to define the needed
+ * specializations to be able to use boost::geometry. This can be done with the
+ * already provided model.
+ */
+namespace boost {
+namespace geometry {
+namespace traits {
+
+/* ************************************************************************** */
+/* Point concept adaptaion ************************************************** */
+/* ************************************************************************** */
+
+template<> struct tag<bp2d::PointImpl> {
+    using type = point_tag;
+};
+
+template<> struct coordinate_type<bp2d::PointImpl> {
+    using type = bp2d::Coord;
+};
+
+template<> struct coordinate_system<bp2d::PointImpl> {
+    using type = cs::cartesian;
+};
+
+template<> struct dimension<bp2d::PointImpl>: boost::mpl::int_<2> {};
+
+template<>
+struct access<bp2d::PointImpl, 0 > {
+    static inline bp2d::Coord get(bp2d::PointImpl const& a) {
+        return libnest2d::getX(a);
+    }
+
+    static inline void set(bp2d::PointImpl& a,
+                           bp2d::Coord const& value) {
+        libnest2d::setX(a, value);
+    }
+};
+
+template<>
+struct access<bp2d::PointImpl, 1 > {
+    static inline bp2d::Coord get(bp2d::PointImpl const& a) {
+        return libnest2d::getY(a);
+    }
+
+    static inline void set(bp2d::PointImpl& a,
+                           bp2d::Coord const& value) {
+        libnest2d::setY(a, value);
+    }
+};
+
+
+/* ************************************************************************** */
+/* Box concept adaptaion **************************************************** */
+/* ************************************************************************** */
+
+template<> struct tag<bp2d::Box> {
+    using type = box_tag;
+};
+
+template<> struct point_type<bp2d::Box> {
+    using type = bp2d::PointImpl;
+};
+
+template<> struct indexed_access<bp2d::Box, min_corner, 0> {
+    static inline bp2d::Coord get(bp2d::Box const& box) {
+        return bp2d::getX(box.minCorner());
+    }
+    static inline void set(bp2d::Box &box, bp2d::Coord const& coord) {
+        bp2d::setX(box.minCorner(), coord);
+    }
+};
+
+template<> struct indexed_access<bp2d::Box, min_corner, 1> {
+    static inline bp2d::Coord get(bp2d::Box const& box) {
+        return bp2d::getY(box.minCorner());
+    }
+    static inline void set(bp2d::Box &box, bp2d::Coord const& coord) {
+        bp2d::setY(box.minCorner(), coord);
+    }
+};
+
+template<> struct indexed_access<bp2d::Box, max_corner, 0> {
+    static inline bp2d::Coord get(bp2d::Box const& box) {
+        return bp2d::getX(box.maxCorner());
+    }
+    static inline void set(bp2d::Box &box, bp2d::Coord const& coord) {
+        bp2d::setX(box.maxCorner(), coord);
+    }
+};
+
+template<> struct indexed_access<bp2d::Box, max_corner, 1> {
+    static inline bp2d::Coord get(bp2d::Box const& box) {
+        return bp2d::getY(box.maxCorner());
+    }
+    static inline void set(bp2d::Box &box, bp2d::Coord const& coord) {
+        bp2d::setY(box.maxCorner(), coord);
+    }
+};
+
+/* ************************************************************************** */
+/* Segment concept adaptaion ************************************************ */
+/* ************************************************************************** */
+
+template<> struct tag<bp2d::Segment> {
+    using type = segment_tag;
+};
+
+template<> struct point_type<bp2d::Segment> {
+    using type = bp2d::PointImpl;
+};
+
+template<> struct indexed_access<bp2d::Segment, 0, 0> {
+    static inline bp2d::Coord get(bp2d::Segment const& seg) {
+        return bp2d::getX(seg.first());
+    }
+    static inline void set(bp2d::Segment &seg, bp2d::Coord const& coord) {
+        auto p = seg.first(); bp2d::setX(p, coord); seg.first(p);
+    }
+};
+
+template<> struct indexed_access<bp2d::Segment, 0, 1> {
+    static inline bp2d::Coord get(bp2d::Segment const& seg) {
+        return bp2d::getY(seg.first());
+    }
+    static inline void set(bp2d::Segment &seg, bp2d::Coord const& coord) {
+        auto p = seg.first(); bp2d::setY(p, coord); seg.first(p);
+    }
+};
+
+template<> struct indexed_access<bp2d::Segment, 1, 0> {
+    static inline bp2d::Coord get(bp2d::Segment const& seg) {
+        return bp2d::getX(seg.second());
+    }
+    static inline void set(bp2d::Segment &seg, bp2d::Coord const& coord) {
+        auto p = seg.second();  bp2d::setX(p, coord); seg.second(p);
+    }
+};
+
+template<> struct indexed_access<bp2d::Segment, 1, 1> {
+    static inline bp2d::Coord get(bp2d::Segment const& seg) {
+        return bp2d::getY(seg.second());
+    }
+    static inline void set(bp2d::Segment &seg, bp2d::Coord const& coord) {
+        auto p = seg.second(); bp2d::setY(p, coord); seg.second(p);
+    }
+};
+
+
+/* ************************************************************************** */
+/* Polygon concept adaptation *********************************************** */
+/* ************************************************************************** */
+
+// Connversion between libnest2d::Orientation and order_selector ///////////////
+
+template<bp2d::Orientation> struct ToBoostOrienation {};
+
+template<>
+struct ToBoostOrienation<bp2d::Orientation::CLOCKWISE> {
+    static const order_selector Value = clockwise;
+};
+
+template<>
+struct ToBoostOrienation<bp2d::Orientation::COUNTER_CLOCKWISE> {
+    static const order_selector Value = counterclockwise;
+};
+
+static const bp2d::Orientation RealOrientation =
+        bp2d::OrientationType<bp2d::PolygonImpl>::Value;
+
+// Ring implementation /////////////////////////////////////////////////////////
+
+// Boost would refer to ClipperLib::Path (alias bp2d::PolygonImpl) as a ring
+template<> struct tag<bp2d::PathImpl> {
+    using type = ring_tag;
+};
+
+template<> struct point_order<bp2d::PathImpl> {
+    static const order_selector value =
+            ToBoostOrienation<RealOrientation>::Value;
+};
+
+// All our Paths should be closed for the bin packing application
+template<> struct closure<bp2d::PathImpl> {
+    static const closure_selector value = closed;
+};
+
+// Polygon implementation //////////////////////////////////////////////////////
+
+template<> struct tag<bp2d::PolygonImpl> {
+    using type = polygon_tag;
+};
+
+template<> struct exterior_ring<bp2d::PolygonImpl> {
+    static inline bp2d::PathImpl& get(bp2d::PolygonImpl& p) {
+        return libnest2d::shapelike::contour(p);
+    }
+
+    static inline bp2d::PathImpl const& get(bp2d::PolygonImpl const& p) {
+        return libnest2d::shapelike::contour(p);
+    }
+};
+
+template<> struct ring_const_type<bp2d::PolygonImpl> {
+    using type = const bp2d::PathImpl&;
+};
+
+template<> struct ring_mutable_type<bp2d::PolygonImpl> {
+    using type = bp2d::PathImpl&;
+};
+
+template<> struct interior_const_type<bp2d::PolygonImpl> {
+   using type = const libnest2d::THolesContainer<bp2d::PolygonImpl>&;
+};
+
+template<> struct interior_mutable_type<bp2d::PolygonImpl> {
+   using type = libnest2d::THolesContainer<bp2d::PolygonImpl>&;
+};
+
+template<>
+struct interior_rings<bp2d::PolygonImpl> {
+
+    static inline libnest2d::THolesContainer<bp2d::PolygonImpl>& get(
+            bp2d::PolygonImpl& p)
+    {
+        return libnest2d::shapelike::holes(p);
+    }
+
+    static inline const libnest2d::THolesContainer<bp2d::PolygonImpl>& get(
+            bp2d::PolygonImpl const& p)
+    {
+        return libnest2d::shapelike::holes(p);
+    }
+};
+
+/* ************************************************************************** */
+/* MultiPolygon concept adaptation ****************************************** */
+/* ************************************************************************** */
+
+template<> struct tag<bp2d::Shapes> {
+    using type = multi_polygon_tag;
+};
+
+}   // traits
+}   // geometry
+
+// This is an addition to the ring implementation of Polygon concept
+template<>
+struct range_value<bp2d::PathImpl> {
+    using type = bp2d::PointImpl;
+};
+
+template<>
+struct range_value<bp2d::Shapes> {
+    using type = bp2d::PolygonImpl;
+};
+
+}   // boost
+
+/* ************************************************************************** */
+/* Algorithms *************************************************************** */
+/* ************************************************************************** */
+
+namespace libnest2d { // Now the algorithms that boost can provide...
+
+namespace pointlike {
+template<>
+inline double distance(const PointImpl& p1, const PointImpl& p2 )
+{
+    return boost::geometry::distance(p1, p2);
+}
+
+template<>
+inline double distance(const PointImpl& p, const bp2d::Segment& seg )
+{
+    return boost::geometry::distance(p, seg);
+}
+}
+
+namespace shapelike {
+// Tell libnest2d how to make string out of a ClipperPolygon object
+template<>
+inline bool intersects(const PathImpl& sh1, const PathImpl& sh2)
+{
+    return boost::geometry::intersects(sh1, sh2);
+}
+
+// Tell libnest2d how to make string out of a ClipperPolygon object
+template<>
+inline bool intersects(const PolygonImpl& sh1, const PolygonImpl& sh2)
+{
+    return boost::geometry::intersects(sh1, sh2);
+}
+
+// Tell libnest2d how to make string out of a ClipperPolygon object
+template<>
+inline bool intersects(const bp2d::Segment& s1, const bp2d::Segment& s2)
+{
+    return boost::geometry::intersects(s1, s2);
+}
+
+#ifndef DISABLE_BOOST_AREA
+template<>
+inline double area(const PolygonImpl& shape, const PolygonTag&)
+{
+    return boost::geometry::area(shape);
+}
+#endif
+
+template<>
+inline bool isInside(const PointImpl& point, const PolygonImpl& shape)
+{
+    return boost::geometry::within(point, shape);
+}
+
+template<>
+inline bool isInside(const PolygonImpl& sh1, const PolygonImpl& sh2)
+{
+    return boost::geometry::within(sh1, sh2);
+}
+
+template<>
+inline bool touches(const PolygonImpl& sh1, const PolygonImpl& sh2)
+{
+    return boost::geometry::touches(sh1, sh2);
+}
+
+template<>
+inline bool touches( const PointImpl& point, const PolygonImpl& shape)
+{
+    return boost::geometry::touches(point, shape);
+}
+
+#ifndef DISABLE_BOOST_BOUNDING_BOX
+template<>
+inline bp2d::Box boundingBox(const PolygonImpl& sh, const PolygonTag&)
+{
+    bp2d::Box b;
+    boost::geometry::envelope(sh, b);
+    return b;
+}
+
+template<>
+inline bp2d::Box boundingBox(const PathImpl& sh, const PolygonTag&)
+{
+    bp2d::Box b;
+    boost::geometry::envelope(sh, b);
+    return b;
+}
+
+template<>
+inline bp2d::Box boundingBox<bp2d::Shapes>(const bp2d::Shapes& shapes,
+                                           const MultiPolygonTag&)
+{
+    bp2d::Box b;
+    boost::geometry::envelope(shapes, b);
+    return b;
+}
+#endif
+
+#ifndef DISABLE_BOOST_CONVEX_HULL
+template<>
+inline PolygonImpl convexHull(const PolygonImpl& sh, const PolygonTag&)
+{
+    PolygonImpl ret;
+    boost::geometry::convex_hull(sh, ret);
+    return ret;
+}
+
+template<>
+inline PolygonImpl convexHull(const TMultiShape<PolygonImpl>& shapes,
+                              const MultiPolygonTag&)
+{
+    PolygonImpl ret;
+    boost::geometry::convex_hull(shapes, ret);
+    return ret;
+}
+#endif
+
+#ifndef DISABLE_BOOST_OFFSET
+template<>
+inline void offset(PolygonImpl& sh, bp2d::Coord distance)
+{
+    PolygonImpl cpy = sh;
+    boost::geometry::buffer(cpy, sh, distance);
+}
+#endif
+
+#ifndef DISABLE_BOOST_SERIALIZE
+template<> inline std::string serialize<libnest2d::Formats::SVG>(
+        const PolygonImpl& sh, double scale)
+{
+    std::stringstream ss;
+    std::string style = "fill: none; stroke: black; stroke-width: 1px;";
+
+    using namespace boost::geometry;
+    using Pointf = model::point<double, 2, cs::cartesian>;
+    using Polygonf = model::polygon<Pointf>;
+
+    Polygonf::ring_type ring;
+    Polygonf::inner_container_type holes;
+    ring.reserve(shapelike::contourVertexCount(sh));
+
+    for(auto it = shapelike::cbegin(sh); it != shapelike::cend(sh); it++) {
+        auto& v = *it;
+        ring.emplace_back(getX(v)*scale, getY(v)*scale);
+    };
+
+    auto H = shapelike::holes(sh);
+    for(PathImpl& h : H ) {
+        Polygonf::ring_type hf;
+        for(auto it = h.begin(); it != h.end(); it++) {
+            auto& v = *it;
+            hf.emplace_back(getX(v)*scale, getY(v)*scale);
+        };
+        holes.push_back(hf);
+    }
+
+    Polygonf poly;
+    poly.outer() = ring;
+    poly.inners() = holes;
+    auto svg_data = boost::geometry::svg(poly, style);
+
+    ss << svg_data << std::endl;
+
+    return ss.str();
+}
+#endif
+
+#ifndef DISABLE_BOOST_UNSERIALIZE
+template<>
+inline void unserialize<libnest2d::Formats::SVG>(
+        PolygonImpl& sh,
+        const std::string& str)
+{
+}
+#endif
+
+template<> inline std::pair<bool, std::string> isValid(const PolygonImpl& sh)
+{
+    std::string message;
+    bool ret = boost::geometry::is_valid(sh, message);
+
+    return {ret, message};
+}
+}
+
+namespace nfp {
+
+#ifndef DISABLE_BOOST_NFP_MERGE
+
+// Warning: I could not get boost union_ to work. Geometries will overlap.
+template<>
+inline bp2d::Shapes nfp::merge(const bp2d::Shapes& shapes,
+                               const PolygonImpl& sh)
+{
+    bp2d::Shapes retv;
+    boost::geometry::union_(shapes, sh, retv);
+    return retv;
+}
+
+template<>
+inline bp2d::Shapes nfp::merge(const bp2d::Shapes& shapes)
+{
+    bp2d::Shapes retv;
+    boost::geometry::union_(shapes, shapes.back(), retv);
+    return retv;
+}
+#endif
+
+}
+
+
+}
+
+
+
+#endif // BOOST_ALG_HPP

src/libnest2d/include/libnest2d/utils/metaloop.hpp

@@ -0,0 +1,227 @@
+#ifndef METALOOP_HPP
+#define METALOOP_HPP
+
+#include <libnest2d/common.hpp>
+#include <tuple>
+#include <functional>
+
+namespace libnest2d {
+
+/* ************************************************************************** */
+/* C++14 std::index_sequence implementation:                                  */
+/* ************************************************************************** */
+
+/**
+ * \brief C++11 compatible implementation of the index_sequence type from C++14
+ */
+template<size_t...Ints> struct index_sequence {
+    using value_type = size_t;
+    BP2D_CONSTEXPR value_type size() const { return sizeof...(Ints); }
+};
+
+// A Help structure to generate the integer list
+template<size_t...Nseq> struct genSeq;
+
+// Recursive template to generate the list
+template<size_t I, size_t...Nseq> struct genSeq<I, Nseq...> {
+    // Type will contain a genSeq with Nseq appended by one element
+    using Type = typename genSeq< I - 1, I - 1, Nseq...>::Type;
+};
+
+// Terminating recursion
+template <size_t ... Nseq> struct genSeq<0, Nseq...> {
+    // If I is zero, Type will contain index_sequence with the fuly generated
+    // integer list.
+    using Type = index_sequence<Nseq...>;
+};
+
+/// Helper alias to make an index sequence from 0 to N
+template<size_t N> using make_index_sequence = typename genSeq<N>::Type;
+
+/// Helper alias to make an index sequence for a parameter pack
+template<class...Args>
+using index_sequence_for = make_index_sequence<sizeof...(Args)>;
+
+
+/* ************************************************************************** */
+
+namespace opt {
+
+using std::forward;
+using std::tuple;
+using std::get;
+using std::tuple_element;
+
+/**
+ * @brief Helper class to be able to loop over a parameter pack's elements.
+ */
+class metaloop {
+
+// The implementation is based on partial struct template specializations.
+// Basically we need a template type that is callable and takes an integer
+// non-type template parameter which can be used to implement recursive calls.
+//
+// C++11 will not allow the usage of a plain template function that is why we
+// use struct with overloaded call operator. At the same time C++11 prohibits
+// partial template specialization with a non type parameter such as int. We
+// need to wrap that in a type (see metaloop::Int).
+
+/*
+ * A helper alias to create integer values wrapped as a type. It is necessary
+ * 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 necessary 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 lambda 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 meta-programming pattern.
+ */
+template<class...Args>
+using MetaLoop = _MetaLoop<Int<sizeof...(Args)-1>, Args...>;
+
+public:
+
+/**
+ * \brief The final usable function template.
+ *
+ * This is similar to what varags was on C but in compile time C++11.
+ * You can call:
+ * apply(<the mapping function>, <arbitrary number of arguments of any type>);
+ * For example:
+ *
+ *      struct mapfunc {
+ *          template<class T> void operator()(int N, T&& element) {
+ *              std::cout << "The value of the parameter "<< N <<": "
+ *                        << element << std::endl;
+ *          }
+ *      };
+ *
+ *      apply(mapfunc(), 'a', 10, 151.545);
+ *
+ * C++14:
+ *      apply([](int N, auto&& element){
+ *          std::cout << "The value of the parameter "<< N <<": "
+ *                        << element << std::endl;
+ *      }, 'a', 10, 151.545);
+ *
+ * This yields the output:
+ * The value of the parameter 0: a
+ * The value of the parameter 1: 10
+ * The value of the parameter 2: 151.545
+ *
+ * As an addition, the function can be called with a tuple as the second
+ * parameter holding the arguments instead of a parameter pack.
+ *
+ */
+template<class...Args, class Fn>
+inline static void apply(Fn&& fn, Args&&...args) {
+    MetaLoop<Args...>().run(tuple<Args&&...>(forward<Args>(args)...),
+                            forward<Fn>(fn));
+}
+
+/// The version of apply with a tuple rvalue reference.
+template<class...Args, class Fn>
+inline static void apply(Fn&& fn, tuple<Args...>&& tup) {
+    MetaLoop<Args...>().run(std::move(tup), forward<Fn>(fn));
+}
+
+/// The version of apply with a tuple lvalue reference.
+template<class...Args, class Fn>
+inline static void apply(Fn&& fn, tuple<Args...>& tup) {
+    MetaLoop<Args...>().run(tup, forward<Fn>(fn));
+}
+
+/// The version of apply with a tuple const reference.
+template<class...Args, class Fn>
+inline static void apply(Fn&& fn, const tuple<Args...>& tup) {
+    MetaLoop<Args...>().run(tup, forward<Fn>(fn));
+}
+
+/**
+ * Call a function with its arguments encapsualted in a tuple.
+ */
+template<class Fn, class Tup, std::size_t...Is>
+inline static auto
+callFunWithTuple(Fn&& fn, Tup&& tup, index_sequence<Is...>) ->
+    decltype(fn(std::get<Is>(tup)...))
+{
+    return fn(std::get<Is>(tup)...);
+}
+
+};
+}
+}
+
+#endif // METALOOP_HPP

src/libnest2d/include/libnest2d/utils/rotfinder.hpp

@@ -0,0 +1,41 @@
+#ifndef ROTFINDER_HPP
+#define ROTFINDER_HPP
+
+#include <libnest2d/libnest2d.hpp>
+#include <libnest2d/optimizer.hpp>
+#include <iterator>
+
+namespace libnest2d {
+
+template<class RawShape>
+Radians findBestRotation(_Item<RawShape>& item) {
+    opt::StopCriteria stopcr;
+    stopcr.absolute_score_difference = 0.01;
+    stopcr.max_iterations = 10000;
+    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 Iterator>
+void findMinimumBoundingBoxRotations(Iterator from, Iterator to) {
+    using V = typename std::iterator_traits<Iterator>::value_type;
+    std::for_each(from, to, [](V& item){
+        Radians rot = findBestRotation(item);
+        item.rotate(rot);
+    });
+}
+
+}
+
+#endif // ROTFINDER_HPP