Merge remote-tracking branch 'remotes/prusa/master'

WIP, in particular, have to review the custom beds svg/stl

8 files changed, 1272 insertions, 536 deletions

diff --git a/src/libslic3r/SLA/SLAAutoSupports.cpp b/src/libslic3r/SLA/SLAAutoSupports.cpp
index 97c6d61f5..cbd48796d 100644
--- a/src/libslic3r/SLA/SLAAutoSupports.cpp
+++ b/src/libslic3r/SLA/SLAAutoSupports.cpp
@@ -59,8 +59,6 @@ SLAAutoSupports::SLAAutoSupports(const TriangleMesh& mesh, const sla::EigenMesh3
 void SLAAutoSupports::project_onto_mesh(std::vector<sla::SupportPoint>& points) const
 {
     // The function  makes sure that all the points are really exactly placed on the mesh.
-    igl::Hit hit_up{0, 0, 0.f, 0.f, 0.f};
-    igl::Hit hit_down{0, 0, 0.f, 0.f, 0.f};
 
     // Use a reasonable granularity to account for the worker thread synchronization cost.
     tbb::parallel_for(tbb::blocked_range<size_t>(0, points.size(), 64),
@@ -140,7 +138,6 @@ static std::vector<SLAAutoSupports::MyLayer> make_layers(
 			SLAAutoSupports::MyLayer &layer_above = layers[layer_id];
 			SLAAutoSupports::MyLayer &layer_below = layers[layer_id - 1];
             //FIXME WTF?
-            const float height = (layer_id>2 ? heights[layer_id-3] : heights[0]-(heights[1]-heights[0]));
             const float layer_height = (layer_id!=0 ? heights[layer_id]-heights[layer_id-1] : heights[0]);
             const float safe_angle = 5.f * (float(M_PI)/180.f); // smaller number - less supports
             const double between_layers_offset =  double(scale_(layer_height / std::tan(safe_angle)));
@@ -212,7 +209,7 @@ void SLAAutoSupports::process(const std::vector<ExPolygons>& slices, const std::
         for (Structure &top : layer_top->islands)
 			for (Structure::Link &bottom_link : top.islands_below) {
                 Structure &bottom = *bottom_link.island;
-                float centroids_dist = (bottom.centroid - top.centroid).norm();
+                //float centroids_dist = (bottom.centroid - top.centroid).norm();
                 // Penalization resulting from centroid offset:
 //                  bottom.supports_force *= std::min(1.f, 1.f - std::min(1.f, (1600.f * layer_height) * centroids_dist * centroids_dist / bottom.area));
                 float &support_force = support_force_bottom[&bottom - layer_bottom->islands.data()];
@@ -239,7 +236,7 @@ void SLAAutoSupports::process(const std::vector<ExPolygons>& slices, const std::
 //            s.supports_force_inherited /= std::max(1.f, (layer_height / 0.3f) * e_area / s.area);
             s.supports_force_inherited /= std::max(1.f, 0.17f * (s.overhangs_area) / s.area);
 
-            float force_deficit = s.support_force_deficit(m_config.tear_pressure());
+            //float force_deficit = s.support_force_deficit(m_config.tear_pressure());
             if (s.islands_below.empty()) { // completely new island - needs support no doubt
                 uniformly_cover({ *s.polygon }, s, point_grid, true);
             } else if (! s.dangling_areas.empty()) {
@@ -380,7 +377,7 @@ static inline std::vector<Vec2f> poisson_disk_from_samples(const std::vector<Vec
     {
         typename Cells::iterator last_cell_id_it;
         Vec2i32           last_cell_id(-1, -1);
-        for (int i = 0; i < raw_samples_sorted.size(); ++ i) {
+        for (size_t i = 0; i < raw_samples_sorted.size(); ++ i) {
             const RawSample &sample = raw_samples_sorted[i];
             if (sample.cell_id == last_cell_id) {
                 // This sample is in the same cell as the previous, so just increase the count.  Cells are
@@ -389,7 +386,7 @@ static inline std::vector<Vec2f> poisson_disk_from_samples(const std::vector<Vec
             } else {
                 // This is a new cell.
                 PoissonDiskGridEntry data;
-                data.first_sample_idx = i;
+                data.first_sample_idx = int(i);
                 data.sample_cnt       = 1;
                 auto result     = cells.insert({sample.cell_id, data});
                 last_cell_id    = sample.cell_id;
diff --git a/src/libslic3r/SLA/SLABasePool.cpp b/src/libslic3r/SLA/SLABasePool.cpp
index 14b4a2181..c337a5ac5 100644
--- a/src/libslic3r/SLA/SLABasePool.cpp
+++ b/src/libslic3r/SLA/SLABasePool.cpp
@@ -8,9 +8,9 @@
 #include "MTUtils.hpp"
 
 // For debugging:
-//#include <fstream>
-//#include <libnest2d/tools/benchmark.h>
-//#include "SVG.hpp"
+// #include <fstream>
+// #include <libnest2d/tools/benchmark.h>
+// #include "SVG.hpp"
 
 namespace Slic3r { namespace sla {
 
@@ -184,9 +184,10 @@ Contour3D walls(const Polygon& lower, const Polygon& upper,
 }
 
 /// Offsetting with clipper and smoothing the edges into a curvature.
-void offset(ExPolygon& sh, coord_t distance) {
+void offset(ExPolygon& sh, coord_t distance, bool edgerounding = true) {
     using ClipperLib::ClipperOffset;
     using ClipperLib::jtRound;
+    using ClipperLib::jtMiter;
     using ClipperLib::etClosedPolygon;
     using ClipperLib::Paths;
     using ClipperLib::Path;
@@ -203,11 +204,13 @@ void offset(ExPolygon& sh, coord_t distance) {
         return;
     }
 
+    auto jointype = edgerounding? jtRound : jtMiter;
+    
     ClipperOffset offs;
     offs.ArcTolerance = scaled<double>(0.01);
     Paths result;
-    offs.AddPath(ctour, jtRound, etClosedPolygon);
-    offs.AddPaths(holes, jtRound, etClosedPolygon);
+    offs.AddPath(ctour, jointype, etClosedPolygon);
+    offs.AddPaths(holes, jointype, etClosedPolygon);
     offs.Execute(result, static_cast<double>(distance));
 
     // Offsetting reverts the orientation and also removes the last vertex
@@ -237,6 +240,50 @@ void offset(ExPolygon& sh, coord_t distance) {
     }
 }
 
+void offset(Polygon &sh, coord_t distance, bool edgerounding = true)
+{
+    using ClipperLib::ClipperOffset;
+    using ClipperLib::jtRound;
+    using ClipperLib::jtMiter;
+    using ClipperLib::etClosedPolygon;
+    using ClipperLib::Paths;
+    using ClipperLib::Path;
+
+    auto &&ctour = Slic3rMultiPoint_to_ClipperPath(sh);
+
+    // If the input is not at least a triangle, we can not do this algorithm
+    if (ctour.size() < 3) {
+        BOOST_LOG_TRIVIAL(error) << "Invalid geometry for offsetting!";
+        return;
+    }
+
+    ClipperOffset offs;
+    offs.ArcTolerance = 0.01 * scaled(1.);
+    Paths result;
+    offs.AddPath(ctour, edgerounding ? jtRound : 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) {
+                auto rr = ClipperPath_to_Slic3rPolygon(r);
+                sh.points.swap(rr.points);
+                found_the_contour = true;
+            } else {
+                BOOST_LOG_TRIVIAL(warning)
+                    << "Warning: offsetting result is invalid!";
+            }
+        }
+    }
+}
+
 /// Unification of polygons (with clipper) preserving holes as well.
 ExPolygons unify(const ExPolygons& shapes) {
     using ClipperLib::ptSubject;
@@ -307,6 +354,116 @@ ExPolygons unify(const ExPolygons& shapes) {
     return retv;
 }
 
+Polygons unify(const Polygons& shapes) {
+    using ClipperLib::ptSubject;
+    
+    bool closed = true;
+    bool valid = true;
+
+    ClipperLib::Clipper clipper;
+
+    for(auto& path : shapes) {
+        auto clipperpath = Slic3rMultiPoint_to_ClipperPath(path);
+
+        if(!clipperpath.empty())
+            valid &= clipper.AddPath(clipperpath, ptSubject, closed);
+    }
+
+    if(!valid) BOOST_LOG_TRIVIAL(warning) << "Unification of invalid shapes!";
+
+    ClipperLib::Paths result;
+    clipper.Execute(ClipperLib::ctUnion, result, ClipperLib::pftNonZero);
+
+    Polygons ret;
+    for (ClipperLib::Path &p : result) {
+        Polygon pp = ClipperPath_to_Slic3rPolygon(p);
+        if (!pp.is_clockwise()) ret.emplace_back(std::move(pp));
+    }
+    
+    return ret;
+}
+
+// Function to cut tiny connector cavities for a given polygon. The input poly
+// will be offsetted by "padding" and small rectangle shaped cavities will be
+// inserted along the perimeter in every "stride" distance. The stick rectangles
+// will have a with about "stick_width". The input dimensions are in world 
+// measure, not the scaled clipper units.
+void breakstick_holes(ExPolygon& poly,
+                      double padding,
+                      double stride,
+                      double stick_width,
+                      double penetration)
+{   
+    // SVG svg("bridgestick_plate.svg");
+    // svg.draw(poly);
+
+    auto transf = [stick_width, penetration, padding, stride](Points &pts) {
+        // The connector stick will be a small rectangle with dimensions
+        // stick_width x (penetration + padding) to have some penetration
+        // into the input polygon.
+
+        Points out;
+        out.reserve(2 * pts.size()); // output polygon points
+
+        // stick bottom and right edge dimensions
+        double sbottom = scaled(stick_width);
+        double sright  = scaled(penetration + padding);
+
+        // scaled stride distance
+        double sstride = scaled(stride);
+        double t       = 0;
+
+        // process pairs of vertices as an edge, start with the last and
+        // first point
+        for (size_t i = pts.size() - 1, j = 0; j < pts.size(); i = j, ++j) {
+            // Get vertices and the direction vectors
+            const Point &a = pts[i], &b = pts[j];
+            Vec2d        dir = b.cast<double>() - a.cast<double>();
+            double       nrm = dir.norm();
+            dir /= nrm;
+            Vec2d dirp(-dir(Y), dir(X));
+
+            // Insert start point
+            out.emplace_back(a);
+
+            // dodge the start point, do not make sticks on the joins
+            while (t < sbottom) t += sbottom;
+            double tend = nrm - sbottom;
+
+            while (t < tend) { // insert the stick on the polygon perimeter
+
+                // calculate the stick rectangle vertices and insert them
+                // into the output.
+                Point p1 = a + (t * dir).cast<coord_t>();
+                Point p2 = p1 + (sright * dirp).cast<coord_t>();
+                Point p3 = p2 + (sbottom * dir).cast<coord_t>();
+                Point p4 = p3 + (sright * -dirp).cast<coord_t>();
+                out.insert(out.end(), {p1, p2, p3, p4});
+
+                // continue along the perimeter
+                t += sstride;
+            }
+
+            t = t - nrm;
+
+            // Insert edge endpoint
+            out.emplace_back(b);
+        }
+        
+        // move the new points
+        out.shrink_to_fit();
+        pts.swap(out);
+    };
+    
+    if(stride > 0.0 && stick_width > 0.0 && padding > 0.0) {
+        transf(poly.contour.points);
+        for (auto &h : poly.holes) transf(h.points);
+    }
+    
+    // svg.draw(poly);
+    // svg.Close();
+}
+
 /// This method will create a rounded edge around a flat polygon in 3d space.
 /// 'base_plate' parameter is the target plate.
 /// 'radius' is the radius of the edges.
@@ -426,41 +583,37 @@ inline Point centroid(Points& pp) {
     return c;
 }
 
-inline Point centroid(const ExPolygon& poly) {
-    return poly.contour.centroid();
+inline Point centroid(const Polygon& poly) {
+    return poly.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,
-                        ThrowOnCancel throw_on_cancel = [](){})
+Polygons concave_hull(const Polygons& polys, double maxd_mm, ThrowOnCancel thr)
 {
     namespace bgi = boost::geometry::index;
-    using SpatElement = std::pair<BoundingBox, unsigned>;
+    using SpatElement = std::pair<Point, unsigned>;
     using SpatIndex = bgi::rtree< SpatElement, bgi::rstar<16, 4> >;
 
-    if(polys.empty()) return ExPolygons();
+    if(polys.empty()) return Polygons();
+    
+    const double max_dist = scaled(maxd_mm);
 
-    ExPolygons punion = unify(polys);   // could be redundant
+    Polygons 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++));
-    });
-
+                   [](const Polygon& poly) { return centroid(poly); });
 
+    SpatIndex ctrindex;
+    unsigned  idx = 0;
+    for(const Point &ct : centroids) ctrindex.insert(std::make_pair(ct, idx++));
+    
     // Centroid of the centroids of islands. This is where the additional
     // connector sticks are routed.
     Point cc = centroid(centroids);
@@ -470,25 +623,32 @@ ExPolygons concave_hull(const ExPolygons& polys, double max_dist_mm = 50,
     idx = 0;
     std::transform(centroids.begin(), centroids.end(),
                    std::back_inserter(punion),
-                   [&punion, &boxindex, cc, max_dist_mm, &idx, throw_on_cancel]
+                   [&centroids, &ctrindex, cc, max_dist, &idx, thr]
                    (const Point& c)
     {
-        throw_on_cancel();
+        thr();
         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 = scaled<double>(max_dist_mm);
-
-        ExPolygon& expo = punion[idx++];
-        BoundingBox querybb(expo);
-
-        querybb.offset(max_dist);
+        
+        Point& ct = centroids[idx];
+        
         std::vector<SpatElement> result;
-        boxindex.query(bgi::intersects(querybb), std::back_inserter(result));
-        if(result.size() <= 1) return ExPolygon();
+        ctrindex.query(bgi::nearest(ct, 2), std::back_inserter(result));
 
-        ExPolygon r;
-        auto& ctour = r.contour.points;
+        double dist = max_dist;
+        for (const SpatElement &el : result)
+            if (el.second != idx) {
+                dist = Line(el.first, ct).length();
+                break;
+            }
+        
+        idx++;
+        
+        if (dist >= max_dist) return Polygon();
+        
+        Polygon r;
+        auto& ctour = r.points;
 
         ctour.reserve(3);
         ctour.emplace_back(cc);
@@ -507,24 +667,20 @@ ExPolygons concave_hull(const ExPolygons& polys, double max_dist_mm = 50,
     return punion;
 }
 
-void base_plate(const TriangleMesh &mesh, ExPolygons &output, float h,
-                float layerh, ThrowOnCancel thrfn)
+void base_plate(const TriangleMesh &      mesh,
+                ExPolygons &              output,
+                const std::vector<float> &heights,
+                ThrowOnCancel             thrfn)
 {
-    TriangleMesh m = mesh;
-    m.require_shared_vertices(); // TriangleMeshSlicer needs this
-    TriangleMeshSlicer slicer(&m);
-
-    auto bb = mesh.bounding_box();
-    float gnd = float(bb.min(Z));
-    std::vector<float> heights = {float(bb.min(Z))};
-    for(float hi = gnd + layerh; hi <= gnd + h; hi += layerh)
-        heights.emplace_back(hi);
-
-    std::vector<ExPolygons> out; out.reserve(size_t(std::ceil(h/layerh)));
+    if (mesh.empty()) return;
+    //    m.require_shared_vertices(); // TriangleMeshSlicer needs this    
+    TriangleMeshSlicer slicer(&mesh);
+    
+    std::vector<ExPolygons> out; out.reserve(heights.size());
     slicer.slice(heights, &out, thrfn);
-
+    
     size_t count = 0; for(auto& o : out) count += o.size();
-
+    
     // Now we have to unify all slice layers which can be an expensive operation
     // so we will try to simplify the polygons
     ExPolygons tmp; tmp.reserve(count);
@@ -533,16 +689,33 @@ void base_plate(const TriangleMesh &mesh, ExPolygons &output, float h,
             auto&& exss = e.simplify(scaled<double>(0.1));
             for(ExPolygon& ep : exss) tmp.emplace_back(std::move(ep));
         }
-
+    
     ExPolygons utmp = unify(tmp);
-
+    
     for(auto& o : utmp) {
         auto&& smp = o.simplify(scaled<double>(0.1));
         output.insert(output.end(), smp.begin(), smp.end());
     }
 }
 
-Contour3D create_base_pool(const ExPolygons &ground_layer, 
+void base_plate(const TriangleMesh &mesh,
+                ExPolygons &        output,
+                float               h,
+                float               layerh,
+                ThrowOnCancel       thrfn)
+{
+    auto bb = mesh.bounding_box();
+    float gnd = float(bb.min(Z));
+    std::vector<float> heights = {float(bb.min(Z))};
+    
+    for(float hi = gnd + layerh; hi <= gnd + h; hi += layerh)
+        heights.emplace_back(hi);
+    
+    base_plate(mesh, output, heights, thrfn);
+}
+
+Contour3D create_base_pool(const Polygons &ground_layer, 
+                           const ExPolygons &obj_self_pad = {},
                            const PoolConfig& cfg = PoolConfig()) 
 {
     // for debugging:
@@ -557,7 +730,7 @@ Contour3D create_base_pool(const ExPolygons &ground_layer,
     // serve as the bottom plate of the pad. We will offset this concave hull
     // and then offset back the result with clipper with rounding edges ON. This
     // trick will create a nice rounded pad shape.
-    ExPolygons concavehs = concave_hull(ground_layer, mergedist, cfg.throw_on_cancel);
+    Polygons concavehs = concave_hull(ground_layer, mergedist, cfg.throw_on_cancel);
 
     const double thickness      = cfg.min_wall_thickness_mm;
     const double wingheight     = cfg.min_wall_height_mm;
@@ -577,42 +750,37 @@ Contour3D create_base_pool(const ExPolygons &ground_layer,
 
     Contour3D pool;
 
-    for(ExPolygon& concaveh : concavehs) {
-        if(concaveh.contour.points.empty()) return pool;
-
-        // Get rid of any holes in the concave hull output.
-        concaveh.holes.clear();
+    for(Polygon& concaveh : concavehs) {
+        if(concaveh.points.empty()) return pool;
 
         // Here lies the trick that does the smoothing only with clipper offset
         // calls. The offset is configured to round edges. Inner edges will
         // be rounded because we offset twice: ones to get the outer (top) plate
         // and again to get the inner (bottom) plate
         auto outer_base = concaveh;
-        outer_base.holes.clear();
         offset(outer_base, s_safety_dist + s_wingdist + s_thickness);
 
-        ExPolygon bottom_poly = outer_base;
-        bottom_poly.holes.clear();
+        ExPolygon bottom_poly; bottom_poly.contour = outer_base;
         offset(bottom_poly, -s_bottom_offs);
 
         // Punching a hole in the top plate for the cavity
         ExPolygon top_poly;
         ExPolygon middle_base;
         ExPolygon inner_base;
-        top_poly.contour = outer_base.contour;
+        top_poly.contour = outer_base;
 
         if(wingheight > 0) {
-            inner_base = outer_base;
+            inner_base.contour = outer_base;
             offset(inner_base, -(s_thickness + s_wingdist + s_eradius));
 
-            middle_base = outer_base;
+            middle_base.contour = outer_base;
             offset(middle_base, -s_thickness);
             top_poly.holes.emplace_back(middle_base.contour);
             auto& tph = top_poly.holes.back().points;
             std::reverse(tph.begin(), tph.end());
         }
 
-        ExPolygon ob = outer_base; double wh = 0;
+        ExPolygon ob; ob.contour = 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.
@@ -659,6 +827,7 @@ Contour3D create_base_pool(const ExPolygons &ground_layer,
         if(wingheight > 0) {
             // Generate the smoothed edge geometry
             wh = 0;
+            ob = middle_base;
             if(s_eradius) pool.merge(round_edges(middle_base,
                                    r,
                                    phi - 90, // from tangent lines
@@ -673,11 +842,59 @@ Contour3D create_base_pool(const ExPolygons &ground_layer,
                              wh, -wingdist, thrcl));
         }
 
-        // Now we need to triangulate the top and bottom plates as well as the
-        // cavity bottom plate which is the same as the bottom plate but it is
-        // elevated by the thickness.
+        if (cfg.embed_object) {
+            ExPolygons bttms = diff_ex(to_polygons(bottom_poly),
+                                       to_polygons(obj_self_pad));
+            
+            assert(!bttms.empty());
+            
+            std::sort(bttms.begin(), bttms.end(),
+                      [](const ExPolygon& e1, const ExPolygon& e2) {
+                          return e1.contour.area() > e2.contour.area();
+                      });
+            
+            if(wingheight > 0) inner_base.holes = bttms.front().holes;
+            else top_poly.holes = bttms.front().holes;
+
+            auto straight_walls =
+                [&pool](const Polygon &cntr, coord_t z_low, coord_t z_high) {
+                    
+                auto lines = cntr.lines();
+                
+                for (auto &l : lines) {
+                    auto s = coord_t(pool.points.size());
+                    auto& pts = pool.points;
+                    pts.emplace_back(unscale(l.a.x(), l.a.y(), z_low));
+                    pts.emplace_back(unscale(l.b.x(), l.b.y(), z_low));
+                    pts.emplace_back(unscale(l.a.x(), l.a.y(), z_high));
+                    pts.emplace_back(unscale(l.b.x(), l.b.y(), z_high));
+                    
+                    pool.indices.emplace_back(s, s + 1, s + 3);
+                    pool.indices.emplace_back(s, s + 3, s + 2);
+                }
+            };
+            
+            coord_t z_lo = -scaled(fullheight), z_hi = -scaled(wingheight);
+            for (ExPolygon &ep : bttms) {
+                pool.merge(triangulate_expolygon_3d(ep, -fullheight, true));
+                for (auto &h : ep.holes) straight_walls(h, z_lo, z_hi);
+            }
+            
+            // Skip the outer contour, triangulate the holes
+            for (auto it = std::next(bttms.begin()); it != bttms.end(); ++it) {
+                pool.merge(triangulate_expolygon_3d(*it, -wingheight));
+                straight_walls(it->contour, z_lo, z_hi);
+            }
+            
+        } else {
+            // Now we need to triangulate the top and bottom plates as well as
+            // the cavity bottom plate which is the same as the bottom plate
+            // but it is elevated by the thickness.
+            
+            pool.merge(triangulate_expolygon_3d(bottom_poly, -fullheight, true));
+        }
+        
         pool.merge(triangulate_expolygon_3d(top_poly));
-        pool.merge(triangulate_expolygon_3d(bottom_poly, -fullheight, true));
 
         if(wingheight > 0)
             pool.merge(triangulate_expolygon_3d(inner_base, -wingheight));
@@ -687,8 +904,8 @@ Contour3D create_base_pool(const ExPolygons &ground_layer,
     return pool;
 }
 
-void create_base_pool(const ExPolygons &ground_layer, TriangleMesh& out,
-                      const PoolConfig& cfg)
+void create_base_pool(const Polygons &ground_layer, TriangleMesh& out,
+                      const ExPolygons &holes, const PoolConfig& cfg)
 {
     
 
@@ -698,7 +915,7 @@ void create_base_pool(const ExPolygons &ground_layer, TriangleMesh& out,
     // std::fstream fout("pad_debug.obj", std::fstream::out);
     // if(fout.good()) pool.to_obj(fout);
 
-    out.merge(mesh(create_base_pool(ground_layer, cfg)));
+    out.merge(mesh(create_base_pool(ground_layer, holes, cfg)));
 }
 
 }
diff --git a/src/libslic3r/SLA/SLABasePool.hpp b/src/libslic3r/SLA/SLABasePool.hpp
index 3c88e58c8..eec426bbf 100644
--- a/src/libslic3r/SLA/SLABasePool.hpp
+++ b/src/libslic3r/SLA/SLABasePool.hpp
@@ -8,7 +8,9 @@
 namespace Slic3r {
 
 class ExPolygon;
+class Polygon;
 using ExPolygons = std::vector<ExPolygon>;
+using Polygons = std::vector<Polygon>;
 
 class TriangleMesh;
 
@@ -19,16 +21,43 @@ using ThrowOnCancel = std::function<void(void)>;
 /// Calculate the polygon representing the silhouette from the specified height
 void base_plate(const TriangleMesh& mesh,       // input mesh
                 ExPolygons& output,             // Output will be merged with
-                float zlevel = 0.1f,            // Plate creation level
+                float samplingheight = 0.1f,    // The height range to sample
                 float layerheight = 0.05f,      // The sampling height
                 ThrowOnCancel thrfn = [](){});  // Will be called frequently
 
+void base_plate(const TriangleMesh& mesh,       // input mesh
+                ExPolygons& output,             // Output will be merged with
+                const std::vector<float>&,      // Exact Z levels to sample
+                ThrowOnCancel thrfn = [](){});  // Will be called frequently
+
+// Function to cut tiny connector cavities for a given polygon. The input poly
+// will be offsetted by "padding" and small rectangle shaped cavities will be
+// inserted along the perimeter in every "stride" distance. The stick rectangles
+// will have a with about "stick_width". The input dimensions are in world 
+// measure, not the scaled clipper units.
+void breakstick_holes(ExPolygon &poly,
+                      double     padding,
+                      double     stride,
+                      double     stick_width,
+                      double     penetration = 0.0);
+
+Polygons concave_hull(const Polygons& polys, double max_dist_mm = 50,
+                      ThrowOnCancel throw_on_cancel = [](){});
+
 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;
     double wall_slope = std::atan(1.0);          // Universal constant for Pi/4
+    struct EmbedObject {
+        double object_gap_mm = 0.5;
+        double stick_stride_mm = 10;
+        double stick_width_mm = 0.3;
+        double stick_penetration_mm = 0.1;
+        bool enabled = false;
+        operator bool() const { return enabled; }
+    } embed_object;
 
     ThrowOnCancel throw_on_cancel = [](){};
 
@@ -42,15 +71,12 @@ struct PoolConfig {
 };
 
 /// Calculate the pool for the mesh for SLA printing
-void create_base_pool(const ExPolygons& base_plate,
+void create_base_pool(const Polygons& base_plate,
                       TriangleMesh& output_mesh,
+                      const ExPolygons& holes,
                       const PoolConfig& = PoolConfig());
 
-/// TODO: Currently the base plate of the pool will have half the height of the
-/// whole pool. So the carved out space has also half the height. This is not
-/// a particularly elegant solution, the thickness should be exactly
-/// min_wall_thickness and it should be corrected in the future. This method
-/// will return the correct value for further processing.
+/// Returns the elevation needed for compensating the pad.
 inline double get_pad_elevation(const PoolConfig& cfg) {
     return cfg.min_wall_thickness_mm;
 }
diff --git a/src/libslic3r/SLA/SLACommon.hpp b/src/libslic3r/SLA/SLACommon.hpp
index 855802759..d8e258035 100644
--- a/src/libslic3r/SLA/SLACommon.hpp
+++ b/src/libslic3r/SLA/SLACommon.hpp
@@ -18,49 +18,65 @@ namespace sla {
     
 // An enum to keep track of where the current points on the ModelObject came from.
 enum class PointsStatus {
-    None,           // No points were generated so far.
+    NoPoints,           // No points were generated so far.
     Generating,     // The autogeneration algorithm triggered, but not yet finished.
     AutoGenerated,  // Points were autogenerated (i.e. copied from the backend).
     UserModified    // User has done some edits.
 };
 
-struct SupportPoint {
+struct SupportPoint
+{
     Vec3f pos;
     float head_front_radius;
-    bool is_new_island;
-
-    SupportPoint() :
-    pos(Vec3f::Zero()), head_front_radius(0.f), is_new_island(false) {}
-
-    SupportPoint(float pos_x, float pos_y, float pos_z, float head_radius, bool new_island) :
-    pos(pos_x, pos_y, pos_z), head_front_radius(head_radius), is_new_island(new_island) {}
-
-    SupportPoint(Vec3f position, float head_radius, bool new_island) :
-    pos(position), head_front_radius(head_radius), is_new_island(new_island) {}
-
-    SupportPoint(Eigen::Matrix<float, 5, 1, Eigen::DontAlign> data) :
-    pos(data(0), data(1), data(2)), head_front_radius(data(3)), is_new_island(data(4) != 0.f) {}
-
-    bool operator==(const SupportPoint& sp) const { return (pos==sp.pos) && head_front_radius==sp.head_front_radius && is_new_island==sp.is_new_island; }
-    bool operator!=(const SupportPoint& sp) const { return !(sp == (*this)); }
+    bool  is_new_island;
+
+    SupportPoint()
+        : pos(Vec3f::Zero()), head_front_radius(0.f), is_new_island(false)
+    {}
+
+    SupportPoint(float pos_x,
+                 float pos_y,
+                 float pos_z,
+                 float head_radius,
+                 bool  new_island)
+        : pos(pos_x, pos_y, pos_z)
+        , head_front_radius(head_radius)
+        , is_new_island(new_island)
+    {}
+
+    SupportPoint(Vec3f position, float head_radius, bool new_island)
+        : pos(position)
+        , head_front_radius(head_radius)
+        , is_new_island(new_island)
+    {}
+
+    SupportPoint(Eigen::Matrix<float, 5, 1, Eigen::DontAlign> data)
+        : pos(data(0), data(1), data(2))
+        , head_front_radius(data(3))
+        , is_new_island(data(4) != 0.f)
+    {}
+
+    bool operator==(const SupportPoint &sp) const
+    {
+        return (pos == sp.pos) && head_front_radius == sp.head_front_radius &&
+               is_new_island == sp.is_new_island;
+    }
+    bool operator!=(const SupportPoint &sp) const { return !(sp == (*this)); }
+
+    template<class Archive> void serialize(Archive &ar)
+    {
+        ar(pos, head_front_radius, is_new_island);
+    }
 };
 
 /// An index-triangle structure for libIGL functions. Also serves as an
 /// alternative (raw) input format for the SLASupportTree
-/*struct EigenMesh3D {
-    Eigen::MatrixXd V;
-    Eigen::MatrixXi F;
-    double ground_level = 0;
-};*/
-
-/// An index-triangle structure for libIGL functions. Also serves as an
-/// alternative (raw) input format for the SLASupportTree
 class EigenMesh3D {
     class AABBImpl;
 
     Eigen::MatrixXd m_V;
     Eigen::MatrixXi m_F;
-    double m_ground_level = 0;
+    double m_ground_level = 0, m_gnd_offset = 0;
 
     std::unique_ptr<AABBImpl> m_aabb;
 public:
@@ -71,7 +87,9 @@ public:
 
     ~EigenMesh3D();
 
-    inline double ground_level() const { return m_ground_level; }
+    inline double ground_level() const { return m_ground_level + m_gnd_offset; }
+    inline void ground_level_offset(double o) { m_gnd_offset = o; }
+    inline double ground_level_offset() const { return m_gnd_offset; }
 
     inline const Eigen::MatrixXd& V() const { return m_V; }
     inline const Eigen::MatrixXi& F() const { return m_F; }
@@ -149,11 +167,14 @@ public:
 #endif /* SLIC3R_SLA_NEEDS_WINDTREE */
 
     double squared_distance(const Vec3d& p, int& i, Vec3d& c) const;
+    inline double squared_distance(const Vec3d &p) const
+    {
+        int   i;
+        Vec3d c;
+        return squared_distance(p, i, c);
+    }
 };
 
-
-
-
 } // namespace sla
 } // namespace Slic3r
 
diff --git a/src/libslic3r/SLA/SLASpatIndex.hpp b/src/libslic3r/SLA/SLASpatIndex.hpp
index e5fbfa7d4..90dcdc362 100644
--- a/src/libslic3r/SLA/SLASpatIndex.hpp
+++ b/src/libslic3r/SLA/SLASpatIndex.hpp
@@ -7,13 +7,15 @@
 
 #include <Eigen/Geometry>
 
+#include <libslic3r/BoundingBox.hpp>
+
 namespace Slic3r {
 namespace sla {
 
 typedef Eigen::Matrix<double,   3, 1, Eigen::DontAlign> Vec3d;
-using SpatElement = std::pair<Vec3d, unsigned>;
+using PointIndexEl = std::pair<Vec3d, unsigned>;
 
-class SpatIndex {
+class PointIndex {
     class Impl;
 
     // We use Pimpl because it takes a long time to compile boost headers which
@@ -21,30 +23,67 @@ class SpatIndex {
     std::unique_ptr<Impl> m_impl;
 public:
 
-    SpatIndex();
-    ~SpatIndex();
+    PointIndex();
+    ~PointIndex();
 
-    SpatIndex(const SpatIndex&);
-    SpatIndex(SpatIndex&&);
-    SpatIndex& operator=(const SpatIndex&);
-    SpatIndex& operator=(SpatIndex&&);
+    PointIndex(const PointIndex&);
+    PointIndex(PointIndex&&);
+    PointIndex& operator=(const PointIndex&);
+    PointIndex& operator=(PointIndex&&);
 
-    void insert(const SpatElement&);
-    bool remove(const SpatElement&);
+    void insert(const PointIndexEl&);
+    bool remove(const PointIndexEl&);
 
     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);
+    std::vector<PointIndexEl> query(std::function<bool(const PointIndexEl&)>);
+    std::vector<PointIndexEl> nearest(const Vec3d&, unsigned k);
 
     // For testing
     size_t size() const;
     bool empty() const { return size() == 0; }
 
-    void foreach(std::function<void(const SpatElement& el)> fn);
+    void foreach(std::function<void(const PointIndexEl& el)> fn);
+};
+
+using BoxIndexEl = std::pair<Slic3r::BoundingBox, unsigned>;
+
+class BoxIndex {
+    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:
+    
+    BoxIndex();
+    ~BoxIndex();
+    
+    BoxIndex(const BoxIndex&);
+    BoxIndex(BoxIndex&&);
+    BoxIndex& operator=(const BoxIndex&);
+    BoxIndex& operator=(BoxIndex&&);
+    
+    void insert(const BoxIndexEl&);
+    inline void insert(const BoundingBox& bb, unsigned idx)
+    {
+        insert(std::make_pair(bb, unsigned(idx)));
+    }
+    
+    bool remove(const BoxIndexEl&);
+
+    enum QueryType { qtIntersects, qtWithin };
+
+    std::vector<BoxIndexEl> query(const BoundingBox&, QueryType qt);
+    
+    // For testing
+    size_t size() const;
+    bool empty() const { return size() == 0; }
+    
+    void foreach(std::function<void(const BoxIndexEl& el)> fn);
 };
 
 }
diff --git a/src/libslic3r/SLA/SLASupportTree.cpp b/src/libslic3r/SLA/SLASupportTree.cpp
index dbabfbfa7..48b78fcce 100644
--- a/src/libslic3r/SLA/SLASupportTree.cpp
+++ b/src/libslic3r/SLA/SLASupportTree.cpp
@@ -9,10 +9,12 @@
 #include "SLASpatIndex.hpp"
 #include "SLABasePool.hpp"
 
+#include <libslic3r/MTUtils.hpp>
 #include <libslic3r/ClipperUtils.hpp>
 #include <libslic3r/Model.hpp>
 
 #include <libnest2d/optimizers/nlopt/genetic.hpp>
+#include <libnest2d/optimizers/nlopt/subplex.hpp>
 #include <boost/log/trivial.hpp>
 #include <tbb/parallel_for.h>
 #include <libslic3r/I18N.hpp>
@@ -81,6 +83,42 @@ const unsigned SupportConfig::max_bridges_on_pillar = 3;
 using Coordf = double;
 using Portion = std::tuple<double, double>;
 
+// Set this to true to enable full parallelism in this module.
+// Only the well tested parts will be concurrent if this is set to false.
+const constexpr bool USE_FULL_CONCURRENCY = false;
+
+template<bool> struct _ccr {};
+
+template<> struct _ccr<true>
+{
+    using Mutex = SpinMutex;
+
+    template<class It, class Fn>
+    static inline void enumerate(It from, It to, Fn fn)
+    {
+        using TN = size_t;
+        auto iN  = to - from;
+        TN   N   = iN < 0 ? 0 : TN(iN);
+
+        tbb::parallel_for(TN(0), N, [from, fn](TN n) { fn(*(from + n), n); });
+    }
+};
+
+template<> struct _ccr<false>
+{
+    struct Mutex { inline void lock() {} inline void unlock() {} };
+
+    template<class It, class Fn>
+    static inline void enumerate(It from, It to, Fn fn)
+    {
+        for (auto it = from; it != to; ++it) fn(*it, it - from);
+    }
+};
+
+using ccr = _ccr<USE_FULL_CONCURRENCY>;
+using ccr_seq = _ccr<false>;
+using ccr_par = _ccr<true>;
+
 inline Portion make_portion(double a, double b) {
     return std::make_tuple(a, b);
 }
@@ -413,7 +451,7 @@ struct Pillar {
         assert(steps > 0);
 
         height = jp(Z) - endp(Z);
-        if(height > 0) { // Endpoint is below the starting point
+        if(height > EPSILON) { // Endpoint is below the starting point
 
             // We just create a bridge geometry with the pillar parameters and
             // move the data.
@@ -528,6 +566,7 @@ struct CompactBridge {
                   const Vec3d& ep,
                   const Vec3d& n,
                   double r,
+                  bool endball = true,
                   size_t steps = 45)
     {
         Vec3d startp = sp + r * n;
@@ -541,12 +580,14 @@ struct CompactBridge {
         double fa = 2*PI/steps;
         auto upperball = sphere(r, Portion{PI / 2 - fa, PI}, fa);
         for(auto& p : upperball.points) p += startp;
-
-        auto lowerball = sphere(r, Portion{0, PI/2 + 2*fa}, fa);
-        for(auto& p : lowerball.points) p += endp;
-
+        
+        if(endball) {
+            auto lowerball = sphere(r, Portion{0, PI/2 + 2*fa}, fa);
+            for(auto& p : lowerball.points) p += endp;
+            mesh.merge(lowerball);
+        }
+        
         mesh.merge(upperball);
-        mesh.merge(lowerball);
     }
 };
 
@@ -556,29 +597,123 @@ struct Pad {
     PoolConfig cfg;
     double zlevel = 0;
 
-    Pad() {}
+    Pad() = default;
 
-    Pad(const TriangleMesh& object_support_mesh,
-        const ExPolygons& baseplate,
+    Pad(const TriangleMesh& support_mesh,
+        const ExPolygons& modelbase,
         double ground_level,
         const PoolConfig& pcfg) :
         cfg(pcfg),
-        zlevel(ground_level +
-               (sla::get_pad_fullheight(pcfg) - sla::get_pad_elevation(pcfg)) )
+        zlevel(ground_level + 
+               sla::get_pad_fullheight(pcfg) -
+               sla::get_pad_elevation(pcfg))
     {
-        ExPolygons basep;
-        cfg.throw_on_cancel();
+        Polygons basep;
+        auto &thr = cfg.throw_on_cancel;
+        
+        thr();
+        
+        // Get a sample for the pad from the support mesh
+        {
+            ExPolygons platetmp;
 
-        // The 0.1f is the layer height with which the mesh is sampled and then
-        // the layers are unified into one vector of polygons.
-        base_plate(object_support_mesh, basep,
-                   float(cfg.min_wall_height_mm + cfg.min_wall_thickness_mm),
-                   0.1f, pcfg.throw_on_cancel);
+            float zstart = float(zlevel);
+            float zend   = zstart + float(get_pad_fullheight(pcfg) + EPSILON);
 
-        for(auto& bp : baseplate) basep.emplace_back(bp);
+            base_plate(support_mesh, platetmp, grid(zstart, zend, 0.1f), thr);
+
+            // We don't need no... holes control...
+            for (const ExPolygon &bp : platetmp)
+                basep.emplace_back(std::move(bp.contour));
+        }
+        
+        if(pcfg.embed_object) {
+            
+            // If the zero elevation mode is ON, we need to process the model
+            // base silhouette. Create the offsetted version and punch the
+            // breaksticks across its perimeter.
+            
+            ExPolygons modelbase_offs = modelbase;
+            
+            if (pcfg.embed_object.object_gap_mm > 0.0)
+                modelbase_offs
+                    = offset_ex(modelbase_offs,
+                                float(scaled(pcfg.embed_object.object_gap_mm)));
+            
+            // Create a spatial index of the support silhouette polygons.
+            // This will be used to check for intersections with the model
+            // silhouette polygons. If there is no intersection, then a certain
+            // part of the pad is redundant as it does not host any supports.
+            BoxIndex bindex;
+            {
+                unsigned idx = 0;
+                for(auto &bp : basep) {
+                    auto bb = bp.bounding_box();
+                    bb.offset(float(scaled(pcfg.min_wall_thickness_mm)));
+                    bindex.insert(bb, idx++);
+                }
+            }
+            
+            ExPolygons concaveh = offset_ex(
+                concave_hull(basep, pcfg.max_merge_distance_mm, thr),
+                scaled<float>(pcfg.min_wall_thickness_mm));
+            
+            // Punching the breaksticks across the offsetted polygon perimeters
+            auto pad_stickholes = reserve_vector<ExPolygon>(modelbase.size());
+            for(auto& poly : modelbase_offs) {
+                
+                bool overlap = false;
+                for (const ExPolygon &p : concaveh)
+                    overlap = overlap || poly.overlaps(p);
+                
+                auto bb = poly.contour.bounding_box();
+                bb.offset(scaled<float>(pcfg.min_wall_thickness_mm));
+                
+                std::vector<BoxIndexEl> qres =
+                    bindex.query(bb, BoxIndex::qtIntersects);
+                
+                if (!qres.empty() || overlap) {
+                    
+                    // The model silhouette polygon 'poly' HAS an intersection
+                    // with the support silhouettes. Include this polygon
+                    // in the pad holes with the breaksticks and merge the
+                    // original (offsetted) version with the rest of the pad
+                    // base plate.
+                    
+                    basep.emplace_back(poly.contour);
+                    
+                    // The holes of 'poly' will become positive parts of the
+                    // pad, so they has to be checked for intersections as well
+                    // and erased if there is no intersection with the supports
+                    auto it = poly.holes.begin();
+                    while(it != poly.holes.end()) {
+                        if (bindex.query(it->bounding_box(),
+                                         BoxIndex::qtIntersects).empty())
+                            it = poly.holes.erase(it);
+                        else
+                            ++it;
+                    }
+                    
+                    // Punch the breaksticks
+                    sla::breakstick_holes(
+                        poly,
+                        pcfg.embed_object.object_gap_mm,   // padding
+                        pcfg.embed_object.stick_stride_mm,
+                        pcfg.embed_object.stick_width_mm,
+                        pcfg.embed_object.stick_penetration_mm);
+                    
+                    pad_stickholes.emplace_back(poly);
+                }
+            }
+            
+            create_base_pool(basep, tmesh, pad_stickholes, cfg);
+        } else {
+            for (const ExPolygon &bp : modelbase) basep.emplace_back(bp.contour);
+            create_base_pool(basep, tmesh, {}, cfg);
+        }
 
-        create_base_pool(basep, tmesh, cfg);
         tmesh.translate(0, 0, float(zlevel));
+        tmesh.require_shared_vertices();
     }
 
     bool empty() const { return tmesh.facets_count() == 0; }
@@ -603,7 +738,7 @@ inline Vec2d to_vec2(const Vec3d& v3) {
     return {v3(X), v3(Y)};
 }
 
-bool operator==(const SpatElement& e1, const SpatElement& e2) {
+bool operator==(const PointIndexEl& e1, const PointIndexEl& e2) {
     return e1.second == e2.second;
 }
 
@@ -620,7 +755,7 @@ ClusteredPoints cluster(const PointSet& points,
 ClusteredPoints cluster(
         const std::vector<unsigned>& indices,
         std::function<Vec3d(unsigned)> pointfn,
-        std::function<bool(const SpatElement&, const SpatElement&)> predicate,
+        std::function<bool(const PointIndexEl&, const PointIndexEl&)> predicate,
         unsigned max_points);
 
 // This class will hold the support tree meshes with some additional bookkeeping
@@ -637,63 +772,84 @@ ClusteredPoints cluster(
 // The support pad is considered an auxiliary geometry and is not part of the
 // merged mesh. It can be retrieved using a dedicated method (pad())
 class SLASupportTree::Impl {
-    std::map<unsigned, Head> m_heads;
+    // For heads it is beneficial to use the same IDs as for the support points.
+    std::vector<Head> m_heads;
+    std::vector<size_t> m_head_indices;
+    
     std::vector<Pillar> m_pillars;
     std::vector<Junction> m_junctions;
     std::vector<Bridge> m_bridges;
     std::vector<CompactBridge> m_compact_bridges;
     Controller m_ctl;
-
+    
     Pad m_pad;
+    
+    using Mutex = ccr::Mutex;
+    
+    mutable Mutex m_mutex;
     mutable TriangleMesh meshcache; mutable bool meshcache_valid = false;
     mutable double model_height = 0; // the full height of the model
+    
 public:
     double ground_level = 0;
-
+    
     Impl() = default;
     inline Impl(const Controller& ctl): m_ctl(ctl) {}
-
+    
     const Controller& ctl() const { return m_ctl; }
-
-    template<class...Args> Head& add_head(unsigned id, Args&&... args) {
-        auto el = m_heads.emplace(std::piecewise_construct,
-                            std::forward_as_tuple(id),
-                            std::forward_as_tuple(std::forward<Args>(args)...));
-        el.first->second.id = id;
+    
+    template<class...Args> Head& add_head(unsigned id, Args&&... args)
+    {
+        std::lock_guard<Mutex> lk(m_mutex);
+        m_heads.emplace_back(std::forward<Args>(args)...);
+        m_heads.back().id = id;
+        
+        if (id >= m_head_indices.size()) m_head_indices.resize(id + 1);
+        m_head_indices[id] = m_heads.size() - 1;
+        
         meshcache_valid = false;
-        return el.first->second;
+        return m_heads.back();
     }
-
-    template<class...Args> Pillar& add_pillar(unsigned headid, Args&&... args) {
-        auto it = m_heads.find(headid);
-        assert(it != m_heads.end());
-        Head& head = it->second;
+    
+    template<class...Args> Pillar& add_pillar(unsigned headid, Args&&... args)
+    {
+        std::lock_guard<Mutex> lk(m_mutex);
+        
+        assert(headid < m_head_indices.size());
+        Head &head = m_heads[m_head_indices[headid]];
+        
         m_pillars.emplace_back(head, std::forward<Args>(args)...);
         Pillar& pillar = m_pillars.back();
         pillar.id = long(m_pillars.size() - 1);
         head.pillar_id = pillar.id;
         pillar.start_junction_id = head.id;
         pillar.starts_from_head = true;
+        
         meshcache_valid = false;
         return m_pillars.back();
     }
-
-    void increment_bridges(const Pillar& pillar) {
+    
+    void increment_bridges(const Pillar& pillar)
+    {
+        std::lock_guard<Mutex> lk(m_mutex);
         assert(pillar.id >= 0 && size_t(pillar.id) < m_pillars.size());
-
-        if(pillar.id >= 0 && size_t(pillar.id) < m_pillars.size())
+        
+        if(pillar.id >= 0 && size_t(pillar.id) < m_pillars.size()) 
             m_pillars[size_t(pillar.id)].bridges++;
     }
-
-    void increment_links(const Pillar& pillar) {
+    
+    void increment_links(const Pillar& pillar)
+    {
+        std::lock_guard<Mutex> lk(m_mutex);
         assert(pillar.id >= 0 && size_t(pillar.id) < m_pillars.size());
-
+        
         if(pillar.id >= 0 && size_t(pillar.id) < m_pillars.size())
-            m_pillars[size_t(pillar.id)].links++;
+            m_pillars[size_t(pillar.id)].links++;        
     }
-
+    
     template<class...Args> Pillar& add_pillar(Args&&...args)
     {
+        std::lock_guard<Mutex> lk(m_mutex);
         m_pillars.emplace_back(std::forward<Args>(args)...);
         Pillar& pillar = m_pillars.back();
         pillar.id = long(m_pillars.size() - 1);
@@ -701,162 +857,175 @@ public:
         meshcache_valid = false;
         return m_pillars.back();
     }
-
-    const Head& pillar_head(long pillar_id) const {
+    
+    const Head& pillar_head(long pillar_id) const
+    {
+        std::lock_guard<Mutex> lk(m_mutex);
         assert(pillar_id >= 0 && pillar_id < long(m_pillars.size()));
+        
         const Pillar& p = m_pillars[size_t(pillar_id)];
         assert(p.starts_from_head && p.start_junction_id >= 0);
-        auto it = m_heads.find(unsigned(p.start_junction_id));
-        assert(it != m_heads.end());
-        return it->second;
+        assert(size_t(p.start_junction_id) < m_head_indices.size());
+        
+        return m_heads[m_head_indices[p.start_junction_id]];
     }
-
-    const Pillar& head_pillar(unsigned headid) const {
-        auto it = m_heads.find(headid);
-        assert(it != m_heads.end());
-        const Head& h = it->second;
+    
+    const Pillar& head_pillar(unsigned headid) const
+    {
+        std::lock_guard<Mutex> lk(m_mutex);
+        assert(headid < m_head_indices.size());
+        
+        const Head& h = m_heads[m_head_indices[headid]];
         assert(h.pillar_id >= 0 && h.pillar_id < long(m_pillars.size()));
-        return pillar(h.pillar_id);
+        
+        return m_pillars[size_t(h.pillar_id)];
     }
-
-    template<class...Args> const Junction& add_junction(Args&&... args) {
+    
+    template<class...Args> const Junction& add_junction(Args&&... args)
+    {
+        std::lock_guard<Mutex> lk(m_mutex);
         m_junctions.emplace_back(std::forward<Args>(args)...);
         m_junctions.back().id = long(m_junctions.size() - 1);
         meshcache_valid = false;
         return m_junctions.back();
     }
-
-    template<class...Args> const Bridge& add_bridge(Args&&... args) {
+    
+    template<class...Args> const Bridge& add_bridge(Args&&... args)
+    {
+        std::lock_guard<Mutex> lk(m_mutex);
         m_bridges.emplace_back(std::forward<Args>(args)...);
         m_bridges.back().id = long(m_bridges.size() - 1);
         meshcache_valid = false;
         return m_bridges.back();
     }
-
-    template<class...Args>
-    const CompactBridge& add_compact_bridge(Args&&...args) {
+    
+    template<class...Args> const CompactBridge& add_compact_bridge(Args&&...args)
+    {
+        std::lock_guard<Mutex> lk(m_mutex);
         m_compact_bridges.emplace_back(std::forward<Args>(args)...);
         m_compact_bridges.back().id = long(m_compact_bridges.size() - 1);
         meshcache_valid = false;
         return m_compact_bridges.back();
     }
-
-    const std::map<unsigned, Head>& heads() const { return m_heads; }
-    Head& head(unsigned idx) {
+    
+    Head &head(unsigned id)
+    {
+        std::lock_guard<Mutex> lk(m_mutex);
+        assert(id < m_head_indices.size());
+        
         meshcache_valid = false;
-        auto it = m_heads.find(idx);
-        assert(it != m_heads.end());
-        return it->second;
+        return m_heads[m_head_indices[id]];
     }
-    const std::vector<Pillar>& pillars() const { return m_pillars; }
-    const std::vector<Bridge>& bridges() const { return m_bridges; }
-    const std::vector<Junction>& junctions() const { return m_junctions; }
-    const std::vector<CompactBridge>& compact_bridges() const {
-        return m_compact_bridges;
+    
+    inline size_t pillarcount() const {
+        std::lock_guard<Mutex> lk(m_mutex);
+        return m_pillars.size();
     }
-
-    template<class T> inline const Pillar& pillar(T id) const {
-        static_assert(std::is_integral<T>::value, "Invalid index type");
+    
+    template<class T> inline IntegerOnly<T, const Pillar&> pillar(T id) const
+    {
+        std::lock_guard<Mutex> lk(m_mutex);
         assert(id >= 0 && size_t(id) < m_pillars.size() &&
                size_t(id) < std::numeric_limits<size_t>::max());
+        
         return m_pillars[size_t(id)];
     }
-
-    const Pad& create_pad(const TriangleMesh& object_supports,
-                          const ExPolygons& baseplate,
-                          const PoolConfig& cfg) {
-        m_pad = Pad(object_supports, baseplate, ground_level, cfg);
+    
+    const Pad &create_pad(const TriangleMesh &object_supports,
+                          const ExPolygons &  modelbase,
+                          const PoolConfig &  cfg)
+    {
+        m_pad = Pad(object_supports, modelbase, ground_level, cfg);
         return m_pad;
     }
-
-    void remove_pad() {
-        m_pad = Pad();
-    }
-
+    
+    void remove_pad() { m_pad = Pad(); }
+    
     const Pad& pad() const { return m_pad; }
-
+    
     // WITHOUT THE PAD!!!
-    const TriangleMesh& merged_mesh() const {
-        if(meshcache_valid) return meshcache;
-
+    const TriangleMesh &merged_mesh() const
+    {
+        if (meshcache_valid) return meshcache;
+        
         Contour3D merged;
-
-        for(auto& headel : heads()) {
-            if(m_ctl.stopcondition()) break;
-            if(headel.second.is_valid())
-                merged.merge(headel.second.mesh);
+        
+        for (auto &head : m_heads) {
+            if (m_ctl.stopcondition()) break;
+            if (head.is_valid()) merged.merge(head.mesh);
         }
-
-        for(auto& stick : pillars()) {
-            if(m_ctl.stopcondition()) break;
+        
+        for (auto &stick : m_pillars) {
+            if (m_ctl.stopcondition()) break;
             merged.merge(stick.mesh);
             merged.merge(stick.base);
         }
-
-        for(auto& j : junctions()) {
-            if(m_ctl.stopcondition()) break;
+        
+        for (auto &j : m_junctions) {
+            if (m_ctl.stopcondition()) break;
             merged.merge(j.mesh);
         }
-
-        for(auto& cb : compact_bridges()) {
-            if(m_ctl.stopcondition()) break;
+        
+        for (auto &cb : m_compact_bridges) {
+            if (m_ctl.stopcondition()) break;
             merged.merge(cb.mesh);
         }
-
-        for(auto& bs : bridges()) {
-            if(m_ctl.stopcondition()) break;
+        
+        for (auto &bs : m_bridges) {
+            if (m_ctl.stopcondition()) break;
             merged.merge(bs.mesh);
         }
-
-
-        if(m_ctl.stopcondition()) {
+        
+        if (m_ctl.stopcondition()) {
             // In case of failure we have to return an empty mesh
             meshcache = TriangleMesh();
             return meshcache;
         }
-
+        
         meshcache = mesh(merged);
-
+        
         // The mesh will be passed by const-pointer to TriangleMeshSlicer,
         // which will need this.
-        meshcache.require_shared_vertices();
-
-        // TODO: Is this necessary?
-        //meshcache.repair();
-
-        BoundingBoxf3&& bb = meshcache.bounding_box();
-        model_height = bb.max(Z) - bb.min(Z);
-
+        if (!meshcache.empty()) meshcache.require_shared_vertices();
+        
+        BoundingBoxf3 &&bb = meshcache.bounding_box();
+        model_height       = bb.max(Z) - bb.min(Z);
+        
         meshcache_valid = true;
         return meshcache;
     }
-
+    
     // WITH THE PAD
-    double full_height() const {
-        if(merged_mesh().empty() && !pad().empty())
+    double full_height() const
+    {
+        if (merged_mesh().empty() && !pad().empty())
             return get_pad_fullheight(pad().cfg);
-
+        
         double h = mesh_height();
-        if(!pad().empty()) h += sla::get_pad_elevation(pad().cfg);
+        if (!pad().empty()) h += sla::get_pad_elevation(pad().cfg);
         return h;
     }
-
+    
     // WITHOUT THE PAD!!!
-    double mesh_height() const {
-        if(!meshcache_valid) merged_mesh();
+    double mesh_height() const
+    {
+        if (!meshcache_valid) merged_mesh();
         return model_height;
     }
     
     // Intended to be called after the generation is fully complete
-    void clear_support_data() {
+    void merge_and_cleanup()
+    {
         merged_mesh(); // in case the mesh is not generated, it should be...
-        m_heads.clear();
-        m_pillars.clear();
-        m_junctions.clear();
-        m_bridges.clear();
-        m_compact_bridges.clear();
+        
+        // Doing clear() does not garantee to release the memory.
+        m_heads = {};
+        m_head_indices = {};
+        m_pillars = {};
+        m_junctions = {};
+        m_bridges = {};
+        m_compact_bridges = {};
     }
-
 };
 
 // This function returns the position of the centroid in the input 'clust'
@@ -947,7 +1116,7 @@ class SLASupportTree::Algorithm {
     ThrowOnCancel m_thr;
 
     // A spatial index to easily find strong pillars to connect to.
-    SpatIndex m_pillar_index;
+    PointIndex m_pillar_index;
 
     inline double ray_mesh_intersect(const Vec3d& s,
                                      const Vec3d& dir)
@@ -1025,11 +1194,10 @@ class SLASupportTree::Algorithm {
         // Now a and b vectors are perpendicular to v and to each other.
         // Together they define the plane where we have to iterate with the
         // given angles in the 'phis' vector
-        tbb::parallel_for(size_t(0), phis.size(),
-                          [&phis, &hits, &m, sd, r_pin, r_back, s, a, b, c]
-                          (size_t i)
+        ccr_par::enumerate(phis.begin(), phis.end(),
+                           [&hits, &m, sd, r_pin, r_back, s, a, b, c]
+                           (double phi, size_t i)
         {
-            double& phi = phis[i];
             double sinphi = std::sin(phi);
             double cosphi = std::cos(phi);
 
@@ -1128,12 +1296,11 @@ class SLASupportTree::Algorithm {
 
         // Hit results
         std::array<HitResult, SAMPLES> hits;
-
-        tbb::parallel_for(size_t(0), phis.size(),
-                          [&m, &phis, a, b, sd, dir, r, s, ins_check, &hits]
-                          (size_t i)
+        
+        ccr_par::enumerate(phis.begin(), phis.end(),
+                           [&m, a, b, sd, dir, r, s, ins_check, &hits]
+                           (double phi, size_t i)
         {
-            double& phi = phis[i];
             double sinphi = std::sin(phi);
             double cosphi = std::cos(phi);
 
@@ -1149,7 +1316,7 @@ class SLASupportTree::Algorithm {
             auto hr = m.query_ray_hit(p + sd*dir, dir);
 
             if(ins_check && hr.is_inside()) {
-                if(hr.distance() > r + sd) hits[i] = HitResult(0.0);
+                if(hr.distance() > 2 * r + sd) hits[i] = HitResult(0.0);
                 else {
                     // re-cast the ray from the outside of the object
                     auto hr2 =
@@ -1264,9 +1431,12 @@ class SLASupportTree::Algorithm {
 
     // For connecting a head to a nearby pillar.
     bool connect_to_nearpillar(const Head& head, long nearpillar_id) {
-
-        auto nearpillar = [this, nearpillar_id]() { return m_result.pillar(nearpillar_id); };
-        if(nearpillar().bridges > m_cfg.max_bridges_on_pillar) return false;
+        
+        auto nearpillar = [this, nearpillar_id]() {
+            return m_result.pillar(nearpillar_id);
+        };
+        
+        if (nearpillar().bridges > m_cfg.max_bridges_on_pillar) return false;
 
         Vec3d headjp = head.junction_point();
         Vec3d nearjp_u = nearpillar().startpoint();
@@ -1337,7 +1507,7 @@ class SLASupportTree::Algorithm {
     }
 
     bool search_pillar_and_connect(const Head& head) {
-        SpatIndex spindex = m_pillar_index;
+        PointIndex spindex = m_pillar_index;
 
         long nearest_id = -1;
 
@@ -1358,7 +1528,7 @@ class SLASupportTree::Algorithm {
 
             if(nearest_id >= 0) {
                 auto nearpillarID = unsigned(nearest_id);
-                if(nearpillarID < m_result.pillars().size()) {
+                if(nearpillarID < m_result.pillarcount()) {
                     if(!connect_to_nearpillar(head, nearpillarID)) {
                         nearest_id = -1;    // continue searching
                         spindex.remove(ne); // without the current pillar
@@ -1369,6 +1539,120 @@ class SLASupportTree::Algorithm {
 
         return nearest_id >= 0;
     }
+    
+    // This is a proxy function for pillar creation which will mind the gap
+    // between the pad and the model bottom in zero elevation mode.
+    void create_ground_pillar(const Vec3d &jp,
+                              const Vec3d &sourcedir,
+                              double       radius,
+                              int          head_id = -1)
+    {
+        // People were killed for this number (seriously)
+        static const double SQR2 = std::sqrt(2.0);
+        static const Vec3d  DOWN = {0.0, 0.0, -1.0};
+
+        double gndlvl       = m_result.ground_level;
+        Vec3d  endp         = {jp(X), jp(Y), gndlvl};
+        double sd           = m_cfg.pillar_base_safety_distance_mm;
+        int    pillar_id    = -1;
+        double min_dist     = sd + m_cfg.base_radius_mm + EPSILON;
+        double dist         = 0;
+        bool   can_add_base = true;
+        bool   normal_mode  = true;
+
+        if (m_cfg.object_elevation_mm < EPSILON
+            && (dist = std::sqrt(m_mesh.squared_distance(endp))) < min_dist) {
+            // Get the distance from the mesh. This can be later optimized
+            // to get the distance in 2D plane because we are dealing with
+            // the ground level only.
+
+            normal_mode     = false;
+            double mv       = min_dist - dist;
+            double azimuth  = std::atan2(sourcedir(Y), sourcedir(X));
+            double sinpolar = std::sin(PI - m_cfg.bridge_slope);
+            double cospolar = std::cos(PI - m_cfg.bridge_slope);
+            double cosazm   = std::cos(azimuth);
+            double sinazm   = std::sin(azimuth);
+
+            auto dir = Vec3d(cosazm * sinpolar, sinazm * sinpolar, cospolar)
+                           .normalized();
+
+            using namespace libnest2d::opt;
+            StopCriteria scr;
+            scr.stop_score = min_dist;
+            SubplexOptimizer solver(scr);
+
+            auto result = solver.optimize_max(
+                [this, dir, jp, gndlvl](double mv) {
+                    Vec3d endp = jp + SQR2 * mv * dir;
+                    endp(Z)    = gndlvl;
+                    return std::sqrt(m_mesh.squared_distance(endp));
+                },
+                initvals(mv), bound(0.0, 2 * min_dist));
+
+            mv           = std::get<0>(result.optimum);
+            endp         = jp + SQR2 * mv * dir;
+            Vec3d pgnd   = {endp(X), endp(Y), gndlvl};
+            can_add_base = result.score > min_dist;
+            
+            double gnd_offs = m_mesh.ground_level_offset();
+            auto abort_in_shame =
+                [gnd_offs, &normal_mode, &can_add_base, &endp, jp, gndlvl]()
+            {
+                normal_mode  = true;
+                can_add_base = false;   // Nothing left to do, hope for the best
+                endp         = {jp(X), jp(Y), gndlvl - gnd_offs };
+            };
+
+            // We have to check if the bridge is feasible.
+            if (bridge_mesh_intersect(jp, dir, radius) < (endp - jp).norm())
+                abort_in_shame();
+            else {
+                // If the new endpoint is below ground, do not make a pillar
+                if (endp(Z) < gndlvl)
+                    endp = endp - SQR2 * (gndlvl - endp(Z)) * dir; // back off
+                else {
+                    
+                    auto hit = bridge_mesh_intersect(endp, DOWN, radius);
+                    if (!std::isinf(hit.distance())) abort_in_shame();
+
+                    Pillar &plr = m_result.add_pillar(endp, pgnd, radius);
+
+                    if (can_add_base)
+                        plr.add_base(m_cfg.base_height_mm,
+                                     m_cfg.base_radius_mm);
+
+                    pillar_id = plr.id;
+                }
+
+                m_result.add_bridge(jp, endp, radius);
+                m_result.add_junction(endp, radius);
+
+                // Add a degenerated pillar and the bridge.
+                // The degenerate pillar will have zero length and it will
+                // prevent from queries of head_pillar() to have non-existing
+                // pillar when the head should have one.
+                if (head_id >= 0)
+                    m_result.add_pillar(unsigned(head_id), jp, radius);
+            }
+        }
+        
+        if (normal_mode) {
+            Pillar &plr = head_id >= 0
+                              ? m_result.add_pillar(unsigned(head_id),
+                                                    endp,
+                                                    radius)
+                              : m_result.add_pillar(jp, endp, radius);
+
+            if (can_add_base)
+                plr.add_base(m_cfg.base_height_mm, m_cfg.base_radius_mm);
+
+            pillar_id = plr.id;
+        } 
+            
+        if(pillar_id >= 0) // Save the pillar endpoint in the spatial index
+            m_pillar_index.insert(endp, pillar_id);
+    }
 
 public:
 
@@ -1432,84 +1716,90 @@ public:
         using libnest2d::opt::initvals;
         using libnest2d::opt::GeneticOptimizer;
         using libnest2d::opt::StopCriteria;
-
-        for(unsigned i = 0, fidx = 0; i < filtered_indices.size(); ++i)
+        
+        ccr::Mutex mutex;
+        auto addfn = [&mutex](PtIndices &container, unsigned val) {
+            std::lock_guard<ccr::Mutex> lk(mutex);
+            container.emplace_back(val);
+        };
+        
+        ccr::enumerate(filtered_indices.begin(), filtered_indices.end(),
+                       [this, &nmls, addfn](unsigned fidx, size_t i)
         {
             m_thr();
-
-            fidx = filtered_indices[i];
+            
             auto n = nmls.row(i);
-
+            
             // for all normals we generate the spherical coordinates and
             // saturate the polar angle to 45 degrees from the bottom then
             // convert back to standard coordinates to get the new normal.
             // Then we just create a quaternion from the two normals
             // (Quaternion::FromTwoVectors) and apply the rotation to the
             // arrow head.
-
-            double z = n(2);
-            double r = 1.0;     // for normalized vector
-            double polar = std::acos(z / r);
+            
+            double z       = n(2);
+            double r       = 1.0; // for normalized vector
+            double polar   = std::acos(z / r);
             double azimuth = std::atan2(n(1), n(0));
-
+            
             // skip if the tilt is not sane
             if(polar >= PI - m_cfg.normal_cutoff_angle) {
-
+                
                 // We saturate the polar angle to 3pi/4
                 polar = std::max(polar, 3*PI / 4);
-
+                
                 // save the head (pinpoint) position
                 Vec3d hp = m_points.row(fidx);
-
+                
                 double w = m_cfg.head_width_mm +
                            m_cfg.head_back_radius_mm +
                            2*m_cfg.head_front_radius_mm;
-
+                
                 double pin_r = double(m_support_pts[fidx].head_front_radius);
-
+                
                 // Reassemble the now corrected normal
                 auto nn = Vec3d(std::cos(azimuth) * std::sin(polar),
                                 std::sin(azimuth) * std::sin(polar),
                                 std::cos(polar)).normalized();
-
+                
                 // check available distance
-                double t = pinhead_mesh_intersect(
-                                  hp, // touching point
-                                  nn, // normal
-                                  pin_r,
-                                  m_cfg.head_back_radius_mm,
-                                  w);
-
-                if(t <= w) {
-
+                EigenMesh3D::hit_result t
+                    = pinhead_mesh_intersect(hp, // touching point
+                                             nn, // normal
+                                             pin_r,
+                                             m_cfg.head_back_radius_mm,
+                                             w);
+                
+                if(t.distance() <= w) {
+                    
                     // Let's try to optimize this angle, there might be a
                     // viable normal that doesn't collide with the model
                     // geometry and its very close to the default.
-
+                    
                     StopCriteria stc;
                     stc.max_iterations = m_cfg.optimizer_max_iterations;
                     stc.relative_score_difference = m_cfg.optimizer_rel_score_diff;
                     stc.stop_score = w; // space greater than w is enough
                     GeneticOptimizer solver(stc);
                     solver.seed(0); // we want deterministic behavior
-
+                    
                     auto oresult = solver.optimize_max(
                         [this, pin_r, w, hp](double plr, double azm)
-                    {
-                        auto n = Vec3d(std::cos(azm) * std::sin(plr),
-                                       std::sin(azm) * std::sin(plr),
-                                       std::cos(plr)).normalized();
-
-                        double score = pinhead_mesh_intersect( hp, n, pin_r,
-                                          m_cfg.head_back_radius_mm, w);
-
-                        return score;
-                    },
-                    initvals(polar, azimuth), // start with what we have
-                    bound(3*PI/4, PI),  // Must not exceed the tilt limit
-                    bound(-PI, PI)      // azimuth can be a full search
-                    );
-
+                        {
+                            auto n = Vec3d(std::cos(azm) * std::sin(plr),
+                                           std::sin(azm) * std::sin(plr),
+                                           std::cos(plr)).normalized();
+
+                            double score = pinhead_mesh_intersect(
+                                hp, n, pin_r, m_cfg.head_back_radius_mm, w);
+
+                            return score;
+                        },
+                        initvals(polar, azimuth), // start with what we have
+                        bound(3*PI/4, PI),  // Must not exceed the tilt limit
+                        bound(-PI, PI)      // azimuth can be a full search
+                        );
+                    
                     if(oresult.score > w) {
                         polar = std::get<0>(oresult.optimum);
                         azimuth = std::get<1>(oresult.optimum);
@@ -1519,20 +1809,25 @@ public:
                         t = oresult.score;
                     }
                 }
-
+                
                 // save the verified and corrected normal
                 m_support_nmls.row(fidx) = nn;
-
-                if(t > w) {
-                    // mark the point for needing a head.
-                    m_iheads.emplace_back(fidx);
-                } else if( polar >= 3*PI/4 ) {
-                    // Headless supports do not tilt like the headed ones so
-                    // the normal should point almost to the ground.
-                    m_iheadless.emplace_back(fidx);
+                
+                if (t.distance() > w) {
+                    // Check distance from ground, we might have zero elevation.
+                    if (hp(Z) + w * nn(Z) < m_result.ground_level) {
+                        addfn(m_iheadless, fidx);
+                    } else {
+                        // mark the point for needing a head.
+                        addfn(m_iheads, fidx);
+                    }
+                } else if (polar >= 3 * PI / 4) {
+                    // Headless supports do not tilt like the headed ones
+                    // so the normal should point almost to the ground.
+                    addfn(m_iheadless, fidx);
                 }
             }
-        }
+        });
 
         m_thr();
     }
@@ -1594,16 +1889,22 @@ public:
         // from each other in the XY plane to not cross their pillar bases
         // These clusters of support points will join in one pillar,
         // possibly in their centroid support point.
+        
         auto pointfn = [this](unsigned i) {
             return m_result.head(i).junction_point();
         };
-        auto predicate = [this](const SpatElement& e1, const SpatElement& e2) {
+
+        auto predicate = [this](const PointIndexEl &e1,
+                                const PointIndexEl &e2) {
             double d2d = distance(to_2d(e1.first), to_2d(e2.first));
             double d3d = distance(e1.first, e2.first);
-            return d2d < 2 * m_cfg.base_radius_mm &&
-                   d3d < m_cfg.max_bridge_length_mm;
+            return d2d < 2 * m_cfg.base_radius_mm
+                   && d3d < m_cfg.max_bridge_length_mm;
         };
-        m_pillar_clusters = cluster(ground_head_indices, pointfn, predicate,
+
+        m_pillar_clusters = cluster(ground_head_indices,
+                                    pointfn,
+                                    predicate,
                                     m_cfg.max_bridges_on_pillar);
     }
 
@@ -1615,7 +1916,7 @@ public:
     void routing_to_ground()
     {
         const double pradius = m_cfg.head_back_radius_mm;
-        const double gndlvl = m_result.ground_level;
+        // const double gndlvl = m_result.ground_level;
 
         ClusterEl cl_centroids;
         cl_centroids.reserve(m_pillar_clusters.size());
@@ -1648,13 +1949,8 @@ public:
 
             Head& h = m_result.head(hid);
             h.transform();
-            Vec3d p = h.junction_point(); p(Z) = gndlvl;
-            auto& plr = m_result.add_pillar(hid, p, h.r_back_mm)
-                                .add_base(m_cfg.base_height_mm,
-                                          m_cfg.base_radius_mm);
 
-            // Save the pillar endpoint and the pillar id in the spatial index
-            m_pillar_index.insert(plr.endpoint(), unsigned(plr.id));
+            create_ground_pillar(h.junction_point(), h.dir, h.r_back_mm, h.id);
         }
 
         // now we will go through the clusters ones again and connect the
@@ -1681,15 +1977,12 @@ public:
                    !search_pillar_and_connect(sidehead))
                 {
                     Vec3d pstart = sidehead.junction_point();
-                    Vec3d pend = Vec3d{pstart(X), pstart(Y), gndlvl};
+                    //Vec3d pend = Vec3d{pstart(X), pstart(Y), gndlvl};
                     // Could not find a pillar, create one
-                    auto& pillar = m_result.add_pillar(unsigned(sidehead.id),
-                                                       pend, pradius)
-                                           .add_base(m_cfg.base_height_mm,
-                                                     m_cfg.base_radius_mm);
-
-                    // connects to ground, eligible for bridging
-                    m_pillar_index.insert(pend, unsigned(pillar.id));
+                    create_ground_pillar(pstart,
+                                         sidehead.dir,
+                                         pradius,
+                                         sidehead.id);
                 }
             }
         }
@@ -1718,20 +2011,21 @@ public:
             m_result.add_bridge(hjp, endp, head.r_back_mm);
             m_result.add_junction(endp, head.r_back_mm);
 
-            auto groundp = endp;
-            groundp(Z) = m_result.ground_level;
-            auto& newpillar = m_result.add_pillar(endp, groundp, head.r_back_mm)
-                                      .add_base(m_cfg.base_height_mm,
-                                                m_cfg.base_radius_mm);
-            m_pillar_index.insert(groundp, unsigned(newpillar.id));
+            this->create_ground_pillar(endp, dir, head.r_back_mm);
         };
 
         std::vector<unsigned> modelpillars;
+        ccr::Mutex mutex;
 
         // TODO: connect these to the ground pillars if possible
-        for(auto item : m_iheads_onmodel) { m_thr();
-            unsigned idx = item.first;
-            EigenMesh3D::hit_result hit = item.second;
+        ccr::enumerate(m_iheads_onmodel.begin(), m_iheads_onmodel.end(),
+                       [this, routedown, &modelpillars, &mutex]
+                       (const std::pair<unsigned, EigenMesh3D::hit_result> &el,
+                        size_t)
+        {
+            m_thr();
+            unsigned idx = el.first;
+            EigenMesh3D::hit_result hit = el.second;
 
             auto& head = m_result.head(idx);
             Vec3d hjp = head.junction_point();
@@ -1740,7 +2034,7 @@ public:
             // Search nearby pillar
             // /////////////////////////////////////////////////////////////////
 
-            if(search_pillar_and_connect(head)) { head.transform(); continue; }
+            if(search_pillar_and_connect(head)) { head.transform(); return; }
 
             // /////////////////////////////////////////////////////////////////
             // Try straight path
@@ -1762,7 +2056,7 @@ public:
             }
 
             if(std::isinf(tdown)) { // we heave found a route to the ground
-                routedown(head, head.dir, d); continue;
+                routedown(head, head.dir, d); return;
             }
 
             // /////////////////////////////////////////////////////////////////
@@ -1824,7 +2118,7 @@ public:
             }
 
             if(std::isinf(tdown)) { // we heave found a route to the ground
-                routedown(head, bridgedir, d); continue;
+                routedown(head, bridgedir, d); return;
             }
 
             // /////////////////////////////////////////////////////////////////
@@ -1867,8 +2161,9 @@ public:
                 pill.base = tailhead.mesh;
 
                 // Experimental: add the pillar to the index for cascading
+                std::lock_guard<ccr::Mutex> lk(mutex);
                 modelpillars.emplace_back(unsigned(pill.id));
-                continue;
+                return;
             }
 
             // We have failed to route this head.
@@ -1876,13 +2171,35 @@ public:
                     << "Failed to route model facing support point."
                     << " ID: " << idx;
             head.invalidate();
-        }
+        });
 
         for(auto pillid : modelpillars) {
             auto& pillar = m_result.pillar(pillid);
             m_pillar_index.insert(pillar.endpoint(), pillid);
         }
     }
+    
+    // Helper function for interconnect_pillars where pairs of already connected
+    // pillars should be checked for not to be processed again. This can be done
+    // in O(log) or even constant time with a set or an unordered set of hash
+    // values uniquely representing a pair of integers. The order of numbers
+    // within the pair should not matter, it has the same unique hash.
+    template<class I> static I pairhash(I a, I b)
+    {
+        using std::ceil; using std::log2; using std::max; using std::min;
+        
+        static_assert(std::is_integral<I>::value,
+                      "This function works only for integral types.");
+
+        I g = min(a, b), l = max(a, b);
+        
+        auto bits_g = g ? int(ceil(log2(g))) : 0;
+
+        // Assume the hash will fit into the output variable
+        assert((l ? (ceil(log2(l))) : 0) + bits_g < int(sizeof(I) * CHAR_BIT));
+        
+        return (l << bits_g) + g;
+    }
 
     void interconnect_pillars() {
         // Now comes the algorithm that connects pillars with each other.
@@ -1900,45 +2217,51 @@ public:
         double min_height_ratio = 0.5;
 
         std::set<unsigned long> pairs;
-
+        
+        // A function to connect one pillar with its neighbors. THe number of
+        // neighbors is given in the configuration. This function if called
+        // for every pillar in the pillar index. A pair of pillar will not
+        // be connected multiple times this is ensured by the 'pairs' set which
+        // remembers the processed pillar pairs
         auto cascadefn =
-                [this, d, &pairs, min_height_ratio, H1] (const SpatElement& el)
+                [this, d, &pairs, min_height_ratio, H1] (const PointIndexEl& el)
         {
-            Vec3d qp = el.first;
-
-            const Pillar& pillar = m_result.pillar(el.second);
+            Vec3d qp = el.first;    // endpoint of the pillar
 
+            const Pillar& pillar = m_result.pillar(el.second); // actual pillar
+            
+            // Get the max number of neighbors a pillar should connect to
             unsigned neighbors = m_cfg.pillar_cascade_neighbors;
 
-            // connections are enough for one pillar
+            // connections are already enough for the pillar
             if(pillar.links >= neighbors) return;
 
             // Query all remaining points within reach
-            auto qres = m_pillar_index.query([qp, d](const SpatElement& e){
+            auto qres = m_pillar_index.query([qp, d](const PointIndexEl& e){
                 return distance(e.first, qp) < d;
             });
 
             // sort the result by distance (have to check if this is needed)
             std::sort(qres.begin(), qres.end(),
-                      [qp](const SpatElement& e1, const SpatElement& e2){
+                      [qp](const PointIndexEl& e1, const PointIndexEl& e2){
                 return distance(e1.first, qp) < distance(e2.first, qp);
             });
 
-            for(auto& re : qres) {
+            for(auto& re : qres) { // process the queried neighbors
 
-                if(re.second == el.second) continue;
+                if(re.second == el.second) continue; // Skip self
 
                 auto a = el.second, b = re.second;
 
-                // I hope that the area of a square is never equal to its
-                // circumference
-                auto hashval = 2 * (a + b) + a * b;
-
+                // Get unique hash for the given pair (order doesn't matter)
+                auto hashval = pairhash(a, b);
+                
+                // Search for the pair amongst the remembered pairs
                 if(pairs.find(hashval) != pairs.end()) continue;
+                
+                const Pillar& neighborpillar = m_result.pillar(re.second);
 
-                const Pillar& neighborpillar = m_result.pillars()[re.second];
-
-                // this neighbor is occupied
+                // this neighbor is occupied, skip
                 if(neighborpillar.links >= neighbors) continue;
 
                 if(interconnect(pillar, neighborpillar)) {
@@ -1960,47 +2283,79 @@ public:
                 if(pillar.links >= neighbors) break;
             }
         };
-
+        
+        // Run the cascade for the pillars in the index
         m_pillar_index.foreach(cascadefn);
-
-        size_t pillarcount = m_result.pillars().size();
-
+       
+        // We would be done here if we could allow some pillars to not be
+        // connected with any neighbors. But this might leave the support tree
+        // unprintable.
+        //
+        // The current solution is to insert additional pillars next to these
+        // lonely pillars. One or even two additional pillar might get inserted
+        // depending on the length of the lonely pillar.
+        
+        size_t pillarcount = m_result.pillarcount();
+        
+        // Again, go through all pillars, this time in the whole support tree
+        // not just the index.
         for(size_t pid = 0; pid < pillarcount; pid++) {
             auto pillar = [this, pid]() { return m_result.pillar(pid); };
-
+           
+            // Decide how many additional pillars will be needed:
+            
             unsigned needpillars = 0;
-            if(pillar().bridges > m_cfg.max_bridges_on_pillar) needpillars = 3;
-            else if(pillar().links < 2 && pillar().height > H2) {
+            if (pillar().bridges > m_cfg.max_bridges_on_pillar)
+                needpillars = 3;
+            else if (pillar().links < 2 && pillar().height > H2) {
                 // Not enough neighbors to support this pillar
                 needpillars = 2 - pillar().links;
-            }
-            else if(pillar().links < 1 && pillar().height > H1) {
+            } else if (pillar().links < 1 && pillar().height > H1) {
                 // No neighbors could be found and the pillar is too long.
                 needpillars = 1;
             }
 
-            // Search for new pillar locations
-            bool found = false;
-            double alpha = 0; // goes to 2Pi
-            double r = 2 * m_cfg.base_radius_mm;
-            Vec3d pillarsp = pillar().startpoint();
+            // Search for new pillar locations:
+
+            bool   found    = false;
+            double alpha    = 0; // goes to 2Pi
+            double r        = 2 * m_cfg.base_radius_mm;
+            Vec3d  pillarsp = pillar().startpoint();
+
+            // temp value for starting point detection
             Vec3d sp(pillarsp(X), pillarsp(Y), pillarsp(Z) - r);
-            std::vector<bool> tv(needpillars, false);
-            std::vector<Vec3d> spts(needpillars);
 
-            while(!found && alpha < 2*PI) {
+            // A vector of bool for placement feasbility
+            std::vector<bool>  canplace(needpillars, false);
+            std::vector<Vec3d> spts(needpillars); // vector of starting points
 
-                for(unsigned n = 0; n < needpillars; n++) {
-                    double a = alpha + n * PI/3;
-                    Vec3d s = sp;
+            double gnd      = m_result.ground_level;
+            double min_dist = m_cfg.pillar_base_safety_distance_mm +
+                              m_cfg.base_radius_mm + EPSILON;
+            
+            while(!found && alpha < 2*PI) {
+                for (unsigned n = 0;
+                     n < needpillars && (!n || canplace[n - 1]);
+                     n++)
+                {
+                    double a = alpha + n * PI / 3;
+                    Vec3d  s = sp;
                     s(X) += std::cos(a) * r;
                     s(Y) += std::sin(a) * r;
                     spts[n] = s;
+                    
+                    // Check the path vertically down                    
                     auto hr = bridge_mesh_intersect(s, {0, 0, -1}, pillar().r);
-                    tv[n] = std::isinf(hr.distance());
+                    Vec3d gndsp{s(X), s(Y), gnd};
+                    
+                    // If the path is clear, check for pillar base collisions
+                    canplace[n] = std::isinf(hr.distance()) &&
+                                  std::sqrt(m_mesh.squared_distance(gndsp)) >
+                                      min_dist;
                 }
 
-                found = std::all_of(tv.begin(), tv.end(), [](bool v){return v;});
+                found = std::all_of(canplace.begin(), canplace.end(),
+                                    [](bool v) { return v; });
 
                 // 20 angles will be tried...
                 alpha += 0.1 * PI;
@@ -2010,7 +2365,7 @@ public:
             newpills.reserve(needpillars);
 
             if(found) for(unsigned n = 0; n < needpillars; n++) {
-                Vec3d s = spts[n]; double gnd = m_result.ground_level;
+                Vec3d s = spts[n]; 
                 Pillar p(s, Vec3d(s(X), s(Y), gnd), pillar().r);
                 p.add_base(m_cfg.base_height_mm, m_cfg.base_radius_mm);
 
@@ -2075,9 +2430,13 @@ public:
             // This is only for checking
             double idist = bridge_mesh_intersect(sph, dir, R, true);
             double dist = ray_mesh_intersect(sj, dir);
+            if (std::isinf(dist))
+                dist = sph(Z) - m_mesh.ground_level()
+                       + m_mesh.ground_level_offset();
 
-            if(std::isinf(idist) || std::isnan(idist) || idist < 2*R ||
-               std::isinf(dist)  || std::isnan(dist)   || dist < 2*R) {
+            if(std::isnan(idist) || idist < 2*R ||
+               std::isnan(dist)  || dist  < 2*R)
+            {
                 BOOST_LOG_TRIVIAL(warning) << "Can not find route for headless"
                                            << " support stick at: "
                                            << sj.transpose();
@@ -2085,9 +2444,11 @@ public:
             }
 
             Vec3d ej = sj + (dist + HWIDTH_MM)* dir;
-            m_result.add_compact_bridge(sp, ej, n, R);
+            m_result.add_compact_bridge(sp, ej, n, R, !std::isinf(dist));
         }
     }
+    
+    void merge_result() { m_result.merge_and_cleanup(); }
 };
 
 bool SLASupportTree::generate(const std::vector<SupportPoint> &support_points,
@@ -2096,9 +2457,9 @@ bool SLASupportTree::generate(const std::vector<SupportPoint> &support_points,
                               const Controller &ctl)
 {
     if(support_points.empty()) return false;
-
+    
     Algorithm alg(cfg, mesh, support_points, *m_impl, ctl.cancelfn);
-
+    
     // Let's define the individual steps of the processing. We can experiment
     // later with the ordering and the dependencies between them.
     enum Steps {
@@ -2110,55 +2471,58 @@ bool SLASupportTree::generate(const std::vector<SupportPoint> &support_points,
         ROUTING_NONGROUND,
         CASCADE_PILLARS,
         HEADLESS,
+        MERGE_RESULT,
         DONE,
         ABORT,
         NUM_STEPS
         //...
     };
-
+    
     // Collect the algorithm steps into a nice sequence
     std::array<std::function<void()>, NUM_STEPS> program = {
         [] () {
             // Begin...
             // Potentially clear up the shared data (not needed for now)
         },
-
+        
         std::bind(&Algorithm::filter, &alg),
-
+        
         std::bind(&Algorithm::add_pinheads, &alg),
-
+        
         std::bind(&Algorithm::classify, &alg),
-
+        
         std::bind(&Algorithm::routing_to_ground, &alg),
-
+        
         std::bind(&Algorithm::routing_to_model, &alg),
-
+        
         std::bind(&Algorithm::interconnect_pillars, &alg),
-
+        
         std::bind(&Algorithm::routing_headless, &alg),
-
+        
+        std::bind(&Algorithm::merge_result, &alg),
+        
         [] () {
             // Done
         },
-
+        
         [] () {
             // Abort
         }
     };
-
+    
     Steps pc = BEGIN;
-
+    
     if(cfg.ground_facing_only) {
         program[ROUTING_NONGROUND] = []() {
             BOOST_LOG_TRIVIAL(info)
-                    << "Skipping model-facing supports as requested.";
+                << "Skipping model-facing supports as requested.";
         };
         program[HEADLESS] = []() {
             BOOST_LOG_TRIVIAL(info) << "Skipping headless stick generation as"
                                        " requested.";
         };
     }
-
+    
     // Let's define a simple automaton that will run our program.
     auto progress = [&ctl, &pc] () {
         static const std::array<std::string, NUM_STEPS> stepstr {
@@ -2170,10 +2534,11 @@ bool SLASupportTree::generate(const std::vector<SupportPoint> &support_points,
             "Routing supports to model surface",
             "Interconnecting pillars",
             "Processing small holes",
+            "Merging support mesh",
             "Done",
             "Abort"
         };
-
+        
         static const std::array<unsigned, NUM_STEPS> stepstate {
             0,
             10,
@@ -2182,13 +2547,14 @@ bool SLASupportTree::generate(const std::vector<SupportPoint> &support_points,
             60,
             70,
             80,
-            90,
+            85,
+            99,
             100,
             0
         };
-
+        
         if(ctl.stopcondition()) pc = ABORT;
-
+        
         switch(pc) {
         case BEGIN: pc = FILTER; break;
         case FILTER: pc = PINHEADS; break;
@@ -2197,24 +2563,28 @@ bool SLASupportTree::generate(const std::vector<SupportPoint> &support_points,
         case ROUTING_GROUND: pc = ROUTING_NONGROUND; break;
         case ROUTING_NONGROUND: pc = CASCADE_PILLARS; break;
         case CASCADE_PILLARS: pc = HEADLESS; break;
-        case HEADLESS: pc = DONE; break;
+        case HEADLESS: pc = MERGE_RESULT; break;
+        case MERGE_RESULT: pc = DONE; break;
         case DONE:
         case ABORT: break;
         default: ;
         }
+        
         ctl.statuscb(stepstate[pc], stepstr[pc]);
     };
-
+    
     // Just here we run the computation...
     while(pc < DONE) {
         progress();
         program[pc]();
     }
-
+    
     return pc == ABORT;
 }
 
-SLASupportTree::SLASupportTree(): m_impl(new Impl()) {}
+SLASupportTree::SLASupportTree(double gnd_lvl): m_impl(new Impl()) {
+    m_impl->ground_level = gnd_lvl;
+}
 
 const TriangleMesh &SLASupportTree::merged_mesh() const
 {
@@ -2226,57 +2596,48 @@ void SLASupportTree::merged_mesh_with_pad(TriangleMesh &outmesh) const {
     outmesh.merge(get_pad());
 }
 
-SlicedSupports SLASupportTree::slice(float layerh, float init_layerh) const
+std::vector<ExPolygons> SLASupportTree::slice(
+    const std::vector<float> &heights, float cr) const
 {
-    if(init_layerh < 0) init_layerh = layerh;
-    auto& stree = get();
-
-    const auto modelh = float(stree.full_height());
-    auto gndlvl = float(this->m_impl->ground_level);
-    const Pad& pad = m_impl->pad();
-    if(!pad.empty()) gndlvl -= float(get_pad_elevation(pad.cfg));
-
-    std::vector<float> heights;
-    heights.reserve(size_t(modelh/layerh) + 1);
-
-    for(float h = gndlvl + init_layerh; h < gndlvl + modelh; h += layerh) {
-        heights.emplace_back(h);
+    const TriangleMesh &sup_mesh = m_impl->merged_mesh();
+    const TriangleMesh &pad_mesh = get_pad();
+    
+    std::vector<ExPolygons> sup_slices;
+    if (!sup_mesh.empty()) { 
+        TriangleMeshSlicer sup_slicer(&sup_mesh);
+        sup_slicer.closing_radius = cr;
+        sup_slicer.slice(heights, &sup_slices, m_impl->ctl().cancelfn);
     }
-
-    TriangleMesh fullmesh = m_impl->merged_mesh();
-    fullmesh.merge(get_pad());
-    fullmesh.require_shared_vertices(); // TriangleMeshSlicer needs this
-    TriangleMeshSlicer slicer(&fullmesh);
-    SlicedSupports ret;
-    slicer.slice(heights, &ret, get().ctl().cancelfn);
-
-    return ret;
-}
-
-SlicedSupports SLASupportTree::slice(const std::vector<float> &heights,
-                                     float cr) const
-{
-    TriangleMesh fullmesh = m_impl->merged_mesh();
-    fullmesh.merge(get_pad());
-    fullmesh.require_shared_vertices(); // TriangleMeshSlicer needs this
-    TriangleMeshSlicer slicer(cr, 0);
-    //TriangleMeshSlicer slicer(&fullmesh);
-    slicer.init(&fullmesh, []() {});
-    SlicedSupports ret;
-    slicer.slice(heights, &ret, get().ctl().cancelfn);
-
-    return ret;
+    
+    auto bb = pad_mesh.bounding_box();
+    auto maxzit = std::upper_bound(heights.begin(), heights.end(), bb.max.z());
+    
+    auto padgrid = reserve_vector<float>(heights.end() - maxzit);
+    std::copy(heights.begin(), maxzit, std::back_inserter(padgrid));
+    
+    std::vector<ExPolygons> pad_slices;
+    if (!pad_mesh.empty()) { 
+        TriangleMeshSlicer pad_slicer(&pad_mesh);
+        pad_slicer.closing_radius = cr;
+        pad_slicer.slice(padgrid, &pad_slices, m_impl->ctl().cancelfn);
+    }
+    
+    size_t len = std::min(heights.size(), pad_slices.size());
+    len = std::min(len, sup_slices.size());
+    
+    for (size_t i = 0; i < len; ++i) {
+        std::copy(pad_slices[i].begin(), pad_slices[i].end(),
+                  std::back_inserter(sup_slices[i]));
+        pad_slices[i] = {}; 
+    }
+    
+    return sup_slices;
 }
 
-const TriangleMesh &SLASupportTree::add_pad(const SliceLayer& baseplate,
+const TriangleMesh &SLASupportTree::add_pad(const ExPolygons& modelbase,
                                             const PoolConfig& pcfg) const
 {
-//    PoolConfig pcfg;
-//    pcfg.min_wall_thickness_mm = min_wall_thickness_mm;
-//    pcfg.min_wall_height_mm    = min_wall_height_mm;
-//    pcfg.max_merge_distance_mm = max_merge_distance_mm;
-//    pcfg.edge_radius_mm        = edge_radius_mm;
-    return m_impl->create_pad(merged_mesh(), baseplate, pcfg).tmesh;
+    return m_impl->create_pad(merged_mesh(), modelbase, pcfg).tmesh;
 }
 
 const TriangleMesh &SLASupportTree::get_pad() const
@@ -2297,16 +2658,6 @@ SLASupportTree::SLASupportTree(const std::vector<SupportPoint> &points,
 {
     m_impl->ground_level = emesh.ground_level() - cfg.object_elevation_mm;
     generate(points, emesh, cfg, ctl);
-    m_impl->clear_support_data();
-}
-
-SLASupportTree::SLASupportTree(const SLASupportTree &c):
-    m_impl(new Impl(*c.m_impl)) {}
-
-SLASupportTree &SLASupportTree::operator=(const SLASupportTree &c)
-{
-    m_impl = make_unique<Impl>(*c.m_impl);
-    return *this;
 }
 
 SLASupportTree::~SLASupportTree() {}
diff --git a/src/libslic3r/SLA/SLASupportTree.hpp b/src/libslic3r/SLA/SLASupportTree.hpp
index bb821534e..59485ba77 100644
--- a/src/libslic3r/SLA/SLASupportTree.hpp
+++ b/src/libslic3r/SLA/SLASupportTree.hpp
@@ -17,10 +17,11 @@ class TriangleMesh;
 class Model;
 class ModelInstance;
 class ModelObject;
+class Polygon;
 class ExPolygon;
 
-using SliceLayer = std::vector<ExPolygon>;
-using SlicedSupports = std::vector<SliceLayer>;
+using Polygons = std::vector<Polygon>;
+using ExPolygons = std::vector<ExPolygon>;
 
 namespace sla {
 
@@ -73,6 +74,10 @@ struct SupportConfig {
     // The elevation in Z direction upwards. This is the space between the pad
     // and the model object's bounding box bottom.
     double object_elevation_mm = 10;
+    
+    // The shortest distance between a pillar base perimeter from the model
+    // body. This is only useful when elevation is set to zero.
+    double pillar_base_safety_distance_mm = 0.5;
 
     // /////////////////////////////////////////////////////////////////////////
     // Compile time configuration values (candidates for runtime)
@@ -153,15 +158,15 @@ class SLASupportTree {
 
 public:
 
-    SLASupportTree();
+    SLASupportTree(double ground_level = 0.0);
 
     SLASupportTree(const std::vector<SupportPoint>& pts,
                    const EigenMesh3D& em,
                    const SupportConfig& cfg = {},
                    const Controller& ctl = {});
-
-    SLASupportTree(const SLASupportTree&);
-    SLASupportTree& operator=(const SLASupportTree&);
+    
+    SLASupportTree(const SLASupportTree&) = delete;
+    SLASupportTree& operator=(const SLASupportTree&) = delete;
 
     ~SLASupportTree();
 
@@ -171,13 +176,15 @@ public:
 
     void merged_mesh_with_pad(TriangleMesh&) const;
 
-    /// Get the sliced 2d layers of the support geometry.
-    SlicedSupports slice(float layerh, float init_layerh = -1.0) const;
-
-    SlicedSupports slice(const std::vector<float>&, float closing_radius) const;
+    std::vector<ExPolygons> slice(const std::vector<float> &,
+                                  float closing_radius) const;
 
     /// Adding the "pad" (base pool) under the supports
-    const TriangleMesh& add_pad(const SliceLayer& baseplate,
+    /// modelbase will be used according to the embed_object flag in PoolConfig.
+    /// If set, the plate will interpreted as the model's intrinsic pad. 
+    /// Otherwise, the modelbase will be unified with the base plate calculated
+    /// from the supports.
+    const TriangleMesh& add_pad(const ExPolygons& modelbase,
                                 const PoolConfig& pcfg) const;
 
     /// Get the pad geometry
diff --git a/src/libslic3r/SLA/SLASupportTreeIGL.cpp b/src/libslic3r/SLA/SLASupportTreeIGL.cpp
index fbbd645d1..4016b31f8 100644
--- a/src/libslic3r/SLA/SLASupportTreeIGL.cpp
+++ b/src/libslic3r/SLA/SLASupportTreeIGL.cpp
@@ -12,10 +12,18 @@
 #include "SLABoostAdapter.hpp"
 #include "boost/geometry/index/rtree.hpp"
 
+#ifdef _MSC_VER
+#pragma warning(push)
+#pragma warning(disable: 4244)
+#pragma warning(disable: 4267)
+#endif
 #include <igl/ray_mesh_intersect.h>
 #include <igl/point_mesh_squared_distance.h>
 #include <igl/remove_duplicate_vertices.h>
 #include <igl/signed_distance.h>
+#ifdef _MSC_VER
+#pragma warning(pop)
+#endif
 
 #include <tbb/parallel_for.h>
 
@@ -29,69 +37,137 @@ namespace sla {
 using igl::PI;
 
 /* **************************************************************************
- * SpatIndex implementation
+ * PointIndex implementation
  * ************************************************************************** */
 
-class SpatIndex::Impl {
+class PointIndex::Impl {
 public:
-    using BoostIndex = boost::geometry::index::rtree< SpatElement,
+    using BoostIndex = boost::geometry::index::rtree< PointIndexEl,
                        boost::geometry::index::rstar<16, 4> /* ? */ >;
 
     BoostIndex m_store;
 };
 
-SpatIndex::SpatIndex(): m_impl(new Impl()) {}
-SpatIndex::~SpatIndex() {}
+PointIndex::PointIndex(): m_impl(new Impl()) {}
+PointIndex::~PointIndex() {}
 
-SpatIndex::SpatIndex(const SpatIndex &cpy): m_impl(new Impl(*cpy.m_impl)) {}
-SpatIndex::SpatIndex(SpatIndex&& cpy): m_impl(std::move(cpy.m_impl)) {}
+PointIndex::PointIndex(const PointIndex &cpy): m_impl(new Impl(*cpy.m_impl)) {}
+PointIndex::PointIndex(PointIndex&& cpy): m_impl(std::move(cpy.m_impl)) {}
 
-SpatIndex& SpatIndex::operator=(const SpatIndex &cpy)
+PointIndex& PointIndex::operator=(const PointIndex &cpy)
 {
     m_impl.reset(new Impl(*cpy.m_impl));
     return *this;
 }
 
-SpatIndex& SpatIndex::operator=(SpatIndex &&cpy)
+PointIndex& PointIndex::operator=(PointIndex &&cpy)
 {
     m_impl.swap(cpy.m_impl);
     return *this;
 }
 
-void SpatIndex::insert(const SpatElement &el)
+void PointIndex::insert(const PointIndexEl &el)
 {
     m_impl->m_store.insert(el);
 }
 
-bool SpatIndex::remove(const SpatElement& el)
+bool PointIndex::remove(const PointIndexEl& el)
 {
     return m_impl->m_store.remove(el) == 1;
 }
 
-std::vector<SpatElement>
-SpatIndex::query(std::function<bool(const SpatElement &)> fn)
+std::vector<PointIndexEl>
+PointIndex::query(std::function<bool(const PointIndexEl &)> fn)
 {
     namespace bgi = boost::geometry::index;
 
-    std::vector<SpatElement> ret;
+    std::vector<PointIndexEl> 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)
+std::vector<PointIndexEl> PointIndex::nearest(const Vec3d &el, unsigned k = 1)
 {
     namespace bgi = boost::geometry::index;
-    std::vector<SpatElement> ret; ret.reserve(k);
+    std::vector<PointIndexEl> ret; ret.reserve(k);
     m_impl->m_store.query(bgi::nearest(el, k), std::back_inserter(ret));
     return ret;
 }
 
-size_t SpatIndex::size() const
+size_t PointIndex::size() const
 {
     return m_impl->m_store.size();
 }
 
-void SpatIndex::foreach(std::function<void (const SpatElement &)> fn)
+void PointIndex::foreach(std::function<void (const PointIndexEl &)> fn)
+{
+    for(auto& el : m_impl->m_store) fn(el);
+}
+
+/* **************************************************************************
+ * BoxIndex implementation
+ * ************************************************************************** */
+
+class BoxIndex::Impl {
+public:
+    using BoostIndex = boost::geometry::index::
+        rtree<BoxIndexEl, boost::geometry::index::rstar<16, 4> /* ? */>;
+
+    BoostIndex m_store;
+};
+
+BoxIndex::BoxIndex(): m_impl(new Impl()) {}
+BoxIndex::~BoxIndex() {}
+
+BoxIndex::BoxIndex(const BoxIndex &cpy): m_impl(new Impl(*cpy.m_impl)) {}
+BoxIndex::BoxIndex(BoxIndex&& cpy): m_impl(std::move(cpy.m_impl)) {}
+
+BoxIndex& BoxIndex::operator=(const BoxIndex &cpy)
+{
+    m_impl.reset(new Impl(*cpy.m_impl));
+    return *this;
+}
+
+BoxIndex& BoxIndex::operator=(BoxIndex &&cpy)
+{
+    m_impl.swap(cpy.m_impl);
+    return *this;
+}
+
+void BoxIndex::insert(const BoxIndexEl &el)
+{
+    m_impl->m_store.insert(el);
+}
+
+bool BoxIndex::remove(const BoxIndexEl& el)
+{
+    return m_impl->m_store.remove(el) == 1;
+}
+
+std::vector<BoxIndexEl> BoxIndex::query(const BoundingBox &qrbb,
+                                        BoxIndex::QueryType qt)
+{
+    namespace bgi = boost::geometry::index;
+    
+    std::vector<BoxIndexEl> ret; ret.reserve(m_impl->m_store.size());
+    
+    switch (qt) {
+    case qtIntersects:
+        m_impl->m_store.query(bgi::intersects(qrbb), std::back_inserter(ret));
+        break;
+    case qtWithin:
+        m_impl->m_store.query(bgi::within(qrbb), std::back_inserter(ret));
+    }
+    
+    return ret;
+}
+
+size_t BoxIndex::size() const
+{
+    return m_impl->m_store.size();
+}
+
+void BoxIndex::foreach(std::function<void (const BoxIndexEl &)> fn)
 {
     for(auto& el : m_impl->m_store) fn(el);
 }
@@ -343,12 +419,14 @@ PointSet normals(const PointSet& points,
     return ret;
 }
 namespace bgi = boost::geometry::index;
-using Index3D = bgi::rtree< SpatElement, bgi::rstar<16, 4> /* ? */ >;
+using Index3D = bgi::rtree< PointIndexEl, bgi::rstar<16, 4> /* ? */ >;
 
-ClusteredPoints cluster(Index3D& sindex, unsigned max_points,
-                        std::function<std::vector<SpatElement>(const Index3D&, const SpatElement&)> qfn)
+ClusteredPoints cluster(Index3D &sindex,
+                        unsigned max_points,
+                        std::function<std::vector<PointIndexEl>(
+                            const Index3D &, const PointIndexEl &)> qfn)
 {
-    using Elems = std::vector<SpatElement>;
+    using Elems = std::vector<PointIndexEl>;
 
     // Recursive function for visiting all the points in a given distance to
     // each other
@@ -356,8 +434,8 @@ ClusteredPoints cluster(Index3D& sindex, unsigned max_points,
     [&sindex, &group, max_points, qfn](Elems& pts, Elems& cluster)
     {
         for(auto& p : pts) {
-            std::vector<SpatElement> tmp = qfn(sindex, p);
-            auto cmp = [](const SpatElement& e1, const SpatElement& e2){
+            std::vector<PointIndexEl> tmp = qfn(sindex, p);
+            auto cmp = [](const PointIndexEl& e1, const PointIndexEl& e2){
                 return e1.second < e2.second;
             };
 
@@ -401,12 +479,12 @@ ClusteredPoints cluster(Index3D& sindex, unsigned max_points,
 }
 
 namespace {
-std::vector<SpatElement> distance_queryfn(const Index3D& sindex,
-                                          const SpatElement& p,
+std::vector<PointIndexEl> distance_queryfn(const Index3D& sindex,
+                                          const PointIndexEl& p,
                                           double dist,
                                           unsigned max_points)
 {
-    std::vector<SpatElement> tmp; tmp.reserve(max_points);
+    std::vector<PointIndexEl> tmp; tmp.reserve(max_points);
     sindex.query(
         bgi::nearest(p.first, max_points),
         std::back_inserter(tmp)
@@ -433,7 +511,7 @@ ClusteredPoints cluster(
     for(auto idx : indices) sindex.insert( std::make_pair(pointfn(idx), idx));
 
     return cluster(sindex, max_points,
-                   [dist, max_points](const Index3D& sidx, const SpatElement& p)
+                   [dist, max_points](const Index3D& sidx, const PointIndexEl& p)
     {
         return distance_queryfn(sidx, p, dist, max_points);
     });
@@ -443,7 +521,7 @@ ClusteredPoints cluster(
 ClusteredPoints cluster(
         const std::vector<unsigned>& indices,
         std::function<Vec3d(unsigned)> pointfn,
-        std::function<bool(const SpatElement&, const SpatElement&)> predicate,
+        std::function<bool(const PointIndexEl&, const PointIndexEl&)> predicate,
         unsigned max_points)
 {
     // A spatial index for querying the nearest points
@@ -453,10 +531,10 @@ ClusteredPoints cluster(
     for(auto idx : indices) sindex.insert( std::make_pair(pointfn(idx), idx));
 
     return cluster(sindex, max_points,
-        [max_points, predicate](const Index3D& sidx, const SpatElement& p)
+        [max_points, predicate](const Index3D& sidx, const PointIndexEl& p)
     {
-        std::vector<SpatElement> tmp; tmp.reserve(max_points);
-        sidx.query(bgi::satisfies([p, predicate](const SpatElement& e){
+        std::vector<PointIndexEl> tmp; tmp.reserve(max_points);
+        sidx.query(bgi::satisfies([p, predicate](const PointIndexEl& e){
             return predicate(p, e);
         }), std::back_inserter(tmp));
         return tmp;
@@ -473,7 +551,7 @@ ClusteredPoints cluster(const PointSet& pts, double dist, unsigned max_points)
         sindex.insert(std::make_pair(Vec3d(pts.row(i)), unsigned(i)));
 
     return cluster(sindex, max_points,
-                   [dist, max_points](const Index3D& sidx, const SpatElement& p)
+                   [dist, max_points](const Index3D& sidx, const PointIndexEl& p)
     {
         return distance_queryfn(sidx, p, dist, max_points);
     });