Add the full source of BambuStudio

using version 1.0.10

src/libslic3r/SLA/AGGRaster.hpp

@@ -0,0 +1,222 @@
+#ifndef AGGRASTER_HPP
+#define AGGRASTER_HPP
+
+#include <libslic3r/SLA/RasterBase.hpp>
+#include "libslic3r/ExPolygon.hpp"
+#include "libslic3r/MTUtils.hpp"
+
+// For rasterizing
+#include <agg/agg_basics.h>
+#include <agg/agg_rendering_buffer.h>
+#include <agg/agg_pixfmt_gray.h>
+#include <agg/agg_pixfmt_rgb.h>
+#include <agg/agg_renderer_base.h>
+#include <agg/agg_renderer_scanline.h>
+
+#include <agg/agg_scanline_p.h>
+#include <agg/agg_rasterizer_scanline_aa.h>
+#include <agg/agg_path_storage.h>
+
+namespace Slic3r {
+
+inline const Polygon& contour(const ExPolygon& p) { return p.contour; }
+inline const Polygons& holes(const ExPolygon& p) { return p.holes; }
+
+namespace sla {
+
+template<class Color> struct Colors {
+    static const Color White;
+    static const Color Black;
+};
+
+template<class Color> const Color Colors<Color>::White = Color{255};
+template<class Color> const Color Colors<Color>::Black = Color{0};
+
+template<class PixelRenderer,
+         template<class /*agg::renderer_base<PixelRenderer>*/> class Renderer,
+         class Rasterizer = agg::rasterizer_scanline_aa<>,
+         class Scanline   = agg::scanline_p8>
+class AGGRaster: public RasterBase {
+public:
+    using TColor = typename PixelRenderer::color_type;
+    using TValue = typename TColor::value_type;
+    using TPixel = typename PixelRenderer::pixel_type;
+    using TRawBuffer = agg::rendering_buffer;
+    
+protected:
+    
+    Resolution m_resolution;
+    PixelDim m_pxdim_scaled;    // used for scaled coordinate polygons
+    
+    std::vector<TPixel> m_buf;
+    agg::rendering_buffer m_rbuf;
+    
+    PixelRenderer m_pixrenderer;
+    
+    agg::renderer_base<PixelRenderer> m_raw_renderer;
+    Renderer<agg::renderer_base<PixelRenderer>> m_renderer;
+    
+    Trafo m_trafo;
+    Scanline m_scanlines;
+    Rasterizer m_rasterizer;
+    
+    void flipy(agg::path_storage &path) const
+    {
+        path.flip_y(0, double(m_resolution.height_px));
+    }
+    
+    void flipx(agg::path_storage &path) const
+    {
+        path.flip_x(0, double(m_resolution.width_px));
+    }
+    
+    double getPx(const Point &p) { return p(0) * m_pxdim_scaled.w_mm; }
+    double getPy(const Point &p) { return p(1) * m_pxdim_scaled.h_mm; }
+    agg::path_storage to_path(const Polygon &poly) { return to_path(poly.points); }
+    
+    template<class PointVec> agg::path_storage _to_path(const PointVec& v)
+    {
+        agg::path_storage path;
+        
+        auto it = v.begin();
+        path.move_to(getPx(*it), getPy(*it));
+        while(++it != v.end()) path.line_to(getPx(*it), getPy(*it));
+        path.line_to(getPx(v.front()), getPy(v.front()));
+        
+        return path;
+    }
+    
+    template<class PointVec> agg::path_storage _to_path_flpxy(const PointVec& v)
+    {
+        agg::path_storage path;
+        
+        auto it = v.begin();
+        path.move_to(getPy(*it), getPx(*it));
+        while(++it != v.end()) path.line_to(getPy(*it), getPx(*it));
+        path.line_to(getPy(v.front()), getPx(v.front()));
+        
+        return path;
+    }
+    
+    template<class PointVec> agg::path_storage to_path(const PointVec &v)
+    {
+        auto path = m_trafo.flipXY ? _to_path_flpxy(v) : _to_path(v);
+        
+        path.translate_all_paths(m_trafo.center_x * m_pxdim_scaled.w_mm,
+                                 m_trafo.center_y * m_pxdim_scaled.h_mm);
+        
+        if(m_trafo.mirror_x) flipx(path);
+        if(m_trafo.mirror_y) flipy(path);
+        
+        return path;
+    }
+    
+    template<class P> void _draw(const P &poly)
+    {
+        m_rasterizer.reset();
+        
+        m_rasterizer.add_path(to_path(contour(poly)));
+        for(auto& h : holes(poly)) m_rasterizer.add_path(to_path(h));
+        
+        agg::render_scanlines(m_rasterizer, m_scanlines, m_renderer);
+    }
+    
+public:
+    template<class GammaFn>
+    AGGRaster(const Resolution &res,
+              const PixelDim &  pd,
+              const Trafo &     trafo,
+              const TColor &    foreground,
+              const TColor &    background,
+              GammaFn &&        gammafn)
+        : m_resolution(res)
+        , m_pxdim_scaled(SCALING_FACTOR, SCALING_FACTOR)
+        , m_buf(res.pixels())
+        , m_rbuf(reinterpret_cast<TValue *>(m_buf.data()),
+                 unsigned(res.width_px),
+                 unsigned(res.height_px),
+                 int(res.width_px *PixelRenderer::num_components))
+        , m_pixrenderer(m_rbuf)
+        , m_raw_renderer(m_pixrenderer)
+        , m_renderer(m_raw_renderer)
+        , m_trafo(trafo)
+    {
+        // Visual Studio compiler gives warnings about possible division by zero.
+        assert(pd.w_mm != 0 && pd.h_mm != 0);
+        if (pd.w_mm != 0 && pd.h_mm != 0) {
+            m_pxdim_scaled.w_mm /= pd.w_mm;
+            m_pxdim_scaled.h_mm /= pd.h_mm;
+        }
+        m_renderer.color(foreground);
+        clear(background);
+        
+        m_rasterizer.gamma(gammafn);
+    }
+    
+    Trafo trafo() const override { return m_trafo; }
+    Resolution resolution() const override { return m_resolution; }
+    PixelDim   pixel_dimensions() const override
+    {
+        return {SCALING_FACTOR / m_pxdim_scaled.w_mm,
+                SCALING_FACTOR / m_pxdim_scaled.h_mm};
+    }
+    
+    void draw(const ExPolygon &poly) override { _draw(poly); }
+    
+    EncodedRaster encode(RasterEncoder encoder) const override
+    {
+        return encoder(m_buf.data(), m_resolution.width_px, m_resolution.height_px, 1);    
+    }
+    
+    void clear(const TColor color) { m_raw_renderer.clear(color); }
+};
+
+/*
+ * Captures an anti-aliased monochrome canvas where vectorial
+ * polygons can be rasterized. Fill color is always white and the background is
+ * black. Contours are anti-aliased.
+ * 
+ * A gamma function can be specified at compile time to make it more flexible.
+ */
+using _RasterGrayscaleAA =
+    AGGRaster<agg::pixfmt_gray8, agg::renderer_scanline_aa_solid>;
+
+class RasterGrayscaleAA : public _RasterGrayscaleAA {
+    using Base = _RasterGrayscaleAA;
+    using typename Base::TColor;
+    using typename Base::TValue;
+public:
+    template<class GammaFn>
+    RasterGrayscaleAA(const RasterBase::Resolution &res,
+                      const RasterBase::PixelDim &  pd,
+                      const RasterBase::Trafo &     trafo,
+                      GammaFn &&                    fn)
+        : Base(res, pd, trafo, Colors<TColor>::White, Colors<TColor>::Black,
+               std::forward<GammaFn>(fn))
+    {}
+    
+    uint8_t read_pixel(size_t col, size_t row) const
+    {
+        static_assert(std::is_same<TValue, uint8_t>::value, "Not grayscale pix");
+        
+        uint8_t px;
+        Base::m_buf[row * Base::resolution().width_px + col].get(px);
+        return px;
+    }
+    
+    void clear() { Base::clear(Colors<TColor>::Black); }
+};
+
+class RasterGrayscaleAAGammaPower: public RasterGrayscaleAA {
+public:
+    RasterGrayscaleAAGammaPower(const RasterBase::Resolution &res,
+                                const RasterBase::PixelDim &  pd,
+                                const RasterBase::Trafo &     trafo,
+                                double                        gamma = 1.)
+        : RasterGrayscaleAA(res, pd, trafo, agg::gamma_power(gamma))
+    {}
+};
+
+}} // namespace Slic3r::sla
+
+#endif // AGGRASTER_HPP

src/libslic3r/SLA/BoostAdapter.hpp

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

src/libslic3r/SLA/Clustering.cpp

@@ -0,0 +1,152 @@
+#include "Clustering.hpp"
+#include "boost/geometry/index/rtree.hpp"
+
+#include <libslic3r/SLA/SpatIndex.hpp>
+#include <libslic3r/SLA/BoostAdapter.hpp>
+
+namespace Slic3r { namespace sla {
+
+namespace bgi = boost::geometry::index;
+using Index3D = bgi::rtree< PointIndexEl, bgi::rstar<16, 4> /* ? */ >;
+
+namespace {
+
+bool cmp_ptidx_elements(const PointIndexEl& e1, const PointIndexEl& e2)
+{
+    return e1.second < e2.second;
+};
+
+ClusteredPoints cluster(Index3D &sindex,
+                        unsigned max_points,
+                        std::function<std::vector<PointIndexEl>(
+                            const Index3D &, const PointIndexEl &)> qfn)
+{
+    using Elems = std::vector<PointIndexEl>;
+
+    // Recursive function for visiting all the points in a given distance to
+    // each other
+    std::function<void(Elems&, Elems&)> group =
+        [&sindex, &group, max_points, qfn](Elems& pts, Elems& cluster)
+    {
+        for(auto& p : pts) {
+            std::vector<PointIndexEl> tmp = qfn(sindex, p);
+
+            std::sort(tmp.begin(), tmp.end(), cmp_ptidx_elements);
+
+            Elems newpts;
+            std::set_difference(tmp.begin(), tmp.end(),
+                                cluster.begin(), cluster.end(),
+                                std::back_inserter(newpts), cmp_ptidx_elements);
+
+            int c = max_points && newpts.size() + cluster.size() > max_points?
+                        int(max_points - cluster.size()) : int(newpts.size());
+
+            cluster.insert(cluster.end(), newpts.begin(), newpts.begin() + c);
+            std::sort(cluster.begin(), cluster.end(), cmp_ptidx_elements);
+
+            if(!newpts.empty() && (!max_points || cluster.size() < max_points))
+                group(newpts, cluster);
+        }
+    };
+
+    std::vector<Elems> clusters;
+    for(auto it = sindex.begin(); it != sindex.end();) {
+        Elems cluster = {};
+        Elems pts = {*it};
+        group(pts, cluster);
+
+        for(auto& c : cluster) sindex.remove(c);
+        it = sindex.begin();
+
+        clusters.emplace_back(cluster);
+    }
+
+    ClusteredPoints result;
+    for(auto& cluster : clusters) {
+        result.emplace_back();
+        for(auto c : cluster) result.back().emplace_back(c.second);
+    }
+
+    return result;
+}
+
+std::vector<PointIndexEl> distance_queryfn(const Index3D& sindex,
+                                           const PointIndexEl& p,
+                                           double dist,
+                                           unsigned max_points)
+{
+    std::vector<PointIndexEl> tmp; tmp.reserve(max_points);
+    sindex.query(
+        bgi::nearest(p.first, max_points),
+        std::back_inserter(tmp)
+        );
+
+    for(auto it = tmp.begin(); it < tmp.end(); ++it)
+        if((p.first - it->first).norm() > dist) it = tmp.erase(it);
+
+    return tmp;
+}
+
+} // namespace
+
+// Clustering a set of points by the given criteria
+ClusteredPoints cluster(
+    const std::vector<unsigned>& indices,
+    std::function<Vec3d(unsigned)> pointfn,
+    double dist,
+    unsigned max_points)
+{
+    // A spatial index for querying the nearest points
+    Index3D sindex;
+
+    // Build the index
+    for(auto idx : indices) sindex.insert( std::make_pair(pointfn(idx), idx));
+
+    return cluster(sindex, max_points,
+                   [dist, max_points](const Index3D& sidx, const PointIndexEl& p)
+                   {
+                       return distance_queryfn(sidx, p, dist, max_points);
+                   });
+}
+
+// Clustering a set of points by the given criteria
+ClusteredPoints cluster(
+    const std::vector<unsigned>& indices,
+    std::function<Vec3d(unsigned)> pointfn,
+    std::function<bool(const PointIndexEl&, const PointIndexEl&)> predicate,
+    unsigned max_points)
+{
+    // A spatial index for querying the nearest points
+    Index3D sindex;
+
+    // Build the index
+    for(auto idx : indices) sindex.insert( std::make_pair(pointfn(idx), idx));
+
+    return cluster(sindex, max_points,
+                   [max_points, predicate](const Index3D& sidx, const PointIndexEl& p)
+                   {
+                       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;
+                   });
+}
+
+ClusteredPoints cluster(const Eigen::MatrixXd& pts, double dist, unsigned max_points)
+{
+    // A spatial index for querying the nearest points
+    Index3D sindex;
+
+    // Build the index
+    for(Eigen::Index i = 0; i < pts.rows(); i++)
+        sindex.insert(std::make_pair(Vec3d(pts.row(i)), unsigned(i)));
+
+    return cluster(sindex, max_points,
+                   [dist, max_points](const Index3D& sidx, const PointIndexEl& p)
+                   {
+                       return distance_queryfn(sidx, p, dist, max_points);
+                   });
+}
+
+}} // namespace Slic3r::sla

src/libslic3r/SLA/Clustering.hpp

@@ -0,0 +1,82 @@
+#ifndef SLA_CLUSTERING_HPP
+#define SLA_CLUSTERING_HPP
+
+#include <vector>
+
+#include <libslic3r/Point.hpp>
+#include <libslic3r/SLA/SpatIndex.hpp>
+
+namespace Slic3r { namespace sla {
+
+using ClusterEl = std::vector<unsigned>;
+using ClusteredPoints = std::vector<ClusterEl>;
+
+// Clustering a set of points by the given distance.
+ClusteredPoints cluster(const std::vector<unsigned>& indices,
+                        std::function<Vec3d(unsigned)> pointfn,
+                        double dist,
+                        unsigned max_points);
+
+ClusteredPoints cluster(const Eigen::MatrixXd& points,
+                        double dist,
+                        unsigned max_points);
+
+ClusteredPoints cluster(
+    const std::vector<unsigned>& indices,
+    std::function<Vec3d(unsigned)> pointfn,
+    std::function<bool(const PointIndexEl&, const PointIndexEl&)> predicate,
+    unsigned max_points);
+
+// This function returns the position of the centroid in the input 'clust'
+// vector of point indices.
+template<class DistFn, class PointFn>
+long cluster_centroid(const ClusterEl &clust, PointFn pointfn, DistFn df)
+{
+    switch(clust.size()) {
+    case 0: /* empty cluster */ return -1;
+    case 1: /* only one element */ return 0;
+    case 2: /* if two elements, there is no center */ return 0;
+    default: ;
+    }
+
+    // The function works by calculating for each point the average distance
+    // from all the other points in the cluster. We create a selector bitmask of
+    // the same size as the cluster. The bitmask will have two true bits and
+    // false bits for the rest of items and we will loop through all the
+    // permutations of the bitmask (combinations of two points). Get the
+    // distance for the two points and add the distance to the averages.
+    // The point with the smallest average than wins.
+
+    // The complexity should be O(n^2) but we will mostly apply this function
+    // for small clusters only (cca 3 elements)
+
+    std::vector<bool> sel(clust.size(), false);   // create full zero bitmask
+    std::fill(sel.end() - 2, sel.end(), true);    // insert the two ones
+    std::vector<double> avgs(clust.size(), 0.0);  // store the average distances
+
+    do {
+        std::array<size_t, 2> idx;
+        for(size_t i = 0, j = 0; i < clust.size(); i++)
+            if(sel[i]) idx[j++] = i;
+
+        double d = df(pointfn(clust[idx[0]]),
+                      pointfn(clust[idx[1]]));
+
+        // add the distance to the sums for both associated points
+        for(auto i : idx) avgs[i] += d;
+
+        // now continue with the next permutation of the bitmask with two 1s
+    } while(std::next_permutation(sel.begin(), sel.end()));
+
+    // Divide by point size in the cluster to get the average (may be redundant)
+    for(auto& a : avgs) a /= clust.size();
+
+    // get the lowest average distance and return the index
+    auto minit = std::min_element(avgs.begin(), avgs.end());
+    return long(minit - avgs.begin());
+}
+
+
+}} // namespace Slic3r::sla
+
+#endif // CLUSTERING_HPP

src/libslic3r/SLA/ConcaveHull.cpp

@@ -0,0 +1,145 @@
+#include <libslic3r/SLA/ConcaveHull.hpp>
+#include <libslic3r/SLA/SpatIndex.hpp>
+
+#include <libslic3r/MTUtils.hpp>
+#include <libslic3r/ClipperUtils.hpp>
+
+#include <boost/log/trivial.hpp>
+
+namespace Slic3r {
+namespace sla {
+
+inline Vec3d to_vec3(const Vec2crd &v2) { return {double(v2(X)), double(v2(Y)), 0.}; }
+inline Vec3d to_vec3(const Vec2d &v2) { return {v2(X), v2(Y), 0.}; }
+inline Vec2crd to_vec2(const Vec3d &v3) { return {coord_t(v3(X)), coord_t(v3(Y))}; }
+
+Point ConcaveHull::centroid(const Points &pp)
+{
+    Point c;
+    switch(pp.size()) {
+    case 0: break;
+    case 1: c = pp.front(); break;
+    case 2: c = (pp[0] + pp[1]) / 2; break;
+    default: {
+        auto MAX = std::numeric_limits<Point::coord_type>::max();
+        auto MIN = std::numeric_limits<Point::coord_type>::min();
+        Point min = {MAX, MAX}, max = {MIN, MIN};
+
+        for(auto& p : pp) {
+            if(p(0) < min(0)) min(0) = p(0);
+            if(p(1) < min(1)) min(1) = p(1);
+            if(p(0) > max(0)) max(0) = p(0);
+            if(p(1) > max(1)) max(1) = p(1);
+        }
+        c(0) = min(0) + (max(0) - min(0)) / 2;
+        c(1) = min(1) + (max(1) - min(1)) / 2;
+        break;
+    }
+    }
+
+    return c;
+}
+
+Points ConcaveHull::calculate_centroids() const
+{
+    // We get the centroids of all the islands in the 2D slice
+    Points centroids = reserve_vector<Point>(m_polys.size());
+    std::transform(m_polys.begin(), m_polys.end(),
+                   std::back_inserter(centroids),
+                   [](const Polygon &poly) { return centroid(poly); });
+
+    return centroids;
+}
+
+void ConcaveHull::merge_polygons() { m_polys = get_contours(union_ex(m_polys)); }
+
+void ConcaveHull::add_connector_rectangles(const Points &centroids,
+                                           coord_t       max_dist,
+                                           ThrowOnCancel thr)
+{
+    // Centroid of the centroids of islands. This is where the additional
+    // connector sticks are routed.
+    Point cc = centroid(centroids);
+
+    PointIndex ctrindex;
+    unsigned  idx = 0;
+    for(const Point &ct : centroids) ctrindex.insert(to_vec3(ct), idx++);
+
+    m_polys.reserve(m_polys.size() + centroids.size());
+
+    idx = 0;
+    for (const Point &c : centroids) {
+        thr();
+
+        double dx = c.x() - cc.x(), dy = c.y() - cc.y();
+        double l  = std::sqrt(dx * dx + dy * dy);
+        double nx = dx / l, ny = dy / l;
+
+        const Point &ct = centroids[idx];
+
+        std::vector<PointIndexEl> result = ctrindex.nearest(to_vec3(ct), 2);
+
+        double dist = max_dist;
+        for (const PointIndexEl &el : result)
+            if (el.second != idx) {
+                dist = Line(to_vec2(el.first), ct).length();
+                break;
+            }
+
+        idx++;
+
+        if (dist >= max_dist) return;
+
+        Polygon r;
+        r.points.reserve(3);
+        r.points.emplace_back(cc);
+
+        Point n(scaled(nx), scaled(ny));
+        r.points.emplace_back(c + Point(n.y(), -n.x()));
+        r.points.emplace_back(c + Point(-n.y(), n.x()));
+        offset(r, scaled<float>(1.));
+
+        m_polys.emplace_back(r);
+    }
+}
+
+ConcaveHull::ConcaveHull(const Polygons &polys, double mergedist, ThrowOnCancel thr)
+{
+    if(polys.empty()) return;
+
+    m_polys = polys;
+    merge_polygons();
+
+    if(m_polys.size() == 1) return;
+
+    Points centroids = calculate_centroids();
+
+    add_connector_rectangles(centroids, scaled(mergedist), thr);
+
+    merge_polygons();
+}
+
+ExPolygons ConcaveHull::to_expolygons() const
+{
+    auto ret = reserve_vector<ExPolygon>(m_polys.size());
+    for (const Polygon &p : m_polys) ret.emplace_back(ExPolygon(p));
+    return ret;
+}
+
+ExPolygons offset_waffle_style_ex(const ConcaveHull &hull, coord_t delta)
+{
+    return to_expolygons(offset_waffle_style(hull, delta));
+}
+
+Polygons offset_waffle_style(const ConcaveHull &hull, coord_t delta)
+{
+    auto arc_tolerance = scaled<double>(0.01);
+    Polygons res = closing(hull.polygons(), 2 * delta, delta, ClipperLib::jtRound, arc_tolerance);
+
+    auto it = std::remove_if(res.begin(), res.end(), [](Polygon &p) { return p.is_clockwise(); });
+    res.erase(it, res.end());
+
+    return res;
+}
+
+}} // namespace Slic3r::sla

src/libslic3r/SLA/ConcaveHull.hpp

@@ -0,0 +1,53 @@
+#ifndef SLA_CONCAVEHULL_HPP
+#define SLA_CONCAVEHULL_HPP
+
+#include <libslic3r/ExPolygon.hpp>
+
+namespace Slic3r {
+namespace sla {
+
+inline Polygons get_contours(const ExPolygons &poly)
+{
+    Polygons ret; ret.reserve(poly.size());
+    for (const ExPolygon &p : poly) ret.emplace_back(p.contour);
+
+    return ret;
+}
+
+using ThrowOnCancel = std::function<void()>;
+
+/// 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...)
+class ConcaveHull {
+    Polygons m_polys;
+
+    static Point centroid(const Points& pp);
+
+    static inline Point centroid(const Polygon &poly) { return poly.centroid(); }
+
+    Points calculate_centroids() const;
+
+    void merge_polygons();
+
+    void add_connector_rectangles(const Points &centroids,
+                                  coord_t       max_dist,
+                                  ThrowOnCancel thr);
+public:
+
+    ConcaveHull(const ExPolygons& polys, double merge_dist, ThrowOnCancel thr)
+        : ConcaveHull{to_polygons(polys), merge_dist, thr} {}
+
+    ConcaveHull(const Polygons& polys, double mergedist, ThrowOnCancel thr);
+
+    const Polygons & polygons() const { return m_polys; }
+
+    ExPolygons to_expolygons() const;
+};
+
+ExPolygons offset_waffle_style_ex(const ConcaveHull &ccvhull, coord_t delta);
+Polygons   offset_waffle_style(const ConcaveHull &polys, coord_t delta);
+
+}}     // namespace Slic3r::sla
+#endif // CONCAVEHULL_HPP

src/libslic3r/SLA/Concurrency.hpp

@@ -0,0 +1,70 @@
+#ifndef SLA_CONCURRENCY_H
+#define SLA_CONCURRENCY_H
+
+// FIXME: Deprecated
+
+#include <libslic3r/Execution/ExecutionSeq.hpp>
+#include <libslic3r/Execution/ExecutionTBB.hpp>
+
+namespace Slic3r {
+namespace sla {
+
+// 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 = true;
+
+template<bool> struct _ccr {};
+
+template<> struct _ccr<true>
+{
+    using SpinningMutex = execution::SpinningMutex<ExecutionTBB>;
+    using BlockingMutex = execution::BlockingMutex<ExecutionTBB>;
+
+    template<class It, class Fn>
+    static void for_each(It from, It to, Fn &&fn, size_t granularity = 1)
+    {
+        execution::for_each(ex_tbb, from, to, std::forward<Fn>(fn), granularity);
+    }
+
+    template<class...Args>
+    static auto reduce(Args&&...args)
+    {
+        return execution::reduce(ex_tbb, std::forward<Args>(args)...);
+    }
+
+    static size_t max_concurreny()
+    {
+        return execution::max_concurrency(ex_tbb);
+    }
+};
+
+template<> struct _ccr<false>
+{
+    using SpinningMutex = execution::SpinningMutex<ExecutionSeq>;
+    using BlockingMutex = execution::BlockingMutex<ExecutionSeq>;
+
+    template<class It, class Fn>
+    static void for_each(It from, It to, Fn &&fn, size_t granularity = 1)
+    {
+        execution::for_each(ex_seq, from, to, std::forward<Fn>(fn), granularity);
+    }
+
+    template<class...Args>
+    static auto reduce(Args&&...args)
+    {
+        return execution::reduce(ex_seq, std::forward<Args>(args)...);
+    }
+
+    static size_t max_concurreny()
+    {
+        return execution::max_concurrency(ex_seq);
+    }
+};
+
+using ccr = _ccr<USE_FULL_CONCURRENCY>;
+using ccr_seq = _ccr<false>;
+using ccr_par = _ccr<true>;
+
+}} // namespace Slic3r::sla
+
+#endif // SLACONCURRENCY_H

src/libslic3r/SLA/Hollowing.cpp

@@ -0,0 +1,563 @@
+#include <functional>
+#include <optional>
+
+#include <libslic3r/OpenVDBUtils.hpp>
+#include <libslic3r/TriangleMesh.hpp>
+#include <libslic3r/TriangleMeshSlicer.hpp>
+#include <libslic3r/SLA/Hollowing.hpp>
+#include <libslic3r/SLA/IndexedMesh.hpp>
+#include <libslic3r/ClipperUtils.hpp>
+#include <libslic3r/QuadricEdgeCollapse.hpp>
+#include <libslic3r/SLA/SupportTreeMesher.hpp>
+
+#include <boost/log/trivial.hpp>
+
+#include <libslic3r/MTUtils.hpp>
+#include <libslic3r/I18N.hpp>
+
+//! macro used to mark string used at localization,
+//! return same string
+#define L(s) Slic3r::I18N::translate(s)
+
+namespace Slic3r {
+namespace sla {
+
+struct Interior {
+    indexed_triangle_set mesh;
+    openvdb::FloatGrid::Ptr gridptr;
+    mutable std::optional<openvdb::FloatGrid::ConstAccessor> accessor;
+
+    double closing_distance = 0.;
+    double thickness = 0.;
+    double voxel_scale = 1.;
+    double nb_in = 3.;  // narrow band width inwards
+    double nb_out = 3.; // narrow band width outwards
+    // Full narrow band is the sum of the two above values.
+
+    void reset_accessor() const  // This resets the accessor and its cache
+    // Not a thread safe call!
+    {
+        if (gridptr)
+            accessor = gridptr->getConstAccessor();
+    }
+};
+
+void InteriorDeleter::operator()(Interior *p)
+{
+    delete p;
+}
+
+indexed_triangle_set &get_mesh(Interior &interior)
+{
+    return interior.mesh;
+}
+
+const indexed_triangle_set &get_mesh(const Interior &interior)
+{
+    return interior.mesh;
+}
+
+static InteriorPtr generate_interior_verbose(const TriangleMesh & mesh,
+                                             const JobController &ctl,
+                                             double min_thickness,
+                                             double voxel_scale,
+                                             double closing_dist)
+{
+    double offset = voxel_scale * min_thickness;
+    double D = voxel_scale * closing_dist;
+    float  out_range = 0.1f * float(offset);
+    float  in_range = 1.1f * float(offset + D);
+
+    if (ctl.stopcondition()) return {};
+    else ctl.statuscb(0, L("Hollowing"));
+
+    auto gridptr = mesh_to_grid(mesh.its, {}, voxel_scale, out_range, in_range);
+
+    assert(gridptr);
+
+    if (!gridptr) {
+        BOOST_LOG_TRIVIAL(error) << "Returned OpenVDB grid is NULL";
+        return {};
+    }
+
+    if (ctl.stopcondition()) return {};
+    else ctl.statuscb(30, L("Hollowing"));
+
+    double iso_surface = D;
+    auto   narrowb = double(in_range);
+    gridptr = redistance_grid(*gridptr, -(offset + D), narrowb, narrowb);
+
+    if (ctl.stopcondition()) return {};
+    else ctl.statuscb(70, L("Hollowing"));
+
+    double adaptivity = 0.;
+    InteriorPtr interior = InteriorPtr{new Interior{}};
+
+    interior->mesh = grid_to_mesh(*gridptr, iso_surface, adaptivity);
+    interior->gridptr = gridptr;
+
+    if (ctl.stopcondition()) return {};
+    else ctl.statuscb(100, L("Hollowing"));
+
+    interior->closing_distance = D;
+    interior->thickness = offset;
+    interior->voxel_scale = voxel_scale;
+    interior->nb_in = narrowb;
+    interior->nb_out = narrowb;
+
+    return interior;
+}
+
+InteriorPtr generate_interior(const TriangleMesh &   mesh,
+                              const HollowingConfig &hc,
+                              const JobController &  ctl)
+{
+    static const double MIN_OVERSAMPL = 3.5;
+    static const double MAX_OVERSAMPL = 8.;
+
+    // I can't figure out how to increase the grid resolution through openvdb
+    // API so the model will be scaled up before conversion and the result
+    // scaled down. Voxels have a unit size. If I set voxelSize smaller, it
+    // scales the whole geometry down, and doesn't increase the number of
+    // voxels.
+    //
+    // max 8x upscale, min is native voxel size
+    auto voxel_scale = MIN_OVERSAMPL + (MAX_OVERSAMPL - MIN_OVERSAMPL) * hc.quality;
+
+    InteriorPtr interior =
+        generate_interior_verbose(mesh, ctl, hc.min_thickness, voxel_scale,
+                                  hc.closing_distance);
+
+    if (interior && !interior->mesh.empty()) {
+
+        // flip normals back...
+        swap_normals(interior->mesh);
+
+        // simplify mesh lossless
+        float loss_less_max_error = 2*std::numeric_limits<float>::epsilon();
+        its_quadric_edge_collapse(interior->mesh, 0U, &loss_less_max_error);
+
+        its_compactify_vertices(interior->mesh);
+        its_merge_vertices(interior->mesh);
+
+        // flip normals back...
+        swap_normals(interior->mesh);
+    }
+
+    return interior;
+}
+
+indexed_triangle_set DrainHole::to_mesh() const
+{
+    auto r = double(radius);
+    auto h = double(height);
+    indexed_triangle_set hole = sla::cylinder(r, h, steps);
+    Eigen::Quaternionf q;
+    q.setFromTwoVectors(Vec3f{0.f, 0.f, 1.f}, normal);
+    for(auto& p : hole.vertices) p = q * p + pos;
+    
+    return hole;
+}
+
+bool DrainHole::operator==(const DrainHole &sp) const
+{
+    return (pos == sp.pos) && (normal == sp.normal) &&
+            is_approx(radius, sp.radius) &&
+            is_approx(height, sp.height);
+}
+
+bool DrainHole::is_inside(const Vec3f& pt) const
+{
+    Eigen::Hyperplane<float, 3> plane(normal, pos);
+    float dist = plane.signedDistance(pt);
+    if (dist < float(EPSILON) || dist > height)
+        return false;
+
+    Eigen::ParametrizedLine<float, 3> axis(pos, normal);
+    if ( axis.squaredDistance(pt) < pow(radius, 2.f))
+        return true;
+
+    return false;
+}
+
+
+// Given a line s+dir*t, find parameter t of intersections with the hole
+// and the normal (points inside the hole). Outputs through out reference,
+// returns true if two intersections were found.
+bool DrainHole::get_intersections(const Vec3f& s, const Vec3f& dir,
+                                  std::array<std::pair<float, Vec3d>, 2>& out)
+                                  const
+{
+    assert(is_approx(normal.norm(), 1.f));
+    const Eigen::ParametrizedLine<float, 3> ray(s, dir.normalized());
+
+    for (size_t i=0; i<2; ++i)
+        out[i] = std::make_pair(sla::IndexedMesh::hit_result::infty(), Vec3d::Zero());
+
+    const float sqr_radius = pow(radius, 2.f);
+
+    // first check a bounding sphere of the hole:
+    Vec3f center = pos+normal*height/2.f;
+    float sqr_dist_limit = pow(height/2.f, 2.f) + sqr_radius ;
+    if (ray.squaredDistance(center) > sqr_dist_limit)
+        return false;
+
+    // The line intersects the bounding sphere, look for intersections with
+    // bases of the cylinder.
+
+    size_t found = 0; // counts how many intersections were found
+    Eigen::Hyperplane<float, 3> base;
+    if (! is_approx(ray.direction().dot(normal), 0.f)) {
+        for (size_t i=1; i<=1; --i) {
+            Vec3f cylinder_center = pos+i*height*normal;
+            if (i == 0) {
+                // The hole base can be identical to mesh surface if it is flat
+                // let's better move the base outward a bit
+                cylinder_center -= EPSILON*normal;
+            }
+            base = Eigen::Hyperplane<float, 3>(normal, cylinder_center);
+            Vec3f intersection = ray.intersectionPoint(base);
+            // Only accept the point if it is inside the cylinder base.
+            if ((cylinder_center-intersection).squaredNorm() < sqr_radius) {
+                out[found].first = ray.intersectionParameter(base);
+                out[found].second = (i==0 ? 1. : -1.) * normal.cast<double>();
+                ++found;
+            }
+        }
+    }
+    else
+    {
+        // In case the line was perpendicular to the cylinder axis, previous
+        // block was skipped, but base will later be assumed to be valid.
+        base = Eigen::Hyperplane<float, 3>(normal, pos-EPSILON*normal);
+    }
+
+    // In case there is still an intersection to be found, check the wall
+    if (found != 2 && ! is_approx(std::abs(ray.direction().dot(normal)), 1.f)) {
+        // Project the ray onto the base plane
+        Vec3f proj_origin = base.projection(ray.origin());
+        Vec3f proj_dir = base.projection(ray.origin()+ray.direction())-proj_origin;
+        // save how the parameter scales and normalize the projected direction
+        float par_scale = proj_dir.norm();
+        proj_dir = proj_dir/par_scale;
+        Eigen::ParametrizedLine<float, 3> projected_ray(proj_origin, proj_dir);
+        // Calculate point on the secant that's closest to the center
+        // and its distance to the circle along the projected line
+        Vec3f closest = projected_ray.projection(pos);
+        float dist = sqrt((sqr_radius - (closest-pos).squaredNorm()));
+        // Unproject both intersections on the original line and check
+        // they are on the cylinder and not past it:
+        for (int i=-1; i<=1 && found !=2; i+=2) {
+            Vec3f isect = closest + i*dist * projected_ray.direction();
+            Vec3f to_isect = isect-proj_origin;
+            float par = to_isect.norm() / par_scale;
+            if (to_isect.normalized().dot(proj_dir.normalized()) < 0.f)
+                par *= -1.f;
+            Vec3d hit_normal = (pos-isect).normalized().cast<double>();
+            isect = ray.pointAt(par);
+            // check that the intersection is between the base planes:
+            float vert_dist = base.signedDistance(isect);
+            if (vert_dist > 0.f && vert_dist < height) {
+                out[found].first = par;
+                out[found].second = hit_normal;
+                ++found;
+            }
+        }
+    }
+
+    // If only one intersection was found, it is some corner case,
+    // no intersection will be returned:
+    if (found != 2)
+        return false;
+
+    // Sort the intersections:
+    if (out[0].first > out[1].first)
+        std::swap(out[0], out[1]);
+
+    return true;
+}
+
+void cut_drainholes(std::vector<ExPolygons> & obj_slices,
+                    const std::vector<float> &slicegrid,
+                    float                     closing_radius,
+                    const sla::DrainHoles &   holes,
+                    std::function<void(void)> thr)
+{
+    TriangleMesh mesh;
+    for (const sla::DrainHole &holept : holes)
+        mesh.merge(TriangleMesh{holept.to_mesh()});
+    
+    if (mesh.empty()) return;
+    
+    std::vector<ExPolygons> hole_slices = slice_mesh_ex(mesh.its, slicegrid, closing_radius, thr);
+    
+    if (obj_slices.size() != hole_slices.size())
+        BOOST_LOG_TRIVIAL(warning)
+            << "Sliced object and drain-holes layer count does not match!";
+
+    size_t until = std::min(obj_slices.size(), hole_slices.size());
+    
+    for (size_t i = 0; i < until; ++i)
+        obj_slices[i] = diff_ex(obj_slices[i], hole_slices[i]);
+}
+
+void hollow_mesh(TriangleMesh &mesh, const HollowingConfig &cfg, int flags)
+{
+    InteriorPtr interior = generate_interior(mesh, cfg, JobController{});
+    if (!interior) return;
+
+    hollow_mesh(mesh, *interior, flags);
+}
+
+void hollow_mesh(TriangleMesh &mesh, const Interior &interior, int flags)
+{
+    if (mesh.empty() || interior.mesh.empty()) return;
+
+    if (flags & hfRemoveInsideTriangles && interior.gridptr)
+        remove_inside_triangles(mesh, interior);
+
+    mesh.merge(TriangleMesh{interior.mesh});
+}
+
+// Get the distance of p to the interior's zero iso_surface. Interior should
+// have its zero isosurface positioned at offset + closing_distance inwards form
+// the model surface.
+static double get_distance_raw(const Vec3f &p, const Interior &interior)
+{
+    assert(interior.gridptr);
+
+    if (!interior.accessor) interior.reset_accessor();
+
+    auto v       = (p * interior.voxel_scale).cast<double>();
+    auto grididx = interior.gridptr->transform().worldToIndexCellCentered(
+        {v.x(), v.y(), v.z()});
+
+    return interior.accessor->getValue(grididx) ;
+}
+
+struct TriangleBubble { Vec3f center; double R; };
+
+// Return the distance of bubble center to the interior boundary or NaN if the
+// triangle is too big to be measured.
+static double get_distance(const TriangleBubble &b, const Interior &interior)
+{
+    double R = b.R * interior.voxel_scale;
+    double D = get_distance_raw(b.center, interior);
+
+    return (D > 0. && R >= interior.nb_out) ||
+           (D < 0. && R >= interior.nb_in)  ||
+           ((D - R) < 0. && 2 * R > interior.thickness) ?
+                std::nan("") :
+                // FIXME: Adding interior.voxel_scale is a compromise supposed
+                // to prevent the deletion of the triangles forming the interior
+                // itself. This has a side effect that a small portion of the
+                // bad triangles will still be visible.
+                D - interior.closing_distance /*+ 2 * interior.voxel_scale*/;
+}
+
+double get_distance(const Vec3f &p, const Interior &interior)
+{
+    double d = get_distance_raw(p, interior) - interior.closing_distance;
+    return d / interior.voxel_scale;
+}
+
+// A face that can be divided. Stores the indices into the original mesh if its
+// part of that mesh and the vertices it consists of.
+enum { NEW_FACE = -1};
+struct DivFace {
+    Vec3i indx;
+    std::array<Vec3f, 3> verts;
+    long faceid = NEW_FACE;
+    long parent = NEW_FACE;
+};
+
+// Divide a face recursively and call visitor on all the sub-faces.
+template<class Fn>
+void divide_triangle(const DivFace &face, Fn &&visitor)
+{
+    std::array<Vec3f, 3> edges = {(face.verts[0] - face.verts[1]),
+                                  (face.verts[1] - face.verts[2]),
+                                  (face.verts[2] - face.verts[0])};
+
+    std::array<size_t, 3> edgeidx = {0, 1, 2};
+
+    std::sort(edgeidx.begin(), edgeidx.end(), [&edges](size_t e1, size_t e2) {
+        return edges[e1].squaredNorm() > edges[e2].squaredNorm();
+    });
+
+    DivFace child1, child2;
+
+    child1.parent   = face.faceid == NEW_FACE ? face.parent : face.faceid;
+    child1.indx(0)  = -1;
+    child1.indx(1)  = face.indx(edgeidx[1]);
+    child1.indx(2)  = face.indx((edgeidx[1] + 1) % 3);
+    child1.verts[0] = (face.verts[edgeidx[0]] + face.verts[(edgeidx[0] + 1) % 3]) / 2.;
+    child1.verts[1] = face.verts[edgeidx[1]];
+    child1.verts[2] = face.verts[(edgeidx[1] + 1) % 3];
+
+    if (visitor(child1))
+        divide_triangle(child1, std::forward<Fn>(visitor));
+
+    child2.parent   = face.faceid == NEW_FACE ? face.parent : face.faceid;
+    child2.indx(0)  = -1;
+    child2.indx(1)  = face.indx(edgeidx[2]);
+    child2.indx(2)  = face.indx((edgeidx[2] + 1) % 3);
+    child2.verts[0] = child1.verts[0];
+    child2.verts[1] = face.verts[edgeidx[2]];
+    child2.verts[2] = face.verts[(edgeidx[2] + 1) % 3];
+
+    if (visitor(child2))
+        divide_triangle(child2, std::forward<Fn>(visitor));
+}
+
+void remove_inside_triangles(TriangleMesh &mesh, const Interior &interior,
+                             const std::vector<bool> &exclude_mask)
+{
+    enum TrPos { posInside, posTouch, posOutside };
+
+    auto &faces       = mesh.its.indices;
+    auto &vertices    = mesh.its.vertices;
+    auto bb           = mesh.bounding_box();
+
+    bool use_exclude_mask = faces.size() == exclude_mask.size();
+    auto is_excluded = [&exclude_mask, use_exclude_mask](size_t face_id) {
+        return use_exclude_mask && exclude_mask[face_id];
+    };
+
+    // TODO: Parallel mode not working yet
+    using exec_policy = ccr_seq;
+
+    // Info about the needed modifications on the input mesh.
+    struct MeshMods {
+
+        // Just a thread safe wrapper for a vector of triangles.
+        struct {
+            std::vector<std::array<Vec3f, 3>> data;
+            exec_policy::SpinningMutex        mutex;
+
+            void emplace_back(const std::array<Vec3f, 3> &pts)
+            {
+                std::lock_guard lk{mutex};
+                data.emplace_back(pts);
+            }
+
+            size_t size() const { return data.size(); }
+            const std::array<Vec3f, 3>& operator[](size_t idx) const
+            {
+                return data[idx];
+            }
+
+        } new_triangles;
+
+        // A vector of bool for all faces signaling if it needs to be removed
+        // or not.
+        std::vector<bool> to_remove;
+
+        MeshMods(const TriangleMesh &mesh):
+            to_remove(mesh.its.indices.size(), false) {}
+
+        // Number of triangles that need to be removed.
+        size_t to_remove_cnt() const
+        {
+            return std::accumulate(to_remove.begin(), to_remove.end(), size_t(0));
+        }
+
+    } mesh_mods{mesh};
+
+    // Must return true if further division of the face is needed.
+    auto divfn = [&interior, bb, &mesh_mods](const DivFace &f) {
+        BoundingBoxf3 facebb { f.verts.begin(), f.verts.end() };
+
+        // Face is certainly outside the cavity
+        if (! facebb.intersects(bb) && f.faceid != NEW_FACE) {
+            return false;
+        }
+
+        TriangleBubble bubble{facebb.center().cast<float>(), facebb.radius()};
+
+        double D = get_distance(bubble, interior);
+        double R = bubble.R * interior.voxel_scale;
+
+        if (std::isnan(D)) // The distance cannot be measured, triangle too big
+            return true;
+
+        // Distance of the bubble wall to the interior wall. Negative if the
+        // bubble is overlapping with the interior
+        double bubble_distance = D - R;
+
+        // The face is crossing the interior or inside, it must be removed and
+        // parts of it re-added, that are outside the interior
+        if (bubble_distance < 0.) {
+            if (f.faceid != NEW_FACE)
+                mesh_mods.to_remove[f.faceid] = true;
+
+            if (f.parent != NEW_FACE) // Top parent needs to be removed as well
+                mesh_mods.to_remove[f.parent] = true;
+
+            // If the outside part is between the interior end the exterior
+            // (inside the wall being invisible), no further division is needed.
+            if ((R + D) < interior.thickness)
+                return false;
+
+            return true;
+        } else if (f.faceid == NEW_FACE) {
+            // New face completely outside needs to be re-added.
+            mesh_mods.new_triangles.emplace_back(f.verts);
+        }
+
+        return false;
+    };
+
+    interior.reset_accessor();
+
+    exec_policy::for_each(size_t(0), faces.size(), [&] (size_t face_idx) {
+        const Vec3i &face = faces[face_idx];
+
+        // If the triangle is excluded, we need to keep it.
+        if (is_excluded(face_idx))
+            return;
+
+        std::array<Vec3f, 3> pts =
+            { vertices[face(0)], vertices[face(1)], vertices[face(2)] };
+
+        BoundingBoxf3 facebb { pts.begin(), pts.end() };
+
+        // Face is certainly outside the cavity
+        if (! facebb.intersects(bb)) return;
+
+        DivFace df{face, pts, long(face_idx)};
+
+        if (divfn(df))
+            divide_triangle(df, divfn);
+
+    }, exec_policy::max_concurreny());
+
+    auto new_faces = reserve_vector<Vec3i>(faces.size() +
+                                           mesh_mods.new_triangles.size());
+
+    for (size_t face_idx = 0; face_idx < faces.size(); ++face_idx) {
+        if (!mesh_mods.to_remove[face_idx])
+            new_faces.emplace_back(faces[face_idx]);
+    }
+
+    for(size_t i = 0; i < mesh_mods.new_triangles.size(); ++i) {
+        size_t o = vertices.size();
+        vertices.emplace_back(mesh_mods.new_triangles[i][0]);
+        vertices.emplace_back(mesh_mods.new_triangles[i][1]);
+        vertices.emplace_back(mesh_mods.new_triangles[i][2]);
+        new_faces.emplace_back(int(o), int(o + 1), int(o + 2));
+    }
+
+    BOOST_LOG_TRIVIAL(info)
+            << "Trimming: " << mesh_mods.to_remove_cnt() << " triangles removed";
+    BOOST_LOG_TRIVIAL(info)
+            << "Trimming: " << mesh_mods.new_triangles.size() << " triangles added";
+
+    faces.swap(new_faces);
+    new_faces = {};
+
+    mesh = TriangleMesh{mesh.its};
+    //FIXME do we want to repair the mesh? Are there duplicate vertices or flipped triangles?
+}
+
+}} // namespace Slic3r::sla

src/libslic3r/SLA/Hollowing.hpp

@@ -0,0 +1,109 @@
+#ifndef SLA_HOLLOWING_HPP
+#define SLA_HOLLOWING_HPP
+
+#include <memory>
+#include <libslic3r/TriangleMesh.hpp>
+#include <libslic3r/SLA/JobController.hpp>
+
+namespace Slic3r {
+
+namespace sla {
+
+struct HollowingConfig
+{
+    double min_thickness    = 2.;
+    double quality          = 0.5;
+    double closing_distance = 0.5;
+    bool enabled = true;
+};
+
+enum HollowingFlags { hfRemoveInsideTriangles = 0x1 };
+
+// All data related to a generated mesh interior. Includes the 3D grid and mesh
+// and various metadata. No need to manipulate from outside.
+struct Interior;
+struct InteriorDeleter { void operator()(Interior *p); };
+using  InteriorPtr = std::unique_ptr<Interior, InteriorDeleter>;
+
+indexed_triangle_set &      get_mesh(Interior &interior);
+const indexed_triangle_set &get_mesh(const Interior &interior);
+
+struct DrainHole
+{
+    Vec3f pos;
+    Vec3f normal;
+    float radius;
+    float height;
+    bool  failed = false;
+
+    DrainHole()
+        : pos(Vec3f::Zero()), normal(Vec3f::UnitZ()), radius(5.f), height(10.f)
+    {}
+
+    DrainHole(Vec3f p, Vec3f n, float r, float h, bool fl = false)
+        : pos(p), normal(n), radius(r), height(h), failed(fl)
+    {}
+
+    DrainHole(const DrainHole& rhs) :
+        DrainHole(rhs.pos, rhs.normal, rhs.radius, rhs.height, rhs.failed) {}
+    
+    bool operator==(const DrainHole &sp) const;
+    
+    bool operator!=(const DrainHole &sp) const { return !(sp == (*this)); }
+
+    bool is_inside(const Vec3f& pt) const;
+
+    bool get_intersections(const Vec3f& s, const Vec3f& dir,
+                           std::array<std::pair<float, Vec3d>, 2>& out) const;
+    
+    indexed_triangle_set to_mesh() const;
+    
+    template<class Archive> inline void serialize(Archive &ar)
+    {
+        ar(pos, normal, radius, height, failed);
+    }
+
+    static constexpr size_t steps = 32;
+};
+
+using DrainHoles = std::vector<DrainHole>;
+
+constexpr float HoleStickOutLength = 1.f;
+
+InteriorPtr generate_interior(const TriangleMesh &mesh,
+                              const HollowingConfig &  = {},
+                              const JobController &ctl = {});
+
+// Will do the hollowing
+void hollow_mesh(TriangleMesh &mesh, const HollowingConfig &cfg, int flags = 0);
+
+// Hollowing prepared in "interior", merge with original mesh
+void hollow_mesh(TriangleMesh &mesh, const Interior &interior, int flags = 0);
+
+void remove_inside_triangles(TriangleMesh &mesh, const Interior &interior,
+                             const std::vector<bool> &exclude_mask = {});
+
+double get_distance(const Vec3f &p, const Interior &interior);
+
+template<class T>
+FloatingOnly<T> get_distance(const Vec<3, T> &p, const Interior &interior)
+{
+    return get_distance(Vec3f(p.template cast<float>()), interior);
+}
+
+void cut_drainholes(std::vector<ExPolygons> & obj_slices,
+                    const std::vector<float> &slicegrid,
+                    float                     closing_radius,
+                    const sla::DrainHoles &   holes,
+                    std::function<void(void)> thr);
+
+inline void swap_normals(indexed_triangle_set &its)
+{
+    for (auto &face : its.indices)
+        std::swap(face(0), face(2));
+}
+
+} // namespace sla
+} // namespace Slic3r
+
+#endif // HOLLOWINGFILTER_H

src/libslic3r/SLA/IndexedMesh.cpp

@@ -0,0 +1,455 @@
+#include "IndexedMesh.hpp"
+#include "Concurrency.hpp"
+
+#include <libslic3r/AABBTreeIndirect.hpp>
+#include <libslic3r/TriangleMesh.hpp>
+
+#include <numeric>
+
+#ifdef SLIC3R_HOLE_RAYCASTER
+#include <libslic3r/SLA/Hollowing.hpp>
+#endif
+
+namespace Slic3r {
+
+namespace sla {
+
+class IndexedMesh::AABBImpl {
+private:
+    AABBTreeIndirect::Tree3f m_tree;
+    double                   m_triangle_ray_epsilon;
+
+public:
+    void init(const indexed_triangle_set &its, bool calculate_epsilon)
+    {
+        m_triangle_ray_epsilon = 0.000001;
+        if (calculate_epsilon) {
+            // Calculate epsilon from average triangle edge length.
+            double l = its_average_edge_length(its);
+            if (l > 0)
+                m_triangle_ray_epsilon = 0.000001 * l * l;
+        }
+        m_tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(
+            its.vertices, its.indices);
+    }
+
+    void intersect_ray(const indexed_triangle_set &its,
+                       const Vec3d &               s,
+                       const Vec3d &               dir,
+                       igl::Hit &                  hit)
+    {
+        AABBTreeIndirect::intersect_ray_first_hit(its.vertices, its.indices,
+                                                  m_tree, s, dir, hit, m_triangle_ray_epsilon);
+    }
+
+    void intersect_ray(const indexed_triangle_set &its,
+                       const Vec3d &               s,
+                       const Vec3d &               dir,
+                       std::vector<igl::Hit> &     hits)
+    {
+        AABBTreeIndirect::intersect_ray_all_hits(its.vertices, its.indices,
+                                                 m_tree, s, dir, hits, m_triangle_ray_epsilon);
+    }
+
+    double squared_distance(const indexed_triangle_set & its,
+                            const Vec3d &                point,
+                            int &                        i,
+                            Eigen::Matrix<double, 1, 3> &closest)
+    {
+        size_t idx_unsigned = 0;
+        Vec3d  closest_vec3d(closest);
+        double dist =
+            AABBTreeIndirect::squared_distance_to_indexed_triangle_set(
+                its.vertices, its.indices, m_tree, point, idx_unsigned,
+                closest_vec3d);
+        i       = int(idx_unsigned);
+        closest = closest_vec3d;
+        return dist;
+    }
+};
+
+template<class M> void IndexedMesh::init(const M &mesh, bool calculate_epsilon)
+{
+    BoundingBoxf3 bb = bounding_box(mesh);
+    m_ground_level += bb.min(Z);
+
+    // Build the AABB accelaration tree
+    m_aabb->init(*m_tm, calculate_epsilon);
+}
+
+IndexedMesh::IndexedMesh(const indexed_triangle_set& tmesh, bool calculate_epsilon)
+    : m_aabb(new AABBImpl()), m_tm(&tmesh)
+{
+    init(tmesh, calculate_epsilon);
+}
+
+IndexedMesh::IndexedMesh(const TriangleMesh &mesh, bool calculate_epsilon)
+    : m_aabb(new AABBImpl()), m_tm(&mesh.its)
+{
+    init(mesh, calculate_epsilon);
+}
+
+IndexedMesh::~IndexedMesh() {}
+
+IndexedMesh::IndexedMesh(const IndexedMesh &other):
+    m_tm(other.m_tm), m_ground_level(other.m_ground_level),
+    m_aabb( new AABBImpl(*other.m_aabb) ) {}
+
+
+IndexedMesh &IndexedMesh::operator=(const IndexedMesh &other)
+{
+    m_tm = other.m_tm;
+    m_ground_level = other.m_ground_level;
+    m_aabb.reset(new AABBImpl(*other.m_aabb)); return *this;
+}
+
+IndexedMesh &IndexedMesh::operator=(IndexedMesh &&other) = default;
+
+IndexedMesh::IndexedMesh(IndexedMesh &&other) = default;
+
+
+
+const std::vector<Vec3f>& IndexedMesh::vertices() const
+{
+    return m_tm->vertices;
+}
+
+
+
+const std::vector<Vec3i>& IndexedMesh::indices()  const
+{
+    return m_tm->indices;
+}
+
+
+
+const Vec3f& IndexedMesh::vertices(size_t idx) const
+{
+    return m_tm->vertices[idx];
+}
+
+
+
+const Vec3i& IndexedMesh::indices(size_t idx) const
+{
+    return m_tm->indices[idx];
+}
+
+
+Vec3d IndexedMesh::normal_by_face_id(int face_id) const {
+
+    return its_unnormalized_normal(*m_tm, face_id).cast<double>().normalized();
+}
+
+
+IndexedMesh::hit_result
+IndexedMesh::query_ray_hit(const Vec3d &s, const Vec3d &dir) const
+{
+    assert(is_approx(dir.norm(), 1.));
+    igl::Hit hit{-1, -1, 0.f, 0.f, 0.f};
+    hit.t = std::numeric_limits<float>::infinity();
+
+#ifdef SLIC3R_HOLE_RAYCASTER
+    if (! m_holes.empty()) {
+
+        // If there are holes, the hit_results will be made by
+        // query_ray_hits (object) and filter_hits (holes):
+        return filter_hits(query_ray_hits(s, dir));
+    }
+#endif
+
+    m_aabb->intersect_ray(*m_tm, s, dir, hit);
+    hit_result ret(*this);
+    ret.m_t = double(hit.t);
+    ret.m_dir = dir;
+    ret.m_source = s;
+    if(!std::isinf(hit.t) && !std::isnan(hit.t)) {
+        ret.m_normal = this->normal_by_face_id(hit.id);
+        ret.m_face_id = hit.id;
+    }
+
+    return ret;
+}
+
+std::vector<IndexedMesh::hit_result>
+IndexedMesh::query_ray_hits(const Vec3d &s, const Vec3d &dir) const
+{
+    std::vector<IndexedMesh::hit_result> outs;
+    std::vector<igl::Hit> hits;
+    m_aabb->intersect_ray(*m_tm, s, dir, hits);
+
+    // The sort is necessary, the hits are not always sorted.
+    std::sort(hits.begin(), hits.end(),
+              [](const igl::Hit& a, const igl::Hit& b) { return a.t < b.t; });
+
+    // Remove duplicates. They sometimes appear, for example when the ray is cast
+    // along an axis of a cube due to floating-point approximations in igl (?)
+    hits.erase(std::unique(hits.begin(), hits.end(),
+                           [](const igl::Hit& a, const igl::Hit& b)
+                           { return a.t == b.t; }),
+               hits.end());
+
+    //  Convert the igl::Hit into hit_result
+    outs.reserve(hits.size());
+    for (const igl::Hit& hit : hits) {
+        outs.emplace_back(IndexedMesh::hit_result(*this));
+        outs.back().m_t = double(hit.t);
+        outs.back().m_dir = dir;
+        outs.back().m_source = s;
+        if(!std::isinf(hit.t) && !std::isnan(hit.t)) {
+            outs.back().m_normal = this->normal_by_face_id(hit.id);
+            outs.back().m_face_id = hit.id;
+        }
+    }
+
+    return outs;
+}
+
+
+#ifdef SLIC3R_HOLE_RAYCASTER
+IndexedMesh::hit_result IndexedMesh::filter_hits(
+    const std::vector<IndexedMesh::hit_result>& object_hits) const
+{
+    assert(! m_holes.empty());
+    hit_result out(*this);
+
+    if (object_hits.empty())
+        return out;
+
+    const Vec3d& s = object_hits.front().source();
+    const Vec3d& dir = object_hits.front().direction();
+
+    // A helper struct to save an intersetion with a hole
+    struct HoleHit {
+        HoleHit(float t_p, const Vec3d& normal_p, bool entry_p) :
+            t(t_p), normal(normal_p), entry(entry_p) {}
+        float t;
+        Vec3d normal;
+        bool entry;
+    };
+    std::vector<HoleHit> hole_isects;
+    hole_isects.reserve(m_holes.size());
+
+    auto sf = s.cast<float>();
+    auto dirf = dir.cast<float>();
+
+    // Collect hits on all holes, preserve information about entry/exit
+    for (const sla::DrainHole& hole : m_holes) {
+        std::array<std::pair<float, Vec3d>, 2> isects;
+        if (hole.get_intersections(sf, dirf, isects)) {
+            // Ignore hole hits behind the source
+            if (isects[0].first > 0.f) hole_isects.emplace_back(isects[0].first, isects[0].second, true);
+            if (isects[1].first > 0.f) hole_isects.emplace_back(isects[1].first, isects[1].second, false);
+        }
+    }
+
+    // Holes can intersect each other, sort the hits by t
+    std::sort(hole_isects.begin(), hole_isects.end(),
+              [](const HoleHit& a, const HoleHit& b) { return a.t < b.t; });
+
+    // Now inspect the intersections with object and holes, in the order of
+    // increasing distance. Keep track how deep are we nested in mesh/holes and
+    // pick the correct intersection.
+    // This needs to be done twice - first to find out how deep in the structure
+    // the source is, then to pick the correct intersection.
+    int hole_nested = 0;
+    int object_nested = 0;
+    for (int dry_run=1; dry_run>=0; --dry_run) {
+        hole_nested = -hole_nested;
+        object_nested = -object_nested;
+
+        bool is_hole = false;
+        bool is_entry = false;
+        const HoleHit* next_hole_hit = hole_isects.empty() ? nullptr : &hole_isects.front();
+        const hit_result* next_mesh_hit = &object_hits.front();
+
+        while (next_hole_hit || next_mesh_hit) {
+            if (next_hole_hit && next_mesh_hit) // still have hole and obj hits
+                is_hole = (next_hole_hit->t < next_mesh_hit->m_t);
+            else
+                is_hole = next_hole_hit; // one or the other ran out
+
+            // Is this entry or exit hit?
+            is_entry = is_hole ? next_hole_hit->entry : ! next_mesh_hit->is_inside();
+
+            if (! dry_run) {
+                if (! is_hole && hole_nested == 0) {
+                    // This is a valid object hit
+                    return *next_mesh_hit;
+                }
+                if (is_hole && ! is_entry && object_nested != 0) {
+                    // This holehit is the one we seek
+                    out.m_t = next_hole_hit->t;
+                    out.m_normal = next_hole_hit->normal;
+                    out.m_source = s;
+                    out.m_dir = dir;
+                    return out;
+                }
+            }
+
+            // Increase/decrease the counter
+            (is_hole ? hole_nested : object_nested) += (is_entry ? 1 : -1);
+
+            // Advance the respective pointer
+            if (is_hole && next_hole_hit++ == &hole_isects.back())
+                next_hole_hit = nullptr;
+            if (! is_hole && next_mesh_hit++ == &object_hits.back())
+                next_mesh_hit = nullptr;
+        }
+    }
+
+    // if we got here, the ray ended up in infinity
+    return out;
+}
+#endif
+
+
+double IndexedMesh::squared_distance(const Vec3d &p, int& i, Vec3d& c) const {
+    double sqdst = 0;
+    Eigen::Matrix<double, 1, 3> pp = p;
+    Eigen::Matrix<double, 1, 3> cc;
+    sqdst = m_aabb->squared_distance(*m_tm, pp, i, cc);
+    c = cc;
+    return sqdst;
+}
+
+
+static bool point_on_edge(const Vec3d& p, const Vec3d& e1, const Vec3d& e2,
+                          double eps = 0.05)
+{
+    using Line3D = Eigen::ParametrizedLine<double, 3>;
+
+    auto line = Line3D::Through(e1, e2);
+    double d = line.distance(p);
+    return std::abs(d) < eps;
+}
+
+PointSet normals(const PointSet& points,
+                 const IndexedMesh& mesh,
+                 double eps,
+                 std::function<void()> thr, // throw on cancel
+                 const std::vector<unsigned>& pt_indices)
+{
+    if (points.rows() == 0 || mesh.vertices().empty() || mesh.indices().empty())
+        return {};
+
+    std::vector<unsigned> range = pt_indices;
+    if (range.empty()) {
+        range.resize(size_t(points.rows()), 0);
+        std::iota(range.begin(), range.end(), 0);
+    }
+
+    PointSet ret(range.size(), 3);
+
+    //    for (size_t ridx = 0; ridx < range.size(); ++ridx)
+    ccr::for_each(size_t(0), range.size(),
+        [&ret, &mesh, &points, thr, eps, &range](size_t ridx) {
+            thr();
+            unsigned el = range[ridx];
+            auto  eidx   = Eigen::Index(el);
+            int   faceid = 0;
+            Vec3d p;
+
+            mesh.squared_distance(points.row(eidx), faceid, p);
+
+            auto trindex = mesh.indices(faceid);
+
+            const Vec3d &p1 = mesh.vertices(trindex(0)).cast<double>();
+            const Vec3d &p2 = mesh.vertices(trindex(1)).cast<double>();
+            const Vec3d &p3 = mesh.vertices(trindex(2)).cast<double>();
+
+            // We should check if the point lies on an edge of the hosting
+            // triangle. If it does then all the other triangles using the
+            // same two points have to be searched and the final normal should
+            // be some kind of aggregation of the participating triangle
+            // normals. We should also consider the cases where the support
+            // point lies right on a vertex of its triangle. The procedure is
+            // the same, get the neighbor triangles and calculate an average
+            // normal.
+
+            // mark the vertex indices of the edge. ia and ib marks and edge
+            // ic will mark a single vertex.
+            int ia = -1, ib = -1, ic = -1;
+
+            if (std::abs((p - p1).norm()) < eps) {
+                ic = trindex(0);
+            } else if (std::abs((p - p2).norm()) < eps) {
+                ic = trindex(1);
+            } else if (std::abs((p - p3).norm()) < eps) {
+                ic = trindex(2);
+            } else if (point_on_edge(p, p1, p2, eps)) {
+                ia = trindex(0);
+                ib = trindex(1);
+            } else if (point_on_edge(p, p2, p3, eps)) {
+                ia = trindex(1);
+                ib = trindex(2);
+            } else if (point_on_edge(p, p1, p3, eps)) {
+                ia = trindex(0);
+                ib = trindex(2);
+            }
+
+            // vector for the neigboring triangles including the detected one.
+            std::vector<size_t> neigh;
+            if (ic >= 0) { // The point is right on a vertex of the triangle
+                for (size_t n = 0; n < mesh.indices().size(); ++n) {
+                    thr();
+                    Vec3i ni = mesh.indices(n);
+                    if ((ni(X) == ic || ni(Y) == ic || ni(Z) == ic))
+                        neigh.emplace_back(n);
+                }
+            } else if (ia >= 0 && ib >= 0) { // the point is on and edge
+                // now get all the neigboring triangles
+                for (size_t n = 0; n < mesh.indices().size(); ++n) {
+                    thr();
+                    Vec3i ni = mesh.indices(n);
+                    if ((ni(X) == ia || ni(Y) == ia || ni(Z) == ia) &&
+                        (ni(X) == ib || ni(Y) == ib || ni(Z) == ib))
+                        neigh.emplace_back(n);
+                }
+            }
+
+            // Calculate the normals for the neighboring triangles
+            std::vector<Vec3d> neighnorms;
+            neighnorms.reserve(neigh.size());
+            for (size_t &tri_id : neigh)
+                neighnorms.emplace_back(mesh.normal_by_face_id(tri_id));
+
+            // Throw out duplicates. They would cause trouble with summing. We
+            // will use std::unique which works on sorted ranges. We will sort
+            // by the coefficient-wise sum of the normals. It should force the
+            // same elements to be consecutive.
+            std::sort(neighnorms.begin(), neighnorms.end(),
+                      [](const Vec3d &v1, const Vec3d &v2) {
+                          return v1.sum() < v2.sum();
+                      });
+
+            auto lend = std::unique(neighnorms.begin(), neighnorms.end(),
+                                    [](const Vec3d &n1, const Vec3d &n2) {
+                                        // Compare normals for equivalence.
+                                        // This is controvers stuff.
+                                        auto deq = [](double a, double b) {
+                                            return std::abs(a - b) < 1e-3;
+                                        };
+                                        return deq(n1(X), n2(X)) &&
+                                               deq(n1(Y), n2(Y)) &&
+                                               deq(n1(Z), n2(Z));
+                                    });
+
+            if (!neighnorms.empty()) { // there were neighbors to count with
+                // sum up the normals and then normalize the result again.
+                // This unification seems to be enough.
+                Vec3d sumnorm(0, 0, 0);
+                sumnorm = std::accumulate(neighnorms.begin(), lend, sumnorm);
+                sumnorm.normalize();
+                ret.row(long(ridx)) = sumnorm;
+            } else { // point lies safely within its triangle
+                Eigen::Vector3d U   = p2 - p1;
+                Eigen::Vector3d V   = p3 - p1;
+                ret.row(long(ridx)) = U.cross(V).normalized();
+            }
+        });
+
+    return ret;
+}
+
+}} // namespace Slic3r::sla

src/libslic3r/SLA/IndexedMesh.hpp

@@ -0,0 +1,153 @@
+#ifndef SLA_INDEXEDMESH_H
+#define SLA_INDEXEDMESH_H
+
+#include <memory>
+#include <vector>
+
+#include <libslic3r/Point.hpp>
+
+// There is an implementation of a hole-aware raycaster that was eventually
+// not used in production version. It is now hidden under following define
+// for possible future use.
+// #define SLIC3R_HOLE_RAYCASTER
+
+#ifdef SLIC3R_HOLE_RAYCASTER
+  #include "libslic3r/SLA/Hollowing.hpp"
+#endif
+
+struct indexed_triangle_set;
+
+namespace Slic3r {
+
+class TriangleMesh;
+
+namespace sla {
+
+using PointSet = Eigen::MatrixXd;
+
+/// An index-triangle structure for libIGL functions. Also serves as an
+/// alternative (raw) input format for the SLASupportTree.
+//  Implemented in libslic3r/SLA/Common.cpp
+class IndexedMesh {
+    class AABBImpl;
+    
+    const indexed_triangle_set* m_tm;
+    double m_ground_level = 0, m_gnd_offset = 0;
+    
+    std::unique_ptr<AABBImpl> m_aabb;
+
+#ifdef SLIC3R_HOLE_RAYCASTER
+    // This holds a copy of holes in the mesh. Initialized externally
+    // by load_mesh setter.
+    std::vector<DrainHole> m_holes;
+#endif
+
+    template<class M> void init(const M &mesh, bool calculate_epsilon);
+
+public:
+    
+    // calculate_epsilon ... calculate epsilon for triangle-ray intersection from an average triangle edge length.
+    // If set to false, a default epsilon is used, which works for "reasonable" meshes.
+    explicit IndexedMesh(const indexed_triangle_set &tmesh, bool calculate_epsilon = false);
+    explicit IndexedMesh(const TriangleMesh &mesh, bool calculate_epsilon = false);
+    
+    IndexedMesh(const IndexedMesh& other);
+    IndexedMesh& operator=(const IndexedMesh&);
+    
+    IndexedMesh(IndexedMesh &&other);
+    IndexedMesh& operator=(IndexedMesh &&other);
+    
+    ~IndexedMesh();
+    
+    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; }
+    
+    const std::vector<Vec3f>& vertices() const;
+    const std::vector<Vec3i>& indices()  const;
+    const Vec3f& vertices(size_t idx) const;
+    const Vec3i& indices(size_t idx) const;
+    
+    // Result of a raycast
+    class hit_result {
+        // m_t holds a distance from m_source to the intersection.
+        double m_t = infty();
+        int m_face_id = -1;
+        const IndexedMesh *m_mesh = nullptr;
+        Vec3d m_dir;
+        Vec3d m_source;
+        Vec3d m_normal;
+        friend class IndexedMesh;
+        
+        // A valid object of this class can only be obtained from
+        // IndexedMesh::query_ray_hit method.
+        explicit inline hit_result(const IndexedMesh& em): m_mesh(&em) {}
+    public:
+        // This denotes no hit on the mesh.
+        static inline constexpr double infty() { return std::numeric_limits<double>::infinity(); }
+        
+        explicit inline hit_result(double val = infty()) : m_t(val) {}
+        
+        inline double distance() const { return m_t; }
+        inline const Vec3d& direction() const { return m_dir; }
+        inline const Vec3d& source() const { return m_source; }
+        inline Vec3d position() const { return m_source + m_dir * m_t; }
+        inline int face() const { return m_face_id; }
+        inline bool is_valid() const { return m_mesh != nullptr; }
+        inline bool is_hit() const { return m_face_id >= 0 && !std::isinf(m_t); }
+
+        inline const Vec3d& normal() const {
+            assert(is_valid());
+            return m_normal;
+        }
+
+        inline bool is_inside() const {
+            return is_hit() && normal().dot(m_dir) > 0;
+        }
+    };
+
+#ifdef SLIC3R_HOLE_RAYCASTER
+    // Inform the object about location of holes
+    // creates internal copy of the vector
+    void load_holes(const std::vector<DrainHole>& holes) {
+        m_holes = holes;
+    }
+
+    // Iterates over hits and holes and returns the true hit, possibly
+    // on the inside of a hole.
+    // This function is currently not used anywhere, it was written when the
+    // holes were subtracted on slices, that is, before we started using CGAL
+    // to actually cut the holes into the mesh.
+    hit_result filter_hits(const std::vector<IndexedMesh::hit_result>& obj_hits) const;
+#endif
+
+    // Casting a ray on the mesh, returns the distance where the hit occures.
+    hit_result query_ray_hit(const Vec3d &s, const Vec3d &dir) const;
+    
+    // Casts a ray on the mesh and returns all hits
+    std::vector<hit_result> query_ray_hits(const Vec3d &s, const Vec3d &dir) const;
+
+    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);
+    }
+
+    Vec3d normal_by_face_id(int face_id) const;
+
+    const indexed_triangle_set * get_triangle_mesh() const { return m_tm; }
+};
+
+// Calculate the normals for the selected points (from 'points' set) on the
+// mesh. This will call squared distance for each point.
+PointSet normals(const PointSet& points,
+    const IndexedMesh& convert_mesh,
+    double eps = 0.05,  // min distance from edges
+    std::function<void()> throw_on_cancel = [](){},
+    const std::vector<unsigned>& selected_points = {});
+
+}} // namespace Slic3r::sla
+
+#endif // INDEXEDMESH_H

src/libslic3r/SLA/JobController.hpp

@@ -0,0 +1,32 @@
+#ifndef SLA_JOBCONTROLLER_HPP
+#define SLA_JOBCONTROLLER_HPP
+
+#include <functional>
+#include <string>
+
+namespace Slic3r { namespace sla {
+
+/// A Control structure for the support calculation. Consists of the status
+/// indicator callback and the stop condition predicate.
+struct JobController
+{
+    using StatusFn = std::function<void(unsigned, const std::string&)>;
+    using StopCond = std::function<bool(void)>;
+    using CancelFn = std::function<void(void)>;
+    
+    // This will signal the status of the calculation to the front-end
+    StatusFn statuscb = [](unsigned, const std::string&){};
+    
+    // Returns true if the calculation should be aborted.
+    StopCond stopcondition = [](){ return false; };
+    
+    // Similar to cancel callback. This should check the stop condition and
+    // if true, throw an appropriate exception. (TriangleMeshSlicer needs this)
+    // consider it a hard abort. stopcondition is permits the algorithm to
+    // terminate itself
+    CancelFn cancelfn = [](){};
+};
+
+}} // namespace Slic3r::sla
+
+#endif // JOBCONTROLLER_HPP

src/libslic3r/SLA/Pad.cpp

@@ -0,0 +1,538 @@
+#include <libslic3r/SLA/Pad.hpp>
+#include <libslic3r/SLA/SpatIndex.hpp>
+#include <libslic3r/SLA/BoostAdapter.hpp>
+//#include <libslic3r/SLA/Contour3D.hpp>
+#include <libslic3r/TriangleMeshSlicer.hpp>
+
+#include "ConcaveHull.hpp"
+
+#include "boost/log/trivial.hpp"
+#include "ClipperUtils.hpp"
+#include "Tesselate.hpp"
+#include "MTUtils.hpp"
+
+#include "TriangulateWall.hpp"
+
+// For debugging:
+// #include <fstream>
+// #include <libnest2d/tools/benchmark.h>
+#include "SVG.hpp"
+
+#include "I18N.hpp"
+#include <boost/log/trivial.hpp>
+
+//! macro used to mark string used at localization,
+//! return same string
+#define L(s) Slic3r::I18N::translate(s)
+
+namespace Slic3r { namespace sla {
+
+namespace {
+
+indexed_triangle_set walls(
+    const Polygon &lower,
+    const Polygon &upper,
+    double         lower_z_mm,
+    double         upper_z_mm)
+{
+    indexed_triangle_set w;
+    triangulate_wall(w.vertices, w.indices, lower, upper, lower_z_mm,
+                     upper_z_mm);
+    
+    return w;
+}
+
+// Same as walls() but with identical higher and lower polygons.
+inline indexed_triangle_set straight_walls(const Polygon &plate,
+                                           double         lo_z,
+                                           double         hi_z)
+{
+    return walls(plate, plate, lo_z, hi_z);
+}
+
+// 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(Points& pts,
+                      double padding,
+                      double stride,
+                      double stick_width,
+                      double penetration)
+{
+    if(stride <= EPSILON || stick_width <= EPSILON || padding <= EPSILON)
+        return;
+
+    // SVG svg("bridgestick_plate.svg");
+    // svg.draw(poly);
+
+    // 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);
+}
+
+template<class...Args>
+ExPolygons breakstick_holes(const ExPolygons &input, Args...args)
+{
+    ExPolygons ret = input;
+    for (ExPolygon &p : ret) {
+        breakstick_holes(p.contour.points, args...);
+        for (auto &h : p.holes) breakstick_holes(h.points, args...);
+    }
+
+    return ret;
+}
+
+static inline coord_t get_waffle_offset(const PadConfig &c)
+{
+    return scaled(c.brim_size_mm + c.wing_distance());
+}
+
+static inline double get_merge_distance(const PadConfig &c)
+{
+    return 2. * (1.8 * c.wall_thickness_mm) + c.max_merge_dist_mm;
+}
+
+// Part of the pad configuration that is used for 3D geometry generation
+struct PadConfig3D {
+    double thickness, height, wing_height, slope;
+
+    explicit PadConfig3D(const PadConfig &cfg2d)
+        : thickness{cfg2d.wall_thickness_mm}
+        , height{cfg2d.full_height()}
+        , wing_height{cfg2d.wall_height_mm}
+        , slope{cfg2d.wall_slope}
+    {}
+
+    inline double bottom_offset() const
+    {
+        return (thickness + wing_height) / std::tan(slope);
+    }
+};
+
+// Outer part of the skeleton is used to generate the waffled edges of the pad.
+// Inner parts will not be waffled or offsetted. Inner parts are only used if
+// pad is generated around the object and correspond to holes and inner polygons
+// in the model blueprint.
+struct PadSkeleton { ExPolygons inner, outer; };
+
+PadSkeleton divide_blueprint(const ExPolygons &bp)
+{
+    ClipperLib::PolyTree ptree = union_pt(bp);
+
+    PadSkeleton ret;
+    ret.inner.reserve(size_t(ptree.Total()));
+    ret.outer.reserve(size_t(ptree.Total()));
+
+    for (ClipperLib::PolyTree::PolyNode *node : ptree.Childs) {
+        ExPolygon poly;
+        poly.contour.points = std::move(node->Contour);
+        for (ClipperLib::PolyTree::PolyNode *child : node->Childs) {
+            poly.holes.emplace_back(std::move(child->Contour));
+
+            traverse_pt(child->Childs, &ret.inner);
+        }
+
+        ret.outer.emplace_back(poly);
+    }
+    
+    return ret;
+}
+
+// A helper class for storing polygons and maintaining a spatial index of their
+// bounding boxes.
+class Intersector {
+    BoxIndex       m_index;
+    ExPolygons     m_polys;
+
+public:
+
+    // Add a new polygon to the index
+    void add(const ExPolygon &ep)
+    {
+        m_polys.emplace_back(ep);
+        m_index.insert(BoundingBox{ep}, unsigned(m_index.size()));
+    }
+
+    // Check an arbitrary polygon for intersection with the indexed polygons
+    bool intersects(const ExPolygon &poly)
+    {
+        // Create a suitable query bounding box.
+        auto bb = poly.contour.bounding_box();
+
+        std::vector<BoxIndexEl> qres = m_index.query(bb, BoxIndex::qtIntersects);
+
+        // Now check intersections on the actual polygons (not just the boxes)
+        bool is_overlap = false;
+        auto qit        = qres.begin();
+        while (!is_overlap && qit != qres.end())
+            is_overlap = is_overlap || poly.overlaps(m_polys[(qit++)->second]);
+
+        return is_overlap;
+    }
+};
+
+// This dummy intersector to implement the "force pad everywhere" feature
+struct DummyIntersector
+{
+    inline void add(const ExPolygon &) {}
+    inline bool intersects(const ExPolygon &) { return true; }
+};
+
+template<class _Intersector>
+class _AroundPadSkeleton : public PadSkeleton
+{
+    // A spatial index used to be able to efficiently find intersections of
+    // support polygons with the model polygons.
+    _Intersector m_intersector;
+
+public:
+    _AroundPadSkeleton(const ExPolygons &support_blueprint,
+                       const ExPolygons &model_blueprint,
+                       const PadConfig & cfg,
+                       ThrowOnCancel     thr)
+    {
+        // We need to merge the support and the model contours in a special
+        // way in which the model contours have to be substracted from the
+        // support contours. The pad has to have a hole in which the model can
+        // fit perfectly (thus the substraction -- diff_ex). Also, the pad has
+        // to be eliminated from areas where there is no need for a pad, due
+        // to missing supports.
+
+        add_supports_to_index(support_blueprint);
+
+        auto model_bp_offs =
+            offset_ex(model_blueprint,
+                      scaled<float>(cfg.embed_object.object_gap_mm),
+                      ClipperLib::jtMiter, 1);
+
+        ExPolygons fullcvh =
+            wafflized_concave_hull(support_blueprint, model_bp_offs, cfg, thr);
+
+        auto model_bp_sticks =
+            breakstick_holes(model_bp_offs, cfg.embed_object.object_gap_mm,
+                             cfg.embed_object.stick_stride_mm,
+                             cfg.embed_object.stick_width_mm,
+                             cfg.embed_object.stick_penetration_mm);
+
+        ExPolygons fullpad = diff_ex(fullcvh, model_bp_sticks);
+
+        PadSkeleton divided = divide_blueprint(fullpad);
+        
+        remove_redundant_parts(divided.outer);
+        remove_redundant_parts(divided.inner);
+
+        outer = std::move(divided.outer);
+        inner = std::move(divided.inner);
+    }
+
+private:
+
+    // Add the support blueprint to the search index to be queried later
+    void add_supports_to_index(const ExPolygons &supp_bp)
+    {
+        for (auto &ep : supp_bp) m_intersector.add(ep);
+    }
+
+    // Create the wafflized pad around all object in the scene. This pad doesnt
+    // have any holes yet.
+    ExPolygons wafflized_concave_hull(const ExPolygons &supp_bp,
+                                       const ExPolygons &model_bp,
+                                       const PadConfig  &cfg,
+                                       ThrowOnCancel     thr)
+    {
+        auto allin = reserve_vector<ExPolygon>(supp_bp.size() + model_bp.size());
+
+        for (auto &ep : supp_bp) allin.emplace_back(ep.contour);
+        for (auto &ep : model_bp) allin.emplace_back(ep.contour);
+
+        ConcaveHull cchull{allin, get_merge_distance(cfg), thr};
+        return offset_waffle_style_ex(cchull, get_waffle_offset(cfg));
+    }
+
+    // To remove parts of the pad skeleton which do not host any supports
+    void remove_redundant_parts(ExPolygons &parts)
+    {
+        auto endit = std::remove_if(parts.begin(), parts.end(),
+                                    [this](const ExPolygon &p) {
+                                        return !m_intersector.intersects(p);
+                                    });
+
+        parts.erase(endit, parts.end());
+    }
+};
+
+using AroundPadSkeleton = _AroundPadSkeleton<Intersector>;
+using BrimPadSkeleton   = _AroundPadSkeleton<DummyIntersector>;
+
+class BelowPadSkeleton : public PadSkeleton
+{
+public:
+    BelowPadSkeleton(const ExPolygons &support_blueprint,
+                     const ExPolygons &model_blueprint,
+                     const PadConfig & cfg,
+                     ThrowOnCancel     thr)
+    {
+        outer.reserve(support_blueprint.size() + model_blueprint.size());
+
+        for (auto &ep : support_blueprint) outer.emplace_back(ep.contour);
+        for (auto &ep : model_blueprint) outer.emplace_back(ep.contour);
+
+        ConcaveHull ochull{outer, get_merge_distance(cfg), thr};
+
+        outer = offset_waffle_style_ex(ochull, get_waffle_offset(cfg));
+    }
+};
+
+// Offset the contour only, leave the holes untouched
+template<class...Args>
+ExPolygon offset_contour_only(const ExPolygon &poly, coord_t delta, Args...args)
+{
+    Polygons tmp = offset(poly.contour, float(delta), args...);
+
+    if (tmp.empty()) return {};
+
+    Polygons holes = poly.holes;
+    for (auto &h : holes) h.reverse();
+
+    ExPolygons tmp2 = diff_ex(tmp, holes);
+
+    if (tmp2.empty()) return {};
+
+    return std::move(tmp2.front());
+}
+
+bool add_cavity(indexed_triangle_set &pad,
+                ExPolygon &           top_poly,
+                const PadConfig3D &   cfg,
+                ThrowOnCancel         thr)
+{
+    auto logerr = []{BOOST_LOG_TRIVIAL(error)<<"Could not create pad cavity";};
+
+    double    wing_distance = cfg.wing_height / std::tan(cfg.slope);
+    coord_t   delta_inner   = -scaled(cfg.thickness + wing_distance);
+    coord_t   delta_middle  = -scaled(cfg.thickness);
+    ExPolygon inner_base    = offset_contour_only(top_poly, delta_inner);
+    ExPolygon middle_base   = offset_contour_only(top_poly, delta_middle);
+
+    if (inner_base.empty() || middle_base.empty()) { logerr(); return false; }
+
+    ExPolygons pdiff = diff_ex(top_poly, middle_base.contour);
+
+    if (pdiff.size() != 1) { logerr(); return false; }
+
+    top_poly = pdiff.front();
+
+    double z_min = -cfg.wing_height, z_max = 0;
+    its_merge(pad, walls(inner_base.contour, middle_base.contour, z_min, z_max));
+    thr();
+    its_merge(pad, triangulate_expolygon_3d(inner_base, z_min, NORMALS_UP));
+
+    return true;
+}
+
+indexed_triangle_set create_outer_pad_geometry(const ExPolygons & skeleton,
+                                               const PadConfig3D &cfg,
+                                               ThrowOnCancel      thr)
+{
+    indexed_triangle_set ret;
+
+    for (const ExPolygon &pad_part : skeleton) {
+        ExPolygon top_poly{pad_part};
+        ExPolygon bottom_poly =
+            offset_contour_only(pad_part, -scaled(cfg.bottom_offset()));
+
+        if (bottom_poly.empty()) continue;
+        thr();
+        
+        double z_min = -cfg.height, z_max = 0;
+        its_merge(ret, walls(top_poly.contour, bottom_poly.contour, z_max, z_min));
+
+        if (cfg.wing_height > 0. && add_cavity(ret, top_poly, cfg, thr))
+            z_max = -cfg.wing_height;
+
+        for (auto &h : bottom_poly.holes)
+            its_merge(ret, straight_walls(h, z_max, z_min));
+        
+        its_merge(ret, triangulate_expolygon_3d(bottom_poly, z_min, NORMALS_DOWN));
+        its_merge(ret, triangulate_expolygon_3d(top_poly, NORMALS_UP));
+    }
+
+    return ret;
+}
+
+indexed_triangle_set create_inner_pad_geometry(const ExPolygons & skeleton,
+                                               const PadConfig3D &cfg,
+                                               ThrowOnCancel      thr)
+{
+    indexed_triangle_set ret;
+
+    double z_max = 0., z_min = -cfg.height;
+    for (const ExPolygon &pad_part : skeleton) {
+        thr();
+        its_merge(ret, straight_walls(pad_part.contour, z_max, z_min));
+
+        for (auto &h : pad_part.holes)
+            its_merge(ret, straight_walls(h, z_max, z_min));
+    
+        its_merge(ret, triangulate_expolygon_3d(pad_part, z_min, NORMALS_DOWN));
+        its_merge(ret, triangulate_expolygon_3d(pad_part, z_max, NORMALS_UP));
+    }
+
+    return ret;
+}
+
+indexed_triangle_set create_pad_geometry(const PadSkeleton &skelet,
+                                         const PadConfig &  cfg,
+                                         ThrowOnCancel      thr)
+{
+#ifndef NDEBUG
+    SVG svg("pad_skeleton.svg");
+    svg.draw(skelet.outer, "green");
+    svg.draw(skelet.inner, "blue");
+    svg.Close();
+#endif
+
+    PadConfig3D cfg3d(cfg);
+    auto pg = create_outer_pad_geometry(skelet.outer, cfg3d, thr);
+    its_merge(pg, create_inner_pad_geometry(skelet.inner, cfg3d, thr));
+
+    return pg;
+}
+
+indexed_triangle_set create_pad_geometry(const ExPolygons &supp_bp,
+                                         const ExPolygons &model_bp,
+                                         const PadConfig & cfg,
+                                         ThrowOnCancel     thr)
+{
+    PadSkeleton skelet;
+
+    if (cfg.embed_object.enabled) {
+        if (cfg.embed_object.everywhere)
+            skelet = BrimPadSkeleton(supp_bp, model_bp, cfg, thr);
+        else
+            skelet = AroundPadSkeleton(supp_bp, model_bp, cfg, thr);
+    } else
+        skelet = BelowPadSkeleton(supp_bp, model_bp, cfg, thr);
+
+    return create_pad_geometry(skelet, cfg, thr);
+}
+
+} // namespace
+
+void pad_blueprint(const indexed_triangle_set &mesh,
+                   ExPolygons &                output,
+                   const std::vector<float> &  heights,
+                   ThrowOnCancel               thrfn)
+{
+    if (mesh.empty()) return;
+
+    std::vector<ExPolygons> out = slice_mesh_ex(mesh, heights, thrfn);
+
+    size_t count = 0;
+    for(auto& o : out) count += o.size();
+
+    // Unification is expensive, a simplify also speeds up the pad generation
+    auto tmp = reserve_vector<ExPolygon>(count);
+    for(ExPolygons& o : out)
+        for(ExPolygon& e : o) {
+            auto&& exss = e.simplify(scaled<double>(0.1));
+            for(ExPolygon& ep : exss) tmp.emplace_back(std::move(ep));
+        }
+
+    ExPolygons utmp = union_ex(tmp);
+
+    for(auto& o : utmp) {
+        auto&& smp = o.simplify(scaled<double>(0.1));
+        output.insert(output.end(), smp.begin(), smp.end());
+    }
+}
+
+void pad_blueprint(const indexed_triangle_set &mesh,
+                   ExPolygons &                output,
+                   float                       h,
+                   float                       layerh,
+                   ThrowOnCancel               thrfn)
+{
+    float gnd = float(bounding_box(mesh).min(Z));
+
+    std::vector<float> slicegrid = grid(gnd, gnd + h, layerh);
+    pad_blueprint(mesh, output, slicegrid, thrfn);
+}
+
+void create_pad(const ExPolygons &    sup_blueprint,
+                const ExPolygons &    model_blueprint,
+                indexed_triangle_set &out,
+                const PadConfig &     cfg,
+                ThrowOnCancel         thr)
+{
+    auto t = create_pad_geometry(sup_blueprint, model_blueprint, cfg, thr);
+    its_merge(out, t);
+}
+
+std::string PadConfig::validate() const
+{
+    static const double constexpr MIN_BRIM_SIZE_MM = .1;
+
+    if (brim_size_mm < MIN_BRIM_SIZE_MM ||
+        bottom_offset() > brim_size_mm + wing_distance() ||
+        get_waffle_offset(*this) <= MIN_BRIM_SIZE_MM)
+        return L("Pad brim size is too small for the current configuration.");
+
+    return "";
+}
+
+}} // namespace Slic3r::sla

src/libslic3r/SLA/Pad.hpp

@@ -0,0 +1,95 @@
+#ifndef SLA_PAD_HPP
+#define SLA_PAD_HPP
+
+#include <vector>
+#include <functional>
+#include <cmath>
+#include <string>
+
+struct indexed_triangle_set;
+
+namespace Slic3r {
+
+class ExPolygon;
+class Polygon;
+using ExPolygons = std::vector<ExPolygon>;
+using Polygons = std::vector<Polygon>;
+
+namespace sla {
+
+using ThrowOnCancel = std::function<void(void)>;
+
+/// Calculate the polygon representing the silhouette.
+void pad_blueprint(
+    const indexed_triangle_set &mesh,       // input mesh
+    ExPolygons &        output,     // Output will be merged with
+    const std::vector<float> &,     // Exact Z levels to sample
+    ThrowOnCancel thrfn = [] {}); // Function that throws if cancel was requested
+
+void pad_blueprint(
+    const indexed_triangle_set &mesh,
+    ExPolygons &                output,
+    float         samplingheight = 0.1f,  // The height range to sample
+    float         layerheight    = 0.05f, // The sampling height
+    ThrowOnCancel thrfn          = [] {});
+
+struct PadConfig {
+    double wall_thickness_mm = 1.;
+    double wall_height_mm = 1.;
+    double max_merge_dist_mm = 50;
+    double wall_slope = std::atan(1.0);          // Universal constant for Pi/4
+    double brim_size_mm = 1.6;
+
+    struct EmbedObject {
+        double object_gap_mm = 1.;
+        double stick_stride_mm = 10.;
+        double stick_width_mm = 0.5;
+        double stick_penetration_mm = 0.1;
+        bool enabled = false;
+        bool everywhere = false;
+        operator bool() const { return enabled; }
+    } embed_object;
+
+    inline PadConfig() = default;
+    inline PadConfig(double thickness,
+                     double height,
+                     double mergedist,
+                     double slope)
+        : wall_thickness_mm(thickness)
+        , wall_height_mm(height)
+        , max_merge_dist_mm(mergedist)
+        , wall_slope(slope)
+    {}
+
+    inline double bottom_offset() const
+    {
+        return (wall_thickness_mm + wall_height_mm) / std::tan(wall_slope);
+    }
+
+    inline double wing_distance() const
+    {
+        return wall_height_mm / std::tan(wall_slope);
+    }
+
+    inline double full_height() const
+    {
+        return wall_height_mm + wall_thickness_mm;
+    }
+
+    /// Returns the elevation needed for compensating the pad.
+    inline double required_elevation() const { return wall_thickness_mm; }
+
+    std::string validate() const;
+};
+
+void create_pad(
+    const ExPolygons &    support_contours,
+    const ExPolygons &    model_contours,
+    indexed_triangle_set &output_mesh,
+    const PadConfig &             = PadConfig(),
+    ThrowOnCancel throw_on_cancel = [] {});
+
+} // namespace sla
+} // namespace Slic3r
+
+#endif // SLABASEPOOL_HPP

src/libslic3r/SLA/RasterBase.cpp

@@ -0,0 +1,89 @@
+#ifndef SLARASTER_CPP
+#define SLARASTER_CPP
+
+#include <functional>
+
+#include <libslic3r/SLA/RasterBase.hpp>
+#include <libslic3r/SLA/AGGRaster.hpp>
+
+// minz image write:
+#include <miniz.h>
+
+namespace Slic3r { namespace sla {
+
+const RasterBase::TMirroring RasterBase::NoMirror = {false, false};
+const RasterBase::TMirroring RasterBase::MirrorX  = {true, false};
+const RasterBase::TMirroring RasterBase::MirrorY  = {false, true};
+const RasterBase::TMirroring RasterBase::MirrorXY = {true, true};
+
+EncodedRaster PNGRasterEncoder::operator()(const void *ptr, size_t w, size_t h,
+                                           size_t      num_components)
+{
+    std::vector<uint8_t> buf;
+    size_t s = 0;
+    
+    void *rawdata = tdefl_write_image_to_png_file_in_memory(
+        ptr, int(w), int(h), int(num_components), &s);
+    
+    // On error, data() will return an empty vector. No other info can be
+    // retrieved from miniz anyway...
+    if (rawdata == nullptr) return EncodedRaster({}, "png");
+    
+    auto pptr = static_cast<std::uint8_t*>(rawdata);
+    
+    buf.reserve(s);
+    std::copy(pptr, pptr + s, std::back_inserter(buf));
+    
+    MZ_FREE(rawdata);
+    return EncodedRaster(std::move(buf), "png");
+}
+
+std::ostream &operator<<(std::ostream &stream, const EncodedRaster &bytes)
+{
+    stream.write(reinterpret_cast<const char *>(bytes.data()),
+                 std::streamsize(bytes.size()));
+    
+    return stream;
+}
+
+EncodedRaster PPMRasterEncoder::operator()(const void *ptr, size_t w, size_t h,
+                                           size_t      num_components)
+{
+    std::vector<uint8_t> buf;
+    
+    auto header = std::string("P5 ") +
+            std::to_string(w) + " " +
+            std::to_string(h) + " " + "255 ";
+    
+    auto sz = w * h * num_components;
+    size_t s = sz + header.size();
+    
+    buf.reserve(s);
+
+    auto buff = reinterpret_cast<const std::uint8_t*>(ptr);
+    std::copy(header.begin(), header.end(), std::back_inserter(buf));
+    std::copy(buff, buff+sz, std::back_inserter(buf));
+    
+    return EncodedRaster(std::move(buf), "ppm");
+}
+
+std::unique_ptr<RasterBase> create_raster_grayscale_aa(
+    const RasterBase::Resolution &res,
+    const RasterBase::PixelDim &  pxdim,
+    double                        gamma,
+    const RasterBase::Trafo &     tr)
+{
+    std::unique_ptr<RasterBase> rst;
+    
+    if (gamma > 0)
+        rst = std::make_unique<RasterGrayscaleAAGammaPower>(res, pxdim, tr, gamma);
+    else
+        rst = std::make_unique<RasterGrayscaleAA>(res, pxdim, tr, agg::gamma_threshold(.5));
+    
+    return rst;
+}
+
+} // namespace sla
+} // namespace Slic3r
+
+#endif // SLARASTER_CPP

src/libslic3r/SLA/RasterBase.hpp

@@ -0,0 +1,119 @@
+#ifndef SLA_RASTERBASE_HPP
+#define SLA_RASTERBASE_HPP
+
+#include <ostream>
+#include <memory>
+#include <vector>
+#include <array>
+#include <utility>
+#include <cstdint>
+
+#include <libslic3r/ExPolygon.hpp>
+#include <libslic3r/SLA/Concurrency.hpp>
+
+namespace Slic3r {
+
+namespace sla {
+
+// Raw byte buffer paired with its size. Suitable for compressed image data.
+class EncodedRaster {
+protected:
+    std::vector<uint8_t> m_buffer;
+    std::string m_ext;
+public:
+    EncodedRaster() = default;
+    explicit EncodedRaster(std::vector<uint8_t> &&buf, std::string ext)
+        : m_buffer(std::move(buf)), m_ext(std::move(ext))
+    {}
+    
+    size_t size() const { return m_buffer.size(); }
+    const void * data() const { return m_buffer.data(); }
+    const char * extension() const { return m_ext.c_str(); }
+};
+
+using RasterEncoder =
+    std::function<EncodedRaster(const void *ptr, size_t w, size_t h, size_t num_components)>;
+
+class RasterBase {
+public:
+    
+    enum Orientation { roLandscape, roPortrait };
+    
+    using TMirroring = std::array<bool, 2>;
+    static const TMirroring NoMirror;
+    static const TMirroring MirrorX;
+    static const TMirroring MirrorY;
+    static const TMirroring MirrorXY;
+    
+    struct Trafo {
+        bool mirror_x = false, mirror_y = false, flipXY = false;
+        coord_t center_x = 0, center_y = 0;
+        
+        // Portrait orientation will make sure the drawed polygons are rotated
+        // by 90 degrees.
+        Trafo(Orientation o = roLandscape, const TMirroring &mirror = NoMirror)
+            // XY flipping implicitly does an X mirror
+            : mirror_x(o == roPortrait ? !mirror[0] : mirror[0])
+            , mirror_y(!mirror[1]) // Makes raster origin to be top left corner
+            , flipXY(o == roPortrait)
+        {}
+        
+        TMirroring get_mirror() const { return { (roPortrait ? !mirror_x : mirror_x), mirror_y}; }
+        Orientation get_orientation() const { return flipXY ? roPortrait : roLandscape; }
+        Point get_center() const { return {center_x, center_y}; }
+    };
+    
+    /// Type that represents a resolution in pixels.
+    struct Resolution {
+        size_t width_px = 0;
+        size_t height_px = 0;
+        
+        Resolution() = default;
+        Resolution(size_t w, size_t h) : width_px(w), height_px(h) {}
+        size_t pixels() const { return width_px * height_px; }
+    };
+    
+    /// Types that represents the dimension of a pixel in millimeters.
+    struct PixelDim {
+        double w_mm = 1.;
+        double h_mm = 1.;
+        
+        PixelDim() = default;
+        PixelDim(double px_width_mm, double px_height_mm)
+            : w_mm(px_width_mm), h_mm(px_height_mm)
+        {}
+    };
+    
+    virtual ~RasterBase() = default;
+    
+    /// Draw a polygon with holes.
+    virtual void draw(const ExPolygon& poly) = 0;
+    
+    /// Get the resolution of the raster.
+    virtual Resolution resolution() const = 0;
+    virtual PixelDim   pixel_dimensions() const = 0;
+    virtual Trafo      trafo() const = 0;
+    
+    virtual EncodedRaster encode(RasterEncoder encoder) const = 0;
+};
+
+struct PNGRasterEncoder {
+    EncodedRaster operator()(const void *ptr, size_t w, size_t h, size_t num_components);
+};
+
+struct PPMRasterEncoder {
+    EncodedRaster operator()(const void *ptr, size_t w, size_t h, size_t num_components);
+};
+
+std::ostream& operator<<(std::ostream &stream, const EncodedRaster &bytes);
+
+// If gamma is zero, thresholding will be performed which disables AA.
+std::unique_ptr<RasterBase> create_raster_grayscale_aa(
+    const RasterBase::Resolution &res,
+    const RasterBase::PixelDim &  pxdim,
+    double                        gamma = 1.0,
+    const RasterBase::Trafo &     tr    = {});
+
+}} // namespace Slic3r::sla
+
+#endif // SLARASTERBASE_HPP

src/libslic3r/SLA/RasterToPolygons.cpp

@@ -0,0 +1,91 @@
+#include "RasterToPolygons.hpp"
+
+#include "AGGRaster.hpp"
+#include "libslic3r/MarchingSquares.hpp"
+#include "MTUtils.hpp"
+#include "ClipperUtils.hpp"
+
+namespace marchsq {
+
+// Specialize this struct to register a raster type for the Marching squares alg
+template<> struct _RasterTraits<Slic3r::sla::RasterGrayscaleAA> {
+    using Rst = Slic3r::sla::RasterGrayscaleAA;
+    
+    // The type of pixel cell in the raster
+    using ValueType = uint8_t;
+    
+    // Value at a given position
+    static uint8_t get(const Rst &rst, size_t row, size_t col) { return rst.read_pixel(col, row); }
+    
+    // Number of rows and cols of the raster
+    static size_t rows(const Rst &rst) { return rst.resolution().height_px; }
+    static size_t cols(const Rst &rst) { return rst.resolution().width_px; }
+};
+
+} // namespace Slic3r::marchsq
+
+namespace Slic3r { namespace sla {
+
+template<class Fn> void foreach_vertex(ExPolygon &poly, Fn &&fn)
+{
+    for (auto &p : poly.contour.points) fn(p);
+    for (auto &h : poly.holes)
+        for (auto &p : h.points) fn(p);
+}
+
+ExPolygons raster_to_polygons(const RasterGrayscaleAA &rst, Vec2i windowsize)
+{    
+    size_t rows = rst.resolution().height_px, cols = rst.resolution().width_px;
+    
+    if (rows < 2 || cols < 2) return {};
+    
+    Polygons polys;
+    long w_rows = std::max(2l, long(windowsize.y()));
+    long w_cols = std::max(2l, long(windowsize.x()));
+    
+    std::vector<marchsq::Ring> rings =
+        marchsq::execute(rst, 128, {w_rows, w_cols});
+    
+    polys.reserve(rings.size());
+    
+    auto pxd = rst.pixel_dimensions();
+    pxd.w_mm = (rst.resolution().width_px * pxd.w_mm) / (rst.resolution().width_px - 1);
+    pxd.h_mm = (rst.resolution().height_px * pxd.h_mm) / (rst.resolution().height_px - 1);
+    
+    for (const marchsq::Ring &ring : rings) {
+        Polygon poly; Points &pts = poly.points;
+        pts.reserve(ring.size());
+        
+        for (const marchsq::Coord &crd : ring)
+            pts.emplace_back(scaled(crd.c * pxd.w_mm), scaled(crd.r * pxd.h_mm));
+        
+        polys.emplace_back(poly);
+    }
+    
+    // reverse the raster transformations
+    ExPolygons unioned = union_ex(polys);
+    coord_t width = scaled(cols * pxd.h_mm), height = scaled(rows * pxd.w_mm);
+    
+    auto tr = rst.trafo();
+    for (ExPolygon &expoly : unioned) {
+        if (tr.mirror_y)
+            foreach_vertex(expoly, [height](Point &p) {p.y() = height - p.y(); });
+        
+        if (tr.mirror_x)
+            foreach_vertex(expoly, [width](Point &p) {p.x() = width - p.x(); });
+        
+        expoly.translate(-tr.center_x, -tr.center_y);
+        
+        if (tr.flipXY)
+            foreach_vertex(expoly, [](Point &p) { std::swap(p.x(), p.y()); });
+        
+        if ((tr.mirror_x + tr.mirror_y + tr.flipXY) % 2) {
+            expoly.contour.reverse();
+            for (auto &h : expoly.holes) h.reverse();
+        }
+    }
+    
+    return unioned;
+}
+
+}} // namespace Slic3r

src/libslic3r/SLA/RasterToPolygons.hpp

@@ -0,0 +1,15 @@
+#ifndef RASTERTOPOLYGONS_HPP
+#define RASTERTOPOLYGONS_HPP
+
+#include "libslic3r/ExPolygon.hpp"
+
+namespace Slic3r {
+namespace sla {
+
+class RasterGrayscaleAA;
+
+ExPolygons raster_to_polygons(const RasterGrayscaleAA &rst, Vec2i windowsize = {2, 2});
+
+}} // namespace Slic3r::sla
+
+#endif // RASTERTOPOLYGONS_HPP

src/libslic3r/SLA/ReprojectPointsOnMesh.hpp

@@ -0,0 +1,46 @@
+#ifndef REPROJECTPOINTSONMESH_HPP
+#define REPROJECTPOINTSONMESH_HPP
+
+#include "libslic3r/Point.hpp"
+#include "SupportPoint.hpp"
+#include "Hollowing.hpp"
+#include "IndexedMesh.hpp"
+#include "libslic3r/Model.hpp"
+
+#include <tbb/parallel_for.h>
+
+namespace Slic3r { namespace sla {
+
+template<class Pt> Vec3d pos(const Pt &p) { return p.pos.template cast<double>(); }
+template<class Pt> void pos(Pt &p, const Vec3d &pp) { p.pos = pp.cast<float>(); }
+
+template<class PointType>
+void reproject_support_points(const IndexedMesh &mesh, std::vector<PointType> &pts)
+{
+    tbb::parallel_for(size_t(0), pts.size(), [&mesh, &pts](size_t idx) {
+        int junk;
+        Vec3d new_pos;
+        mesh.squared_distance(pos(pts[idx]), junk, new_pos);
+        pos(pts[idx], new_pos);
+    });
+}
+
+inline void reproject_points_and_holes(ModelObject *object)
+{
+    bool has_sppoints = !object->sla_support_points.empty();
+    bool has_holes    = !object->sla_drain_holes.empty();
+
+    if (!object || (!has_holes && !has_sppoints)) return;
+
+    TriangleMesh rmsh = object->raw_mesh();
+    IndexedMesh emesh{rmsh};
+
+    if (has_sppoints)
+        reproject_support_points(emesh, object->sla_support_points);
+
+    if (has_holes)
+        reproject_support_points(emesh, object->sla_drain_holes);
+}
+
+}}
+#endif // REPROJECTPOINTSONMESH_HPP

src/libslic3r/SLA/Rotfinder.cpp

@@ -0,0 +1,476 @@
+#include <limits>
+
+#include <libslic3r/SLA/Rotfinder.hpp>
+
+#include <libslic3r/Execution/ExecutionTBB.hpp>
+#include <libslic3r/Execution/ExecutionSeq.hpp>
+
+#include <libslic3r/Optimize/BruteforceOptimizer.hpp>
+#include <libslic3r/Optimize/NLoptOptimizer.hpp>
+
+#include "libslic3r/SLAPrint.hpp"
+#include "libslic3r/PrintConfig.hpp"
+
+#include <libslic3r/Geometry.hpp>
+
+#include <thread>
+
+namespace Slic3r { namespace sla {
+
+namespace {
+
+inline const Vec3f DOWN = {0.f, 0.f, -1.f};
+constexpr double POINTS_PER_UNIT_AREA = 1.f;
+
+// Get the vertices of a triangle directly in an array of 3 points
+std::array<Vec3f, 3> get_triangle_vertices(const TriangleMesh &mesh,
+                                           size_t              faceidx)
+{
+    const auto &face = mesh.its.indices[faceidx];
+    return {mesh.its.vertices[face(0)],
+            mesh.its.vertices[face(1)],
+            mesh.its.vertices[face(2)]};
+}
+
+std::array<Vec3f, 3> get_transformed_triangle(const TriangleMesh &mesh,
+                                              const Transform3f & tr,
+                                              size_t              faceidx)
+{
+    const auto &tri = get_triangle_vertices(mesh, faceidx);
+    return {tr * tri[0], tr * tri[1], tr * tri[2]};
+}
+
+template<class T> Vec<3, T> normal(const std::array<Vec<3, T>, 3> &tri)
+{
+    Vec<3, T> U = tri[1] - tri[0];
+    Vec<3, T> V = tri[2] - tri[0];
+    return U.cross(V).normalized();
+}
+
+template<class T, class AccessFn>
+T sum_score(AccessFn &&accessfn, size_t facecount, size_t Nthreads)
+{
+    T  initv         = 0.;
+    auto   mergefn   = [](T a, T b) { return a + b; };
+    size_t grainsize = facecount / Nthreads;
+    size_t from = 0, to = facecount;
+
+    return execution::reduce(ex_tbb, from, to, initv, mergefn, accessfn, grainsize);
+}
+
+// Get area and normal of a triangle
+struct Facestats {
+    Vec3f  normal;
+    double area;
+
+    explicit Facestats(const std::array<Vec3f, 3> &triangle)
+    {
+        Vec3f U = triangle[1] - triangle[0];
+        Vec3f V = triangle[2] - triangle[0];
+        Vec3f C = U.cross(V);
+        normal = C.normalized();
+        area = 0.5 * C.norm();
+    }
+};
+
+// Try to guess the number of support points needed to support a mesh
+double get_misalginment_score(const TriangleMesh &mesh, const Transform3f &tr)
+{
+    if (mesh.its.vertices.empty()) return std::nan("");
+
+    auto accessfn = [&mesh, &tr](size_t fi) {
+        Facestats fc{get_transformed_triangle(mesh, tr, fi)};
+
+        float score = fc.area
+                      * (std::abs(fc.normal.dot(Vec3f::UnitX()))
+                         + std::abs(fc.normal.dot(Vec3f::UnitY()))
+                         + std::abs(fc.normal.dot(Vec3f::UnitZ())));
+
+        // We should score against the alignment with the reference planes
+        return scaled<int_fast64_t>(score);
+    };
+
+    size_t facecount = mesh.its.indices.size();
+    size_t Nthreads  = std::thread::hardware_concurrency();
+    double S = unscaled(sum_score<int_fast64_t>(accessfn, facecount, Nthreads));
+
+    return S / facecount;
+}
+
+// The score function for a particular face
+inline double get_supportedness_score(const Facestats &fc)
+{
+    // Simply get the angle (acos of dot product) between the face normal and
+    // the DOWN vector.
+    float cosphi = fc.normal.dot(DOWN);
+    float phi = 1.f - std::acos(cosphi) / float(PI);
+
+    // Make the huge slopes more significant than the smaller slopes
+    phi = phi * phi * phi;
+
+    // Multiply with the square root of face area of the current face,
+    // the area is less important as it grows.
+    // This makes many smaller overhangs a bigger impact.
+    return std::sqrt(fc.area) * POINTS_PER_UNIT_AREA * phi;
+}
+
+// Try to guess the number of support points needed to support a mesh
+double get_supportedness_score(const TriangleMesh &mesh, const Transform3f &tr)
+{
+    if (mesh.its.vertices.empty()) return std::nan("");
+
+    auto accessfn = [&mesh, &tr](size_t fi) {
+        Facestats fc{get_transformed_triangle(mesh, tr, fi)};
+        return scaled<int_fast64_t>(get_supportedness_score(fc));
+    };
+
+    size_t facecount = mesh.its.indices.size();
+    size_t Nthreads  = std::thread::hardware_concurrency();
+    double S = unscaled(sum_score<int_fast64_t>(accessfn, facecount, Nthreads));
+
+    return S / facecount;
+}
+
+// Find transformed mesh ground level without copy and with parallel reduce.
+float find_ground_level(const TriangleMesh &mesh,
+                         const Transform3f & tr,
+                         size_t              threads)
+{
+    size_t vsize = mesh.its.vertices.size();
+
+    auto minfn = [](float a, float b) { return std::min(a, b); };
+
+    auto accessfn = [&mesh, &tr] (size_t vi) {
+        return (tr * mesh.its.vertices[vi]).z();
+    };
+
+    auto zmin = std::numeric_limits<float>::max();
+    size_t granularity = vsize / threads;
+    return execution::reduce(ex_tbb, size_t(0), vsize, zmin, minfn, accessfn, granularity);
+}
+
+float get_supportedness_onfloor_score(const TriangleMesh &mesh,
+                                      const Transform3f & tr)
+{
+    if (mesh.its.vertices.empty()) return std::nan("");
+
+    size_t Nthreads = std::thread::hardware_concurrency();
+
+    float zmin = find_ground_level(mesh, tr, Nthreads);
+    float zlvl = zmin + 0.1f; // Set up a slight tolerance from z level
+
+    auto accessfn = [&mesh, &tr, zlvl](size_t fi) {
+        std::array<Vec3f, 3> tri = get_transformed_triangle(mesh, tr, fi);
+        Facestats fc{tri};
+
+        if (tri[0].z() <= zlvl && tri[1].z() <= zlvl && tri[2].z() <= zlvl)
+            return -2 * fc.area * POINTS_PER_UNIT_AREA;
+
+        return get_supportedness_score(fc);
+    };
+
+    size_t facecount = mesh.its.indices.size();
+    double S = unscaled(sum_score<int_fast64_t>(accessfn, facecount, Nthreads));
+
+    return S / facecount;
+}
+
+using XYRotation = std::array<double, 2>;
+
+// prepare the rotation transformation
+Transform3f to_transform3f(const XYRotation &rot)
+{
+    Transform3f rt = Transform3f::Identity();
+    rt.rotate(Eigen::AngleAxisf(float(rot[1]), Vec3f::UnitY()));
+    rt.rotate(Eigen::AngleAxisf(float(rot[0]), Vec3f::UnitX()));
+
+    return rt;
+}
+
+XYRotation from_transform3f(const Transform3f &tr)
+{
+    Vec3d rot3 = Geometry::Transformation{tr.cast<double>()}.get_rotation();
+    return {rot3.x(), rot3.y()};
+}
+
+inline bool is_on_floor(const SLAPrintObjectConfig &cfg)
+{
+    auto opt_elevation = cfg.support_object_elevation.getFloat();
+    auto opt_padaround = cfg.pad_around_object.getBool();
+
+    return opt_elevation < EPSILON || opt_padaround;
+}
+
+// collect the rotations for each face of the convex hull
+std::vector<XYRotation> get_chull_rotations(const TriangleMesh &mesh, size_t max_count)
+{
+    TriangleMesh chull = mesh.convex_hull_3d();
+    double chull2d_area = chull.convex_hull().area();
+    double area_threshold = chull2d_area / (scaled<double>(1e3) * scaled(1.));
+
+    size_t facecount = chull.its.indices.size();
+
+    struct RotArea { XYRotation rot; double area; };
+
+    auto inputs = reserve_vector<RotArea>(facecount);
+
+    auto rotcmp = [](const RotArea &r1, const RotArea &r2) {
+        double xdiff = r1.rot[X] - r2.rot[X], ydiff = r1.rot[Y] - r2.rot[Y];
+        return std::abs(xdiff) < EPSILON ? ydiff < 0. : xdiff < 0.;
+    };
+
+    auto eqcmp = [](const XYRotation &r1, const XYRotation &r2) {
+        double xdiff = r1[X] - r2[X], ydiff = r1[Y] - r2[Y];
+        return std::abs(xdiff) < EPSILON  && std::abs(ydiff) < EPSILON;
+    };
+
+    for (size_t fi = 0; fi < facecount; ++fi) {
+        Facestats fc{get_triangle_vertices(chull, fi)};
+
+        if (fc.area > area_threshold)  {
+            auto q = Eigen::Quaternionf{}.FromTwoVectors(fc.normal, DOWN);
+            XYRotation rot = from_transform3f(Transform3f::Identity() * q);
+            RotArea ra = {rot, fc.area};
+
+            auto it = std::lower_bound(inputs.begin(), inputs.end(), ra, rotcmp);
+
+            if (it == inputs.end() || !eqcmp(it->rot, rot))
+                inputs.insert(it, ra);
+        }
+    }
+
+    inputs.shrink_to_fit();
+    if (!max_count) max_count = inputs.size();
+    std::sort(inputs.begin(), inputs.end(),
+              [](const RotArea &ra, const RotArea &rb) {
+                  return ra.area > rb.area;
+              });
+
+    auto ret = reserve_vector<XYRotation>(std::min(max_count, inputs.size()));
+    for (const RotArea &ra : inputs) ret.emplace_back(ra.rot);
+
+    return ret;
+}
+
+// Find the best score from a set of function inputs. Evaluate for every point.
+template<size_t N, class Fn, class It, class StopCond>
+std::array<double, N> find_min_score(Fn &&fn, It from, It to, StopCond &&stopfn)
+{
+    std::array<double, N> ret = {};
+
+    double score = std::numeric_limits<double>::max();
+
+    size_t Nthreads = std::thread::hardware_concurrency();
+    size_t dist = std::distance(from, to);
+    std::vector<double> scores(dist, score);
+
+    execution::for_each(
+        ex_tbb, size_t(0), dist, [&stopfn, &scores, &fn, &from](size_t i) {
+            if (stopfn()) return;
+
+            scores[i] = fn(*(from + i));
+        },
+        dist / Nthreads);
+
+    auto it = std::min_element(scores.begin(), scores.end());
+
+    if (it != scores.end())
+        ret = *(from + std::distance(scores.begin(), it));
+
+    return ret;
+}
+
+} // namespace
+
+
+
+template<unsigned MAX_ITER>
+struct RotfinderBoilerplate {
+    static constexpr unsigned MAX_TRIES = MAX_ITER;
+
+    int status = 0;
+    TriangleMesh mesh;
+    unsigned max_tries;
+    const RotOptimizeParams &params;
+
+    // Assemble the mesh with the correct transformation to be used in rotation
+    // optimization.
+    static TriangleMesh get_mesh_to_rotate(const ModelObject &mo)
+    {
+        TriangleMesh mesh = mo.raw_mesh();
+
+        ModelInstance *mi = mo.instances[0];
+        auto rotation = Vec3d::Zero();
+        auto offset = Vec3d::Zero();
+        Transform3d trafo_instance =
+            Geometry::assemble_transform(offset, rotation,
+                                         mi->get_scaling_factor(),
+                                         mi->get_mirror());
+
+        mesh.transform(trafo_instance);
+
+        return mesh;
+    }
+
+    RotfinderBoilerplate(const ModelObject &mo, const RotOptimizeParams &p)
+        : mesh{get_mesh_to_rotate(mo)}
+        , params{p}
+        , max_tries(p.accuracy() * MAX_TRIES)
+    {
+
+    }
+
+    void statusfn() { params.statuscb()(++status * 100.0 / max_tries); }
+    bool stopcond() { return ! params.statuscb()(-1); }
+};
+
+Vec2d find_best_misalignment_rotation(const ModelObject &      mo,
+                                      const RotOptimizeParams &params)
+{
+    RotfinderBoilerplate<1000> bp{mo, params};
+
+    // Preparing the optimizer.
+    size_t gridsize = std::sqrt(bp.max_tries);
+    opt::Optimizer<opt::AlgBruteForce> solver(
+        opt::StopCriteria{}.max_iterations(bp.max_tries)
+                           .stop_condition([&bp] { return bp.stopcond(); }),
+        gridsize
+    );
+
+    // We are searching rotations around only two axes x, y. Thus the
+    // problem becomes a 2 dimensional optimization task.
+    // We can specify the bounds for a dimension in the following way:
+    auto bounds = opt::bounds({ {-PI, PI}, {-PI, PI} });
+
+    auto result = solver.to_max().optimize(
+        [&bp] (const XYRotation &rot)
+        {
+            bp.statusfn();
+            return get_misalginment_score(bp.mesh, to_transform3f(rot));
+        }, opt::initvals({0., 0.}), bounds);
+
+    return {result.optimum[0], result.optimum[1]};
+}
+
+Vec2d find_least_supports_rotation(const ModelObject &      mo,
+                                   const RotOptimizeParams &params)
+{
+    RotfinderBoilerplate<1000> bp{mo, params};
+
+    SLAPrintObjectConfig pocfg;
+    if (params.print_config())
+        pocfg.apply(*params.print_config(), true);
+
+    pocfg.apply(mo.config.get());
+
+    XYRotation rot;
+
+    // Different search methods have to be used depending on the model elevation
+    if (is_on_floor(pocfg)) {
+
+        std::vector<XYRotation> inputs = get_chull_rotations(bp.mesh, bp.max_tries);
+        bp.max_tries = inputs.size();
+
+        // If the model can be placed on the bed directly, we only need to
+        // check the 3D convex hull face rotations.
+
+        auto objfn = [&bp](const XYRotation &rot) {
+            bp.statusfn();
+            Transform3f tr = to_transform3f(rot);
+            return get_supportedness_onfloor_score(bp.mesh, tr);
+        };
+
+        rot = find_min_score<2>(objfn, inputs.begin(), inputs.end(), [&bp] {
+            return bp.stopcond();
+        });
+
+    } else {
+        // Preparing the optimizer.
+        size_t gridsize = std::sqrt(bp.max_tries); // 2D grid has gridsize^2 calls
+        opt::Optimizer<opt::AlgBruteForce> solver(
+            opt::StopCriteria{}.max_iterations(bp.max_tries)
+                               .stop_condition([&bp] { return bp.stopcond(); }),
+            gridsize
+        );
+
+        // We are searching rotations around only two axes x, y. Thus the
+        // problem becomes a 2 dimensional optimization task.
+        // We can specify the bounds for a dimension in the following way:
+        auto bounds = opt::bounds({ {-PI, PI}, {-PI, PI} });
+
+        auto result = solver.to_min().optimize(
+            [&bp] (const XYRotation &rot)
+            {
+                bp.statusfn();
+                return get_supportedness_score(bp.mesh, to_transform3f(rot));
+            }, opt::initvals({0., 0.}), bounds);
+
+        // Save the result
+        rot = result.optimum;
+    }
+
+    return {rot[0], rot[1]};
+}
+
+inline BoundingBoxf3 bounding_box_with_tr(const indexed_triangle_set &its,
+                                          const Transform3f &tr)
+{
+    if (its.vertices.empty())
+        return {};
+
+    Vec3f bmin = tr * its.vertices.front(), bmax = tr * its.vertices.front();
+
+    for (const Vec3f &p : its.vertices) {
+        Vec3f pp = tr * p;
+        bmin = pp.cwiseMin(bmin);
+        bmax = pp.cwiseMax(bmax);
+    }
+
+    return {bmin.cast<double>(), bmax.cast<double>()};
+}
+
+Vec2d find_min_z_height_rotation(const ModelObject &mo,
+                                 const RotOptimizeParams &params)
+{
+    RotfinderBoilerplate<1000> bp{mo, params};
+
+    TriangleMesh chull = bp.mesh.convex_hull_3d();
+    auto inputs = reserve_vector<XYRotation>(chull.its.indices.size());
+    auto rotcmp = [](const XYRotation &r1, const XYRotation &r2) {
+        double xdiff = r1[X] - r2[X], ydiff = r1[Y] - r2[Y];
+        return std::abs(xdiff) < EPSILON ? ydiff < 0. : xdiff < 0.;
+    };
+    auto eqcmp = [](const XYRotation &r1, const XYRotation &r2) {
+        double xdiff = r1[X] - r2[X], ydiff = r1[Y] - r2[Y];
+        return std::abs(xdiff) < EPSILON  && std::abs(ydiff) < EPSILON;
+    };
+
+    for (size_t fi = 0; fi < chull.its.indices.size(); ++fi) {
+        Facestats fc{get_triangle_vertices(chull, fi)};
+
+        auto q = Eigen::Quaternionf{}.FromTwoVectors(fc.normal, DOWN);
+        XYRotation rot = from_transform3f(Transform3f::Identity() * q);
+
+        auto it = std::lower_bound(inputs.begin(), inputs.end(), rot, rotcmp);
+
+        if (it == inputs.end() || !eqcmp(*it, rot))
+            inputs.insert(it, rot);
+    }
+
+    inputs.shrink_to_fit();
+    bp.max_tries = inputs.size();
+
+    auto objfn = [&bp, &chull](const XYRotation &rot) {
+        bp.statusfn();
+        Transform3f tr = to_transform3f(rot);
+        return bounding_box_with_tr(chull.its, tr).size().z();
+    };
+
+    XYRotation rot = find_min_score<2>(objfn, inputs.begin(), inputs.end(), [&bp] {
+        return bp.stopcond();
+    });
+
+    return {rot[0], rot[1]};
+}
+
+}} // namespace Slic3r::sla

src/libslic3r/SLA/Rotfinder.hpp

@@ -0,0 +1,72 @@
+#ifndef SLA_ROTFINDER_HPP
+#define SLA_ROTFINDER_HPP
+
+#include <functional>
+#include <array>
+
+#include <libslic3r/Point.hpp>
+
+namespace Slic3r {
+
+class ModelObject;
+class SLAPrintObject;
+class TriangleMesh;
+class DynamicPrintConfig;
+
+namespace sla {
+
+using RotOptimizeStatusCB = std::function<bool(int)>;
+
+class RotOptimizeParams {
+    float m_accuracy = 1.;
+    const DynamicPrintConfig *m_print_config = nullptr;
+    RotOptimizeStatusCB m_statuscb = [](int) { return true; };
+
+public:
+
+    RotOptimizeParams &accuracy(float a) { m_accuracy = a; return *this; }
+    RotOptimizeParams &print_config(const DynamicPrintConfig *c)
+    {
+        m_print_config = c;
+        return *this;
+    }
+    RotOptimizeParams &statucb(RotOptimizeStatusCB cb)
+    {
+        m_statuscb = std::move(cb);
+        return *this;
+    }
+
+    float accuracy() const { return m_accuracy; }
+    const DynamicPrintConfig * print_config() const { return m_print_config; }
+    const RotOptimizeStatusCB &statuscb() const { return m_statuscb; }
+};
+
+/**
+  * The function should find the best rotation for SLA upside down printing.
+  *
+  * @param modelobj The model object representing the 3d mesh.
+  * @param accuracy The optimization accuracy from 0.0f to 1.0f. Currently,
+  * the nlopt genetic optimizer is used and the number of iterations is
+  * accuracy * 100000. This can change in the future.
+  * @param statuscb A status indicator callback called with the int
+  * argument spanning from 0 to 100. May not reach 100 if the optimization finds
+  * an optimum before max iterations are reached. It should return a boolean
+  * signaling if the operation may continue (true) or not (false). A status
+  * value lower than 0 shall not update the status but still return a valid
+  * continuation indicator.
+  *
+  * @return Returns the rotations around each axis (x, y, z)
+  */
+Vec2d find_best_misalignment_rotation(const ModelObject &modelobj,
+                                      const RotOptimizeParams & = {});
+
+Vec2d find_least_supports_rotation(const ModelObject &modelobj,
+                                   const RotOptimizeParams & = {});
+
+Vec2d find_min_z_height_rotation(const ModelObject &mo,
+                                 const RotOptimizeParams &params = {});
+
+} // namespace sla
+} // namespace Slic3r
+
+#endif // SLAROTFINDER_HPP

src/libslic3r/SLA/SpatIndex.cpp

@@ -0,0 +1,161 @@
+#include "SpatIndex.hpp"
+
+// for concave hull merging decisions
+#include <libslic3r/SLA/BoostAdapter.hpp>
+
+#ifdef _MSC_VER
+#pragma warning(push)
+#pragma warning(disable: 4244)
+#pragma warning(disable: 4267)
+#endif
+
+#include "boost/geometry/index/rtree.hpp"
+
+#ifdef _MSC_VER
+#pragma warning(pop)
+#endif
+
+namespace Slic3r { namespace sla {
+
+/* **************************************************************************
+ * PointIndex implementation
+ * ************************************************************************** */
+
+class PointIndex::Impl {
+public:
+    using BoostIndex = boost::geometry::index::rtree< PointIndexEl,
+                                                     boost::geometry::index::rstar<16, 4> /* ? */ >;
+
+    BoostIndex m_store;
+};
+
+PointIndex::PointIndex(): m_impl(new Impl()) {}
+PointIndex::~PointIndex() {}
+
+PointIndex::PointIndex(const PointIndex &cpy): m_impl(new Impl(*cpy.m_impl)) {}
+PointIndex::PointIndex(PointIndex&& cpy): m_impl(std::move(cpy.m_impl)) {}
+
+PointIndex& PointIndex::operator=(const PointIndex &cpy)
+{
+    m_impl.reset(new Impl(*cpy.m_impl));
+    return *this;
+}
+
+PointIndex& PointIndex::operator=(PointIndex &&cpy)
+{
+    m_impl.swap(cpy.m_impl);
+    return *this;
+}
+
+void PointIndex::insert(const PointIndexEl &el)
+{
+    m_impl->m_store.insert(el);
+}
+
+bool PointIndex::remove(const PointIndexEl& el)
+{
+    return m_impl->m_store.remove(el) == 1;
+}
+
+std::vector<PointIndexEl>
+PointIndex::query(std::function<bool(const PointIndexEl &)> fn) const
+{
+    namespace bgi = boost::geometry::index;
+
+    std::vector<PointIndexEl> ret;
+    m_impl->m_store.query(bgi::satisfies(fn), std::back_inserter(ret));
+    return ret;
+}
+
+std::vector<PointIndexEl> PointIndex::nearest(const Vec3d &el, unsigned k = 1) const
+{
+    namespace bgi = boost::geometry::index;
+    std::vector<PointIndexEl> ret; ret.reserve(k);
+    m_impl->m_store.query(bgi::nearest(el, k), std::back_inserter(ret));
+    return ret;
+}
+
+size_t PointIndex::size() const
+{
+    return m_impl->m_store.size();
+}
+
+void PointIndex::foreach(std::function<void (const PointIndexEl &)> fn)
+{
+    for(auto& el : m_impl->m_store) fn(el);
+}
+
+void PointIndex::foreach(std::function<void (const PointIndexEl &)> fn) const
+{
+    for(const 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);
+}
+
+}} // namespace Slic3r::sla

src/libslic3r/SLA/SpatIndex.hpp

@@ -0,0 +1,97 @@
+#ifndef SLA_SPATINDEX_HPP
+#define SLA_SPATINDEX_HPP
+
+#include <memory>
+#include <utility>
+#include <vector>
+
+#include <Eigen/Geometry>
+
+#include <libslic3r/BoundingBox.hpp>
+
+namespace Slic3r {
+namespace sla {
+
+typedef Eigen::Matrix<double,   3, 1, Eigen::DontAlign> Vec3d;
+using PointIndexEl = std::pair<Vec3d, unsigned>;
+
+class PointIndex {
+    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:
+
+    PointIndex();
+    ~PointIndex();
+
+    PointIndex(const PointIndex&);
+    PointIndex(PointIndex&&);
+    PointIndex& operator=(const PointIndex&);
+    PointIndex& operator=(PointIndex&&);
+
+    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<PointIndexEl> query(std::function<bool(const PointIndexEl&)>) const;
+    std::vector<PointIndexEl> nearest(const Vec3d&, unsigned k) const;
+    std::vector<PointIndexEl> query(const Vec3d &v, unsigned k) const // wrapper
+    {
+        return nearest(v, k);
+    }
+
+    // For testing
+    size_t size() const;
+    bool empty() const { return size() == 0; }
+
+    void foreach(std::function<void(const PointIndexEl& el)> fn);
+    void foreach(std::function<void(const PointIndexEl& el)> fn) const;
+};
+
+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&);
+    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);
+};
+
+}
+}
+
+#endif // SPATINDEX_HPP

src/libslic3r/SLA/SupportPoint.hpp

@@ -0,0 +1,67 @@
+#ifndef SLA_SUPPORTPOINT_HPP
+#define SLA_SUPPORTPOINT_HPP
+
+#include <libslic3r/Point.hpp>
+
+namespace Slic3r { namespace sla {
+
+// An enum to keep track of where the current points on the ModelObject came from.
+enum class PointsStatus {
+    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
+{
+    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 = false)
+        : 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 = false)
+        : 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
+    {
+        float rdiff = std::abs(head_front_radius - sp.head_front_radius);
+        return (pos == sp.pos) && rdiff < float(EPSILON) &&
+               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);
+    }
+};
+
+using SupportPoints = std::vector<SupportPoint>;
+
+}} // namespace Slic3r::sla
+
+#endif // SUPPORTPOINT_HPP

src/libslic3r/SLA/SupportPointGenerator.cpp

@@ -0,0 +1,668 @@
+//#include "igl/random_points_on_mesh.h"
+//#include "igl/AABB.h"
+
+#include <tbb/parallel_for.h>
+
+#include "SupportPointGenerator.hpp"
+#include "Concurrency.hpp"
+#include "Model.hpp"
+#include "ExPolygon.hpp"
+#include "SVG.hpp"
+#include "Point.hpp"
+#include "ClipperUtils.hpp"
+#include "Tesselate.hpp"
+#include "ExPolygonCollection.hpp"
+#include "MinAreaBoundingBox.hpp"
+#include "libslic3r.h"
+
+#include <iostream>
+#include <random>
+
+namespace Slic3r {
+namespace sla {
+
+/*float SupportPointGenerator::approximate_geodesic_distance(const Vec3d& p1, const Vec3d& p2, Vec3d& n1, Vec3d& n2)
+{
+    n1.normalize();
+    n2.normalize();
+
+    Vec3d v = (p2-p1);
+    v.normalize();
+
+    float c1 = n1.dot(v);
+    float c2 = n2.dot(v);
+    float result = pow(p1(0)-p2(0), 2) + pow(p1(1)-p2(1), 2) + pow(p1(2)-p2(2), 2);
+    // Check for division by zero:
+    if(fabs(c1 - c2) > 0.0001)
+        result *= (asin(c1) - asin(c2)) / (c1 - c2);
+    return result;
+}
+
+
+float SupportPointGenerator::get_required_density(float angle) const
+{
+    // calculation would be density_0 * cos(angle). To provide one more degree of freedom, we will scale the angle
+    // to get the user-set density for 45 deg. So it ends up as density_0 * cos(K * angle).
+    float K = 4.f * float(acos(m_config.density_at_45/m_config.density_at_horizontal) / M_PI);
+    return std::max(0.f, float(m_config.density_at_horizontal * cos(K*angle)));
+}
+
+float SupportPointGenerator::distance_limit(float angle) const
+{
+    return 1./(2.4*get_required_density(angle));
+}*/
+
+SupportPointGenerator::SupportPointGenerator(
+        const sla::IndexedMesh &emesh,
+        const std::vector<ExPolygons> &slices,
+        const std::vector<float> &     heights,
+        const Config &                 config,
+        std::function<void(void)> throw_on_cancel,
+        std::function<void(int)>  statusfn)
+    : SupportPointGenerator(emesh, config, throw_on_cancel, statusfn)
+{
+    std::random_device rd;
+    m_rng.seed(rd());
+    execute(slices, heights);
+}
+
+SupportPointGenerator::SupportPointGenerator(
+        const IndexedMesh &emesh,
+        const SupportPointGenerator::Config &config,
+        std::function<void ()> throw_on_cancel, 
+        std::function<void (int)> statusfn)
+    : m_config(config)
+    , m_emesh(emesh)
+    , m_throw_on_cancel(throw_on_cancel)
+    , m_statusfn(statusfn)
+{
+}
+
+void SupportPointGenerator::execute(const std::vector<ExPolygons> &slices,
+                                    const std::vector<float> &     heights)
+{
+    process(slices, heights);
+    project_onto_mesh(m_output);
+}
+
+void SupportPointGenerator::project_onto_mesh(std::vector<sla::SupportPoint>& points) const
+{
+    // The function  makes sure that all the points are really exactly placed on the mesh.
+
+    // Use a reasonable granularity to account for the worker thread synchronization cost.
+    static constexpr size_t gransize = 64;
+
+    ccr_par::for_each(size_t(0), points.size(), [this, &points](size_t idx)
+    {
+        if ((idx % 16) == 0)
+            // Don't call the following function too often as it flushes CPU write caches due to synchronization primitves.
+            m_throw_on_cancel();
+
+        Vec3f& p = points[idx].pos;
+        // Project the point upward and downward and choose the closer intersection with the mesh.
+        sla::IndexedMesh::hit_result hit_up   = m_emesh.query_ray_hit(p.cast<double>(), Vec3d(0., 0., 1.));
+        sla::IndexedMesh::hit_result hit_down = m_emesh.query_ray_hit(p.cast<double>(), Vec3d(0., 0., -1.));
+
+        bool up   = hit_up.is_hit();
+        bool down = hit_down.is_hit();
+
+        if (!up && !down)
+            return;
+
+        sla::IndexedMesh::hit_result& hit = (!down || (hit_up.distance() < hit_down.distance())) ? hit_up : hit_down;
+        p = p + (hit.distance() * hit.direction()).cast<float>();
+    }, gransize);
+}
+
+static std::vector<SupportPointGenerator::MyLayer> make_layers(
+    const std::vector<ExPolygons>& slices, const std::vector<float>& heights,
+    std::function<void(void)> throw_on_cancel)
+{
+    assert(slices.size() == heights.size());
+
+    // Allocate empty layers.
+    std::vector<SupportPointGenerator::MyLayer> layers;
+    layers.reserve(slices.size());
+    for (size_t i = 0; i < slices.size(); ++ i)
+        layers.emplace_back(i, heights[i]);
+
+    // FIXME: calculate actual pixel area from printer config:
+    //const float pixel_area = pow(wxGetApp().preset_bundle->project_config.option<ConfigOptionFloat>("display_width") / wxGetApp().preset_bundle->project_config.option<ConfigOptionInt>("display_pixels_x"), 2.f); //
+    const float pixel_area = pow(0.047f, 2.f);
+
+    ccr_par::for_each(size_t(0), layers.size(),
+        [&layers, &slices, &heights, pixel_area, throw_on_cancel](size_t layer_id)
+    {
+        if ((layer_id % 8) == 0)
+            // Don't call the following function too often as it flushes
+            // CPU write caches due to synchronization primitves.
+            throw_on_cancel();
+
+        SupportPointGenerator::MyLayer &layer   = layers[layer_id];
+        const ExPolygons &              islands = slices[layer_id];
+        // FIXME WTF?
+        const float height = (layer_id > 2 ?
+                                  heights[layer_id - 3] :
+                                  heights[0] - (heights[1] - heights[0]));
+        layer.islands.reserve(islands.size());
+        for (const ExPolygon &island : islands) {
+            float area = float(island.area() * SCALING_FACTOR * SCALING_FACTOR);
+            if (area >= pixel_area)
+                // FIXME this is not a correct centroid of a polygon with holes.
+                layer.islands.emplace_back(layer, island, get_extents(island.contour),
+                                           unscaled<float>(island.contour.centroid()), area, height);
+        }
+    }, 32 /*gransize*/);
+
+    // Calculate overlap of successive layers. Link overlapping islands.
+    ccr_par::for_each(size_t(1), layers.size(),
+                      [&layers, &heights, throw_on_cancel] (size_t layer_id)
+    {
+      if ((layer_id % 2) == 0)
+          // Don't call the following function too often as it flushes CPU write caches due to synchronization primitves.
+          throw_on_cancel();
+      SupportPointGenerator::MyLayer &layer_above = layers[layer_id];
+      SupportPointGenerator::MyLayer &layer_below = layers[layer_id - 1];
+      //FIXME WTF?
+      const float layer_height = (layer_id!=0 ? heights[layer_id]-heights[layer_id-1] : heights[0]);
+      const float safe_angle = 35.f * (float(M_PI)/180.f); // smaller number - less supports
+      const float between_layers_offset = scaled<float>(layer_height * std::tan(safe_angle));
+      const float slope_angle = 75.f * (float(M_PI)/180.f); // smaller number - less supports
+      const float slope_offset = scaled<float>(layer_height * std::tan(slope_angle));
+      //FIXME This has a quadratic time complexity, it will be excessively slow for many tiny islands.
+      for (SupportPointGenerator::Structure &top : layer_above.islands) {
+          for (SupportPointGenerator::Structure &bottom : layer_below.islands) {
+              float overlap_area = top.overlap_area(bottom);
+              if (overlap_area > 0) {
+                  top.islands_below.emplace_back(&bottom, overlap_area);
+                  bottom.islands_above.emplace_back(&top, overlap_area);
+              }
+          }
+          if (! top.islands_below.empty()) {
+              Polygons bottom_polygons = top.polygons_below();
+              top.overhangs = diff_ex(*top.polygon, bottom_polygons);
+              if (! top.overhangs.empty()) {
+
+                  // Produce 2 bands around the island, a safe band for dangling overhangs
+                  // and an unsafe band for sloped overhangs.
+                  // These masks include the original island
+                  auto dangl_mask = expand(bottom_polygons, between_layers_offset, ClipperLib::jtSquare);
+                  auto overh_mask = expand(bottom_polygons, slope_offset, ClipperLib::jtSquare);
+
+                  // Absolutely hopeless overhangs are those outside the unsafe band
+                  top.overhangs = diff_ex(*top.polygon, overh_mask);
+
+                  // Now cut out the supported core from the safe band
+                  // and cut the safe band from the unsafe band to get distinct
+                  // zones.
+                  overh_mask = diff(overh_mask, dangl_mask);
+                  dangl_mask = diff(dangl_mask, bottom_polygons);
+
+                  top.dangling_areas = intersection_ex(*top.polygon, dangl_mask);
+                  top.overhangs_slopes = intersection_ex(*top.polygon, overh_mask);
+
+                  top.overhangs_area = 0.f;
+                  std::vector<std::pair<ExPolygon*, float>> expolys_with_areas;
+                  for (ExPolygon &ex : top.overhangs) {
+                      float area = float(ex.area());
+                      expolys_with_areas.emplace_back(&ex, area);
+                      top.overhangs_area += area;
+                  }
+                  std::sort(expolys_with_areas.begin(), expolys_with_areas.end(),
+                            [](const std::pair<ExPolygon*, float> &p1, const std::pair<ExPolygon*, float> &p2)
+                            { return p1.second > p2.second; });
+                  ExPolygons overhangs_sorted;
+                  for (auto &p : expolys_with_areas)
+                      overhangs_sorted.emplace_back(std::move(*p.first));
+                  top.overhangs = std::move(overhangs_sorted);
+                  top.overhangs_area *= float(SCALING_FACTOR * SCALING_FACTOR);
+              }
+          }
+      }
+    }, 8 /* gransize */);
+
+    return layers;
+}
+
+void SupportPointGenerator::process(const std::vector<ExPolygons>& slices, const std::vector<float>& heights)
+{
+#ifdef SLA_SUPPORTPOINTGEN_DEBUG
+    std::vector<std::pair<ExPolygon, coord_t>> islands;
+#endif /* SLA_SUPPORTPOINTGEN_DEBUG */
+
+    std::vector<SupportPointGenerator::MyLayer> layers = make_layers(slices, heights, m_throw_on_cancel);
+
+    PointGrid3D point_grid;
+    point_grid.cell_size = Vec3f(10.f, 10.f, 10.f);
+
+    double increment = 100.0 / layers.size();
+    double status    = 0;
+
+    for (unsigned int layer_id = 0; layer_id < layers.size(); ++ layer_id) {
+        SupportPointGenerator::MyLayer *layer_top     = &layers[layer_id];
+        SupportPointGenerator::MyLayer *layer_bottom  = (layer_id > 0) ? &layers[layer_id - 1] : nullptr;
+        std::vector<float>        support_force_bottom;
+        if (layer_bottom != nullptr) {
+            support_force_bottom.assign(layer_bottom->islands.size(), 0.f);
+            for (size_t i = 0; i < layer_bottom->islands.size(); ++ i)
+                support_force_bottom[i] = layer_bottom->islands[i].supports_force_total();
+        }
+        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();
+                // 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()];
+//FIXME this condition does not reflect a bifurcation into a one large island and one tiny island well, it incorrectly resets the support force to zero.
+// One should rather work with the overlap area vs overhang area.
+//                support_force *= std::min(1.f, 1.f - std::min(1.f, 0.1f * centroids_dist * centroids_dist / bottom.area));
+                // Penalization resulting from increasing polygon area:
+                support_force *= std::min(1.f, 20.f * bottom.area / top.area);
+            }
+        // Let's assign proper support force to each of them:
+        if (layer_id > 0) {
+            for (Structure &below : layer_bottom->islands) {
+                float below_support_force = support_force_bottom[&below - layer_bottom->islands.data()];
+                float above_overlap_area = 0.f;
+                for (Structure::Link &above_link : below.islands_above)
+                    above_overlap_area += above_link.overlap_area;
+                for (Structure::Link &above_link : below.islands_above)
+                    above_link.island->supports_force_inherited += below_support_force * above_link.overlap_area / above_overlap_area;
+            }
+        }
+        // Now iterate over all polygons and append new points if needed.
+        for (Structure &s : layer_top->islands) {
+            // Penalization resulting from large diff from the last layer:
+            s.supports_force_inherited /= std::max(1.f, 0.17f * (s.overhangs_area) / s.area);
+
+            add_support_points(s, point_grid);
+        }
+
+        m_throw_on_cancel();
+
+        status += increment;
+        m_statusfn(int(std::round(status)));
+
+#ifdef SLA_SUPPORTPOINTGEN_DEBUG
+        /*std::string layer_num_str = std::string((i<10 ? "0" : "")) + std::string((i<100 ? "0" : "")) + std::to_string(i);
+        output_expolygons(expolys_top, "top" + layer_num_str + ".svg");
+        output_expolygons(diff, "diff" + layer_num_str + ".svg");
+        if (!islands.empty())
+            output_expolygons(islands, "islands" + layer_num_str + ".svg");*/
+#endif /* SLA_SUPPORTPOINTGEN_DEBUG */
+    }
+}
+
+void SupportPointGenerator::add_support_points(SupportPointGenerator::Structure &s, SupportPointGenerator::PointGrid3D &grid3d)
+{
+    // Select each type of surface (overrhang, dangling, slope), derive the support
+    // force deficit for it and call uniformly conver with the right params
+
+    float tp      = m_config.tear_pressure();
+    float current = s.supports_force_total();
+
+    if (s.islands_below.empty()) {
+        // completely new island - needs support no doubt
+        // deficit is full, there is nothing below that would hold this island
+        uniformly_cover({ *s.polygon }, s, s.area * tp, grid3d, IslandCoverageFlags(icfIsNew | icfWithBoundary) );
+        return;
+    }
+
+    if (! s.overhangs.empty()) {
+        uniformly_cover(s.overhangs, s, s.overhangs_area * tp, grid3d);
+    }
+
+    auto areafn = [](double sum, auto &p) { return sum + p.area() * SCALING_FACTOR * SCALING_FACTOR; };
+
+    current = s.supports_force_total();
+    if (! s.dangling_areas.empty()) {
+        // Let's see if there's anything that overlaps enough to need supports:
+        // What we now have in polygons needs support, regardless of what the forces are, so we can add them.
+
+        double a = std::accumulate(s.dangling_areas.begin(), s.dangling_areas.end(), 0., areafn);
+        uniformly_cover(s.dangling_areas, s, a * tp - a * current * s.area, grid3d, icfWithBoundary);
+    }
+
+    current = s.supports_force_total();
+    if (! s.overhangs_slopes.empty()) {
+        double a = std::accumulate(s.overhangs_slopes.begin(), s.overhangs_slopes.end(), 0., areafn);
+        uniformly_cover(s.overhangs_slopes, s, a * tp - a * current / s.area, grid3d, icfWithBoundary);
+    }
+}
+
+std::vector<Vec2f> sample_expolygon(const ExPolygon &expoly, float samples_per_mm2, std::mt19937 &rng)
+{
+    // Triangulate the polygon with holes into triplets of 3D points.
+    std::vector<Vec2f> triangles = Slic3r::triangulate_expolygon_2f(expoly);
+
+    std::vector<Vec2f> out;
+    if (! triangles.empty())
+    {
+        // Calculate area of each triangle.
+        auto   areas = reserve_vector<float>(triangles.size() / 3);
+        double aback = 0.;
+        for (size_t i = 0; i < triangles.size(); ) {
+            const Vec2f &a  = triangles[i ++];
+            const Vec2f  v1 = triangles[i ++] - a;
+            const Vec2f  v2 = triangles[i ++] - a;
+
+            // Prefix sum of the areas.
+            areas.emplace_back(aback + 0.5f * std::abs(cross2(v1, v2)));
+            aback = areas.back();
+        }
+
+        size_t num_samples = size_t(ceil(areas.back() * samples_per_mm2));
+        std::uniform_real_distribution<> random_triangle(0., double(areas.back()));
+        std::uniform_real_distribution<> random_float(0., 1.);
+        for (size_t i = 0; i < num_samples; ++ i) {
+            double r = random_triangle(rng);
+            size_t idx_triangle = std::min<size_t>(std::upper_bound(areas.begin(), areas.end(), (float)r) - areas.begin(), areas.size() - 1) * 3;
+            // Select a random point on the triangle.
+            const Vec2f &a = triangles[idx_triangle ++];
+            const Vec2f &b = triangles[idx_triangle++];
+            const Vec2f &c = triangles[idx_triangle];
+#if 1
+            // https://www.cs.princeton.edu/~funk/tog02.pdf
+            // page 814, formula 1.
+            double u = float(std::sqrt(random_float(rng)));
+            double v = float(random_float(rng));
+            out.emplace_back(a * (1.f - u) + b * (u * (1.f - v)) + c * (v * u));
+#else
+            // Greg Turk, Graphics Gems
+            // https://devsplorer.wordpress.com/2019/08/07/find-a-random-point-on-a-plane-using-barycentric-coordinates-in-unity/
+            double u = float(random_float(rng));
+            double v = float(random_float(rng));
+            if (u + v >= 1.f) {
+              u = 1.f - u;
+              v = 1.f - v;
+            }
+            out.emplace_back(a + u * (b - a) + v * (c - a));
+#endif
+        }
+    }
+    return out;
+}
+
+
+std::vector<Vec2f> sample_expolygon(const ExPolygons &expolys, float samples_per_mm2, std::mt19937 &rng)
+{
+    std::vector<Vec2f> out;
+    for (const ExPolygon &expoly : expolys)
+        append(out, sample_expolygon(expoly, samples_per_mm2, rng));
+
+    return out;
+}
+
+void sample_expolygon_boundary(const ExPolygon &   expoly,
+                               float               samples_per_mm,
+                               std::vector<Vec2f> &out,
+                               std::mt19937 &      /*rng*/)
+{
+    double  point_stepping_scaled = scale_(1.f) / samples_per_mm;
+    for (size_t i_contour = 0; i_contour <= expoly.holes.size(); ++ i_contour) {
+        const Polygon &contour = (i_contour == 0) ? expoly.contour :
+                                                    expoly.holes[i_contour - 1];
+
+        const Points pts = contour.equally_spaced_points(point_stepping_scaled);
+        for (size_t i = 0; i < pts.size(); ++ i)
+            out.emplace_back(unscale<float>(pts[i].x()),
+                             unscale<float>(pts[i].y()));
+    }
+}
+
+std::vector<Vec2f> sample_expolygon_with_boundary(const ExPolygon &expoly, float samples_per_mm2, float samples_per_mm_boundary, std::mt19937 &rng)
+{
+    std::vector<Vec2f> out = sample_expolygon(expoly, samples_per_mm2, rng);
+    sample_expolygon_boundary(expoly, samples_per_mm_boundary, out, rng);
+    return out;
+}
+
+std::vector<Vec2f> sample_expolygon_with_boundary(const ExPolygons &expolys, float samples_per_mm2, float samples_per_mm_boundary, std::mt19937 &rng)
+{
+    std::vector<Vec2f> out;
+    for (const ExPolygon &expoly : expolys)
+        append(out, sample_expolygon_with_boundary(expoly, samples_per_mm2, samples_per_mm_boundary, rng));
+    return out;
+}
+
+template<typename REFUSE_FUNCTION>
+static inline std::vector<Vec2f> poisson_disk_from_samples(const std::vector<Vec2f> &raw_samples, float radius, REFUSE_FUNCTION refuse_function)
+{
+    Vec2f corner_min(std::numeric_limits<float>::max(), std::numeric_limits<float>::max());
+    for (const Vec2f &pt : raw_samples) {
+        corner_min.x() = std::min(corner_min.x(), pt.x());
+        corner_min.y() = std::min(corner_min.y(), pt.y());
+    }
+
+    // Assign the raw samples to grid cells, sort the grid cells lexicographically.
+    struct RawSample
+    {
+        Vec2f coord;
+        Vec2i cell_id;
+        RawSample(const Vec2f &crd = {}, const Vec2i &id = {}): coord{crd}, cell_id{id} {}
+    };
+
+    auto raw_samples_sorted = reserve_vector<RawSample>(raw_samples.size());
+    for (const Vec2f &pt : raw_samples)
+        raw_samples_sorted.emplace_back(pt, ((pt - corner_min) / radius).cast<int>());
+
+    std::sort(raw_samples_sorted.begin(), raw_samples_sorted.end(), [](const RawSample &lhs, const RawSample &rhs)
+        { return lhs.cell_id.x() < rhs.cell_id.x() || (lhs.cell_id.x() == rhs.cell_id.x() && lhs.cell_id.y() < rhs.cell_id.y()); });
+
+    struct PoissonDiskGridEntry {
+        // Resulting output sample points for this cell:
+        enum {
+            max_positions = 4
+        };
+        Vec2f   poisson_samples[max_positions];
+        int     num_poisson_samples = 0;
+
+        // Index into raw_samples:
+        int     first_sample_idx;
+        int     sample_cnt;
+    };
+
+    struct CellIDHash {
+        std::size_t operator()(const Vec2i &cell_id) const {
+            return std::hash<int>()(cell_id.x()) ^ std::hash<int>()(cell_id.y() * 593);
+        }
+    };
+
+    // Map from cell IDs to hash_data.  Each hash_data points to the range in raw_samples corresponding to that cell.
+    // (We could just store the samples in hash_data.  This implementation is an artifact of the reference paper, which
+    // is optimizing for GPU acceleration that we haven't implemented currently.)
+    typedef std::unordered_map<Vec2i, PoissonDiskGridEntry, CellIDHash> Cells;
+    Cells cells;
+    {
+        typename Cells::iterator last_cell_id_it;
+        Vec2i           last_cell_id(-1, -1);
+        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
+                // always contiguous, since we've sorted raw_samples_sorted by cell ID.
+                ++ last_cell_id_it->second.sample_cnt;
+            } else {
+                // This is a new cell.
+                PoissonDiskGridEntry data;
+                data.first_sample_idx = int(i);
+                data.sample_cnt       = 1;
+                auto result     = cells.insert({sample.cell_id, data});
+                last_cell_id    = sample.cell_id;
+                last_cell_id_it = result.first;
+            }
+        }
+    }
+
+    const int   max_trials = 5;
+    const float radius_squared = radius * radius;
+    for (int trial = 0; trial < max_trials; ++ trial) {
+        // Create sample points for each entry in cells.
+        for (auto &it : cells) {
+            const Vec2i          &cell_id   = it.first;
+            PoissonDiskGridEntry &cell_data = it.second;
+            // This cell's raw sample points start at first_sample_idx.  On trial 0, try the first one. On trial 1, try first_sample_idx + 1.
+            int next_sample_idx = cell_data.first_sample_idx + trial;
+            if (trial >= cell_data.sample_cnt)
+                // There are no more points to try for this cell.
+                continue;
+            const RawSample &candidate = raw_samples_sorted[next_sample_idx];
+            // See if this point conflicts with any other points in this cell, or with any points in
+            // neighboring cells.  Note that it's possible to have more than one point in the same cell.
+            bool conflict = refuse_function(candidate.coord);
+            for (int i = -1; i < 2 && ! conflict; ++ i) {
+                for (int j = -1; j < 2; ++ j) {
+                    const auto &it_neighbor = cells.find(cell_id + Vec2i(i, j));
+                    if (it_neighbor != cells.end()) {
+                        const PoissonDiskGridEntry &neighbor = it_neighbor->second;
+                        for (int i_sample = 0; i_sample < neighbor.num_poisson_samples; ++ i_sample)
+                            if ((neighbor.poisson_samples[i_sample] - candidate.coord).squaredNorm() < radius_squared) {
+                                conflict = true;
+                                break;
+                            }
+                    }
+                }
+            }
+            if (! conflict) {
+                // Store the new sample.
+                assert(cell_data.num_poisson_samples < cell_data.max_positions);
+                if (cell_data.num_poisson_samples < cell_data.max_positions)
+                    cell_data.poisson_samples[cell_data.num_poisson_samples ++] = candidate.coord;
+            }
+        }
+    }
+
+    // Copy the results to the output.
+    std::vector<Vec2f> out;
+    for (const auto& it : cells)
+        for (int i = 0; i < it.second.num_poisson_samples; ++ i)
+            out.emplace_back(it.second.poisson_samples[i]);
+    return out;
+}
+
+
+void SupportPointGenerator::uniformly_cover(const ExPolygons& islands, Structure& structure, float deficit, PointGrid3D &grid3d, IslandCoverageFlags flags)
+{
+    //int num_of_points = std::max(1, (int)((island.area()*pow(SCALING_FACTOR, 2) * m_config.tear_pressure)/m_config.support_force));
+
+    float support_force_deficit = deficit;
+//    auto bb = get_extents(islands);
+
+    if (flags & icfIsNew) {
+        auto chull = ExPolygonCollection{islands}.convex_hull();
+        auto rotbox = MinAreaBoundigBox{chull, MinAreaBoundigBox::pcConvex};
+        Vec2d bbdim = {unscaled(rotbox.width()), unscaled(rotbox.height())};
+
+        if (bbdim.x() > bbdim.y()) std::swap(bbdim.x(), bbdim.y());
+        double aspectr = bbdim.y() / bbdim.x();
+
+        support_force_deficit *= (1 + aspectr / 2.);
+    }
+
+    if (support_force_deficit < 0)
+        return;
+
+    // Number of newly added points.
+    const size_t poisson_samples_target = size_t(ceil(support_force_deficit / m_config.support_force()));
+
+    const float density_horizontal = m_config.tear_pressure() / m_config.support_force();
+    //FIXME why?
+    float poisson_radius		= std::max(m_config.minimal_distance, 1.f / (5.f * density_horizontal));
+//    const float poisson_radius     = 1.f / (15.f * density_horizontal);
+    const float samples_per_mm2 = 30.f / (float(M_PI) * poisson_radius * poisson_radius);
+    // Minimum distance between samples, in 3D space.
+//    float min_spacing			= poisson_radius / 3.f;
+    float min_spacing			= poisson_radius;
+
+    //FIXME share the random generator. The random generator may be not so cheap to initialize, also we don't want the random generator to be restarted for each polygon.
+
+    std::vector<Vec2f> raw_samples =
+        flags & icfWithBoundary ?
+            sample_expolygon_with_boundary(islands, samples_per_mm2,
+                                           5.f / poisson_radius, m_rng) :
+            sample_expolygon(islands, samples_per_mm2, m_rng);
+
+    std::vector<Vec2f>  poisson_samples;
+    for (size_t iter = 0; iter < 4; ++ iter) {
+        poisson_samples = poisson_disk_from_samples(raw_samples, poisson_radius,
+            [&structure, &grid3d, min_spacing](const Vec2f &pos) {
+                return grid3d.collides_with(pos, structure.layer->print_z, min_spacing);
+            });
+        if (poisson_samples.size() >= poisson_samples_target || m_config.minimal_distance > poisson_radius-EPSILON)
+            break;
+        float coeff = 0.5f;
+        if (poisson_samples.size() * 2 > poisson_samples_target)
+            coeff = float(poisson_samples.size()) / float(poisson_samples_target);
+        poisson_radius = std::max(m_config.minimal_distance, poisson_radius * coeff);
+        min_spacing    = std::max(m_config.minimal_distance, min_spacing * coeff);
+    }
+
+#ifdef SLA_SUPPORTPOINTGEN_DEBUG
+    {
+        static int irun = 0;
+        Slic3r::SVG svg(debug_out_path("SLA_supports-uniformly_cover-%d.svg", irun ++), get_extents(islands));
+        for (const ExPolygon &island : islands)
+            svg.draw(island);
+        for (const Vec2f &pt : raw_samples)
+            svg.draw(Point(scale_(pt.x()), scale_(pt.y())), "red");
+        for (const Vec2f &pt : poisson_samples)
+            svg.draw(Point(scale_(pt.x()), scale_(pt.y())), "blue");
+    }
+#endif /* NDEBUG */
+
+//    assert(! poisson_samples.empty());
+    if (poisson_samples_target < poisson_samples.size()) {
+        std::shuffle(poisson_samples.begin(), poisson_samples.end(), m_rng);
+        poisson_samples.erase(poisson_samples.begin() + poisson_samples_target, poisson_samples.end());
+    }
+    for (const Vec2f &pt : poisson_samples) {
+        m_output.emplace_back(float(pt(0)), float(pt(1)), structure.zlevel, m_config.head_diameter/2.f, flags & icfIsNew);
+        structure.supports_force_this_layer += m_config.support_force();
+        grid3d.insert(pt, &structure);
+    }
+}
+
+
+void remove_bottom_points(std::vector<SupportPoint> &pts, float lvl)
+{
+    // get iterator to the reorganized vector end
+    auto endit = std::remove_if(pts.begin(), pts.end(), [lvl]
+                                (const sla::SupportPoint &sp) {
+        return sp.pos.z() <= lvl;
+    });
+
+    // erase all elements after the new end
+    pts.erase(endit, pts.end());
+}
+
+#ifdef SLA_SUPPORTPOINTGEN_DEBUG
+void SupportPointGenerator::output_structures(const std::vector<Structure>& structures)
+{
+    for (unsigned int i=0 ; i<structures.size(); ++i) {
+        std::stringstream ss;
+        ss << structures[i].unique_id.count() << "_" << std::setw(10) << std::setfill('0') << 1000 + (int)structures[i].height/1000 << ".png";
+        output_expolygons(std::vector<ExPolygon>{*structures[i].polygon}, ss.str());
+    }
+}
+
+void SupportPointGenerator::output_expolygons(const ExPolygons& expolys, const std::string &filename)
+{
+    BoundingBox bb(Point(-30000000, -30000000), Point(30000000, 30000000));
+    Slic3r::SVG svg_cummulative(filename, bb);
+    for (size_t i = 0; i < expolys.size(); ++ i) {
+        /*Slic3r::SVG svg("single"+std::to_string(i)+".svg", bb);
+        svg.draw(expolys[i]);
+        svg.draw_outline(expolys[i].contour, "black", scale_(0.05));
+        svg.draw_outline(expolys[i].holes, "blue", scale_(0.05));
+        svg.Close();*/
+
+        svg_cummulative.draw(expolys[i]);
+        svg_cummulative.draw_outline(expolys[i].contour, "black", scale_(0.05));
+        svg_cummulative.draw_outline(expolys[i].holes, "blue", scale_(0.05));
+    }
+}
+#endif
+
+} // namespace sla
+} // namespace Slic3r

src/libslic3r/SLA/SupportPointGenerator.hpp

@@ -0,0 +1,234 @@
+#ifndef SLA_SUPPORTPOINTGENERATOR_HPP
+#define SLA_SUPPORTPOINTGENERATOR_HPP
+
+#include <random>
+
+#include <libslic3r/SLA/SupportPoint.hpp>
+#include <libslic3r/SLA/IndexedMesh.hpp>
+
+#include <libslic3r/BoundingBox.hpp>
+#include <libslic3r/ClipperUtils.hpp>
+#include <libslic3r/Point.hpp>
+
+#include <boost/container/small_vector.hpp>
+
+// #define SLA_SUPPORTPOINTGEN_DEBUG
+
+namespace Slic3r { namespace sla {
+
+class SupportPointGenerator {
+public:
+    struct Config {
+        float density_relative {1.f};
+        float minimal_distance {1.f};
+        float head_diameter {0.4f};
+
+        // Originally calibrated to 7.7f, reduced density by Tamas to 70% which is 11.1 (7.7 / 0.7) to adjust for new algorithm changes in tm_suppt_gen_improve
+        inline float support_force() const { return 11.1f / density_relative; } // a force one point can support       (arbitrary force unit)
+        inline float tear_pressure() const { return 1.f; }  // pressure that the display exerts    (the force unit per mm2)
+    };
+    
+    SupportPointGenerator(const IndexedMesh& emesh, const std::vector<ExPolygons>& slices,
+                    const std::vector<float>& heights, const Config& config, std::function<void(void)> throw_on_cancel, std::function<void(int)> statusfn);
+    
+    SupportPointGenerator(const IndexedMesh& emesh, const Config& config, std::function<void(void)> throw_on_cancel, std::function<void(int)> statusfn);
+    
+    const std::vector<SupportPoint>& output() const { return m_output; }
+    std::vector<SupportPoint>& output() { return m_output; }
+    
+    struct MyLayer;
+    
+    struct Structure {
+        Structure(MyLayer &layer, const ExPolygon& poly, const BoundingBox &bbox, const Vec2f &centroid, float area, float h) :
+            layer(&layer), polygon(&poly), bbox(bbox), centroid(centroid), area(area), zlevel(h)
+#ifdef SLA_SUPPORTPOINTGEN_DEBUG
+            , unique_id(std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()))
+#endif /* SLA_SUPPORTPOINTGEN_DEBUG */
+        {}
+        MyLayer *layer;
+        const ExPolygon* polygon = nullptr;
+        const BoundingBox bbox;
+        const Vec2f centroid = Vec2f::Zero();
+        const float area = 0.f;
+        float zlevel = 0;
+        // How well is this ExPolygon held to the print base?
+        // Positive number, the higher the better.
+        float supports_force_this_layer     = 0.f;
+        float supports_force_inherited      = 0.f;
+        float supports_force_total() const { return this->supports_force_this_layer + this->supports_force_inherited; }
+#ifdef SLA_SUPPORTPOINTGEN_DEBUG
+        std::chrono::milliseconds unique_id;
+#endif /* SLA_SUPPORTPOINTGEN_DEBUG */
+        
+        struct Link {
+            Link(Structure *island, float overlap_area) : island(island), overlap_area(overlap_area) {}
+            Structure   *island;
+            float        overlap_area;
+        };
+
+#ifdef NDEBUG
+        // In release mode, use the optimized container.
+        boost::container::small_vector<Link, 4> islands_above;
+        boost::container::small_vector<Link, 4> islands_below;
+#else
+        // In debug mode, use the standard vector, which is well handled by debugger visualizer.
+        std::vector<Link>					 	islands_above;
+        std::vector<Link>						islands_below;
+#endif
+        // Overhangs, that are dangling considerably.
+        ExPolygons                              dangling_areas;
+        // Complete overhands.
+        ExPolygons                              overhangs;
+        // Overhangs, where the surface must slope.
+        ExPolygons                              overhangs_slopes;
+        float                                   overhangs_area = 0.f;
+        
+        bool overlaps(const Structure &rhs) const { 
+            //FIXME ExPolygon::overlaps() shall be commutative, it is not!
+            return this->bbox.overlap(rhs.bbox) && (this->polygon->overlaps(*rhs.polygon) || rhs.polygon->overlaps(*this->polygon)); 
+        }
+        float overlap_area(const Structure &rhs) const { 
+            double out = 0.;
+            if (this->bbox.overlap(rhs.bbox)) {
+                Polygons polys = intersection(*this->polygon, *rhs.polygon);
+                for (const Polygon &poly : polys)
+                    out += poly.area();
+            }
+            return float(out);
+        }
+        float area_below() const { 
+            float area = 0.f; 
+            for (const Link &below : this->islands_below) 
+                area += below.island->area; 
+            return area;
+        }
+        Polygons polygons_below() const { 
+            size_t cnt = 0;
+            for (const Link &below : this->islands_below)
+                cnt += 1 + below.island->polygon->holes.size();
+            Polygons out;
+            out.reserve(cnt);
+            for (const Link &below : this->islands_below) {
+                out.emplace_back(below.island->polygon->contour);
+                append(out, below.island->polygon->holes);
+            }
+            return out;
+        }
+        ExPolygons expolygons_below() const { 
+            ExPolygons out;
+            out.reserve(this->islands_below.size());
+            for (const Link &below : this->islands_below)
+                out.emplace_back(*below.island->polygon);
+            return out;
+        }
+        // Positive deficit of the supports. If negative, this area is well supported. If positive, more supports need to be added.
+        float support_force_deficit(const float tear_pressure) const { return this->area * tear_pressure - this->supports_force_total(); }
+    };
+    
+    struct MyLayer {
+        MyLayer(const size_t layer_id, coordf_t print_z) : layer_id(layer_id), print_z(print_z) {}
+        size_t                  layer_id;
+        coordf_t                print_z;
+        std::vector<Structure>  islands;
+    };
+    
+    struct RichSupportPoint {
+        Vec3f        position;
+        Structure   *island;
+    };
+    
+    struct PointGrid3D {
+        struct GridHash {
+            std::size_t operator()(const Vec3i &cell_id) const {
+                return std::hash<int>()(cell_id.x()) ^ std::hash<int>()(cell_id.y() * 593) ^ std::hash<int>()(cell_id.z() * 7919);
+            }
+        };
+        typedef std::unordered_multimap<Vec3i, RichSupportPoint, GridHash> Grid;
+        
+        Vec3f   cell_size;
+        Grid    grid;
+        
+        Vec3i cell_id(const Vec3f &pos) {
+            return Vec3i(int(floor(pos.x() / cell_size.x())),
+                         int(floor(pos.y() / cell_size.y())),
+                         int(floor(pos.z() / cell_size.z())));
+        }
+        
+        void insert(const Vec2f &pos, Structure *island) {
+            RichSupportPoint pt;
+            pt.position = Vec3f(pos.x(), pos.y(), float(island->layer->print_z));
+            pt.island   = island;
+            grid.emplace(cell_id(pt.position), pt);
+        }
+        
+        bool collides_with(const Vec2f &pos, float print_z, float radius) {
+            Vec3f pos3d(pos.x(), pos.y(), print_z);
+            Vec3i cell = cell_id(pos3d);
+            std::pair<Grid::const_iterator, Grid::const_iterator> it_pair = grid.equal_range(cell);
+            if (collides_with(pos3d, radius, it_pair.first, it_pair.second))
+                return true;
+            for (int i = -1; i < 2; ++ i)
+                for (int j = -1; j < 2; ++ j)
+                    for (int k = -1; k < 1; ++ k) {
+                        if (i == 0 && j == 0 && k == 0)
+                            continue;
+                        it_pair = grid.equal_range(cell + Vec3i(i, j, k));
+                        if (collides_with(pos3d, radius, it_pair.first, it_pair.second))
+                            return true;
+                    }
+            return false;
+        }
+        
+    private:
+        bool collides_with(const Vec3f &pos, float radius, Grid::const_iterator it_begin, Grid::const_iterator it_end) {
+            for (Grid::const_iterator it = it_begin; it != it_end; ++ it) {
+                float dist2 = (it->second.position - pos).squaredNorm();
+                if (dist2 < radius * radius)
+                    return true;
+            }
+            return false;
+        }
+    };
+    
+    void execute(const std::vector<ExPolygons> &slices,
+                 const std::vector<float> &     heights);
+    
+    void seed(std::mt19937::result_type s) { m_rng.seed(s); }
+private:
+    std::vector<SupportPoint> m_output;
+    
+    SupportPointGenerator::Config m_config;
+    
+    void process(const std::vector<ExPolygons>& slices, const std::vector<float>& heights);
+
+public:
+    enum IslandCoverageFlags : uint8_t { icfNone = 0x0, icfIsNew = 0x1, icfWithBoundary = 0x2 };
+
+private:
+
+    void uniformly_cover(const ExPolygons& islands, Structure& structure, float deficit, PointGrid3D &grid3d, IslandCoverageFlags flags = icfNone);
+
+    void add_support_points(Structure& structure, PointGrid3D &grid3d);
+
+    void project_onto_mesh(std::vector<SupportPoint>& points) const;
+
+#ifdef SLA_SUPPORTPOINTGEN_DEBUG
+    static void output_expolygons(const ExPolygons& expolys, const std::string &filename);
+    static void output_structures(const std::vector<Structure> &structures);
+#endif // SLA_SUPPORTPOINTGEN_DEBUG
+    
+    const IndexedMesh& m_emesh;
+    std::function<void(void)> m_throw_on_cancel;
+    std::function<void(int)>  m_statusfn;
+    
+    std::mt19937 m_rng;
+};
+
+void remove_bottom_points(std::vector<SupportPoint> &pts, float lvl);
+
+std::vector<Vec2f> sample_expolygon(const ExPolygon &expoly, float samples_per_mm2, std::mt19937 &rng);
+void sample_expolygon_boundary(const ExPolygon &expoly, float samples_per_mm, std::vector<Vec2f> &out, std::mt19937 &rng);
+
+}} // namespace Slic3r::sla
+
+#endif // SUPPORTPOINTGENERATOR_HPP

src/libslic3r/SLA/SupportTree.cpp

@@ -0,0 +1,98 @@
+/**
+ * In this file we will implement the automatic SLA support tree generation.
+ *
+ */
+
+#include <numeric>
+#include <libslic3r/SLA/SupportTree.hpp>
+#include <libslic3r/SLA/SpatIndex.hpp>
+#include <libslic3r/SLA/SupportTreeBuilder.hpp>
+#include <libslic3r/SLA/SupportTreeBuildsteps.hpp>
+
+#include <libslic3r/MTUtils.hpp>
+#include <libslic3r/ClipperUtils.hpp>
+#include <libslic3r/Model.hpp>
+#include <libslic3r/TriangleMeshSlicer.hpp>
+
+#include <libnest2d/optimizers/nlopt/genetic.hpp>
+#include <libnest2d/optimizers/nlopt/subplex.hpp>
+#include <boost/log/trivial.hpp>
+#include <libslic3r/I18N.hpp>
+
+//! macro used to mark string used at localization,
+//! return same string
+#define L(s) Slic3r::I18N::translate(s)
+
+namespace Slic3r {
+namespace sla {
+
+void SupportTree::retrieve_full_mesh(indexed_triangle_set &outmesh) const {
+    its_merge(outmesh, retrieve_mesh(MeshType::Support));
+    its_merge(outmesh, retrieve_mesh(MeshType::Pad));
+}
+
+std::vector<ExPolygons> SupportTree::slice(const std::vector<float> &grid,
+                                           float                     cr) const
+{
+    const indexed_triangle_set &sup_mesh = retrieve_mesh(MeshType::Support);
+    const indexed_triangle_set &pad_mesh = retrieve_mesh(MeshType::Pad);
+
+    using Slices = std::vector<ExPolygons>;
+    auto slices = reserve_vector<Slices>(2);
+
+    if (!sup_mesh.empty()) {
+        slices.emplace_back();
+        slices.back() = slice_mesh_ex(sup_mesh, grid, cr, ctl().cancelfn);
+    }
+
+    if (!pad_mesh.empty()) {
+        slices.emplace_back();
+
+        auto bb = bounding_box(pad_mesh);
+        auto maxzit = std::upper_bound(grid.begin(), grid.end(), bb.max.z());
+        
+        auto cap = grid.end() - maxzit;
+        auto padgrid = reserve_vector<float>(size_t(cap > 0 ? cap : 0));
+        std::copy(grid.begin(), maxzit, std::back_inserter(padgrid));
+
+        slices.back() = slice_mesh_ex(pad_mesh, padgrid, cr, ctl().cancelfn);
+    }
+
+    size_t len = grid.size();
+    for (const Slices &slv : slices) { len = std::min(len, slv.size()); }
+
+    // Either the support or the pad or both has to be non empty
+    if (slices.empty()) return {};
+
+    Slices &mrg = slices.front();
+
+    for (auto it = std::next(slices.begin()); it != slices.end(); ++it) {
+        for (size_t i = 0; i < len; ++i) {
+            Slices &slv = *it;
+            std::copy(slv[i].begin(), slv[i].end(), std::back_inserter(mrg[i]));
+            slv[i] = {}; // clear and delete
+        }
+    }
+
+    return mrg;
+}
+
+SupportTree::UPtr SupportTree::create(const SupportableMesh &sm,
+                                      const JobController &  ctl)
+{
+    auto builder = make_unique<SupportTreeBuilder>();
+    builder->m_ctl = ctl;
+    
+    if (sm.cfg.enabled) {
+        // Execute takes care about the ground_level
+        SupportTreeBuildsteps::execute(*builder, sm);
+        builder->merge_and_cleanup();   // clean metadata, leave only the meshes.
+    } else {
+        // If a pad gets added later, it will be in the right Z level
+        builder->ground_level = sm.emesh.ground_level();
+    }
+    
+    return std::move(builder);
+}
+
+}} // namespace Slic3r::sla

src/libslic3r/SLA/SupportTree.hpp

@@ -0,0 +1,176 @@
+#ifndef SLA_SUPPORTTREE_HPP
+#define SLA_SUPPORTTREE_HPP
+
+#include <vector>
+#include <memory>
+#include <Eigen/Geometry>
+
+#include <libslic3r/SLA/Pad.hpp>
+#include <libslic3r/SLA/IndexedMesh.hpp>
+#include <libslic3r/SLA/SupportPoint.hpp>
+#include <libslic3r/SLA/JobController.hpp>
+
+namespace Slic3r {
+
+class TriangleMesh;
+class Model;
+class ModelInstance;
+class ModelObject;
+class Polygon;
+class ExPolygon;
+
+using Polygons = std::vector<Polygon>;
+using ExPolygons = std::vector<ExPolygon>;
+
+namespace sla {
+
+enum class PillarConnectionMode
+{
+    zigzag,
+    cross,
+    dynamic
+};
+
+struct SupportTreeConfig
+{
+    bool   enabled = true;
+    
+    // Radius in mm of the pointing side of the head.
+    double head_front_radius_mm = 0.2;
+
+    // How much the pinhead has to penetrate the model surface
+    double head_penetration_mm = 0.5;
+
+    // Radius of the back side of the 3d arrow.
+    double head_back_radius_mm = 0.5;
+
+    double head_fallback_radius_mm = 0.25;
+
+    // Width in mm from the back sphere center to the front sphere center.
+    double head_width_mm = 1.0;
+
+    // How to connect pillars
+    PillarConnectionMode pillar_connection_mode = PillarConnectionMode::dynamic;
+
+    // Only generate pillars that can be routed to ground
+    bool ground_facing_only = false;
+
+    // TODO: unimplemented at the moment. This coefficient will have an impact
+    // when bridges and pillars are merged. The resulting pillar should be a bit
+    // thicker than the ones merging into it. How much thicker? I don't know
+    // but it will be derived from this value.
+    double pillar_widening_factor = 0.5;
+
+    // Radius in mm of the pillar base.
+    double base_radius_mm = 2.0;
+
+    // The height of the pillar base cone in mm.
+    double base_height_mm = 1.0;
+
+    // The default angle for connecting support sticks and junctions.
+    double bridge_slope = M_PI/4;
+
+    // The max length of a bridge in mm
+    double max_bridge_length_mm = 10.0;
+
+    // The max distance of a pillar to pillar link.
+    double max_pillar_link_distance_mm = 10.0;
+
+    // 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;
+    
+    unsigned max_bridges_on_pillar = 3;
+    
+    double head_fullwidth() const {
+        return 2 * head_front_radius_mm + head_width_mm +
+               2 * head_back_radius_mm - head_penetration_mm;
+    }
+
+    // /////////////////////////////////////////////////////////////////////////
+    // Compile time configuration values (candidates for runtime)
+    // /////////////////////////////////////////////////////////////////////////
+
+    // The max Z angle for a normal at which it will get completely ignored.
+    static const double constexpr normal_cutoff_angle = 150.0 * M_PI / 180.0;
+
+    // The shortest distance of any support structure from the model surface
+    static const double constexpr safety_distance_mm = 0.5;
+
+    static const double constexpr max_solo_pillar_height_mm = 15.0;
+    static const double constexpr max_dual_pillar_height_mm = 35.0;
+    static const double constexpr optimizer_rel_score_diff = 1e-6;
+    static const unsigned constexpr optimizer_max_iterations = 1000;
+    static const unsigned constexpr pillar_cascade_neighbors = 3;
+    
+};
+
+// TODO: Part of future refactor
+//class SupportConfig {
+//    std::optional<SupportTreeConfig> tree_cfg {std::in_place_t{}}; // fill up
+//    std::optional<PadConfig>         pad_cfg;
+//};
+
+enum class MeshType { Support, Pad };
+
+struct SupportableMesh
+{
+    IndexedMesh  emesh;
+    SupportPoints pts;
+    SupportTreeConfig cfg;
+//    PadConfig     pad_cfg;
+
+    explicit SupportableMesh(const indexed_triangle_set & trmsh,
+                             const SupportPoints &sp,
+                             const SupportTreeConfig &c)
+        : emesh{trmsh}, pts{sp}, cfg{c}
+    {}
+    
+    explicit SupportableMesh(const IndexedMesh   &em,
+                             const SupportPoints &sp,
+                             const SupportTreeConfig &c)
+        : emesh{em}, pts{sp}, cfg{c}
+    {}
+};
+
+/// The class containing mesh data for the generated supports.
+class SupportTree
+{
+    JobController m_ctl;
+public:
+    using UPtr = std::unique_ptr<SupportTree>;
+    
+    static UPtr create(const SupportableMesh &input,
+                       const JobController &ctl = {});
+
+    virtual ~SupportTree() = default;
+
+    virtual const indexed_triangle_set &retrieve_mesh(MeshType meshtype) const = 0;
+
+    /// Adding the "pad" under the supports.
+    /// modelbase will be used according to the embed_object flag in PoolConfig.
+    /// If set, the plate will be interpreted as the model's intrinsic pad. 
+    /// Otherwise, the modelbase will be unified with the base plate calculated
+    /// from the supports.
+    virtual const indexed_triangle_set &add_pad(const ExPolygons &modelbase,
+                                                const PadConfig & pcfg) = 0;
+
+    virtual void remove_pad() = 0;
+    
+    std::vector<ExPolygons> slice(const std::vector<float> &,
+                                  float closing_radius) const;
+    
+    void retrieve_full_mesh(indexed_triangle_set &outmesh) const;
+    
+    const JobController &ctl() const { return m_ctl; }
+};
+
+}
+
+}
+
+#endif // SLASUPPORTTREE_HPP

src/libslic3r/SLA/SupportTreeBuilder.cpp

@@ -0,0 +1,225 @@
+#define NOMINMAX
+
+#include <libslic3r/SLA/SupportTreeBuilder.hpp>
+#include <libslic3r/SLA/SupportTreeBuildsteps.hpp>
+#include <libslic3r/SLA/SupportTreeMesher.hpp>
+//#include <libslic3r/SLA/Contour3D.hpp>
+
+namespace Slic3r {
+namespace sla {
+
+Head::Head(double       r_big_mm,
+           double       r_small_mm,
+           double       length_mm,
+           double       penetration,
+           const Vec3d &direction,
+           const Vec3d &offset)
+    : dir(direction)
+    , pos(offset)
+    , r_back_mm(r_big_mm)
+    , r_pin_mm(r_small_mm)
+    , width_mm(length_mm)
+    , penetration_mm(penetration)
+{
+}
+
+Pad::Pad(const indexed_triangle_set &support_mesh,
+         const ExPolygons &          model_contours,
+         double                      ground_level,
+         const PadConfig &           pcfg,
+         ThrowOnCancel               thr)
+    : cfg(pcfg)
+    , zlevel(ground_level + pcfg.full_height() - pcfg.required_elevation())
+{
+    thr();
+    
+    ExPolygons sup_contours;
+    
+    float zstart = float(zlevel);
+    float zend   = zstart + float(pcfg.full_height() + EPSILON);
+    
+    pad_blueprint(support_mesh, sup_contours, grid(zstart, zend, 0.1f), thr);
+    create_pad(sup_contours, model_contours, tmesh, pcfg);
+    
+    Vec3f offs{.0f, .0f, float(zlevel)};
+    for (auto &p : tmesh.vertices) p += offs;
+
+    its_merge_vertices(tmesh);
+}
+
+const indexed_triangle_set &SupportTreeBuilder::add_pad(
+    const ExPolygons &modelbase, const PadConfig &cfg)
+{
+    m_pad = Pad{merged_mesh(), modelbase, ground_level, cfg, ctl().cancelfn};
+    return m_pad.tmesh;
+}
+
+SupportTreeBuilder::SupportTreeBuilder(SupportTreeBuilder &&o)
+    : m_heads(std::move(o.m_heads))
+    , m_head_indices{std::move(o.m_head_indices)}
+    , m_pillars{std::move(o.m_pillars)}
+    , m_bridges{std::move(o.m_bridges)}
+    , m_crossbridges{std::move(o.m_crossbridges)}
+    , m_pad{std::move(o.m_pad)}
+    , m_meshcache{std::move(o.m_meshcache)}
+    , m_meshcache_valid{o.m_meshcache_valid}
+    , m_model_height{o.m_model_height}
+    , ground_level{o.ground_level}
+{}
+
+SupportTreeBuilder::SupportTreeBuilder(const SupportTreeBuilder &o)
+    : m_heads(o.m_heads)
+    , m_head_indices{o.m_head_indices}
+    , m_pillars{o.m_pillars}
+    , m_bridges{o.m_bridges}
+    , m_crossbridges{o.m_crossbridges}
+    , m_pad{o.m_pad}
+    , m_meshcache{o.m_meshcache}
+    , m_meshcache_valid{o.m_meshcache_valid}
+    , m_model_height{o.m_model_height}
+    , ground_level{o.ground_level}
+{}
+
+SupportTreeBuilder &SupportTreeBuilder::operator=(SupportTreeBuilder &&o)
+{
+    m_heads = std::move(o.m_heads);
+    m_head_indices = std::move(o.m_head_indices);
+    m_pillars = std::move(o.m_pillars);
+    m_bridges = std::move(o.m_bridges);
+    m_crossbridges = std::move(o.m_crossbridges);
+    m_pad = std::move(o.m_pad);
+    m_meshcache = std::move(o.m_meshcache);
+    m_meshcache_valid = o.m_meshcache_valid;
+    m_model_height = o.m_model_height;
+    ground_level = o.ground_level;
+    return *this;
+}
+
+SupportTreeBuilder &SupportTreeBuilder::operator=(const SupportTreeBuilder &o)
+{
+    m_heads = o.m_heads;
+    m_head_indices = o.m_head_indices;
+    m_pillars = o.m_pillars;
+    m_bridges = o.m_bridges;
+    m_crossbridges = o.m_crossbridges;
+    m_pad = o.m_pad;
+    m_meshcache = o.m_meshcache;
+    m_meshcache_valid = o.m_meshcache_valid;
+    m_model_height = o.m_model_height;
+    ground_level = o.ground_level;
+    return *this;
+}
+
+void SupportTreeBuilder::add_pillar_base(long pid, double baseheight, double radius)
+{
+    std::lock_guard<Mutex> lk(m_mutex);
+    assert(pid >= 0 && size_t(pid) < m_pillars.size());
+    Pillar& pll = m_pillars[size_t(pid)];
+    m_pedestals.emplace_back(pll.endpt, std::min(baseheight, pll.height),
+                             std::max(radius, pll.r), pll.r);
+
+    m_pedestals.back().id = m_pedestals.size() - 1;
+    m_meshcache_valid = false;
+}
+
+const indexed_triangle_set &SupportTreeBuilder::merged_mesh(size_t steps) const
+{
+    if (m_meshcache_valid) return m_meshcache;
+    
+    indexed_triangle_set merged;
+    
+    for (auto &head : m_heads) {
+        if (ctl().stopcondition()) break;
+        if (head.is_valid()) its_merge(merged, get_mesh(head, steps));
+    }
+    
+    for (auto &pill : m_pillars) {
+        if (ctl().stopcondition()) break;
+        its_merge(merged, get_mesh(pill, steps));
+    }
+
+    for (auto &pedest : m_pedestals) {
+        if (ctl().stopcondition()) break;
+        its_merge(merged, get_mesh(pedest, steps));
+    }
+    
+    for (auto &j : m_junctions) {
+        if (ctl().stopcondition()) break;
+        its_merge(merged, get_mesh(j, steps));
+    }
+
+    for (auto &bs : m_bridges) {
+        if (ctl().stopcondition()) break;
+        its_merge(merged, get_mesh(bs, steps));
+    }
+    
+    for (auto &bs : m_crossbridges) {
+        if (ctl().stopcondition()) break;
+        its_merge(merged, get_mesh(bs, steps));
+    }
+
+    for (auto &bs : m_diffbridges) {
+        if (ctl().stopcondition()) break;
+        its_merge(merged, get_mesh(bs, steps));
+    }
+
+    for (auto &anch : m_anchors) {
+        if (ctl().stopcondition()) break;
+        its_merge(merged, get_mesh(anch, steps));
+    }
+
+    if (ctl().stopcondition()) {
+        // In case of failure we have to return an empty mesh
+        m_meshcache = {};
+        return m_meshcache;
+    }
+    
+    m_meshcache = merged;
+    
+    // The mesh will be passed by const-pointer to TriangleMeshSlicer,
+    // which will need this.
+    its_merge_vertices(m_meshcache);
+    
+    BoundingBoxf3 bb = bounding_box(m_meshcache);
+    m_model_height   = bb.max(Z) - bb.min(Z);
+
+    m_meshcache_valid = true;
+    return m_meshcache;
+}
+
+double SupportTreeBuilder::full_height() const
+{
+    if (merged_mesh().indices.empty() && !pad().empty())
+        return pad().cfg.full_height();
+    
+    double h = mesh_height();
+    if (!pad().empty()) h += pad().cfg.required_elevation();
+    return h;
+}
+
+const indexed_triangle_set &SupportTreeBuilder::merge_and_cleanup()
+{
+    // in case the mesh is not generated, it should be...
+    auto &ret = merged_mesh(); 
+    
+    // Doing clear() does not garantee to release the memory.
+    m_heads = {};
+    m_head_indices = {};
+    m_pillars = {};
+    m_junctions = {};
+    m_bridges = {};
+    
+    return ret;
+}
+
+const indexed_triangle_set &SupportTreeBuilder::retrieve_mesh(MeshType meshtype) const
+{
+    switch(meshtype) {
+    case MeshType::Support: return merged_mesh();
+    case MeshType::Pad:     return pad().tmesh;
+    }
+    
+    return m_meshcache;
+}
+
+}} // namespace Slic3r::sla

src/libslic3r/SLA/SupportTreeBuilder.hpp

@@ -0,0 +1,450 @@
+#ifndef SLA_SUPPORTTREEBUILDER_HPP
+#define SLA_SUPPORTTREEBUILDER_HPP
+
+#include <libslic3r/SLA/Concurrency.hpp>
+#include <libslic3r/SLA/SupportTree.hpp>
+//#include <libslic3r/SLA/Contour3D.hpp>
+#include <libslic3r/TriangleMesh.hpp>
+#include <libslic3r/SLA/Pad.hpp>
+#include <libslic3r/MTUtils.hpp>
+
+namespace Slic3r {
+namespace sla {
+
+/**
+ * Terminology:
+ *
+ * Support point:
+ * The point on the model surface that needs support.
+ *
+ * Pillar:
+ * A thick column that spans from a support point to the ground and has
+ * a thick cone shaped base where it touches the ground.
+ *
+ * Ground facing support point:
+ * A support point that can be directly connected with the ground with a pillar
+ * that does not collide or cut through the model.
+ *
+ * Non ground facing support point:
+ * A support point that cannot be directly connected with the ground (only with
+ * the model surface).
+ *
+ * Head:
+ * The pinhead that connects to the model surface with the sharp end end
+ * to a pillar or bridge stick with the dull end.
+ *
+ * Headless support point:
+ * A support point on the model surface for which there is not enough place for
+ * the head. It is either in a hole or there is some barrier that would collide
+ * with the head geometry. The headless support point can be ground facing and
+ * non ground facing as well.
+ *
+ * Bridge:
+ * A stick that connects two pillars or a head with a pillar.
+ *
+ * Junction:
+ * A small ball in the intersection of two or more sticks (pillar, bridge, ...)
+ *
+ * CompactBridge:
+ * A bridge that connects a headless support point with the model surface or a
+ * nearby pillar.
+ */
+
+template<class Vec> double distance(const Vec& p) {
+    return std::sqrt(p.transpose() * p);
+}
+
+template<class Vec> double distance(const Vec& pp1, const Vec& pp2) {
+    auto p = pp2 - pp1;
+    return distance(p);
+}
+
+const Vec3d DOWN = {0.0, 0.0, -1.0};
+
+struct SupportTreeNode
+{
+    static const constexpr long ID_UNSET = -1;
+
+    long id = ID_UNSET; // For identification withing a tree.
+};
+
+// A pinhead originating from a support point
+struct Head: public SupportTreeNode {
+    Vec3d dir = DOWN;
+    Vec3d pos = {0, 0, 0};
+    
+    double r_back_mm = 1;
+    double r_pin_mm = 0.5;
+    double width_mm = 2;
+    double penetration_mm = 0.5;
+
+    
+    // If there is a pillar connecting to this head, then the id will be set.
+    long pillar_id = ID_UNSET;
+    
+    long bridge_id = ID_UNSET;
+    
+    inline void invalidate() { id = ID_UNSET; }
+    inline bool is_valid() const { return id >= 0; }
+    
+    Head(double r_big_mm,
+         double r_small_mm,
+         double length_mm,
+         double penetration,
+         const Vec3d &direction = DOWN,  // direction (normal to the dull end)
+         const Vec3d &offset = {0, 0, 0}      // displacement
+         );
+
+    inline double real_width() const
+    {
+        return 2 * r_pin_mm + width_mm + 2 * r_back_mm ;
+    }
+
+    inline double fullwidth() const
+    {
+        return real_width() - penetration_mm;
+    }
+    
+    inline Vec3d junction_point() const
+    {
+        return pos + (fullwidth() - r_back_mm) * dir;
+    }
+    
+    inline double request_pillar_radius(double radius) const
+    {
+        const double rmax = r_back_mm;
+        return radius > 0 && radius < rmax ? radius : rmax;
+    }
+};
+
+// A junction connecting bridges and pillars
+struct Junction: public SupportTreeNode {
+    double r = 1;
+    Vec3d pos;
+
+    Junction(const Vec3d &tr, double r_mm) : r(r_mm), pos(tr) {}
+};
+
+struct Pillar: public SupportTreeNode {
+    double height, r;
+    Vec3d endpt;
+    
+    // If the pillar connects to a head, this is the id of that head
+    bool starts_from_head = true; // Could start from a junction as well
+    long start_junction_id = ID_UNSET;
+    
+    // How many bridges are connected to this pillar
+    unsigned bridges = 0;
+    
+    // How many pillars are cascaded with this one
+    unsigned links = 0;
+
+    Pillar(const Vec3d &endp, double h, double radius = 1.):
+        height{h}, r(radius), endpt(endp), starts_from_head(false) {}
+
+    Vec3d startpoint() const
+    {
+        return {endpt.x(), endpt.y(), endpt.z() + height};
+    }
+    
+    const Vec3d& endpoint() const { return endpt; }
+};
+
+// A base for pillars or bridges that end on the ground
+struct Pedestal: public SupportTreeNode {
+    Vec3d pos;
+    double height, r_bottom, r_top;
+
+    Pedestal(const Vec3d &p, double h, double rbottom, double rtop)
+        : pos{p}, height{h}, r_bottom{rbottom}, r_top{rtop}
+    {}
+};
+
+// This is the thing that anchors a pillar or bridge to the model body.
+// It is actually a reverse pinhead.
+struct Anchor: public Head { using Head::Head; };
+
+// A Bridge between two pillars (with junction endpoints)
+struct Bridge: public SupportTreeNode {
+    double r = 0.8;
+    Vec3d startp = Vec3d::Zero(), endp = Vec3d::Zero();
+    
+    Bridge(const Vec3d &j1,
+           const Vec3d &j2,
+           double       r_mm  = 0.8): r{r_mm}, startp{j1}, endp{j2}
+    {}
+
+    double get_length() const { return (endp - startp).norm(); }
+    Vec3d  get_dir() const { return (endp - startp).normalized(); }
+};
+
+struct DiffBridge: public Bridge {
+    double end_r;
+
+    DiffBridge(const Vec3d &p_s, const Vec3d &p_e, double r_s, double r_e)
+        : Bridge{p_s, p_e, r_s}, end_r{r_e}
+    {}
+};
+
+// A wrapper struct around the pad
+struct Pad {
+    indexed_triangle_set tmesh;
+    PadConfig cfg;
+    double zlevel = 0;
+    
+    Pad() = default;
+
+    Pad(const indexed_triangle_set &support_mesh,
+        const ExPolygons &          model_contours,
+        double                      ground_level,
+        const PadConfig &           pcfg,
+        ThrowOnCancel               thr);
+
+    bool empty() const { return tmesh.indices.size() == 0; }
+};
+
+// This class will hold the support tree meshes with some additional
+// bookkeeping as well. Various parts of the support geometry are stored
+// separately and are merged when the caller queries the merged mesh. The
+// merged result is cached for fast subsequent delivery of the merged mesh
+// which can be quite complex. The support tree creation algorithm can use an
+// instance of this class as a somewhat higher level tool for crafting the 3D
+// support mesh. Parts can be added with the appropriate methods such as
+// add_head or add_pillar which forwards the constructor arguments and fills
+// the IDs of these substructures. The IDs are basically indices into the
+// arrays of the appropriate type (heads, pillars, etc...). One can later query
+// e.g. a pillar for a specific head...
+//
+// 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 SupportTreeBuilder: public SupportTree {
+    // 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<Bridge>     m_crossbridges;
+    std::vector<DiffBridge> m_diffbridges;
+    std::vector<Pedestal>   m_pedestals;
+    std::vector<Anchor>     m_anchors;
+
+    Pad m_pad;
+    
+    using Mutex = ccr::SpinningMutex;
+    
+    mutable indexed_triangle_set m_meshcache;
+    mutable Mutex m_mutex;
+    mutable bool m_meshcache_valid = false;
+    mutable double m_model_height = 0; // the full height of the model
+    
+    template<class BridgeT, class...Args>
+    const BridgeT& _add_bridge(std::vector<BridgeT> &br, Args&&... args)
+    {
+        std::lock_guard<Mutex> lk(m_mutex);
+        br.emplace_back(std::forward<Args>(args)...);
+        br.back().id = long(br.size() - 1);
+        m_meshcache_valid = false;
+        return br.back();
+    }
+    
+public:
+    double ground_level = 0;
+    
+    SupportTreeBuilder() = default;
+    SupportTreeBuilder(SupportTreeBuilder &&o);
+    SupportTreeBuilder(const SupportTreeBuilder &o);
+    SupportTreeBuilder& operator=(SupportTreeBuilder &&o);
+    SupportTreeBuilder& operator=(const SupportTreeBuilder &o);
+
+    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;
+        
+        m_meshcache_valid = false;
+        return m_heads.back();
+    }
+    
+    template<class...Args> long add_pillar(long headid, double length)
+    {
+        std::lock_guard<Mutex> lk(m_mutex);
+        if (m_pillars.capacity() < m_heads.size())
+            m_pillars.reserve(m_heads.size() * 10);
+        
+        assert(headid >= 0 && size_t(headid) < m_head_indices.size());
+        Head &head = m_heads[m_head_indices[size_t(headid)]];
+        
+        Vec3d hjp = head.junction_point() - Vec3d{0, 0, length};
+        m_pillars.emplace_back(hjp, length, head.r_back_mm);
+
+        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;
+        
+        m_meshcache_valid = false;
+        return pillar.id;
+    }
+    
+    void add_pillar_base(long pid, double baseheight = 3, double radius = 2);
+
+    template<class...Args> const Anchor& add_anchor(Args&&...args)
+    {
+        std::lock_guard<Mutex> lk(m_mutex);
+        m_anchors.emplace_back(std::forward<Args>(args)...);
+        m_anchors.back().id = long(m_junctions.size() - 1);
+        m_meshcache_valid = false;
+        return m_anchors.back();
+    }
+    
+    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())
+            m_pillars[size_t(pillar.id)].bridges++;
+    }
+    
+    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++;
+    }
+    
+    unsigned bridgecount(const Pillar &pillar) const {
+        std::lock_guard<Mutex> lk(m_mutex);
+        assert(pillar.id >= 0 && size_t(pillar.id) < m_pillars.size());
+        return pillar.bridges;
+    }
+    
+    template<class...Args> long add_pillar(Args&&...args)
+    {
+        std::lock_guard<Mutex> lk(m_mutex);
+        if (m_pillars.capacity() < m_heads.size())
+            m_pillars.reserve(m_heads.size() * 10);
+        
+        m_pillars.emplace_back(std::forward<Args>(args)...);
+        Pillar& pillar = m_pillars.back();
+        pillar.id = long(m_pillars.size() - 1);
+        pillar.starts_from_head = false;
+        m_meshcache_valid = false;
+        return pillar.id;
+    }
+    
+    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);
+        m_meshcache_valid = false;
+        return m_junctions.back();
+    }
+    
+    const Bridge& add_bridge(const Vec3d &s, const Vec3d &e, double r)
+    {
+        return _add_bridge(m_bridges, s, e, r);
+    }
+    
+    const Bridge& add_bridge(long headid, const Vec3d &endp)
+    {
+        std::lock_guard<Mutex> lk(m_mutex);
+        assert(headid >= 0 && size_t(headid) < m_head_indices.size());
+        
+        Head &h = m_heads[m_head_indices[size_t(headid)]];
+        m_bridges.emplace_back(h.junction_point(), endp, h.r_back_mm);
+        m_bridges.back().id = long(m_bridges.size() - 1);
+        
+        h.bridge_id = m_bridges.back().id;
+        m_meshcache_valid = false;
+        return m_bridges.back();
+    }
+    
+    template<class...Args> const Bridge& add_crossbridge(Args&&... args)
+    {
+        return _add_bridge(m_crossbridges, std::forward<Args>(args)...);
+    }
+
+    template<class...Args> const DiffBridge& add_diffbridge(Args&&... args)
+    {
+        return _add_bridge(m_diffbridges, std::forward<Args>(args)...);
+    }
+    
+    Head &head(unsigned id)
+    {
+        std::lock_guard<Mutex> lk(m_mutex);
+        assert(id < m_head_indices.size());
+        
+        m_meshcache_valid = false;
+        return m_heads[m_head_indices[id]];
+    }
+    
+    inline size_t pillarcount() const {
+        std::lock_guard<Mutex> lk(m_mutex);
+        return m_pillars.size();
+    }
+    
+    inline const std::vector<Pillar> &pillars() const { return m_pillars; }
+    inline const std::vector<Head>   &heads() const { return m_heads; }
+    inline const std::vector<Bridge> &bridges() const { return m_bridges; }
+    inline const std::vector<Bridge> &crossbridges() const { return m_crossbridges; }
+    
+    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)];
+    }
+    
+    template<class T> inline IntegerOnly<T, Pillar&> pillar(T id) 
+    {
+        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& pad() const { return m_pad; }
+    
+    // WITHOUT THE PAD!!!
+    const indexed_triangle_set &merged_mesh(size_t steps = 45) const;
+    
+    // WITH THE PAD
+    double full_height() const;
+    
+    // WITHOUT THE PAD!!!
+    inline double mesh_height() const
+    {
+        if (!m_meshcache_valid) merged_mesh();
+        return m_model_height;
+    }
+    
+    // Intended to be called after the generation is fully complete
+    const indexed_triangle_set & merge_and_cleanup();
+    
+    // Implement SupportTree interface:
+
+    const indexed_triangle_set &add_pad(const ExPolygons &modelbase,
+                                        const PadConfig & pcfg) override;
+
+    void remove_pad() override { m_pad = Pad(); }
+
+    virtual const indexed_triangle_set &retrieve_mesh(
+        MeshType meshtype = MeshType::Support) const override;
+};
+
+}} // namespace Slic3r::sla
+
+#endif // SUPPORTTREEBUILDER_HPP

src/libslic3r/SLA/SupportTreeBuildsteps.hpp

@@ -0,0 +1,378 @@
+#ifndef SLASUPPORTTREEALGORITHM_H
+#define SLASUPPORTTREEALGORITHM_H
+
+#include <cstdint>
+#include <optional>
+
+#include <libslic3r/SLA/SupportTreeBuilder.hpp>
+#include <libslic3r/SLA/Clustering.hpp>
+#include <libslic3r/SLA/SpatIndex.hpp>
+
+namespace Slic3r {
+namespace sla {
+
+// The minimum distance for two support points to remain valid.
+const double /*constexpr*/ D_SP = 0.1;
+
+enum { // For indexing Eigen vectors as v(X), v(Y), v(Z) instead of numbers
+    X, Y, Z
+};
+
+inline Vec2d to_vec2(const Vec3d &v3) { return {v3(X), v3(Y)}; }
+
+inline std::pair<double, double> dir_to_spheric(const Vec3d &n, double norm = 1.)
+{
+    double z       = n.z();
+    double r       = norm;
+    double polar   = std::acos(z / r);
+    double azimuth = std::atan2(n(1), n(0));
+    return {polar, azimuth};
+}
+
+inline Vec3d spheric_to_dir(double polar, double azimuth)
+{
+    return {std::cos(azimuth) * std::sin(polar),
+            std::sin(azimuth) * std::sin(polar), std::cos(polar)};
+}
+
+inline Vec3d spheric_to_dir(const std::tuple<double, double> &v)
+{
+    auto [plr, azm] = v;
+    return spheric_to_dir(plr, azm);
+}
+
+inline Vec3d spheric_to_dir(const std::pair<double, double> &v)
+{
+    return spheric_to_dir(v.first, v.second);
+}
+
+inline Vec3d spheric_to_dir(const std::array<double, 2> &v)
+{
+    return spheric_to_dir(v[0], v[1]);
+}
+
+// Give points on a 3D ring with given center, radius and orientation
+// method based on:
+// https://math.stackexchange.com/questions/73237/parametric-equation-of-a-circle-in-3d-space
+template<size_t N>
+class PointRing {
+    std::array<double, N> m_phis;
+
+    // Two vectors that will be perpendicular to each other and to the
+    // axis. Values for a(X) and a(Y) are now arbitrary, a(Z) is just a
+    // placeholder.
+    // a and b vectors are perpendicular to the ring direction and to each other.
+    // Together they define the plane where we have to iterate with the
+    // given angles in the 'm_phis' vector
+    Vec3d a = {0, 1, 0}, b;
+    double m_radius = 0.;
+
+    static inline bool constexpr is_one(double val)
+    {
+        return std::abs(std::abs(val) - 1) < 1e-20;
+    }
+
+public:
+
+    PointRing(const Vec3d &n)
+    {
+        m_phis = linspace_array<N>(0., 2 * PI);
+
+        // We have to address the case when the direction vector v (same as
+        // dir) is coincident with one of the world axes. In this case two of
+        // its components will be completely zero and one is 1.0. Our method
+        // becomes dangerous here due to division with zero. Instead, vector
+        // 'a' can be an element-wise rotated version of 'v'
+        if(is_one(n(X)) || is_one(n(Y)) || is_one(n(Z))) {
+            a = {n(Z), n(X), n(Y)};
+            b = {n(Y), n(Z), n(X)};
+        }
+        else {
+            a(Z) = -(n(Y)*a(Y)) / n(Z); a.normalize();
+            b = a.cross(n);
+        }
+    }
+
+    Vec3d get(size_t idx, const Vec3d src, double r) const
+    {
+        double phi = m_phis[idx];
+        double sinphi = std::sin(phi);
+        double cosphi = std::cos(phi);
+
+        double rpscos = r * cosphi;
+        double rpssin = r * sinphi;
+
+        // Point on the sphere
+        return {src(X) + rpscos * a(X) + rpssin * b(X),
+                src(Y) + rpscos * a(Y) + rpssin * b(Y),
+                src(Z) + rpscos * a(Z) + rpssin * b(Z)};
+    }
+};
+
+//IndexedMesh::hit_result query_hit(const SupportableMesh &msh, const Bridge &br, double safety_d = std::nan(""));
+//IndexedMesh::hit_result query_hit(const SupportableMesh &msh, const Head &br, double safety_d = std::nan(""));
+
+inline Vec3d dirv(const Vec3d& startp, const Vec3d& endp) {
+    return (endp - startp).normalized();
+}
+
+class PillarIndex {
+    PointIndex m_index;
+    using Mutex = ccr::BlockingMutex;
+    mutable Mutex m_mutex;
+
+public:
+
+    template<class...Args> inline void guarded_insert(Args&&...args)
+    {
+        std::lock_guard<Mutex> lck(m_mutex);
+        m_index.insert(std::forward<Args>(args)...);
+    }
+
+    template<class...Args>
+    inline std::vector<PointIndexEl> guarded_query(Args&&...args) const
+    {
+        std::lock_guard<Mutex> lck(m_mutex);
+        return m_index.query(std::forward<Args>(args)...);
+    }
+
+    template<class...Args> inline void insert(Args&&...args)
+    {
+        m_index.insert(std::forward<Args>(args)...);
+    }
+
+    template<class...Args>
+    inline std::vector<PointIndexEl> query(Args&&...args) const
+    {
+        return m_index.query(std::forward<Args>(args)...);
+    }
+
+    template<class Fn> inline void foreach(Fn fn) { m_index.foreach(fn); }
+    template<class Fn> inline void guarded_foreach(Fn fn)
+    {
+        std::lock_guard<Mutex> lck(m_mutex);
+        m_index.foreach(fn);
+    }
+
+    PointIndex guarded_clone()
+    {
+        std::lock_guard<Mutex> lck(m_mutex);
+        return m_index;
+    }
+};
+
+// Helper function for pillar interconnection where pairs of already connected
+// pillars should be checked for not to be processed again. This can be done
+// in constant time with a 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. The hash value has to have twice as many bits as the
+// arguments need. If the same integral type is used for args and return val,
+// make sure the arguments use only the half of the type's bit depth.
+template<class I, class DoubleI = IntegerOnly<I>>
+IntegerOnly<DoubleI> pairhash(I a, I b)
+{
+    using std::ceil; using std::log2; using std::max; using std::min;
+    static const auto constexpr Ibits = int(sizeof(I) * CHAR_BIT);
+    static const auto constexpr DoubleIbits = int(sizeof(DoubleI) * CHAR_BIT);
+    static const auto constexpr shift = DoubleIbits / 2 < Ibits ? Ibits / 2 : Ibits;
+
+    I g = min(a, b), l = max(a, b);
+
+    // Assume the hash will fit into the output variable
+    assert((g ? (ceil(log2(g))) : 0) <= shift);
+    assert((l ? (ceil(log2(l))) : 0) <= shift);
+
+    return (DoubleI(g) << shift) + l;
+}
+
+class SupportTreeBuildsteps {
+    const SupportTreeConfig& m_cfg;
+    const IndexedMesh& m_mesh;
+    const std::vector<SupportPoint>& m_support_pts;
+
+    using PtIndices = std::vector<unsigned>;
+
+    PtIndices m_iheads;            // support points with pinhead
+    PtIndices m_iheads_onmodel;
+    PtIndices m_iheadless;         // headless support points
+    
+    std::map<unsigned, IndexedMesh::hit_result> m_head_to_ground_scans;
+
+    // normals for support points from model faces.
+    PointSet  m_support_nmls;
+
+    // Clusters of points which can reach the ground directly and can be
+    // bridged to one central pillar
+    std::vector<PtIndices> m_pillar_clusters;
+
+    // This algorithm uses the SupportTreeBuilder class to fill gradually
+    // the support elements (heads, pillars, bridges, ...)
+    SupportTreeBuilder& m_builder;
+
+    // support points in Eigen/IGL format
+    PointSet m_points;
+
+    // throw if canceled: It will be called many times so a shorthand will
+    // come in handy.
+    ThrowOnCancel m_thr;
+
+    // A spatial index to easily find strong pillars to connect to.
+    PillarIndex m_pillar_index;
+
+    // When bridging heads to pillars... TODO: find a cleaner solution
+    ccr::BlockingMutex m_bridge_mutex;
+
+    inline IndexedMesh::hit_result ray_mesh_intersect(const Vec3d& s, 
+                                                      const Vec3d& dir)
+    {
+        return m_mesh.query_ray_hit(s, dir);
+    }
+
+    // This function will test if a future pinhead would not collide with the
+    // model geometry. It does not take a 'Head' object because those are
+    // created after this test. Parameters: s: The touching point on the model
+    // surface. dir: This is the direction of the head from the pin to the back
+    // r_pin, r_back: the radiuses of the pin and the back sphere width: This
+    // is the full width from the pin center to the back center m: The object
+    // mesh.
+    // The return value is the hit result from the ray casting. If the starting
+    // point was inside the model, an "invalid" hit_result will be returned
+    // with a zero distance value instead of a NAN. This way the result can
+    // be used safely for comparison with other distances.
+    IndexedMesh::hit_result pinhead_mesh_intersect(
+        const Vec3d& s,
+        const Vec3d& dir,
+        double r_pin,
+        double r_back,
+        double width,
+        double safety_d);
+
+    IndexedMesh::hit_result pinhead_mesh_intersect(
+        const Vec3d& s,
+        const Vec3d& dir,
+        double r_pin,
+        double r_back,
+        double width)
+    {
+        return pinhead_mesh_intersect(s, dir, r_pin, r_back, width,
+                                      r_back * m_cfg.safety_distance_mm /
+                                          m_cfg.head_back_radius_mm);
+    }
+
+    // Checking bridge (pillar and stick as well) intersection with the model.
+    // If the function is used for headless sticks, the ins_check parameter
+    // have to be true as the beginning of the stick might be inside the model
+    // geometry.
+    // The return value is the hit result from the ray casting. If the starting
+    // point was inside the model, an "invalid" hit_result will be returned
+    // with a zero distance value instead of a NAN. This way the result can
+    // be used safely for comparison with other distances.
+    IndexedMesh::hit_result bridge_mesh_intersect(
+        const Vec3d& s,
+        const Vec3d& dir,
+        double r,
+        double safety_d);
+
+    IndexedMesh::hit_result bridge_mesh_intersect(
+        const Vec3d& s,
+        const Vec3d& dir,
+        double r)
+    {
+        return bridge_mesh_intersect(s, dir, r,
+                                     r * m_cfg.safety_distance_mm /
+                                         m_cfg.head_back_radius_mm);
+    }
+    
+    template<class...Args>
+    inline double bridge_mesh_distance(Args&&...args) {
+        return bridge_mesh_intersect(std::forward<Args>(args)...).distance();
+    }
+
+    // Helper function for interconnecting two pillars with zig-zag bridges.
+    bool interconnect(const Pillar& pillar, const Pillar& nextpillar);
+
+    // For connecting a head to a nearby pillar.
+    bool connect_to_nearpillar(const Head& head, long nearpillar_id);
+    
+    // Find route for a head to the ground. Inserts additional bridge from the
+    // head to the pillar if cannot create pillar directly.
+    // The optional dir parameter is the direction of the bridge which is the
+    // direction of the pinhead if omitted.
+    bool connect_to_ground(Head& head, const Vec3d &dir);
+    inline bool connect_to_ground(Head& head);
+    
+    bool connect_to_model_body(Head &head);
+
+    bool search_pillar_and_connect(const Head& source);
+    
+    // This is a proxy function for pillar creation which will mind the gap
+    // between the pad and the model bottom in zero elevation mode.
+    // jp is the starting junction point which needs to be routed down.
+    // sourcedir is the allowed direction of an optional bridge between the
+    // jp junction and the final pillar.
+    bool create_ground_pillar(const Vec3d &jp,
+                              const Vec3d &sourcedir,
+                              double       radius,
+                              long         head_id = SupportTreeNode::ID_UNSET);
+
+    void add_pillar_base(long pid)
+    {
+        m_builder.add_pillar_base(pid, m_cfg.base_height_mm, m_cfg.base_radius_mm);
+    }
+
+    std::optional<DiffBridge> search_widening_path(const Vec3d &jp,
+                                                   const Vec3d &dir,
+                                                   double       radius,
+                                                   double       new_radius);
+
+public:
+    SupportTreeBuildsteps(SupportTreeBuilder & builder, const SupportableMesh &sm);
+
+    // Now let's define the individual steps of the support generation algorithm
+
+    // Filtering step: here we will discard inappropriate support points
+    // and decide the future of the appropriate ones. We will check if a
+    // pinhead is applicable and adjust its angle at each support point. We
+    // will also merge the support points that are just too close and can
+    // be considered as one.
+    void filter();
+
+    // Pinhead creation: based on the filtering results, the Head objects
+    // will be constructed (together with their triangle meshes).
+    void add_pinheads();
+
+    // Further classification of the support points with pinheads. If the
+    // ground is directly reachable through a vertical line parallel to the
+    // Z axis we consider a support point as pillar candidate. If touches
+    // the model geometry, it will be marked as non-ground facing and
+    // further steps will process it. Also, the pillars will be grouped
+    // into clusters that can be interconnected with bridges. Elements of
+    // these groups may or may not be interconnected. Here we only run the
+    // clustering algorithm.
+    void classify();
+
+    // Step: Routing the ground connected pinheads, and interconnecting
+    // them with additional (angled) bridges. Not all of these pinheads
+    // will be a full pillar (ground connected). Some will connect to a
+    // nearby pillar using a bridge. The max number of such side-heads for
+    // a central pillar is limited to avoid bad weight distribution.
+    void routing_to_ground();
+
+    // Step: routing the pinheads that would connect to the model surface
+    // along the Z axis downwards. For now these will actually be connected with
+    // the model surface with a flipped pinhead. In the future here we could use
+    // some smart algorithms to search for a safe path to the ground or to a
+    // nearby pillar that can hold the supported weight.
+    void routing_to_model();
+
+    void interconnect_pillars();
+
+    inline void merge_result() { m_builder.merged_mesh(); }
+
+    static bool execute(SupportTreeBuilder & builder, const SupportableMesh &sm);
+};
+
+}
+}
+
+#endif // SLASUPPORTTREEALGORITHM_H

src/libslic3r/SLA/SupportTreeMesher.cpp

@@ -0,0 +1,270 @@
+#include "SupportTreeMesher.hpp"
+
+namespace Slic3r { namespace sla {
+
+indexed_triangle_set sphere(double rho, Portion portion, double fa) {
+
+    indexed_triangle_set ret;
+
+    // prohibit close to zero radius
+    if(rho <= 1e-6 && rho >= -1e-6) return ret;
+
+    auto& vertices = ret.vertices;
+    auto& facets = ret.indices;
+
+    // Algorithm:
+    // Add points one-by-one to the sphere grid and form facets using relative
+    // coordinates. Sphere is composed effectively of a mesh of stacked circles.
+
+    // adjust via rounding to get an even multiple for any provided angle.
+    double angle = (2 * PI / floor(2*PI / fa) );
+
+    // Ring to be scaled to generate the steps of the sphere
+    std::vector<double> ring;
+
+    for (double i = 0; i < 2*PI; i+=angle) ring.emplace_back(i);
+
+    const auto sbegin = size_t(2*std::get<0>(portion)/angle);
+    const auto send = size_t(2*std::get<1>(portion)/angle);
+
+    const size_t steps = ring.size();
+    const double increment = 1.0 / double(steps);
+
+    // special case: first ring connects to 0,0,0
+    // insert and form facets.
+    if (sbegin == 0)
+        vertices.emplace_back(
+            Vec3f(0.f, 0.f, float(-rho + increment * sbegin * 2. * rho)));
+
+    auto id = coord_t(vertices.size());
+    for (size_t i = 0; i < ring.size(); i++) {
+        // Fixed scaling
+        const double z = -rho + increment*rho*2.0 * (sbegin + 1.0);
+        // radius of the circle for this step.
+        const double r = std::sqrt(std::abs(rho*rho - z*z));
+        Vec2d b = Eigen::Rotation2Dd(ring[i]) * Eigen::Vector2d(0, r);
+        vertices.emplace_back(Vec3d(b(0), b(1), z).cast<float>());
+
+        if (sbegin == 0)
+            (i == 0) ? facets.emplace_back(coord_t(ring.size()), 0, 1) :
+                       facets.emplace_back(id - 1, 0, id);
+        ++id;
+    }
+
+    // General case: insert and form facets for each step,
+    // joining it to the ring below it.
+    for (size_t s = sbegin + 2; s < send - 1; s++) {
+        const double z = -rho + increment * double(s * 2. * rho);
+        const double r = std::sqrt(std::abs(rho*rho - z*z));
+
+        for (size_t i = 0; i < ring.size(); i++) {
+            Vec2d b = Eigen::Rotation2Dd(ring[i]) * Eigen::Vector2d(0, r);
+            vertices.emplace_back(Vec3d(b(0), b(1), z).cast<float>());
+            auto id_ringsize = coord_t(id - int(ring.size()));
+            if (i == 0) {
+                // wrap around
+                facets.emplace_back(id - 1, id, id + coord_t(ring.size() - 1) );
+                facets.emplace_back(id - 1, id_ringsize, id);
+            } else {
+                facets.emplace_back(id_ringsize - 1, id_ringsize, id);
+                facets.emplace_back(id - 1, id_ringsize - 1, id);
+            }
+            id++;
+        }
+    }
+
+    // special case: last ring connects to 0,0,rho*2.0
+    // only form facets.
+    if(send >= size_t(2*PI / angle)) {
+        vertices.emplace_back(0.f, 0.f, float(-rho + increment*send*2.0*rho));
+        for (size_t i = 0; i < ring.size(); i++) {
+            auto id_ringsize = coord_t(id - int(ring.size()));
+            if (i == 0) {
+                // third vertex is on the other side of the ring.
+                facets.emplace_back(id - 1, id_ringsize, id);
+            } else {
+                auto ci = coord_t(id_ringsize + coord_t(i));
+                facets.emplace_back(ci - 1, ci, id);
+            }
+        }
+    }
+    id++;
+
+    return ret;
+}
+
+indexed_triangle_set cylinder(double r, double h, size_t ssteps, const Vec3d &sp)
+{
+    assert(ssteps > 0);
+
+    indexed_triangle_set ret;
+
+    auto steps = int(ssteps);
+    auto& points = ret.vertices;
+    auto& indices = ret.indices;
+    points.reserve(2*ssteps);
+    double a = 2*PI/steps;
+
+    Vec3d jp = sp;
+    Vec3d endp = {sp(X), sp(Y), sp(Z) + h};
+
+    // Upper circle points
+    for(int i = 0; i < steps; ++i) {
+        double phi = i*a;
+        auto ex = float(endp(X) + r*std::cos(phi));
+        auto ey = float(endp(Y) + r*std::sin(phi));
+        points.emplace_back(ex, ey, float(endp(Z)));
+    }
+
+    // Lower circle points
+    for(int i = 0; i < steps; ++i) {
+        double phi = i*a;
+        auto x = float(jp(X) + r*std::cos(phi));
+        auto y = float(jp(Y) + r*std::sin(phi));
+        points.emplace_back(x, y, float(jp(Z)));
+    }
+
+    // Now create long triangles connecting upper and lower circles
+    indices.reserve(2*ssteps);
+    auto offs = steps;
+    for(int i = 0; i < steps - 1; ++i) {
+        indices.emplace_back(i, i + offs, offs + i + 1);
+        indices.emplace_back(i, offs + i + 1, i + 1);
+    }
+
+    // Last triangle connecting the first and last vertices
+    auto last = steps - 1;
+    indices.emplace_back(0, last, offs);
+    indices.emplace_back(last, offs + last, offs);
+
+    // According to the slicing algorithms, we need to aid them with generating
+    // a watertight body. So we create a triangle fan for the upper and lower
+    // ending of the cylinder to close the geometry.
+    points.emplace_back(jp.cast<float>()); int ci = int(points.size() - 1);
+    for(int i = 0; i < steps - 1; ++i)
+        indices.emplace_back(i + offs + 1, i + offs, ci);
+
+    indices.emplace_back(offs, steps + offs - 1, ci);
+
+    points.emplace_back(endp.cast<float>()); ci = int(points.size() - 1);
+    for(int i = 0; i < steps - 1; ++i)
+        indices.emplace_back(ci, i, i + 1);
+
+    indices.emplace_back(steps - 1, 0, ci);
+
+    return ret;
+}
+
+indexed_triangle_set pinhead(double r_pin,
+                             double r_back,
+                             double length,
+                             size_t steps)
+{
+    assert(steps > 0);
+    assert(length >= 0.);
+    assert(r_back > 0.);
+    assert(r_pin > 0.);
+
+    indexed_triangle_set mesh;
+
+    // We create two spheres which will be connected with a robe that fits
+    // both circles perfectly.
+
+    // Set up the model detail level
+    const double detail = 2 * PI / steps;
+
+    // We don't generate whole circles. Instead, we generate only the
+    // portions which are visible (not covered by the robe) To know the
+    // exact portion of the bottom and top circles we need to use some
+    // rules of tangent circles from which we can derive (using simple
+    // triangles the following relations:
+
+    // The height of the whole mesh
+    const double h   = r_back + r_pin + length;
+    double       phi = PI / 2. - std::acos((r_back - r_pin) / h);
+
+    // To generate a whole circle we would pass a portion of (0, Pi)
+    // To generate only a half horizontal circle we can pass (0, Pi/2)
+    // The calculated phi is an offset to the half circles needed to smooth
+    // the transition from the circle to the robe geometry
+
+    auto &&s1 = sphere(r_back, make_portion(0, PI / 2 + phi), detail);
+    auto &&s2 = sphere(r_pin, make_portion(PI / 2 + phi, PI), detail);
+
+    for (auto &p : s2.vertices) p.z() += h;
+
+    its_merge(mesh, s1);
+    its_merge(mesh, s2);
+
+    for (size_t idx1 = s1.vertices.size() - steps, idx2 = s1.vertices.size();
+         idx1 < s1.vertices.size() - 1; idx1++, idx2++) {
+        coord_t i1s1 = coord_t(idx1), i1s2 = coord_t(idx2);
+        coord_t i2s1 = i1s1 + 1, i2s2 = i1s2 + 1;
+
+        mesh.indices.emplace_back(i1s1, i2s1, i2s2);
+        mesh.indices.emplace_back(i1s1, i2s2, i1s2);
+    }
+
+    auto i1s1 = coord_t(s1.vertices.size()) - coord_t(steps);
+    auto i2s1 = coord_t(s1.vertices.size()) - 1;
+    auto i1s2 = coord_t(s1.vertices.size());
+    auto i2s2 = coord_t(s1.vertices.size()) + coord_t(steps) - 1;
+
+    mesh.indices.emplace_back(i2s2, i2s1, i1s1);
+    mesh.indices.emplace_back(i1s2, i2s2, i1s1);
+
+    return mesh;
+}
+
+indexed_triangle_set halfcone(double       baseheight,
+                              double       r_bottom,
+                              double       r_top,
+                              const Vec3d &pos,
+                              size_t       steps)
+{
+    assert(steps > 0);
+
+    if (baseheight <= 0 || steps <= 0) return {};
+
+    indexed_triangle_set base;
+
+    double a    = 2 * PI / steps;
+    auto   last = int(steps - 1);
+    Vec3d  ep{pos.x(), pos.y(), pos.z() + baseheight};
+    for (size_t i = 0; i < steps; ++i) {
+        double phi = i * a;
+        auto x   = float(pos.x() + r_top * std::cos(phi));
+        auto y   = float(pos.y() + r_top * std::sin(phi));
+        base.vertices.emplace_back(x, y, float(ep.z()));
+    }
+
+    for (size_t i = 0; i < steps; ++i) {
+        double phi = i * a;
+        auto x   = float(pos.x() + r_bottom * std::cos(phi));
+        auto y   = float(pos.y() + r_bottom * std::sin(phi));
+        base.vertices.emplace_back(x, y, float(pos.z()));
+    }
+
+    base.vertices.emplace_back(pos.cast<float>());
+    base.vertices.emplace_back(ep.cast<float>());
+
+    auto &indices = base.indices;
+    auto  hcenter = int(base.vertices.size() - 1);
+    auto  lcenter = int(base.vertices.size() - 2);
+    auto  offs    = int(steps);
+    for (int i = 0; i < last; ++i) {
+        indices.emplace_back(i, i + offs, offs + i + 1);
+        indices.emplace_back(i, offs + i + 1, i + 1);
+        indices.emplace_back(i, i + 1, hcenter);
+        indices.emplace_back(lcenter, offs + i + 1, offs + i);
+    }
+
+    indices.emplace_back(0, last, offs);
+    indices.emplace_back(last, offs + last, offs);
+    indices.emplace_back(hcenter, last, 0);
+    indices.emplace_back(offs, offs + last, lcenter);
+
+    return base;
+}
+
+}} // namespace Slic3r::sla

src/libslic3r/SLA/SupportTreeMesher.hpp

@@ -0,0 +1,129 @@
+#ifndef SUPPORTTREEMESHER_HPP
+#define SUPPORTTREEMESHER_HPP
+
+#include "libslic3r/Point.hpp"
+
+#include "libslic3r/SLA/SupportTreeBuilder.hpp"
+#include "libslic3r/TriangleMesh.hpp"
+//#include "libslic3r/SLA/Contour3D.hpp"
+
+namespace Slic3r { namespace sla {
+
+using Portion = std::tuple<double, double>;
+
+inline Portion make_portion(double a, double b)
+{
+    return std::make_tuple(a, b);
+}
+
+indexed_triangle_set sphere(double  rho,
+                            Portion portion = make_portion(0., 2. * PI),
+                            double  fa      = (2. * PI / 360.));
+
+// Down facing cylinder in Z direction with arguments:
+// r: radius
+// h: Height
+// ssteps: how many edges will create the base circle
+// sp: starting point
+indexed_triangle_set cylinder(double       r,
+                              double       h,
+                              size_t       steps = 45,
+                              const Vec3d &sp    = Vec3d::Zero());
+
+indexed_triangle_set pinhead(double r_pin,
+                             double r_back,
+                             double length,
+                             size_t steps = 45);
+
+indexed_triangle_set halfcone(double       baseheight,
+                              double       r_bottom,
+                              double       r_top,
+                              const Vec3d &pt    = Vec3d::Zero(),
+                              size_t       steps = 45);
+
+inline indexed_triangle_set get_mesh(const Head &h, size_t steps)
+{
+    indexed_triangle_set mesh = pinhead(h.r_pin_mm, h.r_back_mm, h.width_mm, steps);
+
+    for (auto& p : mesh.vertices) p.z() -= (h.fullwidth() - h.r_back_mm);
+
+    using Quaternion = Eigen::Quaternion<float>;
+
+    // We rotate the head to the specified direction. The head's pointing
+    // side is facing upwards so this means that it would hold a support
+    // point with a normal pointing straight down. This is the reason of
+    // the -1 z coordinate
+    auto quatern = Quaternion::FromTwoVectors(Vec3f{0.f, 0.f, -1.f},
+                                              h.dir.cast<float>());
+
+    Vec3f pos = h.pos.cast<float>();
+    for (auto& p : mesh.vertices) p = quatern * p + pos;
+
+    return mesh;
+}
+
+inline indexed_triangle_set get_mesh(const Pillar &p, size_t steps)
+{
+    if(p.height > EPSILON) { // Endpoint is below the starting point
+        // We just create a bridge geometry with the pillar parameters and
+        // move the data.
+        return cylinder(p.r, p.height, steps, p.endpoint());
+    }
+
+    return {};
+}
+
+inline indexed_triangle_set get_mesh(const Pedestal &p, size_t steps)
+{
+    return halfcone(p.height, p.r_bottom, p.r_top, p.pos, steps);
+}
+
+inline indexed_triangle_set get_mesh(const Junction &j, size_t steps)
+{
+    indexed_triangle_set mesh = sphere(j.r, make_portion(0, PI), 2 *PI / steps);
+    Vec3f pos = j.pos.cast<float>();
+    for(auto& p : mesh.vertices) p += pos;
+    return mesh;
+}
+
+inline indexed_triangle_set get_mesh(const Bridge &br, size_t steps)
+{
+    using Quaternion = Eigen::Quaternion<float>;
+    Vec3d v = (br.endp - br.startp);
+    Vec3d dir = v.normalized();
+    double d = v.norm();
+
+    indexed_triangle_set mesh = cylinder(br.r, d, steps);
+
+    auto quater = Quaternion::FromTwoVectors(Vec3f{0.f, 0.f, 1.f},
+                                             dir.cast<float>());
+
+    Vec3f startp = br.startp.cast<float>();
+    for(auto& p : mesh.vertices) p = quater * p + startp;
+
+    return mesh;
+}
+
+inline indexed_triangle_set get_mesh(const DiffBridge &br, size_t steps)
+{
+    double h = br.get_length();
+    indexed_triangle_set mesh = halfcone(h, br.r, br.end_r, Vec3d::Zero(), steps);
+
+    using Quaternion = Eigen::Quaternion<float>;
+
+    // We rotate the head to the specified direction. The head's pointing
+    // side is facing upwards so this means that it would hold a support
+    // point with a normal pointing straight down. This is the reason of
+    // the -1 z coordinate
+    auto quatern = Quaternion::FromTwoVectors(Vec3f{0.f, 0.f, 1.f},
+                                              br.get_dir().cast<float>());
+
+    Vec3f startp = br.startp.cast<float>();
+    for(auto& p : mesh.vertices) p = quatern * p + startp;
+
+    return mesh;
+}
+
+}} // namespace Slic3r::sla
+
+#endif // SUPPORTTREEMESHER_HPP

src/libslic3r/SLA/bicubic.h

@@ -0,0 +1,186 @@
+#ifndef BICUBIC_HPP
+#define BICUBIC_HPP
+
+#include <algorithm>
+#include <vector>
+#include <cmath>
+
+#include <Eigen/Dense>
+
+namespace Slic3r {
+
+namespace BicubicInternal {
+	// Linear kernel, to be able to test cubic methods with hat kernels.
+	template<typename T>
+	struct LinearKernel
+	{
+		typedef T	FloatType;
+
+		static T a00() { return T(0.); }
+		static T a01() { return T(0.); }
+		static T a02() { return T(0.); }
+		static T a03() { return T(0.); }
+		static T a10() { return T(1.); }
+		static T a11() { return T(-1.); }
+		static T a12() { return T(0.); }
+		static T a13() { return T(0.); }
+		static T a20() { return T(0.); }
+		static T a21() { return T(1.); }
+		static T a22() { return T(0.); }
+		static T a23() { return T(0.); }
+		static T a30() { return T(0.); }
+		static T a31() { return T(0.); }
+		static T a32() { return T(0.); }
+		static T a33() { return T(0.); }
+	};
+
+	// Interpolation kernel aka Catmul-Rom aka Keyes kernel.
+	template<typename T>
+	struct CubicCatmulRomKernel
+	{
+		typedef T	FloatType;
+
+		static T a00() { return     0;     }
+		static T a01() { return (T)-0.5;   }
+		static T a02() { return (T) 1.;    }
+		static T a03() { return (T)-0.5;   }
+		static T a10() { return (T) 1.;    }
+		static T a11() { return     0;     }
+		static T a12() { return (T)-5./2.; }
+		static T a13() { return (T) 3./2.; }
+		static T a20() { return     0;     }
+		static T a21() { return (T) 0.5;   }
+		static T a22() { return (T) 2.;    }
+		static T a23() { return (T)-3./2.; }
+		static T a30() { return     0;     }
+		static T a31() { return     0;     }
+		static T a32() { return (T)-0.5;   }
+		static T a33() { return (T) 0.5;   }
+	};
+
+	// B-spline kernel
+	template<typename T>
+	struct CubicBSplineKernel
+	{
+		typedef T	FloatType;
+
+		static T a00() { return (T)  1./6.; }
+		static T a01() { return (T) -3./6.; }
+		static T a02() { return (T)  3./6.; }
+		static T a03() { return (T) -1./6.; }
+		static T a10() { return (T)  4./6.; }
+		static T a11() { return      0;     }
+		static T a12() { return (T) -6./6.; }
+		static T a13() { return (T)  3./6.; }
+		static T a20() { return (T)  1./6.; }
+		static T a21() { return (T)  3./6.; }
+		static T a22() { return (T)  3./6.; }
+		static T a23() { return (T)- 3./6.; }
+		static T a30() { return      0;     }
+		static T a31() { return      0;     }
+		static T a32() { return      0;     }
+		static T a33() { return (T)  1./6.; }
+	};
+
+	template<class T>
+	inline T clamp(T a, T lower, T upper)
+	{
+		return (a < lower) ? lower : 
+	   		   (a > upper) ? upper : a;
+	}
+}
+
+template<typename KERNEL>
+struct CubicKernel
+{
+	typedef typename KERNEL					KernelInternal;
+	typedef typename KERNEL::FloatType		FloatType;
+
+	static FloatType kernel(FloatType x)
+	{
+		x = fabs(x);
+		if (x >= (FloatType)2.)
+			return 0.0f;
+		if (x <= (FloatType)1.) {
+			FloatType x2 = x * x;
+			FloatType x3 = x2 * x;
+			return KERNEL::a10() + KERNEL::a11() * x + KERNEL::a12() * x2 + KERNEL::a13() * x3;
+		}
+		assert(x > (FloatType)1. && x < (FloatType)2.);
+		x -= (FloatType)1.;
+		FloatType x2 = x * x;
+		FloatType x3 = x2 * x;
+		return KERNEL::a00() + KERNEL::a01() * x + KERNEL::a02() * x2 + KERNEL::a03() * x3;
+	}
+
+	static FloatType interpolate(FloatType f0, FloatType f1, FloatType f2, FloatType f3, FloatType x)
+	{
+		const FloatType  x2 = x*x;
+		const FloatType  x3 = x*x*x;
+		return f0*(KERNEL::a00() + KERNEL::a01() * x + KERNEL::a02() * x2 + KERNEL::a03() * x3) +
+			   f1*(KERNEL::a10() + KERNEL::a11() * x + KERNEL::a12() * x2 + KERNEL::a13() * x3) +
+			   f2*(KERNEL::a20() + KERNEL::a21() * x + KERNEL::a22() * x2 + KERNEL::a23() * x3) + 
+			   f3*(KERNEL::a30() + KERNEL::a31() * x + KERNEL::a32() * x2 + KERNEL::a33() * x3);
+	}
+};
+
+// Linear splines
+typedef CubicKernel<BicubicInternal::LinearKernel<float>>					LinearKernelf;
+typedef CubicKernel<BicubicInternal::LinearKernel<double>>					LinearKerneld;
+// Catmul-Rom splines
+typedef CubicKernel<BicubicInternal::CubicCatmulRomKernel<float>>			CubicCatmulRomKernelf;
+typedef CubicKernel<BicubicInternal::CubicCatmulRomKernel<double>>			CubicCatmulRomKerneld;
+typedef CubicKernel<BicubicInternal::CubicCatmulRomKernel<float>>			CubicInterpolationKernelf;
+typedef CubicKernel<BicubicInternal::CubicCatmulRomKernel<double>>			CubicInterpolationKerneld;
+// Cubic B-splines
+typedef CubicKernel<BicubicInternal::CubicBSplineKernel<float>>				CubicBSplineKernelf;
+typedef CubicKernel<BicubicInternal::CubicBSplineKernel<double>>			CubicBSplineKerneld;
+
+template<typename KERNEL, typename Derived>
+static float cubic_interpolate(const Eigen::ArrayBase<Derived> &F, const typename KERNEL::FloatType pt, const typename KERNEL::FloatType dx)
+{
+	typedef typename KERNEL::FloatType T;
+	const int w  = int(F.size());
+	const int ix = (int)floor(pt);
+	const T   s  = pt - (T)ix;
+
+	if (ix > 1 && ix + 2 < w) {
+		// Inside the fully interpolated region.
+		return KERNEL::interpolate(F[ix - 1], F[ix], F[ix + 1], F[ix + 2], s);
+	}
+	// Transition region. Extend with a constant function.
+	auto f = [&F, w](x) { return F[BicubicInternal::clamp(x, 0, w - 1)]; }
+	return KERNEL::interpolate(f(ix - 1), f(ix), f(ix + 1), f(ix + 2), s);
+}
+
+template<typename KERNEL, typename Derived>
+static float bicubic_interpolate(const Eigen::MatrixBase<Derived> &F, const Eigen::Matrix<typename KERNEL::FloatType, 2, 1, Eigen::DontAlign> &pt, const typename KERNEL::FloatType dx)
+{
+	typedef typename KERNEL::FloatType T;
+	const int w  = F.cols();
+	const int h  = F.rows();
+	const int ix = (int)floor(pt[0]);
+	const int iy = (int)floor(pt[1]);
+	const T   s  = pt[0] - (T)ix;
+	const T   t  = pt[1] - (T)iy;
+
+	if (ix > 1 && ix + 2 < w && iy > 1 && iy + 2 < h) {
+		// Inside the fully interpolated region.
+		return KERNEL::interpolate(
+			KERNEL::interpolate(F(ix-1,iy-1),F(ix  ,iy-1),F(ix+1,iy-1),F(ix+2,iy-1),s),
+			KERNEL::interpolate(F(ix-1,iy  ),F(ix  ,iy  ),F(ix+1,iy  ),F(ix+2,iy  ),s),
+			KERNEL::interpolate(F(ix-1,iy+1),F(ix  ,iy+1),F(ix+1,iy+1),F(ix+2,iy+1),s),
+			KERNEL::interpolate(F(ix-1,iy+2),F(ix  ,iy+2),F(ix+1,iy+2),F(ix+2,iy+2),s),t);
+	}
+	// Transition region. Extend with a constant function.
+	auto f = [&f, w, h](int x, int y) { return F(BicubicInternal::clamp(x,0,w-1),BicubicInternal::clamp(y,0,h-1)); }
+	return KERNEL::interpolate(
+		KERNEL::interpolate(f(ix-1,iy-1),f(ix  ,iy-1),f(ix+1,iy-1),f(ix+2,iy-1),s),
+		KERNEL::interpolate(f(ix-1,iy  ),f(ix  ,iy  ),f(ix+1,iy  ),f(ix+2,iy  ),s),
+		KERNEL::interpolate(f(ix-1,iy+1),f(ix  ,iy+1),f(ix+1,iy+1),f(ix+2,iy+1),s),
+		KERNEL::interpolate(f(ix-1,iy+2),f(ix  ,iy+2),f(ix+1,iy+2),f(ix+2,iy+2),s),t);
+}
+
+} // namespace Slic3r
+
+#endif /* BICUBIC_HPP */