diff --git a/src/libslic3r/BranchingTree/BranchingTree.cpp b/src/libslic3r/BranchingTree/BranchingTree.cpp
new file mode 100644
index 000000000..6463a30de
--- /dev/null
+++ b/src/libslic3r/BranchingTree/BranchingTree.cpp
@@ -0,0 +1,194 @@
+#include "BranchingTree.hpp"
+#include "PointCloud.hpp"
+
+#include <numeric>
+#include <optional>
+#include <algorithm>
+
+#include "libslic3r/SLA/SupportTreeUtils.hpp"
+
+namespace Slic3r { namespace branchingtree {
+
+void build_tree(PointCloud &nodes, Builder &builder)
+{
+    constexpr size_t initK = 5;
+
+    auto ptsqueue = nodes.start_queue();
+    auto &properties = nodes.properties();
+
+    struct NodeDistance
+    {
+        size_t node_id       = Node::ID_NONE;
+        float  dst_branching = NaNf;
+        float  dst_euql      = NaNf;
+    };
+    auto distances = reserve_vector<NodeDistance>(initK);
+    double prev_dist_max = 0.;
+    size_t K = initK;
+    bool   routed = true;
+    size_t node_id = Node::ID_NONE;
+
+    while ((!ptsqueue.empty() && builder.is_valid()) || !routed) {
+        if (routed) {
+            node_id = ptsqueue.top();
+            ptsqueue.pop();
+        }
+
+        Node node = nodes.get(node_id);
+        nodes.mark_unreachable(node_id);
+
+        distances.clear();
+        distances.reserve(K);
+
+        float dmax = 0.;
+        nodes.foreach_reachable(
+            node.pos,
+            [&distances, &dmax](size_t id, float dst_branching, float dst_euql) {
+                distances.emplace_back(NodeDistance{id, dst_branching, dst_euql});
+                dmax = std::max(dmax, dst_euql);
+            }, K, prev_dist_max);
+
+        std::sort(distances.begin(), distances.end(),
+                  [](auto &a, auto &b) { return a.dst_branching < b.dst_branching; });
+
+        if (distances.empty()) {
+            builder.report_unroutable(node);
+            K = initK;
+            prev_dist_max = 0.;
+            routed = true;
+
+            continue;
+        }
+
+        prev_dist_max = dmax;
+        K *= 2;
+
+        auto closest_it = distances.begin();
+        routed = false;
+        while (closest_it != distances.end() && !routed && builder.is_valid()) {
+            size_t closest_node_id = closest_it->node_id;
+            Node closest_node = nodes.get(closest_node_id);
+
+            auto type = nodes.get_type(closest_node_id);
+            float w = nodes.get(node_id).weight + closest_it->dst_branching;
+            closest_node.Rmin = std::max(node.Rmin, closest_node.Rmin);
+
+            switch (type) {
+            case BED: {
+                closest_node.weight = w;
+                if (closest_it->dst_branching > nodes.properties().max_branch_length()) {
+                    auto hl_br_len = float(nodes.properties().max_branch_length()) / 2.f;
+                    Node new_node {{node.pos.x(), node.pos.y(), node.pos.z() - hl_br_len}, node.Rmin};
+                    new_node.id = int(nodes.next_junction_id());
+                    new_node.weight = nodes.get(node_id).weight + hl_br_len;
+                    new_node.left = node.id;
+                    if ((routed = builder.add_bridge(node, new_node))) {
+                        size_t new_idx = nodes.insert_junction(new_node);
+                        ptsqueue.push(new_idx);
+                    }
+                }
+                else if ((routed = builder.add_ground_bridge(node, closest_node))) {
+                    closest_node.left = closest_node.right = node_id;
+                    nodes.get(closest_node_id) = closest_node;
+                    nodes.mark_unreachable(closest_node_id);
+                }
+
+                break;
+            }
+            case MESH: {
+                closest_node.weight = w;
+                if ((routed = builder.add_mesh_bridge(node, closest_node))) {
+                    closest_node.left = closest_node.right = node_id;
+                    nodes.get(closest_node_id) = closest_node;
+                    nodes.mark_unreachable(closest_node_id);
+                }
+
+                break;
+            }
+            case LEAF:
+            case JUNCTION: {
+                auto max_slope = float(properties.max_slope());
+
+                if (auto mergept = find_merge_pt(node.pos, closest_node.pos, max_slope)) {
+
+                    float mergedist_closest = (*mergept - closest_node.pos).norm();
+                    float mergedist_node = (*mergept - node.pos).norm();
+                    float Wnode = nodes.get(node_id).weight;
+                    float Wclosest = nodes.get(closest_node_id).weight;
+                    float Wsum = std::max(Wnode, Wclosest);
+                    float distsum = std::max(mergedist_closest, mergedist_node);
+                    w = Wsum + distsum;
+
+                    if (mergedist_closest > EPSILON && mergedist_node > EPSILON) {
+                        Node mergenode{*mergept, closest_node.Rmin};
+                        mergenode.weight = w;
+                        mergenode.id = int(nodes.next_junction_id());
+
+                        if ((routed = builder.add_merger(node, closest_node, mergenode))) {
+                            mergenode.left = node_id;
+                            mergenode.right = closest_node_id;
+                            size_t new_idx = nodes.insert_junction(mergenode);
+                            ptsqueue.push(new_idx);
+                            size_t qid = nodes.get_queue_idx(closest_node_id);
+
+                            if (qid != PointCloud::Unqueued)
+                                ptsqueue.remove(nodes.get_queue_idx(closest_node_id));
+
+                            nodes.mark_unreachable(closest_node_id);
+                        }
+                    } else if (closest_node.pos.z() < node.pos.z() &&
+                               (closest_node.left == Node::ID_NONE ||
+                                closest_node.right == Node::ID_NONE)) {
+                        closest_node.weight = w;
+                        if ((routed = builder.add_bridge(node, closest_node))) {
+                            if (closest_node.left == Node::ID_NONE)
+                                closest_node.left = node_id;
+                            else if (closest_node.right == Node::ID_NONE)
+                                closest_node.right = node_id;
+
+                            nodes.get(closest_node_id) = closest_node;
+                        }
+                    }
+                }
+
+                break;
+            }
+            case NONE:;
+            }
+
+            ++closest_it;
+        }
+
+        if (routed) {
+            prev_dist_max = 0.;
+            K = initK;
+        }
+    }
+}
+
+void build_tree(const indexed_triangle_set &its,
+                const std::vector<Node>    &support_roots,
+                Builder                    &builder,
+                const Properties           &properties)
+{
+    PointCloud nodes(its, support_roots, properties);
+
+    build_tree(nodes, builder);
+}
+
+ExPolygon make_bed_poly(const indexed_triangle_set &its)
+{
+    auto bb = bounding_box(its);
+
+    BoundingBox bbcrd{scaled(to_2d(bb.min)), scaled(to_2d(bb.max))};
+    bbcrd.offset(scaled(10.));
+    Point     min = bbcrd.min, max = bbcrd.max;
+    ExPolygon ret = {{min.x(), min.y()},
+                     {max.x(), min.y()},
+                     {max.x(), max.y()},
+                     {min.x(), max.y()}};
+
+    return ret;
+}
+
+}} // namespace Slic3r::branchingtree
diff --git a/src/libslic3r/BranchingTree/BranchingTree.hpp b/src/libslic3r/BranchingTree/BranchingTree.hpp
new file mode 100644
index 000000000..c88410b3a
--- /dev/null
+++ b/src/libslic3r/BranchingTree/BranchingTree.hpp
@@ -0,0 +1,143 @@
+#ifndef SUPPORTTREEBRANCHING_HPP
+#define SUPPORTTREEBRANCHING_HPP
+
+// For indexed_triangle_set
+#include <admesh/stl.h>
+
+#include "libslic3r/ExPolygon.hpp"
+#include "libslic3r/BoundingBox.hpp"
+
+namespace Slic3r { namespace branchingtree {
+
+// Branching tree input parameters. This is an in-line fillable structure with
+// setters returning self references.
+class Properties
+{
+    double m_max_slope       = PI / 4.;
+    double m_ground_level    = 0.;
+    double m_sampling_radius = .5;
+    double m_max_branch_len  = 10.;
+
+    ExPolygons m_bed_shape;
+
+public:
+    // Maximum slope for bridges of the tree
+    Properties &max_slope(double val) noexcept
+    {
+        m_max_slope = val;
+        return *this;
+    }
+    // Z level of the ground
+    Properties &ground_level(double val) noexcept
+    {
+        m_ground_level = val;
+        return *this;
+    }
+    // How far should sample points be in the mesh and the ground
+    Properties &sampling_radius(double val) noexcept
+    {
+        m_sampling_radius = val;
+        return *this;
+    }
+    // Shape of the print bed (ground)
+    Properties &bed_shape(ExPolygons bed) noexcept
+    {
+        m_bed_shape = std::move(bed);
+        return *this;
+    }
+
+    Properties &max_branch_length(double val) noexcept
+    {
+        m_max_branch_len = val;
+        return *this;
+    }
+
+    double max_slope() const noexcept { return m_max_slope; }
+    double ground_level() const noexcept { return m_ground_level; }
+    double sampling_radius() const noexcept { return m_sampling_radius; }
+    double max_branch_length() const noexcept { return m_max_branch_len; }
+    const ExPolygons &bed_shape() const noexcept { return m_bed_shape; }
+};
+
+// A junction of the branching tree with position and radius.
+struct Node
+{
+    static constexpr int ID_NONE = -1;
+
+    int id = ID_NONE, left = ID_NONE, right = ID_NONE;
+
+    Vec3f pos;
+    float Rmin = 0.f;
+
+    // Tracking the weight of each junction, which is essentially the sum of
+    // the lenghts of all branches emanating from this junction.
+    float weight = 0.f;
+
+    Node(const Vec3f &p, float r_min = .0f) : pos{p}, Rmin{r_min}, weight{0.f}
+    {}
+};
+
+// An output interface for the branching tree generator function. Consider each
+// method as a callback and implement the actions that need to be done.
+class Builder
+{
+public:
+    virtual ~Builder() = default;
+
+    // A simple bridge from junction to junction.
+    virtual bool add_bridge(const Node &from, const Node &to) = 0;
+
+    // An Y shaped structure with two starting points and a merge point below
+    // them. The angles will respect the max_slope setting.
+    virtual bool add_merger(const Node &node,
+                            const Node &closest,
+                            const Node &merge_node) = 0;
+
+    // Add an anchor bridge to the ground (print bed)
+    virtual bool add_ground_bridge(const Node &from,
+                                   const Node &to) = 0;
+
+    // Add an anchor bridge to the model body
+    virtual bool add_mesh_bridge(const Node &from, const Node &to) = 0;
+
+    // Report nodes that can not be routed to an endpoint (model or ground)
+    virtual void report_unroutable(const Node &j) = 0;
+
+    // If returns false, the tree building process shall stop
+    virtual bool is_valid() const { return true; }
+};
+
+// Build the actual tree.
+// its:           The input mesh
+// support_leafs: The input support points
+// builder:       The output interface, describes how to build the tree
+// properties:    Parameters of the tree
+//
+// Notes:
+// The original algorithm implicitly ensures that the generated tree avoids
+// the model body. This implementation uses point sampling of the mesh thus an
+// explicit check is needed if the part of the tree being inserted properly
+// avoids the model. This can be done in the builder implementation. Each
+// method can return a boolean indicating whether the given branch can or
+// cannot be inserted. If a particular path is unavailable, the algorithm
+// will try a few other paths as well. If all of them fail, one of the
+// report_unroutable_* methods will be called as a last resort.
+void build_tree(const indexed_triangle_set &its,
+                const std::vector<Node>    &support_leafs,
+                Builder                    &builder,
+                const Properties           &properties = {});
+
+inline void build_tree(const indexed_triangle_set &its,
+                       const std::vector<Node>    &support_leafs,
+                       Builder                   &&builder,
+                       const Properties           &properties = {})
+{
+    build_tree(its, support_leafs, builder, properties);
+}
+
+// Helper function to derive a bed polygon only from the model bounding box.
+ExPolygon make_bed_poly(const indexed_triangle_set &its);
+
+}} // namespace Slic3r::branchingtree
+
+#endif // SUPPORTTREEBRANCHING_HPP
diff --git a/src/libslic3r/BranchingTree/PointCloud.cpp b/src/libslic3r/BranchingTree/PointCloud.cpp
new file mode 100644
index 000000000..f1d7ae521
--- /dev/null
+++ b/src/libslic3r/BranchingTree/PointCloud.cpp
@@ -0,0 +1,252 @@
+#include "PointCloud.hpp"
+
+#include "libslic3r/Geometry.hpp"
+#include "libslic3r/Tesselate.hpp"
+
+#include <igl/random_points_on_mesh.h>
+
+namespace Slic3r { namespace branchingtree {
+
+std::optional<Vec3f> find_merge_pt(const Vec3f &A,
+                                   const Vec3f &B,
+                                   float        critical_angle)
+{
+    // The idea is that A and B both have their support cones. But searching
+    // for the intersection of these support cones is difficult and its enough
+    // to reduce this problem to 2D and search for the intersection of two
+    // rays that merge somewhere between A and B. The 2D plane is a vertical
+    // slice of the 3D scene where the 2D Y axis is equal to the 3D Z axis and
+    // the 2D X axis is determined by the XY direction of the AB vector.
+    //
+    // Z^
+    //  |    A *
+    //  |     . .   B *
+    //  |    .   .   . .
+    //  |   .     . .   .
+    //  |  .       x     .
+    //  -------------------> XY
+
+    // Determine the transformation matrix for the 2D projection:
+    Vec3f diff = {B.x() - A.x(), B.y() - A.y(), 0.f};
+    Vec3f dir  = diff.normalized(); // TODO: avoid normalization
+
+    Eigen::Matrix<float, 2, 3> tr2D;
+    tr2D.row(0) = Vec3f{dir.x(), dir.y(), dir.z()};
+    tr2D.row(1) = Vec3f{0.f, 0.f, 1.f};
+
+    // Transform the 2 vectors A and B into 2D vector 'a' and 'b'. Here we can
+    // omit 'a', pretend that its the origin and use BA as the vector b.
+    Vec2f b = tr2D * (B - A);
+
+    // Get the square sine of the ray emanating from 'a' towards 'b'. This ray might
+    // exceed the allowed angle but that is corrected subsequently.
+    // The sign of the original sine is also needed, hence b.y is multiplied by
+    // abs(b.y)
+    float b_sqn = b.squaredNorm();
+    float sin2sig_a = b_sqn > EPSILON ? (b.y() * std::abs(b.y())) / b_sqn : 0.f;
+
+    // sine2 from 'b' to 'a' is the opposite of sine2 from a to b
+    float sin2sig_b = -sin2sig_a;
+
+    // Derive the allowed angles from the given critical angle.
+    // critical_angle is measured from the horizontal X axis.
+    // The rays need to go downwards which corresponds to negative angles
+
+    float sincrit = std::sin(critical_angle);  // sine of the critical angle
+    float sin2crit = -sincrit * sincrit;       // signed sine squared
+    sin2sig_a = std::min(sin2sig_a, sin2crit); // Do the angle saturation of both rays
+    sin2sig_b = std::min(sin2sig_b, sin2crit); //
+    float sin2_a = std::abs(sin2sig_a);        // Get cosine squared values
+    float sin2_b = std::abs(sin2sig_b);
+    float cos2_a = 1.f - sin2_a;
+    float cos2_b = 1.f - sin2_b;
+
+    // Derive the new direction vectors. This is by square rooting the sin2
+    // and cos2 values and restoring the original signs
+    Vec2f Da = {std::copysign(std::sqrt(cos2_a), b.x()), std::copysign(std::sqrt(sin2_a), sin2sig_a)};
+    Vec2f Db = {-std::copysign(std::sqrt(cos2_b), b.x()), std::copysign(std::sqrt(sin2_b), sin2sig_b)};
+
+    // Determine where two rays ([0, 0], Da), (b, Db) intersect.
+    // Based on
+    // https://stackoverflow.com/questions/27459080/given-two-points-and-two-direction-vectors-find-the-point-where-they-intersect
+    // One ray is emanating from (0, 0) so the formula is simplified
+    double t1 = (Db.y() * b.x() - b.y() * Db.x()) /
+                (Da.x() * Db.y() - Da.y() * Db.x());
+
+    Vec2f mp = t1 * Da;
+    Vec3f Mp = A + tr2D.transpose() * mp;
+
+    return t1 >= 0.f ? Mp : Vec3f{};
+}
+
+void to_eigen_mesh(const indexed_triangle_set &its,
+                   Eigen::MatrixXd            &V,
+                   Eigen::MatrixXi            &F)
+{
+    V.resize(its.vertices.size(), 3);
+    F.resize(its.indices.size(), 3);
+    for (unsigned int i = 0; i < its.indices.size(); ++i)
+        F.row(i) = its.indices[i];
+
+    for (unsigned int i = 0; i < its.vertices.size(); ++i)
+        V.row(i) = its.vertices[i].cast<double>();
+}
+
+std::vector<Node> sample_mesh(const indexed_triangle_set &its, double radius)
+{
+    std::vector<Node> ret;
+
+    double surface_area = 0.;
+    for (const Vec3i &face : its.indices) {
+        std::array<Vec3f, 3> tri = {its.vertices[face(0)],
+                                    its.vertices[face(1)],
+                                    its.vertices[face(2)]};
+
+        auto U = tri[1] - tri[0], V = tri[2] - tri[0];
+        surface_area += 0.5 * U.cross(V).norm();
+    }
+
+    int N = surface_area / (PI * radius * radius);
+
+    Eigen::MatrixXd B;
+    Eigen::MatrixXi FI;
+    Eigen::MatrixXd V;
+    Eigen::MatrixXi F;
+    to_eigen_mesh(its, V, F);
+    igl::random_points_on_mesh(N, V, F, B, FI);
+
+    ret.reserve(size_t(N));
+    for (int i = 0; i < FI.size(); i++) {
+        int face_id = FI(i);
+
+        if (face_id < 0 || face_id >= int(its.indices.size()))
+            continue;
+
+        Vec3i face = its.indices[face_id];
+
+        if (face(0) >= int(its.vertices.size()) ||
+            face(1) >= int(its.vertices.size()) ||
+            face(2) >= int(its.vertices.size()))
+            continue;
+
+        Vec3f c = B.row(i)(0) * its.vertices[face(0)] +
+                  B.row(i)(1) * its.vertices[face(1)] +
+                  B.row(i)(2) * its.vertices[face(2)];
+
+        ret.emplace_back(c);
+    }
+
+    return ret;
+}
+
+std::vector<Node> sample_bed(const ExPolygons &bed, float z, double radius)
+{
+    std::vector<Vec3f> ret;
+
+    auto triangles = triangulate_expolygons_3d(bed, z);
+    indexed_triangle_set its;
+    its.vertices.reserve(triangles.size());
+
+    for (size_t i = 0; i < triangles.size(); i += 3) {
+        its.vertices.emplace_back(triangles[i].cast<float>());
+        its.vertices.emplace_back(triangles[i + 1].cast<float>());
+        its.vertices.emplace_back(triangles[i + 2].cast<float>());
+
+        its.indices.emplace_back(i, i + 1, i + 2);
+    }
+
+    return sample_mesh(its, radius);
+}
+
+PointCloud::PointCloud(const indexed_triangle_set &M,
+                       std::vector<Node>           support_leafs,
+                       const Properties           &props)
+    : PointCloud{sample_mesh(M, props.sampling_radius()),
+                 sample_bed(props.bed_shape(),
+                            props.ground_level(),
+                            props.sampling_radius()),
+                 std::move(support_leafs), props}
+{}
+
+PointCloud::PointCloud(std::vector<Node> meshpts,
+                       std::vector<Node> bedpts,
+                       std::vector<Node> support_leafs,
+                       const Properties &props)
+    : m_leafs{std::move(support_leafs)}
+    , m_meshpoints{std::move(meshpts)}
+    , m_bedpoints{std::move(bedpts)}
+    , m_props{props}
+    , cos2bridge_slope{std::cos(props.max_slope()) *
+                       std::abs(std::cos(props.max_slope()))}
+    , MESHPTS_BEGIN{m_bedpoints.size()}
+    , LEAFS_BEGIN{MESHPTS_BEGIN + m_meshpoints.size()}
+    , JUNCTIONS_BEGIN{LEAFS_BEGIN + m_leafs.size()}
+    , m_searchable_indices(JUNCTIONS_BEGIN + m_junctions.size(), true)
+    , m_queue_indices(JUNCTIONS_BEGIN + m_junctions.size(), Unqueued)
+    , m_reachable_cnt{JUNCTIONS_BEGIN + m_junctions.size()}
+{
+    for (size_t i = 0; i < m_bedpoints.size(); ++i) {
+        m_bedpoints[i].id = int(i);
+        m_ktree.insert({m_bedpoints[i].pos, i});
+    }
+    
+    for (size_t i = 0; i < m_meshpoints.size(); ++i) {
+        Node &n = m_meshpoints[i];
+        n.id = int(MESHPTS_BEGIN + i);
+        m_ktree.insert({n.pos, n.id});
+    }
+    
+    for (size_t i = 0; i < m_leafs.size(); ++i) {
+        Node &n = m_leafs[i];
+        n.id = int(LEAFS_BEGIN + i);
+        m_ktree.insert({n.pos, n.id});
+    }
+}
+
+float PointCloud::get_distance(const Vec3f &p, size_t node_id) const
+{
+    auto t = get_type(node_id);
+    auto ret = std::numeric_limits<float>::infinity();
+    const auto &node = get(node_id);
+    
+    switch (t) {
+    case MESH:
+    case BED: {
+        // Points of mesh or bed which are outside of the support cone of
+        // 'pos' must be discarded.
+        if (is_outside_support_cone(p, node.pos))
+            ret = std::numeric_limits<float>::infinity();
+        else
+            ret  = (node.pos - p).norm();
+        
+        break;
+    }
+    case LEAF:
+    case JUNCTION:{
+        auto mergept = find_merge_pt(p, node.pos, m_props.max_slope());
+        double maxL2 = m_props.max_branch_length() * m_props.max_branch_length();
+
+        if (!mergept || mergept->z() < (m_props.ground_level() + 2 * node.Rmin))
+            ret = std::numeric_limits<float>::infinity();
+        else if (double a = (node.pos - *mergept).squaredNorm(),
+                 b        = (p - *mergept).squaredNorm();
+                 a < maxL2 && b < maxL2)
+            ret = std::sqrt(b);
+
+        break;
+    }
+    case NONE:
+        ;
+    }
+    
+    // Setting the ret val to infinity will effectively discard this
+    // connection of nodes. max_branch_length property is used here
+    // to discard node=>node and node=>mesh connections longer than this
+    // property.
+    if (t != BED && ret > m_props.max_branch_length())
+        ret = std::numeric_limits<float>::infinity();
+    
+    return ret;
+}
+
+}} // namespace Slic3r::branchingtree
diff --git a/src/libslic3r/BranchingTree/PointCloud.hpp b/src/libslic3r/BranchingTree/PointCloud.hpp
new file mode 100644
index 000000000..f99b17990
--- /dev/null
+++ b/src/libslic3r/BranchingTree/PointCloud.hpp
@@ -0,0 +1,272 @@
+#ifndef POINTCLOUD_HPP
+#define POINTCLOUD_HPP
+
+#include <optional>
+
+#include "BranchingTree.hpp"
+
+#include "libslic3r/Execution/Execution.hpp"
+#include "libslic3r/MutablePriorityQueue.hpp"
+
+#include "libslic3r/BoostAdapter.hpp"
+#include "boost/geometry/index/rtree.hpp"
+
+namespace Slic3r { namespace branchingtree {
+
+std::optional<Vec3f> find_merge_pt(const Vec3f &A,
+                                   const Vec3f &B,
+                                   float        max_slope);
+
+void to_eigen_mesh(const indexed_triangle_set &its,
+                   Eigen::MatrixXd            &V,
+                   Eigen::MatrixXi            &F);
+
+std::vector<Node> sample_mesh(const indexed_triangle_set &its, double radius);
+
+std::vector<Node> sample_bed(const ExPolygons &bed,
+                                 float             z,
+                                 double            radius = 1.);
+
+enum PtType { LEAF, MESH, BED, JUNCTION, NONE };
+
+inline BoundingBox3Base<Vec3f> get_support_cone_bb(const Vec3f &p, const Properties &props)
+{
+    double gnd = props.ground_level() - EPSILON;
+    double h   = p.z() - gnd;
+    double phi = PI / 2 - props.max_slope();
+    auto   r   = float(std::min(h * std::tan(phi), props.max_branch_length() * std::sin(phi)));
+
+    Vec3f bb_min = {p.x() - r, p.y() - r, float(gnd)};
+    Vec3f bb_max = {p.x() + r, p.y() + r, p.z()};
+
+    return {bb_min, bb_max};
+}
+
+// A cloud of points including support points, mesh points, junction points
+// and anchor points on the bed. Junction points can be added or removed, all
+// the other point types are established on creation and remain unchangeable.
+class PointCloud {
+    std::vector<Node> m_leafs, m_junctions, m_meshpoints, m_bedpoints;
+
+    const branchingtree::Properties &m_props;
+
+    const double cos2bridge_slope;
+    const size_t MESHPTS_BEGIN, LEAFS_BEGIN, JUNCTIONS_BEGIN;
+
+private:
+
+    // These vectors have the same size as there are indices for nodes to keep
+    // access complexity constant. WARN: there might be cache non-locality costs
+    std::vector<bool> m_searchable_indices; // searchable flag value of a node
+    std::vector<size_t> m_queue_indices;    // queue id of a node if queued
+
+    size_t m_reachable_cnt;
+
+    struct CoordFn
+    {
+        const PointCloud *self;
+        CoordFn(const PointCloud *s) : self{s} {}
+        float operator()(size_t nodeid, size_t dim) const
+        {
+            return self->get(nodeid).pos(int(dim));
+        }
+    };
+
+    using PointIndexEl = std::pair<Vec3f, unsigned>;
+
+    boost::geometry::index::
+        rtree<PointIndexEl, boost::geometry::index::rstar<16, 4> /* ? */>
+            m_ktree;
+
+    bool is_outside_support_cone(const Vec3f &supp, const Vec3f &pt) const
+    {
+        Vec3d D = (pt - supp).cast<double>();
+        double dot_sq = -D.z() * std::abs(-D.z());
+
+        return dot_sq < D.squaredNorm() * cos2bridge_slope;
+    }
+
+    template<class PC>
+    static auto *get_node(PC &&pc, size_t id)
+    {
+        auto *ret = decltype(pc.m_bedpoints.data())(nullptr);
+
+        switch(pc.get_type(id)) {
+        case BED: ret = &pc.m_bedpoints[id]; break;
+        case MESH: ret = &pc.m_meshpoints[id - pc.MESHPTS_BEGIN]; break;
+        case LEAF: ret = &pc.m_leafs [id - pc.LEAFS_BEGIN]; break;
+        case JUNCTION: ret = &pc.m_junctions[id - pc.JUNCTIONS_BEGIN]; break;
+        case NONE: assert(false);
+        }
+
+        return ret;
+    }
+
+public:
+
+    static constexpr auto Unqueued = size_t(-1);
+
+    struct ZCompareFn
+    {
+        const PointCloud *self;
+        ZCompareFn(const PointCloud *s) : self{s} {}
+        bool operator()(size_t node_a, size_t node_b) const
+        {
+            return self->get(node_a).pos.z() > self->get(node_b).pos.z();
+        }
+    };
+
+    PointCloud(const indexed_triangle_set &M,
+               std::vector<Node>           support_leafs,
+               const Properties           &props);
+
+    PointCloud(std::vector<Node> meshpts,
+               std::vector<Node> bedpts,
+               std::vector<Node> support_leafs,
+               const Properties &props);
+
+    PtType get_type(size_t node_id) const
+    {
+        PtType ret = NONE;
+
+        if (node_id < MESHPTS_BEGIN && !m_bedpoints.empty()) ret = BED;
+        else if (node_id < LEAFS_BEGIN && !m_meshpoints.empty()) ret = MESH;
+        else if (node_id < JUNCTIONS_BEGIN && !m_leafs.empty()) ret = LEAF;
+        else if (node_id >= JUNCTIONS_BEGIN && !m_junctions.empty()) ret = JUNCTION;
+
+        return ret;
+    }
+
+    const Node &get(size_t node_id) const
+    {
+        return *get_node(*this, node_id);
+    }
+
+    Node &get(size_t node_id)
+    {
+        return *get_node(*this, node_id);
+    }
+
+    const Node *find(size_t node_id) const { return get_node(*this, node_id); }
+    Node *find(size_t node_id) { return get_node(*this, node_id); }
+
+    // Return the original index of a leaf in the input array, if the given
+    // node id is indeed of type SUPP
+    int get_leaf_id(size_t node_id) const
+    {
+        return node_id >= LEAFS_BEGIN && node_id < JUNCTIONS_BEGIN ?
+                   node_id - LEAFS_BEGIN :
+                   Node::ID_NONE;
+    }
+
+    size_t get_queue_idx(size_t node_id) const { return m_queue_indices[node_id]; }
+
+    float get_distance(const Vec3f &p, size_t node) const;
+
+    size_t next_junction_id() const
+    {
+        return JUNCTIONS_BEGIN + m_junctions.size();
+    }
+
+    size_t insert_junction(const Node &p)
+    {
+        size_t new_id = next_junction_id();
+        m_junctions.emplace_back(p);
+        m_junctions.back().id = int(new_id);
+        m_ktree.insert({m_junctions.back().pos, new_id});
+        m_searchable_indices.emplace_back(true);
+        m_queue_indices.emplace_back(Unqueued);
+        ++m_reachable_cnt;
+
+        return new_id;
+    }
+
+    const std::vector<Node> &get_junctions() const noexcept  { return m_junctions; }
+    const std::vector<Node> &get_bedpoints() const noexcept  { return m_bedpoints; }
+    const std::vector<Node> &get_meshpoints() const noexcept { return m_meshpoints; }
+    const std::vector<Node> &get_leafs() const noexcept      { return m_leafs; }
+    const Properties & properties() const noexcept { return m_props; }
+
+    void mark_unreachable(size_t node_id)
+    {
+        assert(node_id < m_searchable_indices.size());
+
+        m_searchable_indices[node_id] = false;
+        m_queue_indices[node_id] = Unqueued;
+        --m_reachable_cnt;
+    }
+
+    size_t reachable_count() const { return m_reachable_cnt; }
+
+    template<class Fn>
+    void foreach_reachable(const Vec3f &pos,
+                           Fn         &&visitor,
+                           size_t       k,
+                           double       min_dist = 0.)
+    {
+        // Fake output iterator
+        struct Output {
+            const PointCloud *self;
+            Vec3f p;
+            Fn &visitorfn;
+
+            Output& operator *() { return *this; }
+            void operator=(const PointIndexEl &el) {
+                visitorfn(el.second, self->get_distance(p, el.second),
+                             (p - el.first).squaredNorm());
+            }
+            Output& operator++() { return *this; }
+        };
+
+        namespace bgi = boost::geometry::index;
+        float brln = 2 * m_props.max_branch_length();
+        BoundingBox3Base<Vec3f> bb{{pos.x() - brln, pos.y() - brln,
+                                    float(m_props.ground_level() - EPSILON)},
+                                   {pos.x() + brln, pos.y() + brln,
+                                    m_ktree.bounds().max_corner().get<Z>()}};
+
+        // Extend upwards to find mergable junctions and support points
+        bb.max.z() = m_ktree.bounds().max_corner().get<Z>();
+
+        auto filter = bgi::satisfies(
+                          [this, &pos, min_dist](const PointIndexEl &e) {
+                              double D_branching = get_distance(pos, e.second);
+                              double D_euql      = (pos - e.first).squaredNorm() ;
+                              return m_searchable_indices[e.second] &&
+                                     !std::isinf(D_branching) && D_euql > min_dist;
+                          });
+
+        m_ktree.query(bgi::intersects(bb) && filter && bgi::nearest(pos, k),
+                      Output{this, pos, visitor});
+    }
+
+    auto start_queue()
+    {
+        auto ptsqueue = make_mutable_priority_queue<size_t, true>(
+            [this](size_t el, size_t idx) { m_queue_indices[el] = idx; },
+            ZCompareFn{this});
+
+        ptsqueue.reserve(m_leafs.size());
+        size_t iend = LEAFS_BEGIN + m_leafs.size();
+        for (size_t i = LEAFS_BEGIN; i < iend; ++i)
+            ptsqueue.push(i);
+
+        return ptsqueue;
+    }
+};
+
+template<class PC, class Fn> void traverse(PC &&pc, size_t root, Fn &&fn)
+{
+    if (auto nodeptr = pc.find(root); nodeptr != nullptr) {
+        auto &nroot = *nodeptr;
+        fn(nroot);
+        if (nroot.left  >= 0) traverse(pc, nroot.left, fn);
+        if (nroot.right >= 0) traverse(pc, nroot.right, fn);
+    }
+}
+
+void build_tree(PointCloud &pcloud, Builder &builder);
+
+}} // namespace Slic3r::branchingtree
+
+#endif // POINTCLOUD_HPP
diff --git a/src/libslic3r/CMakeLists.txt b/src/libslic3r/CMakeLists.txt
index e637b1402..2c5944fef 100644
--- a/src/libslic3r/CMakeLists.txt
+++ b/src/libslic3r/CMakeLists.txt
@@ -319,6 +319,12 @@ set(SLIC3R_SOURCES
     SLA/ReprojectPointsOnMesh.hpp
     SLA/DefaultSupportTree.hpp
     SLA/DefaultSupportTree.cpp
+    SLA/BranchingTreeSLA.hpp
+    SLA/BranchingTreeSLA.cpp
+    BranchingTree/BranchingTree.cpp
+    BranchingTree/BranchingTree.hpp
+    BranchingTree/PointCloud.cpp
+    BranchingTree/PointCloud.hpp
 
     Arachne/BeadingStrategy/BeadingStrategy.hpp
     Arachne/BeadingStrategy/BeadingStrategy.cpp
diff --git a/src/libslic3r/Geometry.hpp b/src/libslic3r/Geometry.hpp
index c04bdf07a..d0ae84d86 100644
--- a/src/libslic3r/Geometry.hpp
+++ b/src/libslic3r/Geometry.hpp
@@ -560,13 +560,13 @@ inline bool is_rotation_ninety_degrees(const Vec3d &rotation)
     return is_rotation_ninety_degrees(rotation.x()) && is_rotation_ninety_degrees(rotation.y()) && is_rotation_ninety_degrees(rotation.z());
 }
 
-template <class T>
-std::pair<T, T> dir_to_spheric(const Vec<3, T> &n, T norm = 1.)
+template <class Tout = double, class Tin>
+std::pair<Tout, Tout> dir_to_spheric(const Vec<3, Tin> &n, Tout norm = 1.)
 {
-    T z       = n.z();
-    T r       = norm;
-    T polar   = std::acos(z / r);
-    T azimuth = std::atan2(n(1), n(0));
+    Tout z       = n.z();
+    Tout r       = norm;
+    Tout polar   = std::acos(z / r);
+    Tout azimuth = std::atan2(n(1), n(0));
     return {polar, azimuth};
 }
 
diff --git a/src/libslic3r/Preset.cpp b/src/libslic3r/Preset.cpp
index 092c0ed9b..db6ebe62c 100644
--- a/src/libslic3r/Preset.cpp
+++ b/src/libslic3r/Preset.cpp
@@ -493,6 +493,7 @@ static std::vector<std::string> s_Preset_sla_print_options {
     "layer_height",
     "faded_layers",
     "supports_enable",
+    "support_tree_type",
     "support_head_front_diameter",
     "support_head_penetration",
     "support_head_width",
diff --git a/src/libslic3r/PrintConfig.cpp b/src/libslic3r/PrintConfig.cpp
index bd4f98dfd..40a3cb3d5 100644
--- a/src/libslic3r/PrintConfig.cpp
+++ b/src/libslic3r/PrintConfig.cpp
@@ -176,6 +176,12 @@ static const t_config_enum_values s_keys_map_SLAMaterialSpeed = {
 };
 CONFIG_OPTION_ENUM_DEFINE_STATIC_MAPS(SLAMaterialSpeed);
 
+static inline const t_config_enum_values s_keys_map_SLASupportTreeType = {
+    {"default", int(sla::SupportTreeType::Default)},
+    {"branching",   int(sla::SupportTreeType::Branching)}
+};
+CONFIG_OPTION_ENUM_DEFINE_STATIC_MAPS(SLASupportTreeType);
+
 static const t_config_enum_values s_keys_map_BrimType = {
     {"no_brim",         btNoBrim},
     {"outer_only",      btOuterOnly},
@@ -3564,6 +3570,17 @@ void PrintConfigDef::init_sla_params()
     def->mode = comSimple;
     def->set_default_value(new ConfigOptionBool(true));
 
+    def = this->add("support_tree_type", coEnum);
+    def->label = L("Support tree type");
+    def->tooltip = L("Support tree building strategy");
+    def->enum_keys_map = &ConfigOptionEnum<sla::SupportTreeType>::get_enum_values();
+    def->enum_values = ConfigOptionEnum<sla::SupportTreeType>::get_enum_names();
+    def->enum_labels = ConfigOptionEnum<sla::SupportTreeType>::get_enum_names();
+    def->enum_labels[0] = L("Default");
+    def->enum_labels[1] = L("Branching");
+    def->mode = comAdvanced;
+    def->set_default_value(new ConfigOptionEnum(sla::SupportTreeType::Default));
+
     def = this->add("support_head_front_diameter", coFloat);
     def->label = L("Pinhead front diameter");
     def->category = L("Supports");
@@ -3649,13 +3666,17 @@ void PrintConfigDef::init_sla_params()
     def = this->add("support_pillar_widening_factor", coFloat);
     def->label = L("Pillar widening factor");
     def->category = L("Supports");
-    def->tooltip = L("Merging bridges or pillars into another pillars can "
-                     "increase the radius. Zero means no increase, one means "
-                     "full increase.");
+    def->tooltip  = L(
+         "Merging bridges or pillars into another pillars can "
+          "increase the radius. Zero means no increase, one means "
+          "full increase. The exact amount of increase is unspecified and can "
+          "change in the future. What is garanteed is that thickness will not "
+          "exceed \"support_base_diameter\"");
+
     def->min = 0;
     def->max = 1;
     def->mode = comExpert;
-    def->set_default_value(new ConfigOptionFloat(0.0));
+    def->set_default_value(new ConfigOptionFloat(0.15));
 
     def = this->add("support_base_diameter", coFloat);
     def->label = L("Support base diameter");
diff --git a/src/libslic3r/PrintConfig.hpp b/src/libslic3r/PrintConfig.hpp
index 9e1d7989d..5ec163572 100644
--- a/src/libslic3r/PrintConfig.hpp
+++ b/src/libslic3r/PrintConfig.hpp
@@ -155,6 +155,7 @@ CONFIG_OPTION_ENUM_DECLARE_STATIC_MAPS(SupportMaterialInterfacePattern)
 CONFIG_OPTION_ENUM_DECLARE_STATIC_MAPS(SeamPosition)
 CONFIG_OPTION_ENUM_DECLARE_STATIC_MAPS(SLADisplayOrientation)
 CONFIG_OPTION_ENUM_DECLARE_STATIC_MAPS(SLAPillarConnectionMode)
+CONFIG_OPTION_ENUM_DECLARE_STATIC_MAPS(SLASupportTreeType)
 CONFIG_OPTION_ENUM_DECLARE_STATIC_MAPS(BrimType)
 CONFIG_OPTION_ENUM_DECLARE_STATIC_MAPS(DraftShield)
 CONFIG_OPTION_ENUM_DECLARE_STATIC_MAPS(GCodeThumbnailsFormat)
@@ -829,6 +830,8 @@ PRINT_CONFIG_CLASS_DEFINE(
     // Enabling or disabling support creation
     ((ConfigOptionBool,  supports_enable))
 
+    ((ConfigOptionEnum<sla::SupportTreeType>, support_tree_type))
+
     // Diameter in mm of the pointing side of the head.
     ((ConfigOptionFloat, support_head_front_diameter))/*= 0.2*/
 
diff --git a/src/libslic3r/SLA/BranchingTreeSLA.cpp b/src/libslic3r/SLA/BranchingTreeSLA.cpp
new file mode 100644
index 000000000..72314c9d8
--- /dev/null
+++ b/src/libslic3r/SLA/BranchingTreeSLA.cpp
@@ -0,0 +1,257 @@
+#include "BranchingTreeSLA.hpp"
+
+#include "libslic3r/Execution/ExecutionTBB.hpp"
+
+#include "libslic3r/KDTreeIndirect.hpp"
+
+#include "SupportTreeUtils.hpp"
+#include "BranchingTree/PointCloud.hpp"
+
+#include "Pad.hpp"
+
+#include <map>
+
+namespace Slic3r { namespace sla {
+
+class BranchingTreeBuilder: public branchingtree::Builder {
+    SupportTreeBuilder &m_builder;
+    const SupportableMesh  &m_sm;
+    const branchingtree::PointCloud &m_cloud;
+
+    // Scaling of the input value 'widening_factor:<0, 1>' to produce resonable
+    // widening behaviour
+    static constexpr double WIDENING_SCALE = 0.02;
+
+    double get_radius(const branchingtree::Node &j)
+    {
+        double w = WIDENING_SCALE * m_sm.cfg.pillar_widening_factor * j.weight;
+
+        return std::min(m_sm.cfg.base_radius_mm, double(j.Rmin) + w);
+    }
+
+    std::vector<size_t>  m_unroutable_pinheads;
+
+    void build_subtree(size_t root)
+    {
+        traverse(m_cloud, root, [this](const branchingtree::Node &node) {
+            if (node.left >= 0 && node.right >= 0) {
+                auto nparent = m_cloud.get(node.id);
+                auto nleft = m_cloud.get(node.left);
+                auto nright = m_cloud.get(node.right);
+                Vec3d from1d = nleft.pos.cast<double>();
+                Vec3d from2d = nright.pos.cast<double>();
+                Vec3d tod    = nparent.pos.cast<double>();
+                double mergeR = get_radius(nparent);
+                double leftR  = get_radius(nleft);
+                double rightR = get_radius(nright);
+
+                m_builder.add_diffbridge(from1d, tod, leftR, mergeR);
+                m_builder.add_diffbridge(from2d, tod, rightR, mergeR);
+                m_builder.add_junction(tod, mergeR);
+            } else if (int child = node.left + node.right + 1; child >= 0) {
+                auto from = m_cloud.get(child);
+                auto to   = m_cloud.get(node.id);
+                auto tod  = to.pos.cast<double>();
+                double toR = get_radius(to);
+                m_builder.add_diffbridge(from.pos.cast<double>(),
+                                         tod,
+                                         get_radius(from),
+                                         toR);
+                m_builder.add_junction(tod, toR);
+            }
+        });
+    }
+
+    void discard_subtree(size_t root)
+    {
+        // Discard all the support points connecting to this branch.
+        traverse(m_cloud, root, [this](const branchingtree::Node &node) {
+            int suppid_parent = m_cloud.get_leaf_id(node.id);
+            int suppid_left   = m_cloud.get_leaf_id(node.left);
+            int suppid_right  = m_cloud.get_leaf_id(node.right);
+            if (suppid_parent >= 0)
+                m_unroutable_pinheads.emplace_back(suppid_parent);
+            if (suppid_left >= 0)
+                m_unroutable_pinheads.emplace_back(suppid_left);
+            if (suppid_right >= 0)
+                m_unroutable_pinheads.emplace_back(suppid_right);
+        });
+    }
+
+public:
+    BranchingTreeBuilder(SupportTreeBuilder          &builder,
+                     const SupportableMesh       &sm,
+                     const branchingtree::PointCloud &cloud)
+        : m_builder{builder}, m_sm{sm}, m_cloud{cloud}
+    {}
+
+    bool add_bridge(const branchingtree::Node &from,
+                    const branchingtree::Node &to) override;
+
+    bool add_merger(const branchingtree::Node &node,
+                    const branchingtree::Node &closest,
+                    const branchingtree::Node &merge_node) override;
+
+    bool add_ground_bridge(const branchingtree::Node &from,
+                           const branchingtree::Node &/*to*/) override;
+
+    bool add_mesh_bridge(const branchingtree::Node &from,
+                         const branchingtree::Node &to) override;
+
+    void report_unroutable(const branchingtree::Node &j) override
+    {
+        BOOST_LOG_TRIVIAL(error) << "Cannot route junction at " << j.pos.x()
+                                 << " " << j.pos.y() << " " << j.pos.z();
+
+        // Discard all the support points connecting to this branch.
+        discard_subtree(j.id);
+    }
+
+    const std::vector<size_t>& unroutable_pinheads() const
+    {
+        return m_unroutable_pinheads;
+    }
+
+    bool is_valid() const override { return !m_builder.ctl().stopcondition(); }
+};
+
+bool BranchingTreeBuilder::add_bridge(const branchingtree::Node &from,
+                                  const branchingtree::Node &to)
+{
+    Vec3d fromd = from.pos.cast<double>(), tod = to.pos.cast<double>();
+    double fromR = get_radius(from), toR = get_radius(to);
+    Beam beam{Ball{fromd, fromR}, Ball{tod, toR}};
+    auto   hit = beam_mesh_hit(ex_tbb, m_sm.emesh, beam,
+                               m_sm.cfg.safety_distance_mm);
+
+    bool ret = hit.distance() > (tod - fromd).norm();
+
+    return ret;
+}
+
+bool BranchingTreeBuilder::add_merger(const branchingtree::Node &node,
+                                      const branchingtree::Node &closest,
+                                      const branchingtree::Node &merge_node)
+{
+    Vec3d from1d = node.pos.cast<double>(),
+          from2d = closest.pos.cast<double>(),
+          tod    = merge_node.pos.cast<double>();
+
+    double mergeR   = get_radius(merge_node);
+    double nodeR    = get_radius(node);
+    double closestR = get_radius(closest);
+    Beam beam1{Ball{from1d, nodeR}, Ball{tod, mergeR}};
+    Beam beam2{Ball{from2d, closestR}, Ball{tod, mergeR}};
+
+    auto sd = m_sm.cfg.safety_distance_mm ;
+    auto hit1 = beam_mesh_hit(ex_tbb, m_sm.emesh, beam1, sd);
+    auto hit2 = beam_mesh_hit(ex_tbb, m_sm.emesh, beam2, sd);
+
+    bool ret = hit1.distance() > (tod - from1d).norm() &&
+               hit2.distance() > (tod - from2d).norm();
+
+    return ret;
+}
+
+bool BranchingTreeBuilder::add_ground_bridge(const branchingtree::Node &from,
+                                             const branchingtree::Node &to)
+{
+    bool ret = search_ground_route(ex_tbb, m_builder, m_sm,
+                                   sla::Junction{from.pos.cast<double>(),
+                                                 get_radius(from)},
+                                   get_radius(to)).first;
+
+    if (ret) {
+        build_subtree(from.id);
+    }
+
+    return ret;
+}
+
+bool BranchingTreeBuilder::add_mesh_bridge(const branchingtree::Node &from,
+                                           const branchingtree::Node &to)
+{
+    sla::Junction fromj = {from.pos.cast<double>(), get_radius(from)};
+
+    auto anchor = m_sm.cfg.ground_facing_only ?
+                      std::optional<Anchor>{} : // If no mesh connections are allowed
+                      calculate_anchor_placement(ex_tbb, m_sm, fromj,
+                                                 to.pos.cast<double>());
+
+    if (anchor) {
+        sla::Junction toj = {anchor->junction_point(), anchor->r_back_mm};
+
+        auto hit = beam_mesh_hit(ex_tbb, m_sm.emesh,
+                                 Beam{{fromj.pos, fromj.r}, {toj.pos, toj.r}}, 0.);
+
+        if (hit.distance() > distance(fromj.pos, toj.pos)) {
+            m_builder.add_diffbridge(fromj.pos, toj.pos, fromj.r, toj.r);
+            m_builder.add_anchor(*anchor);
+
+            build_subtree(from.id);
+        } else {
+            anchor.reset();
+        }
+    }
+
+    return bool(anchor);
+}
+
+void create_branching_tree(SupportTreeBuilder &builder, const SupportableMesh &sm)
+{
+    auto coordfn = [&sm](size_t id, size_t dim) { return sm.pts[id].pos(dim); };
+    KDTreeIndirect<3, float, decltype (coordfn)> tree{coordfn, sm.pts.size()};
+
+    auto nondup_idx = non_duplicate_suppt_indices(tree, sm.pts, 0.1);
+    std::vector<std::optional<Head>> heads(nondup_idx.size());
+    auto leafs = reserve_vector<branchingtree::Node>(nondup_idx.size());
+
+    execution::for_each(
+        ex_tbb, size_t(0), nondup_idx.size(),
+        [&sm, &heads, &nondup_idx, &builder](size_t i) {
+            if (!builder.ctl().stopcondition())
+                heads[i] = calculate_pinhead_placement(ex_seq, sm, nondup_idx[i]);
+        },
+        execution::max_concurrency(ex_tbb)
+    );
+
+    if (builder.ctl().stopcondition())
+        return;
+
+    for (auto &h : heads)
+        if (h && h->is_valid()) {
+            leafs.emplace_back(h->junction_point().cast<float>(), h->r_back_mm);
+            h->id = leafs.size() - 1;
+            builder.add_head(h->id, *h);
+        }
+
+    auto &its = *sm.emesh.get_triangle_mesh();
+    ExPolygons bedpolys = {branchingtree::make_bed_poly(its)};
+
+    auto props = branchingtree::Properties{}
+                     .bed_shape(bedpolys)
+                     .ground_level(sla::ground_level(sm))
+                     .max_slope(sm.cfg.bridge_slope)
+                     .max_branch_length(sm.cfg.max_bridge_length_mm);
+
+    auto meshpts = sm.cfg.ground_facing_only ?
+                       std::vector<branchingtree::Node>{} :
+                       branchingtree::sample_mesh(its,
+                                                  props.sampling_radius());
+
+    auto bedpts  = branchingtree::sample_bed(props.bed_shape(),
+                                             props.ground_level(),
+                                             props.sampling_radius());
+
+    branchingtree::PointCloud nodes{std::move(meshpts), std::move(bedpts),
+                                    std::move(leafs), props};
+
+    BranchingTreeBuilder vbuilder{builder, sm, nodes};
+    branchingtree::build_tree(nodes, vbuilder);
+
+    for (size_t id : vbuilder.unroutable_pinheads())
+        builder.head(id).invalidate();
+
+}
+
+}} // namespace Slic3r::sla
diff --git a/src/libslic3r/SLA/BranchingTreeSLA.hpp b/src/libslic3r/SLA/BranchingTreeSLA.hpp
new file mode 100644
index 000000000..c101029bd
--- /dev/null
+++ b/src/libslic3r/SLA/BranchingTreeSLA.hpp
@@ -0,0 +1,15 @@
+#ifndef BRANCHINGTREESLA_HPP
+#define BRANCHINGTREESLA_HPP
+
+#include "libslic3r/BranchingTree/BranchingTree.hpp"
+#include "SupportTreeBuilder.hpp"
+
+#include <boost/log/trivial.hpp>
+
+namespace Slic3r { namespace sla {
+
+void create_branching_tree(SupportTreeBuilder& builder, const SupportableMesh &sm);
+
+}} // namespace Slic3r::sla
+
+#endif // BRANCHINGTREESLA_HPP
diff --git a/src/libslic3r/SLA/SupportTree.cpp b/src/libslic3r/SLA/SupportTree.cpp
index 466fade88..b221c4330 100644
--- a/src/libslic3r/SLA/SupportTree.cpp
+++ b/src/libslic3r/SLA/SupportTree.cpp
@@ -8,6 +8,7 @@
 #include <libslic3r/SLA/SpatIndex.hpp>
 #include <libslic3r/SLA/SupportTreeBuilder.hpp>
 #include <libslic3r/SLA/DefaultSupportTree.hpp>
+#include <libslic3r/SLA/BranchingTreeSLA.hpp>
 
 #include <libslic3r/MTUtils.hpp>
 #include <libslic3r/ClipperUtils.hpp>
@@ -34,6 +35,10 @@ indexed_triangle_set create_support_tree(const SupportableMesh &sm,
             create_default_tree(*builder, sm);
             break;
         }
+        case SupportTreeType::Branching: {
+            create_branching_tree(*builder, sm);
+            break;
+        }
         default:;
         }
 
diff --git a/src/libslic3r/SLA/SupportTreeStrategies.hpp b/src/libslic3r/SLA/SupportTreeStrategies.hpp
index d04ceaa2d..487f43575 100644
--- a/src/libslic3r/SLA/SupportTreeStrategies.hpp
+++ b/src/libslic3r/SLA/SupportTreeStrategies.hpp
@@ -5,7 +5,7 @@
 
 namespace Slic3r { namespace sla {
 
-enum class SupportTreeType { Default, Clever };
+enum class SupportTreeType { Default, Branching };
 enum class PillarConnectionMode { zigzag, cross, dynamic };
 
 }} // namespace Slic3r::sla
diff --git a/src/libslic3r/SLAPrint.cpp b/src/libslic3r/SLAPrint.cpp
index 9ad545338..d9b8e33df 100644
--- a/src/libslic3r/SLAPrint.cpp
+++ b/src/libslic3r/SLAPrint.cpp
@@ -43,6 +43,7 @@ sla::SupportTreeConfig make_support_cfg(const SLAPrintObjectConfig& c)
     sla::SupportTreeConfig scfg;
     
     scfg.enabled = c.supports_enable.getBool();
+    scfg.tree_type = c.support_tree_type.value;
     scfg.head_front_radius_mm = 0.5*c.support_head_front_diameter.getFloat();
     double pillar_r = 0.5 * c.support_pillar_diameter.getFloat();
     scfg.head_back_radius_mm = pillar_r;
@@ -824,6 +825,7 @@ bool SLAPrintObject::invalidate_state_by_config_options(const std::vector<t_conf
             || opt_key == "pad_enable"
             || opt_key == "pad_wall_thickness"
             || opt_key == "supports_enable"
+            || opt_key == "support_tree_type"
             || opt_key == "support_object_elevation"
             || opt_key == "pad_around_object"
             || opt_key == "pad_around_object_everywhere"
diff --git a/src/slic3r/GUI/ConfigManipulation.cpp b/src/slic3r/GUI/ConfigManipulation.cpp
index eb6d012dd..dcdadebad 100644
--- a/src/slic3r/GUI/ConfigManipulation.cpp
+++ b/src/slic3r/GUI/ConfigManipulation.cpp
@@ -363,21 +363,26 @@ void ConfigManipulation::update_print_sla_config(DynamicPrintConfig* config, con
 void ConfigManipulation::toggle_print_sla_options(DynamicPrintConfig* config)
 {
     bool supports_en = config->opt_bool("supports_enable");
+    sla::SupportTreeType treetype = config->opt_enum<sla::SupportTreeType>("support_tree_type");
+    bool is_default_tree = treetype == sla::SupportTreeType::Default;
+    bool is_branching_tree = treetype == sla::SupportTreeType::Branching;
 
     toggle_field("support_head_front_diameter", supports_en);
     toggle_field("support_head_penetration", supports_en);
     toggle_field("support_head_width", supports_en);
     toggle_field("support_pillar_diameter", supports_en);
     toggle_field("support_small_pillar_diameter_percent", supports_en);
-    toggle_field("support_max_bridges_on_pillar", supports_en);
-    toggle_field("support_pillar_connection_mode", supports_en);
+    toggle_field("support_max_bridges_on_pillar", supports_en && is_default_tree);
+    toggle_field("support_pillar_connection_mode", supports_en && is_default_tree);
+    toggle_field("support_tree_type", supports_en);
     toggle_field("support_buildplate_only", supports_en);
     toggle_field("support_base_diameter", supports_en);
     toggle_field("support_base_height", supports_en);
     toggle_field("support_base_safety_distance", supports_en);
     toggle_field("support_critical_angle", supports_en);
     toggle_field("support_max_bridge_length", supports_en);
-    toggle_field("support_max_pillar_link_distance", supports_en);
+    toggle_field("support_max_pillar_link_distance", supports_en && is_default_tree);
+    toggle_field("support_pillar_widening_factor", supports_en && is_branching_tree);
     toggle_field("support_points_density_relative", supports_en);
     toggle_field("support_points_minimal_distance", supports_en);
 
diff --git a/src/slic3r/GUI/Tab.cpp b/src/slic3r/GUI/Tab.cpp
index a5b4c7f7d..8d5d4b067 100644
--- a/src/slic3r/GUI/Tab.cpp
+++ b/src/slic3r/GUI/Tab.cpp
@@ -4637,6 +4637,7 @@ void TabSLAPrint::build()
     page = add_options_page(L("Supports"), "support"/*"sla_supports"*/);
     optgroup = page->new_optgroup(L("Supports"));
     optgroup->append_single_option_line("supports_enable");
+    optgroup->append_single_option_line("support_tree_type");
 
     optgroup = page->new_optgroup(L("Support head"));
     optgroup->append_single_option_line("support_head_front_diameter");
@@ -4650,8 +4651,7 @@ void TabSLAPrint::build()
     
     optgroup->append_single_option_line("support_pillar_connection_mode");
     optgroup->append_single_option_line("support_buildplate_only");
-    // TODO: This parameter is not used at the moment.
-    // optgroup->append_single_option_line("support_pillar_widening_factor");
+    optgroup->append_single_option_line("support_pillar_widening_factor");
     optgroup->append_single_option_line("support_base_diameter");
     optgroup->append_single_option_line("support_base_height");
     optgroup->append_single_option_line("support_base_safety_distance");
diff --git a/tests/sla_print/sla_print_tests.cpp b/tests/sla_print/sla_print_tests.cpp
index a5ed0d029..8ea91d57a 100644
--- a/tests/sla_print/sla_print_tests.cpp
+++ b/tests/sla_print/sla_print_tests.cpp
@@ -8,6 +8,7 @@
 
 #include <libslic3r/TriangleMeshSlicer.hpp>
 #include <libslic3r/SLA/SupportTreeMesher.hpp>
+#include <libslic3r/BranchingTree/PointCloud.hpp>
 
 namespace {
 
@@ -156,37 +157,138 @@ TEST_CASE("DefaultSupports::FloorSupportsDoNotPierceModel", "[SLASupportGenerati
         test_support_model_collision(fname, supportcfg);
 }
 
-//TEST_CASE("CleverSupports::ElevatedSupportGeometryIsValid", "[SLASupportGeneration][Clever]") {
+//TEST_CASE("BranchingSupports::ElevatedSupportGeometryIsValid", "[SLASupportGeneration][Branching]") {
 //    sla::SupportTreeConfig supportcfg;
 //    supportcfg.object_elevation_mm = 10.;
-//    supportcfg.tree_type = sla::SupportTreeType::Clever;
+//    supportcfg.tree_type = sla::SupportTreeType::Branching;
 
 //    for (auto fname : SUPPORT_TEST_MODELS) test_supports(fname, supportcfg);
 //}
 
-//TEST_CASE("CleverSupports::FloorSupportGeometryIsValid", "[SLASupportGeneration][Clever]") {
+//TEST_CASE("BranchingSupports::FloorSupportGeometryIsValid", "[SLASupportGeneration][Branching]") {
 //    sla::SupportTreeConfig supportcfg;
 //    supportcfg.object_elevation_mm = 0;
-//    supportcfg.tree_type = sla::SupportTreeType::Clever;
+//    supportcfg.tree_type = sla::SupportTreeType::Branching;
 
 //    for (auto &fname: SUPPORT_TEST_MODELS) test_supports(fname, supportcfg);
 //}
 
-TEST_CASE("CleverSupports::ElevatedSupportsDoNotPierceModel", "[SLASupportGeneration][Clever]") {
+bool is_outside_support_cone(const Vec3f &supp, const Vec3f &pt, float angle)
+{
+    Vec3d D = (pt - supp).cast<double>();
+    double dot_sq = -D.z() * std::abs(-D.z());
+
+    return dot_sq <
+           D.squaredNorm() * std::cos(angle) * std::abs(std::cos(angle));
+}
+
+TEST_CASE("BranchingSupports::MergePointFinder", "[SLASupportGeneration][Branching]") {
+    SECTION("Identical points have the same merge point") {
+        Vec3f a{0.f, 0.f, 0.f}, b = a;
+        auto slope = float(PI / 4.);
+
+        auto mergept = branchingtree::find_merge_pt(a, b, slope);
+
+        REQUIRE(bool(mergept));
+        REQUIRE((*mergept - b).norm() < EPSILON);
+        REQUIRE((*mergept - a).norm() < EPSILON);
+    }
+
+    // ^ Z
+    // | a *
+    // |
+    // | b * <= mergept
+    SECTION("Points at different heights have the lower point as mergepoint") {
+        Vec3f a{0.f, 0.f, 0.f}, b = {0.f, 0.f, -1.f};
+        auto slope = float(PI / 4.);
+
+        auto mergept = branchingtree::find_merge_pt(a, b, slope);
+
+        REQUIRE(bool(mergept));
+        REQUIRE((*mergept - b).squaredNorm() < 2 * EPSILON);
+    }
+
+    // -|---------> X
+    //  a       b
+    //  *       *
+    //      * <= mergept
+    SECTION("Points at different X have mergept in the middle at lower Z") {
+        Vec3f a{0.f, 0.f, 0.f}, b = {1.f, 0.f, 0.f};
+        auto slope = float(PI / 4.);
+
+        auto mergept = branchingtree::find_merge_pt(a, b, slope);
+
+        REQUIRE(bool(mergept));
+
+           // Distance of mergept should be equal from both input points
+        float D = std::abs((*mergept - b).squaredNorm() - (*mergept - a).squaredNorm());
+
+        REQUIRE(D < EPSILON);
+        REQUIRE(!is_outside_support_cone(a, *mergept, slope));
+        REQUIRE(!is_outside_support_cone(b, *mergept, slope));
+    }
+
+    // -|---------> Y
+    //  a       b
+    //  *       *
+    //      * <= mergept
+    SECTION("Points at different Y have mergept in the middle at lower Z") {
+        Vec3f a{0.f, 0.f, 0.f}, b = {0.f, 1.f, 0.f};
+        auto slope = float(PI / 4.);
+
+        auto mergept = branchingtree::find_merge_pt(a, b, slope);
+
+        REQUIRE(bool(mergept));
+
+        // Distance of mergept should be equal from both input points
+        float D = std::abs((*mergept - b).squaredNorm() - (*mergept - a).squaredNorm());
+
+        REQUIRE(D < EPSILON);
+        REQUIRE(!is_outside_support_cone(a, *mergept, slope));
+        REQUIRE(!is_outside_support_cone(b, *mergept, slope));
+    }
+
+    SECTION("Points separated by less than critical angle have the lower point as mergept") {
+        Vec3f a{-1.f, -1.f, -1.f}, b = {-1.5f, -1.5f, -2.f};
+        auto slope = float(PI / 4.);
+
+        auto mergept = branchingtree::find_merge_pt(a, b, slope);
+
+        REQUIRE(bool(mergept));
+        REQUIRE((*mergept - b).norm() < 2 * EPSILON);
+    }
+
+    // -|----------------------------> Y
+    //  a                          b
+    //  *            * <= mergept  *
+    //
+    SECTION("Points at same height have mergepoint in the middle if critical angle is zero ") {
+        Vec3f a{-1.f, -1.f, -1.f}, b = {-1.5f, -1.5f, -1.f};
+        auto slope = EPSILON;
+
+        auto mergept = branchingtree::find_merge_pt(a, b, slope);
+
+        REQUIRE(bool(mergept));
+        Vec3f middle = (b + a) / 2.;
+        REQUIRE((*mergept - middle).norm() < 4 * EPSILON);
+    }
+}
+
+TEST_CASE("BranchingSupports::ElevatedSupportsDoNotPierceModel", "[SLASupportGeneration][Branching]") {
 
     sla::SupportTreeConfig supportcfg;
     supportcfg.object_elevation_mm = 10.;
-    supportcfg.tree_type = sla::SupportTreeType::Clever;
+    supportcfg.tree_type = sla::SupportTreeType::Branching;
 
     for (auto fname : SUPPORT_TEST_MODELS)
         test_support_model_collision(fname, supportcfg);
 }
 
-TEST_CASE("CleverSupports::FloorSupportsDoNotPierceModel", "[SLASupportGeneration][Clever]") {
+TEST_CASE("BranchingSupports::FloorSupportsDoNotPierceModel", "[SLASupportGeneration][Branching]") {
 
     sla::SupportTreeConfig supportcfg;
     supportcfg.object_elevation_mm = 0;
-    supportcfg.tree_type = sla::SupportTreeType::Clever;
+    supportcfg.tree_type = sla::SupportTreeType::Branching;
 
     for (auto fname : SUPPORT_TEST_MODELS)
         test_support_model_collision(fname, supportcfg);
diff --git a/tests/sla_print/sla_test_utils.cpp b/tests/sla_print/sla_test_utils.cpp
index 07f19300c..737659eed 100644
--- a/tests/sla_print/sla_test_utils.cpp
+++ b/tests/sla_print/sla_test_utils.cpp
@@ -2,57 +2,59 @@
 #include "libslic3r/TriangleMeshSlicer.hpp"
 #include "libslic3r/SLA/AGGRaster.hpp"
 #include "libslic3r/SLA/DefaultSupportTree.hpp"
+#include "libslic3r/SLA/BranchingTreeSLA.hpp"
 
 #include <iomanip>
 
-void test_support_model_collision(const std::string          &obj_filename,
-                                  const sla::SupportTreeConfig   &input_supportcfg,
-                                  const sla::HollowingConfig &hollowingcfg,
-                                  const sla::DrainHoles      &drainholes)
+void test_support_model_collision(
+    const std::string            &obj_filename,
+    const sla::SupportTreeConfig &input_supportcfg,
+    const sla::HollowingConfig   &hollowingcfg,
+    const sla::DrainHoles        &drainholes)
 {
     SupportByproducts byproducts;
-    
+
     sla::SupportTreeConfig supportcfg = input_supportcfg;
-    
+
     // Set head penetration to a small negative value which should ensure that
     // the supports will not touch the model body.
     supportcfg.head_penetration_mm = -0.2;
-    
+
     test_supports(obj_filename, supportcfg, hollowingcfg, drainholes, byproducts);
-    
+
     // Slice the support mesh given the slice grid of the model.
     std::vector<ExPolygons> support_slices =
-        sla::slice(byproducts.supporttree.retrieve_mesh(sla::MeshType::Support),
-              byproducts.supporttree.retrieve_mesh(sla::MeshType::Pad),
+        sla::slice(byproducts.suptree_builder.retrieve_mesh(sla::MeshType::Support),
+              byproducts.suptree_builder.retrieve_mesh(sla::MeshType::Pad),
               byproducts.slicegrid, CLOSING_RADIUS, {});
 
     // The slices originate from the same slice grid so the numbers must match
-    
+
     bool support_mesh_is_empty =
-            byproducts.supporttree.retrieve_mesh(sla::MeshType::Pad).empty() &&
-            byproducts.supporttree.retrieve_mesh(sla::MeshType::Support).empty();
-    
+            byproducts.suptree_builder.retrieve_mesh(sla::MeshType::Pad).empty() &&
+            byproducts.suptree_builder.retrieve_mesh(sla::MeshType::Support).empty();
+
     if (support_mesh_is_empty)
         REQUIRE(support_slices.empty());
     else
         REQUIRE(support_slices.size() == byproducts.model_slices.size());
-    
+
     bool notouch = true;
     for (size_t n = 0; notouch && n < support_slices.size(); ++n) {
         const ExPolygons &sup_slice = support_slices[n];
         const ExPolygons &mod_slice = byproducts.model_slices[n];
-        
+
         Polygons intersections = intersection(sup_slice, mod_slice);
-        
+
         double pinhead_r  = scaled(input_supportcfg.head_front_radius_mm);
 
         // TODO:: make it strict without a threshold of PI * pihead_radius ^ 2
         notouch = notouch && area(intersections) < PI * pinhead_r * pinhead_r;
     }
-    
+
     if (!notouch)
         export_failed_case(support_slices, byproducts);
-    
+
     REQUIRE(notouch);
 }
 
@@ -79,7 +81,7 @@ void export_failed_case(const std::vector<ExPolygons> &support_slices, const Sup
 
     if (do_export_stl) {
         indexed_triangle_set its;
-        byproducts.supporttree.retrieve_full_mesh(its);
+        byproducts.suptree_builder.retrieve_full_mesh(its);
         TriangleMesh m{its};
         m.merge(byproducts.input_mesh);
         m.WriteOBJFile((Catch::getResultCapture().getCurrentTestName() + "_" +
@@ -95,49 +97,49 @@ void test_supports(const std::string          &obj_filename,
 {
     using namespace Slic3r;
     TriangleMesh mesh = load_model(obj_filename);
-    
+
     REQUIRE_FALSE(mesh.empty());
-    
+
     if (hollowingcfg.enabled) {
         sla::InteriorPtr interior = sla::generate_interior(mesh, hollowingcfg);
         REQUIRE(interior);
         mesh.merge(TriangleMesh{sla::get_mesh(*interior)});
     }
-    
+
     auto   bb      = mesh.bounding_box();
     double zmin    = bb.min.z();
     double zmax    = bb.max.z();
     double gnd     = zmin - supportcfg.object_elevation_mm;
     auto   layer_h = 0.05f;
-    
+
     out.slicegrid = grid(float(gnd), float(zmax), layer_h);
     out.model_slices = slice_mesh_ex(mesh.its, out.slicegrid, CLOSING_RADIUS);
     sla::cut_drainholes(out.model_slices, out.slicegrid, CLOSING_RADIUS, drainholes, []{});
-    
+
     // Create the special index-triangle mesh with spatial indexing which
     // is the input of the support point and support mesh generators
     AABBMesh emesh{mesh};
 
-#ifdef SLIC3R_HOLE_RAYCASTER
-    if (hollowingcfg.enabled) 
+    #ifdef SLIC3R_HOLE_RAYCASTER
+    if (hollowingcfg.enabled)
         emesh.load_holes(drainholes);
-#endif
+    #endif
 
     // TODO: do the cgal hole cutting...
-    
+
     // Create the support point generator
     sla::SupportPointGenerator::Config autogencfg;
     autogencfg.head_diameter = float(2 * supportcfg.head_front_radius_mm);
     sla::SupportPointGenerator point_gen{emesh, autogencfg, [] {}, [](int) {}};
-    
+
     point_gen.seed(0); // Make the test repeatable
     point_gen.execute(out.model_slices, out.slicegrid);
-    
+
     // Get the calculated support points.
     std::vector<sla::SupportPoint> support_points = point_gen.output();
-    
+
     int validityflags = ASSUME_NO_REPAIR;
-    
+
     // If there is no elevation, support points shall be removed from the
     // bottom of the object.
     if (std::abs(supportcfg.object_elevation_mm) < EPSILON) {
@@ -145,11 +147,11 @@ void test_supports(const std::string          &obj_filename,
     } else {
         // Should be support points at least on the bottom of the model
         REQUIRE_FALSE(support_points.empty());
-        
+
         // Also the support mesh should not be empty.
         validityflags |= ASSUME_NO_EMPTY;
     }
-    
+
     // Generate the actual support tree
     sla::SupportTreeBuilder treebuilder;
     sla::SupportableMesh    sm{emesh, support_points, supportcfg};
@@ -160,29 +162,34 @@ void test_supports(const std::string          &obj_filename,
         check_support_tree_integrity(treebuilder, supportcfg, sla::ground_level(sm));
         break;
     }
+    case sla::SupportTreeType::Branching: {
+        create_branching_tree(treebuilder, sm);
+        // TODO: check_support_tree_integrity(treebuilder, supportcfg);
+        break;
+    }
     default:;
     }
 
     TriangleMesh output_mesh{treebuilder.retrieve_mesh(sla::MeshType::Support)};
-    
+
     check_validity(output_mesh, validityflags);
-    
+
     // Quick check if the dimensions and placement of supports are correct
     auto obb = output_mesh.bounding_box();
-    
+
     double allowed_zmin = zmin - supportcfg.object_elevation_mm;
-    
+
     if (std::abs(supportcfg.object_elevation_mm) < EPSILON)
         allowed_zmin = zmin - 2 * supportcfg.head_back_radius_mm;
-    
+
     REQUIRE(obb.min.z() >= Approx(allowed_zmin));
     REQUIRE(obb.max.z() <= Approx(zmax));
-    
+
     // Move out the support tree into the byproducts, we can examine it further
     // in various tests.
-    out.obj_fname   = std::move(obj_filename);
-    out.supporttree = std::move(treebuilder);
-    out.input_mesh  = std::move(mesh);
+    out.obj_fname       = std::move(obj_filename);
+    out.suptree_builder = std::move(treebuilder);
+    out.input_mesh      = std::move(mesh);
 }
 
 void check_support_tree_integrity(const sla::SupportTreeBuilder &stree, 
diff --git a/tests/sla_print/sla_test_utils.hpp b/tests/sla_print/sla_test_utils.hpp
index c5943bd5c..187a72d54 100644
--- a/tests/sla_print/sla_test_utils.hpp
+++ b/tests/sla_print/sla_test_utils.hpp
@@ -60,7 +60,7 @@ struct SupportByproducts
     std::string             obj_fname;
     std::vector<float>      slicegrid;
     std::vector<ExPolygons> model_slices;
-    sla::SupportTreeBuilder supporttree;
+    sla::SupportTreeBuilder suptree_builder;
     TriangleMesh            input_mesh;
 };