Importing the SLA computing module into the native source tree.

CMakeLists.txt

@@ -27,6 +27,10 @@ option(SLIC3R_MSVC_COMPILE_PARALLEL "Compile on Visual Studio in parallel" 1)
 option(SLIC3R_MSVC_PDB          "Generate PDB files on MSVC in Release mode" 1)
 option(SLIC3R_PERL_XS           "Compile XS Perl module and enable Perl unit and integration tests" 0)
 
+# Proposal for C++ unit tests and sandboxes
+option(SLIC3R_BUILD_SANDBOXES   "Build development sandboxes" OFF)
+option(SLIC3R_BUILD_TESTS       "Build unit tests" OFF)
+
 if (MSVC)
     if (SLIC3R_MSVC_COMPILE_PARALLEL)
 	   add_compile_options(/MP)
@@ -244,5 +248,13 @@ if (SLIC3R_PERL_XS)
     add_subdirectory(xs)
 endif ()
 
+if(SLIC3R_BUILD_SANDBOXES)
+    add_subdirectory(sandboxes)
+endif()
+
+if(SLIC3R_BUILD_TESTS)
+    add_subdirectory(tests)
+endif()
+
 file(GLOB MyVar var/*.png)
 install(FILES ${MyVar} DESTINATION share/slic3r-prusa3d)

sandboxes/CMakeLists.txt

@@ -0,0 +1,2 @@
+add_subdirectory(slabasebed)
+add_subdirectory(slasupporttree)

sandboxes/slabasebed/CMakeLists.txt

@@ -0,0 +1,2 @@
+add_executable(slabasebed EXCLUDE_FROM_ALL slabasebed.cpp)
+target_link_libraries(slabasebed libslic3r)

sandboxes/slasupporttree/CMakeLists.txt

@@ -0,0 +1,2 @@
+add_executable(slasupporttree EXCLUDE_FROM_ALL slasupporttree.cpp)
+target_link_libraries(slasupporttree libslic3r)

sandboxes/slasupporttree/slasupporttree.cpp

@@ -0,0 +1,48 @@
+#include <iostream>
+#include <string>
+
+#include <libslic3r.h>
+#include "TriangleMesh.hpp"
+#include "Model.hpp"
+#include "callback.hpp"
+#include "SLA/SLASupportTree.hpp"
+#include "benchmark.h"
+
+const std::string USAGE_STR = {
+    "Usage: slasupporttree stlfilename.stl"
+};
+
+void confess_at(const char * /*file*/,
+                int /*line*/,
+                const char * /*func*/,
+                const char * /*pat*/,
+                ...) {}
+
+namespace Slic3r {
+
+void PerlCallback::deregister_callback() {}
+}
+
+int main(const int argc, const char *argv[]) {
+    using namespace Slic3r;
+    using std::cout; using std::endl;
+
+    if(argc < 2) {
+        cout << USAGE_STR << endl;
+        return EXIT_SUCCESS;
+    }
+
+    Benchmark bench;
+    TriangleMesh result;
+
+    bench.start();
+    sla::create_head(result, 3, 1, 4);
+    bench.stop();
+
+    cout << "Support tree creation time: " << std::setprecision(10)
+         << bench.getElapsedSec() << " seconds." << endl;
+
+    result.write_ascii("out.stl");
+
+    return EXIT_SUCCESS;
+}

src/igl/AABB.cpp

@@ -885,7 +885,8 @@ igl::AABB<DerivedV,DIM>::intersect_ray(
   const RowVectorDIMS & dir,
   igl::Hit & hit) const
 {
-#if false
+    // FIXME: Needs a proper path
+#if /*false*/ 0
   // BFS
   std::queue<const AABB *> Q;
   // Or DFS

src/libnest2d/tools/benchmark.h

@@ -0,0 +1,58 @@
+/*
+ * Copyright (C) Tamás Mészáros
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ */
+#ifndef INCLUDE_BENCHMARK_H_
+#define INCLUDE_BENCHMARK_H_
+
+#include <chrono>
+#include <ratio>
+
+/**
+ * A class for doing benchmarks.
+ */
+class Benchmark {
+    typedef std::chrono::high_resolution_clock Clock;
+    typedef Clock::duration Duration;
+    typedef Clock::time_point TimePoint;
+
+    TimePoint t1, t2;
+    Duration d;
+
+    inline double to_sec(Duration d) {
+        return d.count() * double(Duration::period::num) / Duration::period::den;
+    }
+
+public:
+
+    /**
+     *  Measure time from the moment of this call.
+     */
+    void start() {  t1 = Clock::now(); }
+
+    /**
+     *  Measure time to the moment of this call.
+     */
+    void stop() {  t2 = Clock::now(); }
+
+    /**
+     * Get the time elapsed between a start() end a stop() call.
+     * @return Returns the elapsed time in seconds.
+     */
+    double getElapsedSec() {  d = t2 - t1; return to_sec(d); }
+};
+
+
+#endif /* INCLUDE_BENCHMARK_H_ */

src/libslic3r/CMakeLists.txt

@@ -153,6 +153,16 @@ add_library(libslic3r STATIC
     SLABasePool.cpp
     utils.cpp
     Utils.hpp
+    SLA/SLABoilerPlate.hpp
+    SLA/SLABasePool.hpp
+    SLA/SLABasePool.cpp
+    SLA/SLASupportTree.hpp
+    SLA/SLASupportTree.cpp
+    SLA/SLASupportTreeIGL.cpp
+    SLA/SLARotfinder.hpp
+    SLA/SLARotfinder.cpp
+    SLA/SLABoostAdapter.hpp
+    SLA/SLASpatIndex.hpp
 )
 
 add_precompiled_header(libslic3r pchheader.hpp FORCEINCLUDE)

src/libslic3r/SLA/SLABasePool.cpp

@@ -0,0 +1,534 @@
+#include "SLABasePool.hpp"
+#include "SLABoilerPlate.hpp"
+
+#include "boost/log/trivial.hpp"
+#include "SLABoostAdapter.hpp"
+#include "ClipperUtils.hpp"
+
+//#include "SVG.hpp"
+//#include "benchmark.h"
+
+namespace Slic3r { namespace sla {
+
+/// Convert the triangulation output to an intermediate mesh.
+Contour3D convert(const Polygons& triangles, coord_t z, bool dir) {
+
+    Pointf3s points;
+    points.reserve(3*triangles.size());
+    Indices indices;
+    indices.reserve(points.size());
+
+    for(auto& tr : triangles) {
+        auto c = coord_t(points.size()), b = c++, a = c++;
+        if(dir) indices.emplace_back(a, b, c);
+        else indices.emplace_back(c, b, a);
+        for(auto& p : tr.points) {
+            points.emplace_back(unscale(x(p), y(p), z));
+        }
+    }
+
+    return {points, indices};
+}
+
+Contour3D walls(const ExPolygon& floor_plate, const ExPolygon& ceiling,
+                       double floor_z_mm, double ceiling_z_mm) {
+    using std::transform; using std::back_inserter;
+
+    ExPolygon poly;
+    poly.contour.points = floor_plate.contour.points;
+    poly.holes.emplace_back(ceiling.contour);
+    auto& h = poly.holes.front();
+    std::reverse(h.points.begin(), h.points.end());
+    Polygons tri = triangulate(poly);
+
+    Contour3D ret;
+    ret.points.reserve(tri.size() * 3);
+
+    double fz = floor_z_mm;
+    double cz = ceiling_z_mm;
+    auto& rp = ret.points;
+    auto& rpi = ret.indices;
+    ret.indices.reserve(tri.size() * 3);
+
+    coord_t idx = 0;
+
+    auto hlines = h.lines();
+    auto is_upper = [&hlines](const Point& p) {
+        return std::any_of(hlines.begin(), hlines.end(),
+                               [&p](const Line& l) {
+            return l.distance_to(p) < mm(1e-6);
+        });
+    };
+
+    std::for_each(tri.begin(), tri.end(),
+                  [&rp, &rpi, &poly, &idx, is_upper, fz, cz](const Polygon& pp)
+    {
+        for(auto& p : pp.points)
+            if(is_upper(p))
+                rp.emplace_back(unscale(x(p), y(p), mm(cz)));
+            else rp.emplace_back(unscale(x(p), y(p), mm(fz)));
+
+        coord_t a = idx++, b = idx++, c = idx++;
+        if(fz > cz) rpi.emplace_back(c, b, a);
+        else rpi.emplace_back(a, b, c);
+    });
+
+    return ret;
+}
+
+/// Offsetting with clipper and smoothing the edges into a curvature.
+void offset(ExPolygon& sh, coord_t distance) {
+    using ClipperLib::ClipperOffset;
+    using ClipperLib::jtRound;
+    using ClipperLib::etClosedPolygon;
+    using ClipperLib::Paths;
+    using ClipperLib::Path;
+
+    auto&& ctour = Slic3rMultiPoint_to_ClipperPath(sh.contour);
+    auto&& holes = Slic3rMultiPoints_to_ClipperPaths(sh.holes);
+
+    // If the input is not at least a triangle, we can not do this algorithm
+    if(ctour.size() < 3 ||
+       std::any_of(holes.begin(), holes.end(),
+                   [](const Path& p) { return p.size() < 3; })
+            ) {
+        BOOST_LOG_TRIVIAL(error) << "Invalid geometry for offsetting!";
+        return;
+    }
+
+    ClipperOffset offs;
+    offs.ArcTolerance = 0.01*mm(1);
+    Paths result;
+    offs.AddPath(ctour, jtRound, etClosedPolygon);
+    offs.AddPaths(holes, jtRound, 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;
+    sh.holes.clear();
+    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) {
+                auto rr = ClipperPath_to_Slic3rPolygon(r);
+                sh.contour.points.swap(rr.points);
+                found_the_contour = true;
+            } else {
+                BOOST_LOG_TRIVIAL(warning)
+                        << "Warning: offsetting result is invalid!";
+            }
+        } 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.emplace_back(ClipperPath_to_Slic3rPolygon(r));
+        }
+    }
+}
+
+/// Unification of polygons (with clipper) preserving holes as well.
+ExPolygons unify(const ExPolygons& shapes) {
+    using ClipperLib::ptSubject;
+
+    ExPolygons retv;
+
+    bool closed = true;
+    bool valid = true;
+
+    ClipperLib::Clipper clipper;
+
+    for(auto& path : shapes) {
+        auto clipperpath = Slic3rMultiPoint_to_ClipperPath(path.contour);
+
+        if(!clipperpath.empty())
+            valid &= clipper.AddPath(clipperpath, ptSubject, closed);
+
+        auto clipperholes = Slic3rMultiPoints_to_ClipperPaths(path.holes);
+
+        for(auto& hole : clipperholes) {
+            if(!hole.empty())
+                valid &= clipper.AddPath(hole, ptSubject, closed);
+        }
+    }
+
+    if(!valid) BOOST_LOG_TRIVIAL(warning) << "Unification of invalid shapes!";
+
+    ClipperLib::PolyTree result;
+    clipper.Execute(ClipperLib::ctUnion, result, ClipperLib::pftNonZero);
+
+    retv.reserve(static_cast<size_t>(result.Total()));
+
+    // Now we will recursively traverse the polygon tree and serialize it
+    // into an ExPolygon with holes. The polygon tree has the clipper-ish
+    // PolyTree structure which alternates its nodes as contours and holes
+
+    // A "declaration" of function for traversing leafs which are holes
+    std::function<void(ClipperLib::PolyNode*, ExPolygon&)> processHole;
+
+    // Process polygon which calls processHoles which than calls processPoly
+    // again until no leafs are left.
+    auto processPoly = [&retv, &processHole](ClipperLib::PolyNode *pptr) {
+        ExPolygon poly;
+        poly.contour.points = ClipperPath_to_Slic3rPolygon(pptr->Contour);
+        for(auto h : pptr->Childs) { processHole(h, poly); }
+        retv.push_back(poly);
+    };
+
+    // Body of the processHole function
+    processHole = [&processPoly](ClipperLib::PolyNode *pptr, ExPolygon& poly)
+    {
+        poly.holes.emplace_back();
+        poly.holes.back().points = ClipperPath_to_Slic3rPolygon(pptr->Contour);
+        for(auto c : pptr->Childs) processPoly(c);
+    };
+
+    // Wrapper for traversing.
+    auto traverse = [&processPoly] (ClipperLib::PolyNode *node)
+    {
+        for(auto ch : node->Childs) {
+            processPoly(ch);
+        }
+    };
+
+    // Here is the actual traverse
+    traverse(&result);
+
+    return retv;
+}
+
+/// Only a debug function to generate top and bottom plates from a 2D shape.
+/// It is not used in the algorithm directly.
+inline Contour3D roofs(const ExPolygon& poly, coord_t z_distance) {
+    Polygons triangles = triangulate(poly);
+
+    auto lower = convert(triangles, 0, false);
+    auto upper = convert(triangles, z_distance, true);
+    lower.merge(upper);
+    return lower;
+}
+
+template<class ExP, class D>
+Contour3D round_edges(const ExPolygon& base_plate,
+                            double radius_mm,
+                            double degrees,
+                            double ceilheight_mm,
+                            bool dir,
+                            ExP&& last_offset = ExP(), D&& last_height = D())
+{
+    auto ob = base_plate;
+    auto ob_prev = ob;
+    double wh = ceilheight_mm, wh_prev = wh;
+    Contour3D curvedwalls;
+
+    int steps = 15; // int(std::ceil(10*std::pow(radius_mm, 1.0/3)));
+    double stepx = radius_mm / steps;
+    coord_t s = dir? 1 : -1;
+    degrees = std::fmod(degrees, 180);
+
+    if(degrees >= 90) {
+        for(int i = 1; i <= steps; ++i) {
+            ob = base_plate;
+
+            double r2 = radius_mm * radius_mm;
+            double xx = i*stepx;
+            double x2 = xx*xx;
+            double stepy = std::sqrt(r2 - x2);
+
+            offset(ob, s*mm(xx));
+            wh = ceilheight_mm - radius_mm + stepy;
+
+            Contour3D pwalls;
+            pwalls = walls(ob, ob_prev, wh, wh_prev);
+
+            curvedwalls.merge(pwalls);
+            ob_prev = ob;
+            wh_prev = wh;
+        }
+    }
+
+    double tox = radius_mm - radius_mm*std::sin(degrees * PI / 180);
+    int tos = int(tox / stepx);
+
+    for(int i = 1; i <= tos; ++i) {
+        ob = base_plate;
+
+        double r2 = radius_mm * radius_mm;
+        double xx = radius_mm - i*stepx;
+        double x2 = xx*xx;
+        double stepy = std::sqrt(r2 - x2);
+        offset(ob, s*mm(xx));
+        wh = ceilheight_mm - radius_mm - stepy;
+
+        Contour3D pwalls;
+        pwalls = walls(ob_prev, ob, wh_prev, wh);
+
+        curvedwalls.merge(pwalls);
+        ob_prev = ob;
+        wh_prev = wh;
+    }
+
+    last_offset = std::move(ob);
+    last_height = wh;
+
+    return curvedwalls;
+}
+
+/// Generating the concave part of the 3D pool with the bottom plate and the
+/// side walls.
+Contour3D inner_bed(const ExPolygon& poly, double depth_mm,
+                           double begin_h_mm = 0) {
+
+    Polygons triangles = triangulate(poly);
+
+    coord_t depth = mm(depth_mm);
+    coord_t begin_h = mm(begin_h_mm);
+
+    auto bottom = convert(triangles, -depth + begin_h, false);
+    auto lines = poly.lines();
+
+    // Generate outer walls
+    auto fp = [](const Point& p, Point::coord_type z) {
+        return unscale(x(p), y(p), z);
+    };
+
+    for(auto& l : lines) {
+        auto s = coord_t(bottom.points.size());
+
+        bottom.points.emplace_back(fp(l.a, -depth + begin_h));
+        bottom.points.emplace_back(fp(l.b, -depth + begin_h));
+        bottom.points.emplace_back(fp(l.a, begin_h));
+        bottom.points.emplace_back(fp(l.b, begin_h));
+
+        bottom.indices.emplace_back(s + 3, s + 1, s);
+        bottom.indices.emplace_back(s + 2, s + 3, s);
+    }
+
+    return bottom;
+}
+
+inline Point centroid(Points& pp) {
+    Point c;
+    switch(pp.size()) {
+    case 0: break;
+    case 1: c = pp.front(); break;
+    case 2: c = (pp[0] + pp[1]) / 2; break;
+    default: {
+        auto MAX = std::numeric_limits<Point::coord_type>::max();
+        auto MIN = std::numeric_limits<Point::coord_type>::min();
+        Point min = {MAX, MAX}, max = {MIN, MIN};
+
+        for(auto& p : pp) {
+            if(p(0) < min(0)) min(0) = p(0);
+            if(p(1) < min(1)) min(1) = p(1);
+            if(p(0) > max(0)) max(0) = p(0);
+            if(p(1) > max(1)) max(1) = p(1);
+        }
+        c(0) = min(0) + (max(0) - min(0)) / 2;
+        c(1) = min(1) + (max(1) - min(1)) / 2;
+
+        // TODO: fails for non convex cluster
+//        c = std::accumulate(pp.begin(), pp.end(), Point{0, 0});
+//        x(c) /= coord_t(pp.size()); y(c) /= coord_t(pp.size());
+        break;
+    }
+    }
+
+    return c;
+}
+
+inline Point centroid(const ExPolygon& poly) {
+    return poly.contour.centroid();
+}
+
+/// A fake concave hull that is constructed by connecting separate shapes
+/// with explicit bridges. Bridges are generated from each shape's centroid
+/// to the center of the "scene" which is the centroid calculated from the shape
+/// centroids (a star is created...)
+ExPolygons concave_hull(const ExPolygons& polys, double max_dist_mm = 50)
+{
+    namespace bgi = boost::geometry::index;
+    using SpatElement = std::pair<BoundingBox, unsigned>;
+    using SpatIndex = bgi::rtree< SpatElement, bgi::rstar<16, 4> >;
+
+    if(polys.empty()) return ExPolygons();
+
+    ExPolygons punion = unify(polys);   // could be redundant
+
+    if(punion.size() == 1) return punion;
+
+    // We get the centroids of all the islands in the 2D slice
+    Points centroids; centroids.reserve(punion.size());
+    std::transform(punion.begin(), punion.end(), std::back_inserter(centroids),
+                   [](const ExPolygon& poly) { return centroid(poly); });
+
+
+    SpatIndex boxindex; unsigned idx = 0;
+    std::for_each(punion.begin(), punion.end(),
+                  [&boxindex, &idx](const ExPolygon& expo) {
+        BoundingBox bb(expo);
+        boxindex.insert(std::make_pair(bb, idx++));
+    });
+
+
+    // Centroid of the centroids of islands. This is where the additional
+    // connector sticks are routed.
+    Point cc = centroid(centroids);
+
+    punion.reserve(punion.size() + centroids.size());
+
+    idx = 0;
+    std::transform(centroids.begin(), centroids.end(),
+                   std::back_inserter(punion),
+                   [&punion, &boxindex, cc, max_dist_mm, &idx](const Point& c)
+    {
+
+        double dx = x(c) - x(cc), dy = y(c) - y(cc);
+        double l = std::sqrt(dx * dx + dy * dy);
+        double nx = dx / l, ny = dy / l;
+        double max_dist = mm(max_dist_mm);
+
+        ExPolygon& expo = punion[idx++];
+        BoundingBox querybb(expo);
+
+        querybb.offset(max_dist);
+        std::vector<SpatElement> result;
+        boxindex.query(bgi::intersects(querybb), std::back_inserter(result));
+        if(result.size() <= 1) return ExPolygon();
+
+        ExPolygon r;
+        auto& ctour = r.contour.points;
+
+        ctour.reserve(3);
+        ctour.emplace_back(cc);
+
+        Point d(coord_t(mm(1)*nx), coord_t(mm(1)*ny));
+        ctour.emplace_back(c + Point( -y(d),  x(d) ));
+        ctour.emplace_back(c + Point(  y(d), -x(d) ));
+        offset(r, mm(1));
+
+        return r;
+    });
+
+    punion = unify(punion);
+
+    return punion;
+}
+
+void base_plate(const TriangleMesh &mesh, ExPolygons &output, float h)
+{
+    TriangleMesh m = mesh;
+    TriangleMeshSlicer slicer(&m);
+
+    TriangleMesh upper, lower;
+    slicer.cut(h, &upper, &lower);
+    output = lower.horizontal_projection();
+
+    for(auto& o : output) o = o.simplify(0.1/SCALING_FACTOR).front();
+}
+
+void create_base_pool(const ExPolygons &ground_layer, TriangleMesh& out,
+                      const PoolConfig& cfg)
+{
+    double mdist = 2*(1.8*cfg.min_wall_thickness_mm + 4*cfg.edge_radius_mm) +
+                   cfg.max_merge_distance_mm;
+
+    auto concavehs = concave_hull(ground_layer, mdist);
+    for(ExPolygon& concaveh : concavehs) {
+        if(concaveh.contour.points.empty()) return;
+        concaveh.holes.clear();
+
+        const coord_t WALL_THICKNESS = mm(cfg.min_wall_thickness_mm);
+
+        const coord_t WALL_DISTANCE = mm(2*cfg.edge_radius_mm) +
+                                      coord_t(0.8*WALL_THICKNESS);
+
+        const coord_t HEIGHT = mm(cfg.min_wall_height_mm);
+
+        auto outer_base = concaveh;
+        offset(outer_base, WALL_THICKNESS+WALL_DISTANCE);
+        auto inner_base = outer_base;
+        offset(inner_base, -WALL_THICKNESS);
+        inner_base.holes.clear(); outer_base.holes.clear();
+
+        ExPolygon top_poly;
+        top_poly.contour = outer_base.contour;
+        top_poly.holes.emplace_back(inner_base.contour);
+        auto& tph = top_poly.holes.back().points;
+        std::reverse(tph.begin(), tph.end());
+
+        Contour3D pool;
+
+        ExPolygon ob = outer_base; double wh = 0;
+
+        // now we will calculate the angle or portion of the circle from
+        // pi/2 that will connect perfectly with the bottom plate.
+        // this is a tangent point calculation problem and the equation can
+        // be found for example here:
+        // http://www.ambrsoft.com/TrigoCalc/Circles2/CirclePoint/CirclePointDistance.htm
+        // the y coordinate would be:
+        // y = cy + (r^2*py - r*px*sqrt(px^2 + py^2 - r^2) / (px^2 + py^2)
+        // where px and py are the coordinates of the point outside the circle
+        // cx and cy are the circle center, r is the radius
+        // to get the angle we use arcsin function and subtract 90 degrees then
+        // flip the sign to get the right input to the round_edge function.
+        double r = cfg.edge_radius_mm;
+        double cy = 0;
+        double cx = 0;
+        double px = cfg.min_wall_thickness_mm;
+        double py = r - cfg.min_wall_height_mm;
+
+        double pxcx = px - cx;
+        double pycy = py - cy;
+        double b_2 = pxcx*pxcx + pycy*pycy;
+        double r_2 = r*r;
+        double D = std::sqrt(b_2 - r_2);
+        double vy = (r_2*pycy - r*pxcx*D) / b_2;
+        double phi = -(std::asin(vy/r) * 180 / PI - 90);
+
+        auto curvedwalls = round_edges(ob,
+                                       r,
+                                       phi,  // 170 degrees
+                                       0,    // z position of the input plane
+                                       true,
+                                       ob, wh);
+
+        pool.merge(curvedwalls);
+
+        ExPolygon ob_contr = ob;
+        ob_contr.holes.clear();
+
+        auto pwalls = walls(ob_contr, inner_base, wh, -cfg.min_wall_height_mm);
+        pool.merge(pwalls);
+
+        Polygons top_triangles, bottom_triangles;
+        triangulate(top_poly, top_triangles);
+        triangulate(inner_base, bottom_triangles);
+        auto top_plate = convert(top_triangles, 0, false);
+        auto bottom_plate = convert(bottom_triangles, -HEIGHT, true);
+
+        ob = inner_base; wh = 0;
+        // rounded edge generation for the inner bed
+        curvedwalls = round_edges(ob,
+                                  cfg.edge_radius_mm,
+                                  90,   // 90 degrees
+                                  0,    // z position of the input plane
+                                  false,
+                                  ob, wh);
+        pool.merge(curvedwalls);
+
+        auto innerbed = inner_bed(ob, cfg.min_wall_height_mm/2 + wh, wh);
+
+        pool.merge(top_plate);
+        pool.merge(bottom_plate);
+        pool.merge(innerbed);
+
+        out.merge(mesh(pool));
+    }
+}
+
+}
+}

src/libslic3r/SLA/SLABasePool.hpp

@@ -0,0 +1,37 @@
+#ifndef SLABASEPOOL_HPP
+#define SLABASEPOOL_HPP
+
+#include <vector>
+
+namespace Slic3r {
+
+class ExPolygon;
+class TriangleMesh;
+
+namespace sla {
+
+using ExPolygons = std::vector<ExPolygon>;
+
+/// Calculate the polygon representing the silhouette from the specified height
+void base_plate(const TriangleMesh& mesh,
+                ExPolygons& output,
+                float height = 0.1f);
+
+struct PoolConfig {
+    double min_wall_thickness_mm = 2;
+    double min_wall_height_mm = 5;
+    double max_merge_distance_mm = 50;
+    double edge_radius_mm = 1;
+};
+
+/// Calculate the pool for the mesh for SLA printing
+void create_base_pool(const ExPolygons& base_plate,
+                      TriangleMesh& output_mesh,
+                      const PoolConfig& = PoolConfig()
+                      );
+
+}
+
+}
+
+#endif // SLABASEPOOL_HPP

src/libslic3r/SLA/SLABoilerPlate.hpp

@@ -0,0 +1,92 @@
+#ifndef SLABOILERPLATE_HPP
+#define SLABOILERPLATE_HPP
+
+#include <iostream>
+#include <functional>
+#include <numeric>
+
+#include "ExPolygon.hpp"
+#include "TriangleMesh.hpp"
+
+namespace Slic3r {
+namespace sla {
+
+using coord_t = Point::coord_type;
+
+/// get the scaled clipper units for a millimeter value
+inline coord_t mm(double v) { return coord_t(v/SCALING_FACTOR); }
+
+/// Get x and y coordinates (because we are eigenizing...)
+inline coord_t x(const Point& p) { return p(0); }
+inline coord_t y(const Point& p) { return p(1); }
+inline coord_t& x(Point& p) { return p(0); }
+inline coord_t& y(Point& p) { return p(1); }
+
+inline coordf_t x(const Vec3d& p) { return p(0); }
+inline coordf_t y(const Vec3d& p) { return p(1); }
+inline coordf_t z(const Vec3d& p) { return p(2); }
+inline coordf_t& x(Vec3d& p) { return p(0); }
+inline coordf_t& y(Vec3d& p) { return p(1); }
+inline coordf_t& z(Vec3d& p) { return p(2); }
+
+inline coord_t& x(Vec3crd& p) { return p(0); }
+inline coord_t& y(Vec3crd& p) { return p(1); }
+inline coord_t& z(Vec3crd& p) { return p(2); }
+inline coord_t x(const Vec3crd& p) { return p(0); }
+inline coord_t y(const Vec3crd& p) { return p(1); }
+inline coord_t z(const Vec3crd& p) { return p(2); }
+
+inline void triangulate(const ExPolygon& expoly, Polygons& triangles) {
+    expoly.triangulate_p2t(&triangles);
+}
+
+inline Polygons triangulate(const ExPolygon& expoly) {
+    Polygons tri; triangulate(expoly, tri); return tri;
+}
+
+using Indices = std::vector<Vec3crd>;
+
+/// Intermediate struct for a 3D mesh
+struct Contour3D {
+    Pointf3s points;
+    Indices indices;
+
+    void merge(const Contour3D& ctr) {
+        auto s3 = coord_t(points.size());
+        auto s = coord_t(indices.size());
+
+        points.insert(points.end(), ctr.points.begin(), ctr.points.end());
+        indices.insert(indices.end(), ctr.indices.begin(), ctr.indices.end());
+
+        for(auto n = s; n < indices.size(); n++) {
+            auto& idx = indices[n]; x(idx) += s3; y(idx) += s3; z(idx) += s3;
+        }
+    }
+};
+
+struct EigenMesh3D {
+    Eigen::MatrixXd V;
+    Eigen::MatrixXi F;
+};
+
+using PointSet = Eigen::MatrixXd;
+using ClusterEl = std::vector<unsigned>;
+using ClusteredPoints = std::vector<ClusterEl>;
+
+/// Convert the triangulation output to an intermediate mesh.
+Contour3D convert(const Polygons& triangles, coord_t z, bool dir);
+
+/// Mesh from an existing contour.
+inline TriangleMesh mesh(const Contour3D& ctour) {
+    return {ctour.points, ctour.indices};
+}
+
+/// Mesh from an evaporating 3D contour
+inline TriangleMesh mesh(Contour3D&& ctour) {
+    return {std::move(ctour.points), std::move(ctour.indices)};
+}
+
+}
+}
+
+#endif // SLABOILERPLATE_HPP

src/libslic3r/SLA/SLABoostAdapter.hpp

@@ -0,0 +1,132 @@
+#ifndef SLABOOSTADAPTER_HPP
+#define SLABOOSTADAPTER_HPP
+
+#include "SLA/SLABoilerPlate.hpp"
+#include <boost/geometry.hpp>
+
+namespace boost {
+namespace geometry {
+namespace traits {
+
+/* ************************************************************************** */
+/* Point concept adaptation ************************************************* */
+/* ************************************************************************** */
+
+template<> struct tag<Slic3r::Point> {
+    using type = point_tag;
+};
+
+template<> struct coordinate_type<Slic3r::Point> {
+    using type = coord_t;
+};
+
+template<> struct coordinate_system<Slic3r::Point> {
+    using type = cs::cartesian;
+};
+
+template<> struct dimension<Slic3r::Point>: boost::mpl::int_<2> {};
+
+template<int d> struct access<Slic3r::Point, d > {
+    static inline coord_t get(Slic3r::Point const& a) {
+        return a(d);
+    }
+
+    static inline void set(Slic3r::Point& a, coord_t const& value) {
+        a(d) = value;
+    }
+};
+
+// For Vec2d ///////////////////////////////////////////////////////////////////
+
+template<> struct tag<Slic3r::Vec2d> {
+    using type = point_tag;
+};
+
+template<> struct coordinate_type<Slic3r::Vec2d> {
+    using type = double;
+};
+
+template<> struct coordinate_system<Slic3r::Vec2d> {
+    using type = cs::cartesian;
+};
+
+template<> struct dimension<Slic3r::Vec2d>: boost::mpl::int_<2> {};
+
+template<int d> struct access<Slic3r::Vec2d, d > {
+    static inline double get(Slic3r::Vec2d const& a) {
+        return a(d);
+    }
+
+    static inline void set(Slic3r::Vec2d& a, double const& value) {
+        a(d) = value;
+    }
+};
+
+// For Vec3d ///////////////////////////////////////////////////////////////////
+
+template<> struct tag<Slic3r::Vec3d> {
+    using type = point_tag;
+};
+
+template<> struct coordinate_type<Slic3r::Vec3d> {
+    using type = double;
+};
+
+template<> struct coordinate_system<Slic3r::Vec3d> {
+    using type = cs::cartesian;
+};
+
+template<> struct dimension<Slic3r::Vec3d>: boost::mpl::int_<3> {};
+
+template<int d> struct access<Slic3r::Vec3d, d > {
+    static inline double get(Slic3r::Vec3d const& a) {
+        return a(d);
+    }
+
+    static inline void set(Slic3r::Vec3d& a, double const& value) {
+        a(d) = value;
+    }
+};
+
+/* ************************************************************************** */
+/* Box concept adaptation *************************************************** */
+/* ************************************************************************** */
+
+template<> struct tag<Slic3r::BoundingBox> {
+    using type = box_tag;
+};
+
+template<> struct point_type<Slic3r::BoundingBox> {
+    using type = Slic3r::Point;
+};
+
+template<std::size_t d>
+struct indexed_access<Slic3r::BoundingBox, 0, d> {
+    static inline coord_t get(Slic3r::BoundingBox const& box) {
+        return box.min(d);
+    }
+    static inline void set(Slic3r::BoundingBox &box, coord_t const& coord) {
+        box.min(d) = coord;
+    }
+};
+
+template<std::size_t d>
+struct indexed_access<Slic3r::BoundingBox, 1, d> {
+    static inline coord_t get(Slic3r::BoundingBox const& box) {
+        return box.max(d);
+    }
+    static inline void set(Slic3r::BoundingBox &box, coord_t const& coord) {
+        box.max(d) = coord;
+    }
+};
+
+}
+}
+
+template<> struct range_value<std::vector<Slic3r::Vec2d>> {
+    using type = Slic3r::Vec2d;
+};
+
+}
+
+#endif // SLABOOSTADAPTER_HPP

src/libslic3r/SLA/SLARotfinder.cpp

@@ -0,0 +1,156 @@
+#include <limits>
+#include <exception>
+
+#include <libnest2d/optimizers/nlopt/genetic.hpp>
+#include "SLABoilerPlate.hpp"
+#include "SLARotfinder.hpp"
+#include "Model.hpp"
+
+namespace Slic3r {
+namespace sla {
+
+EigenMesh3D to_eigenmesh(const TriangleMesh& tmesh) {
+
+    const stl_file& stl = tmesh.stl;
+
+    EigenMesh3D outmesh;
+    auto& V = outmesh.V;
+    auto& F = outmesh.F;
+
+    V.resize(3*stl.stats.number_of_facets, 3);
+    F.resize(stl.stats.number_of_facets, 3);
+    for (unsigned int i=0; i<stl.stats.number_of_facets; ++i) {
+        const stl_facet* facet = stl.facet_start+i;
+        V(3*i+0, 0) = facet->vertex[0](0); V(3*i+0, 1) =
+                facet->vertex[0](1); V(3*i+0, 2) = facet->vertex[0](2);
+        V(3*i+1, 0) = facet->vertex[1](0); V(3*i+1, 1) =
+                facet->vertex[1](1); V(3*i+1, 2) = facet->vertex[1](2);
+        V(3*i+2, 0) = facet->vertex[2](0); V(3*i+2, 1) =
+                facet->vertex[2](1); V(3*i+2, 2) = facet->vertex[2](2);
+
+        F(i, 0) = 3*i+0;
+        F(i, 1) = 3*i+1;
+        F(i, 2) = 3*i+2;
+    }
+
+    return outmesh;
+}
+
+inline EigenMesh3D to_eigenmesh(const ModelObject& modelobj) {
+    return to_eigenmesh(modelobj.raw_mesh());
+}
+
+std::array<double, 3> find_best_rotation(const ModelObject& modelobj,
+                                         float accuracy,
+                                         std::function<void(unsigned)> statuscb,
+                                         std::function<bool()> stopcond)
+{
+    using libnest2d::opt::Method;
+    using libnest2d::opt::bound;
+    using libnest2d::opt::Optimizer;
+    using libnest2d::opt::TOptimizer;
+    using libnest2d::opt::StopCriteria;
+    using Quaternion = Eigen::Quaternion<double>;
+
+    static const unsigned MAX_TRIES = 100000;
+
+    // return value
+    std::array<double, 3> rot;
+
+    // We will use only one instance of this converted mesh to examine different
+    // rotations
+    EigenMesh3D emesh = to_eigenmesh(modelobj);
+
+    // For current iteration number
+    unsigned status = 0;
+
+    // The maximum number of iterations
+    auto max_tries = unsigned(accuracy * MAX_TRIES);
+
+    // call status callback with zero, because we are at the start
+    statuscb(status);
+
+    // So this is the object function which is called by the solver many times
+    // It has to yield a single value representing the current score. We will
+    // call the status callback in each iteration but the actual value may be
+    // the same for subsequent iterations (status goes from 0 to 100 but
+    // iterations can be many more)
+    auto objfunc = [&emesh, &status, &statuscb, max_tries]
+            (double rx, double ry, double rz)
+    {
+        EigenMesh3D& m = emesh;
+
+        // prepare the rotation transformation
+        Transform3d rt = Transform3d::Identity();
+
+        rt.rotate(Eigen::AngleAxisd(rz, Vec3d::UnitZ()));
+        rt.rotate(Eigen::AngleAxisd(ry, Vec3d::UnitY()));
+        rt.rotate(Eigen::AngleAxisd(rx, Vec3d::UnitX()));
+
+        double score = 0;
+
+        // For all triangles we calculate the normal and sum up the dot product
+        // (a scalar indicating how much are two vectors aligned) with each axis
+        // this will result in a value that is greater if a normal is aligned
+        // with all axes. If the normal is aligned than the triangle itself is
+        // orthogonal to the axes and that is good for print quality.
+
+        // TODO: some applications optimize for minimum z-axis cross section
+        // area. The current function is only an example of how to optimize.
+
+        // Later we can add more criteria like the number of overhangs, etc...
+        for(int i = 0; i < m.F.rows(); i++) {
+            auto idx = m.F.row(i);
+
+            Vec3d p1 = m.V.row(idx(0));
+            Vec3d p2 = m.V.row(idx(1));
+            Vec3d p3 = m.V.row(idx(2));
+
+            Eigen::Vector3d U = p2 - p1;
+            Eigen::Vector3d V = p3 - p1;
+
+            // So this is the normal
+            auto n = U.cross(V).normalized();
+
+            // rotate the normal with the current rotation given by the solver
+            n = rt * n;
+
+            // We should score against the alignment with the reference planes
+            score += std::abs(n.dot(Vec3d::UnitX()));
+            score += std::abs(n.dot(Vec3d::UnitY()));
+            score += std::abs(n.dot(Vec3d::UnitZ()));
+        }
+
+        // report status
+        statuscb( unsigned(++status * 100.0/max_tries) );
+
+        return score;
+    };
+
+    // Firing up the genetic optimizer. For now it uses the nlopt library.
+    StopCriteria stc;
+    stc.max_iterations = max_tries;
+    stc.relative_score_difference = 1e-3;
+    stc.stop_condition = stopcond;      // stop when stopcond returns true
+    TOptimizer<Method::G_GENETIC> solver(stc);
+
+    // We are searching rotations around the three axes x, y, z. Thus the
+    // problem becomes a 3 dimensional optimization task.
+    // We can specify the bounds for a dimension in the following way:
+    auto b = bound(-PI/2, PI/2);
+
+    // Now we start the optimization process with initial angles (0, 0, 0)
+    auto result = solver.optimize_max(objfunc,
+                                      libnest2d::opt::initvals(0.0, 0.0, 0.0),
+                                      b, b, b);
+
+    // Save the result and fck off
+    rot[0] = std::get<0>(result.optimum);
+    rot[1] = std::get<1>(result.optimum);
+    rot[2] = std::get<2>(result.optimum);
+
+    return rot;
+}
+
+}
+}

src/libslic3r/SLA/SLARotfinder.hpp

@@ -0,0 +1,38 @@
+#ifndef SLAROTFINDER_HPP
+#define SLAROTFINDER_HPP
+
+#include <functional>
+#include <array>
+
+namespace Slic3r {
+
+class ModelObject;
+
+namespace sla {
+
+/**
+  * The function should find the best rotation for SLA upside down printing.
+  *
+  * @param modelobj The model object representing the 3d mesh.
+  * @param accuracy The optimization accuracy from 0.0f to 1.0f. Currently,
+  * the nlopt genetic optimizer is used and the number of iterations is
+  * accuracy * 100000. This can change in the future.
+  * @param statuscb A status indicator callback called with the unsigned
+  * argument spanning from 0 to 100. May not reach 100 if the optimization finds
+  * an optimum before max iterations are reached.
+  * @param stopcond A function that if returns true, the search process will be
+  * terminated and the best solution found will be returned.
+  *
+  * @return Returns the rotations around each axis (x, y, z)
+  */
+std::array<double, 3> find_best_rotation(
+        const ModelObject& modelobj,
+        float accuracy = .0f,
+        std::function<void(unsigned)> statuscb = [] (unsigned) {},
+        std::function<bool()> stopcond = [] () { return false; }
+        );
+
+}
+}
+
+#endif // SLAROTFINDER_HPP

src/libslic3r/SLA/SLASpatIndex.hpp

@@ -0,0 +1,112 @@
+#ifndef SPATINDEX_HPP
+#define SPATINDEX_HPP
+
+#include <memory>
+#include <utility>
+#include <vector>
+
+#include <Eigen/Geometry>
+
+namespace Slic3r {
+namespace sla {
+
+typedef Eigen::Matrix<double,   3, 1, Eigen::DontAlign> Vec3d;
+using SpatElement = std::pair<Vec3d, unsigned>;
+
+/**
+ * This class is intended for enhancing range based for loops with indexing.
+ * So instead of this:
+ * { int i = 0; for(auto c : container) { process(c, i); ++i; }
+ *
+ * you can use this:
+ * for(auto ic : container) process(ic.value, ic.index);
+ */
+template<class Container> class Enumerable {
+    Container&& m;
+    using C = typename std::remove_reference<Container>::type;
+    using CC = typename std::remove_const<C>::type;
+
+    template<class S> struct get_iterator {};
+    template<> struct get_iterator<CC> { using type = typename CC::iterator; };
+    template<> struct get_iterator<const CC> {
+        using type = typename CC::const_iterator;
+    };
+
+    template<class Vref> struct _Pair {
+        Vref value;
+        size_t index;
+        _Pair(Vref v, size_t i) : value(v), index(i) {}
+        operator Vref() { return value; }
+    };
+
+    template<class Cit>
+    class _iterator {
+        Cit start;
+        Cit it;
+        using Pair = _Pair<typename std::iterator_traits<Cit>::reference>;
+    public:
+        _iterator(Cit b, Cit i): start(b), it(i) {}
+        _iterator& operator++() { ++it; return *this;}
+        _iterator operator++(int) { auto tmp = it; ++it; return tmp;}
+
+        bool operator!=(_iterator other) { return it != other.it; }
+        Pair operator*() { return Pair(*it, it - start); }
+        using value_type = typename Enumerable::value_type;
+    };
+
+public:
+
+    Enumerable(Container&& c): m(c) {}
+
+    using value_type = typename CC::value_type;
+
+    using iterator = _iterator<typename get_iterator<C>::type>;
+    using const_iterator = _iterator<typename CC::const_iterator>;
+
+    iterator begin() { return iterator(m.begin(), m.begin()); }
+    iterator end()   { return iterator(m.begin(), m.end()); }
+    const_iterator begin() const {return const_iterator(m.cbegin(), m.cbegin());}
+    const_iterator end() const { return const_iterator(m.cbegin(), m.cend());}
+
+};
+
+template<class C> inline Enumerable<C> enumerate(C&& c) {
+    return Enumerable<C>(std::forward<C>(c));
+}
+
+class SpatIndex {
+    class Impl;
+
+    // We use Pimpl because it takes a long time to compile boost headers which
+    // is the engine of this class. We include it only in the cpp file.
+    std::unique_ptr<Impl> m_impl;
+public:
+
+    SpatIndex();
+    ~SpatIndex();
+
+    SpatIndex(const SpatIndex&);
+    SpatIndex(SpatIndex&&);
+    SpatIndex& operator=(const SpatIndex&);
+    SpatIndex& operator=(SpatIndex&&);
+
+    void insert(const SpatElement&);
+    bool remove(const SpatElement&);
+
+    inline void insert(const Vec3d& v, unsigned idx)
+    {
+        insert(std::make_pair(v, unsigned(idx)));
+    }
+
+    std::vector<SpatElement> query(std::function<bool(const SpatElement&)>);
+    std::vector<SpatElement> nearest(const Vec3d&, unsigned k);
+
+    // For testing
+    size_t size() const;
+    bool empty() const { return size() == 0; }
+};
+
+}
+}
+
+#endif // SPATINDEX_HPP

src/libslic3r/SLA/SLASupportTree.hpp

@@ -0,0 +1,132 @@
+#ifndef SLASUPPORTTREE_HPP
+#define SLASUPPORTTREE_HPP
+
+#include <vector>
+#include <array>
+#include <cstdint>
+#include <memory>
+#include <Eigen/Geometry>
+
+namespace Slic3r {
+
+// Needed types from Point.hpp
+typedef int32_t coord_t;
+typedef Eigen::Matrix<double,   3, 1, Eigen::DontAlign> Vec3d;
+typedef Eigen::Matrix<float,    3, 1, Eigen::DontAlign> Vec3f;
+typedef Eigen::Matrix<coord_t,  3, 1, Eigen::DontAlign> Vec3crd;
+typedef std::vector<Vec3d>                              Pointf3s;
+typedef std::vector<Vec3crd>                            Points3;
+
+class TriangleMesh;
+class Model;
+class ModelInstance;
+class ExPolygon;
+
+using SliceLayer = std::vector<ExPolygon>;
+using SlicedSupports = std::vector<SliceLayer>;
+
+namespace sla {
+
+struct SupportConfig {
+    // Radius in mm of the pointing side of the head.
+    double head_front_radius_mm = 0.2;
+
+    // Radius of the back side of the 3d arrow.
+    double head_back_radius_mm = 0.5;
+
+    // Width in mm from the back sphere center to the front sphere center.
+    double head_width_mm = 1.0;
+
+    // Radius in mm of the support pillars.
+    // Warning: this value will be at most 65% of head_back_radius_mm
+    double pillar_radius_mm = 0.8;
+
+    // Radius in mm of the pillar base.
+    double base_radius_mm = 2.0;
+
+    // The height of the pillar base cone in mm.
+    double base_height_mm = 1.0;
+
+    // The default angle for connecting support sticks and junctions.
+    double tilt = M_PI/4;
+
+    // The max length of a bridge in mm
+    double max_bridge_length_mm = 15.0;
+};
+
+/// A Control structure for the support calculation. Consists of the status
+/// indicator callback and the stop condition predicate.
+struct Controller {
+    std::function<void(unsigned, const std::string&)> statuscb =
+            [](unsigned, const std::string&){};
+
+    std::function<bool(void)> stopcondition = [](){ return false; };
+};
+
+/* ************************************************************************** */
+/* TODO: May not be needed:                                                   */
+/* ************************************************************************** */
+
+void create_head(TriangleMesh&, double r1_mm, double r2_mm, double width_mm);
+
+/// Add support volumes to the model directly
+void add_sla_supports(Model& model, const SupportConfig& cfg = {},
+                      const Controller& ctl = {});
+
+/* ************************************************************************** */
+
+using PointSet = Eigen::MatrixXd;
+struct EigenMesh3D;
+
+/// Just a wrapper to the runtime error to be recognizable in try blocks
+class SLASupportsStoppedException: public std::runtime_error {
+public:
+    using std::runtime_error::runtime_error;
+    SLASupportsStoppedException(): std::runtime_error("") {}
+};
+
+/// The class containing mesh data for the generated supports.
+class SLASupportTree {
+    class Impl;
+    std::unique_ptr<Impl> m_impl;
+    Impl& get() { return *m_impl; }
+    const Impl& get() const { return *m_impl; }
+
+    friend void add_sla_supports(Model&,
+                                 const SupportConfig&,
+                                 const Controller&);
+
+    /// Generate the 3D supports for a model intended for SLA print.
+    bool generate(const PointSet& pts,
+                  const EigenMesh3D& mesh,
+                  const SupportConfig& cfg = {},
+                  const Controller& ctl = {});
+public:
+
+    // Constructors will throw if the stop condition becomes true.
+    SLASupportTree(const Model& model,
+                   const SupportConfig& cfg = {},
+                   const Controller& ctl = {});
+
+    SLASupportTree(const PointSet& pts,
+                   const EigenMesh3D& em,
+                   const SupportConfig& cfg = {},
+                   const Controller& ctl = {});
+
+    SLASupportTree(const SLASupportTree&);
+    SLASupportTree& operator=(const SLASupportTree&);
+
+    ~SLASupportTree();
+
+    /// Get the whole mesh united into the output TriangleMesh
+    void merged_mesh(TriangleMesh& outmesh) const;
+
+    /// Get the sliced 2d layers of the support geometry.
+    SlicedSupports slice() const;
+};
+
+}
+
+}
+
+#endif // SLASUPPORTTREE_HPP

src/libslic3r/SLA/SLASupportTreeIGL.cpp

@@ -0,0 +1,237 @@
+#include "SLA/SLASupportTree.hpp"
+#include "SLA/SLABoilerPlate.hpp"
+#include "SLA/SLASpatIndex.hpp"
+
+// HEAVY headers... takes eternity to compile
+
+// for concave hull merging decisions
+#include "SLABoostAdapter.hpp"
+#include "boost/geometry/index/rtree.hpp"
+
+#include <igl/ray_mesh_intersect.h>
+#include <igl/point_mesh_squared_distance.h>
+#include "SLASpatIndex.hpp"
+#include "ClipperUtils.hpp"
+
+namespace Slic3r {
+namespace sla {
+
+class SpatIndex::Impl {
+public:
+    using BoostIndex = boost::geometry::index::rtree< SpatElement,
+                       boost::geometry::index::rstar<16, 4> /* ? */ >;
+
+    BoostIndex m_store;
+};
+
+SpatIndex::SpatIndex(): m_impl(new Impl()) {}
+SpatIndex::~SpatIndex() {}
+
+SpatIndex::SpatIndex(const SpatIndex &cpy): m_impl(new Impl(*cpy.m_impl)) {}
+SpatIndex::SpatIndex(SpatIndex&& cpy): m_impl(std::move(cpy.m_impl)) {}
+
+SpatIndex& SpatIndex::operator=(const SpatIndex &cpy)
+{
+    m_impl.reset(new Impl(*cpy.m_impl));
+    return *this;
+}
+
+SpatIndex& SpatIndex::operator=(SpatIndex &&cpy)
+{
+    m_impl.swap(cpy.m_impl);
+    return *this;
+}
+
+void SpatIndex::insert(const SpatElement &el)
+{
+    m_impl->m_store.insert(el);
+}
+
+bool SpatIndex::remove(const SpatElement& el)
+{
+    return m_impl->m_store.remove(el);
+}
+
+std::vector<SpatElement>
+SpatIndex::query(std::function<bool(const SpatElement &)> fn)
+{
+    namespace bgi = boost::geometry::index;
+
+    std::vector<SpatElement> ret;
+    m_impl->m_store.query(bgi::satisfies(fn), std::back_inserter(ret));
+    return ret;
+}
+
+std::vector<SpatElement> SpatIndex::nearest(const Vec3d &el, unsigned k = 1)
+{
+    namespace bgi = boost::geometry::index;
+    std::vector<SpatElement> ret; ret.reserve(k);
+    m_impl->m_store.query(bgi::nearest(el, k), std::back_inserter(ret));
+    return ret;
+}
+
+size_t SpatIndex::size() const
+{
+    return m_impl->m_store.size();
+}
+
+PointSet normals(const PointSet& points, const EigenMesh3D& mesh) {
+    Eigen::VectorXd dists;
+    Eigen::VectorXi I;
+    PointSet C;
+    igl::point_mesh_squared_distance( points, mesh.V, mesh.F, dists, I, C);
+
+    PointSet ret(I.rows(), 3);
+    for(int i = 0; i < I.rows(); i++) {
+        auto idx = I(i);
+        auto trindex = mesh.F.row(idx);
+
+        auto& p1 = mesh.V.row(trindex(0));
+        auto& p2 = mesh.V.row(trindex(1));
+        auto& p3 = mesh.V.row(trindex(2));
+
+        Eigen::Vector3d U = p2 - p1;
+        Eigen::Vector3d V = p3 - p1;
+        ret.row(i) = U.cross(V).normalized();
+    }
+
+    return ret;
+}
+
+double ray_mesh_intersect(const Vec3d& s,
+                          const Vec3d& dir,
+                          const EigenMesh3D& m)
+{
+    igl::Hit hit;
+    hit.t = std::numeric_limits<float>::infinity();
+    igl::ray_mesh_intersect(s, dir, m.V, m.F, hit);
+    return hit.t;
+}
+
+// Clustering a set of points by the given criteria
+ClusteredPoints cluster(
+        const sla::PointSet& points,
+        std::function<bool(const SpatElement&, const SpatElement&)> pred,
+        unsigned max_points = 0)
+{
+
+    namespace bgi = boost::geometry::index;
+    using Index3D = bgi::rtree< SpatElement, bgi::rstar<16, 4> /* ? */ >;
+
+    // A spatial index for querying the nearest points
+    Index3D sindex;
+
+    // Build the index
+    for(unsigned idx = 0; idx < points.rows(); idx++)
+        sindex.insert( std::make_pair(points.row(idx), idx));
+
+    using Elems = std::vector<SpatElement>;
+
+    // Recursive function for visiting all the points in a given distance to
+    // each other
+    std::function<void(Elems&, Elems&)> group =
+    [&sindex, &group, pred, max_points](Elems& pts, Elems& cluster)
+    {
+        for(auto& p : pts) {
+            std::vector<SpatElement> tmp;
+
+            sindex.query(
+                bgi::satisfies([p, pred](const SpatElement& se) {
+                    return pred(p, se);
+                }),
+                std::back_inserter(tmp)
+            );
+
+            auto cmp = [](const SpatElement& e1, const SpatElement& e2){
+                return e1.second < e2.second;
+            };
+
+            std::sort(tmp.begin(), tmp.end(), cmp);
+
+            Elems newpts;
+            std::set_difference(tmp.begin(), tmp.end(),
+                                cluster.begin(), cluster.end(),
+                                std::back_inserter(newpts), cmp);
+
+            int c = max_points && newpts.size() + cluster.size() > max_points?
+                        int(max_points - cluster.size()) : int(newpts.size());
+
+            cluster.insert(cluster.end(), newpts.begin(), newpts.begin() + c);
+            std::sort(cluster.begin(), cluster.end(), cmp);
+
+            if(!newpts.empty() && (!max_points || cluster.size() < max_points))
+                group(newpts, cluster);
+        }
+    };
+
+    std::vector<Elems> clusters;
+    for(auto it = sindex.begin(); it != sindex.end();) {
+        Elems cluster = {};
+        Elems pts = {*it};
+        group(pts, cluster);
+
+        for(auto& c : cluster) sindex.remove(c);
+        it = sindex.begin();
+
+        clusters.emplace_back(cluster);
+    }
+
+    ClusteredPoints result;
+    for(auto& cluster : clusters) {
+        result.emplace_back();
+        for(auto c : cluster) result.back().emplace_back(c.second);
+    }
+
+    return result;
+}
+
+using Segments = std::vector<std::pair<Vec2d, Vec2d>>;
+
+Segments model_boundary(const EigenMesh3D& emesh, double offs)
+{
+    Segments ret;
+    Polygons pp;
+    pp.reserve(emesh.F.rows());
+
+    for (int i = 0; i < emesh.F.rows(); i++) {
+        auto trindex = emesh.F.row(i);
+        auto& p1 = emesh.V.row(trindex(0));
+        auto& p2 = emesh.V.row(trindex(1));
+        auto& p3 = emesh.V.row(trindex(2));
+
+        Polygon p;
+        p.points.resize(3);
+        p.points[0] = Point::new_scale(p1(X), p1(Y));
+        p.points[1] = Point::new_scale(p2(X), p2(Y));
+        p.points[2] = Point::new_scale(p3(X), p3(Y));
+        p.make_counter_clockwise();
+        pp.emplace_back(p);
+    }
+
+    ExPolygons merged = union_ex(offset(pp, float(scale_(offs))), true);
+
+    for(auto& expoly : merged) {
+        auto lines = expoly.lines();
+        for(Line& l : lines) {
+            Vec2d a(l.a(X) * SCALING_FACTOR, l.a(Y) * SCALING_FACTOR);
+            Vec2d b(l.b(X) * SCALING_FACTOR, l.b(Y) * SCALING_FACTOR);
+            ret.emplace_back(std::make_pair(a, b));
+        }
+    }
+
+    return ret;
+}
+
+//struct SegmentIndex {
+
+//};
+
+//using SegmentIndexEl = std::pair<Segment, unsigned>;
+
+//SegmentIndexEl
+
+
+
+
+}
+}

tests/CMakeLists.txt

@@ -0,0 +1,3 @@
+# TODO Add individual tests as executables in separate directories
+
+# add_subirectory(<testcase>)