diff --git a/src/libslic3r/BranchingTree/BranchingTree.cpp b/src/libslic3r/BranchingTree/BranchingTree.cpp
new file mode 100644
index 000000000..4c9ff899c
--- /dev/null
+++ b/src/libslic3r/BranchingTree/BranchingTree.cpp
@@ -0,0 +1,141 @@
+#include "BranchingTree.hpp"
+#include "PointCloud.hpp"
+
+#include <numeric>
+#include <optional>
+#include <algorithm>
+
+#include "libslic3r/SLA/SupportTreeUtils.hpp"
+
+namespace Slic3r { namespace branchingtree {
+
+bool build_tree(PointCloud &nodes, Builder &builder)
+{
+    auto ptsqueue = nodes.start_queue();
+    auto &properties = nodes.properties();
+
+    struct NodeDistance { size_t node_id; float distance; };
+    auto distances = reserve_vector<NodeDistance>(nodes.reachable_count());
+
+    while (!ptsqueue.empty()) {
+        size_t node_id = ptsqueue.top();
+        ptsqueue.pop();
+
+        Node node = nodes.get(node_id);
+        nodes.remove_node(node_id);
+
+        distances.clear();
+        distances.reserve(nodes.reachable_count());
+
+        nodes.foreach_reachable(node.pos, [&distances](size_t id, float distance) {
+            if (!std::isinf(distance))
+                distances.emplace_back(NodeDistance{id, distance});
+        });
+
+        std::sort(distances.begin(), distances.end(),
+                  [](auto &a, auto &b) { return a.distance < b.distance; });
+
+        if (distances.empty()) {
+            builder.report_unroutable(node);
+            continue;
+        }
+
+        auto closest_it = distances.begin();
+        bool routed = false;
+        while (closest_it != distances.end() && !routed) {
+            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->distance;
+            closest_node.Rmin = std::max(node.Rmin, closest_node.Rmin);
+
+            switch (type) {
+            case BED: {
+                closest_node.weight = w;
+                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.remove_node(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.remove_node(closest_node_id);
+                }
+
+                break;
+            }
+            case SUPP:
+            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) {
+                        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);
+                            ptsqueue.remove(nodes.get_queue_idx(closest_node_id));
+                            nodes.remove_node(closest_node_id);
+                        }
+                    } else if (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)
+            builder.report_unroutable(node);
+    }
+
+    return true;
+}
+
+bool build_tree(const indexed_triangle_set & its,
+                const std::vector<Node> &support_roots,
+                Builder &                    builder,
+                const Properties &           properties)
+{
+    PointCloud nodes(its, support_roots, properties);
+
+    return build_tree(nodes, builder);
+}
+
+}} // namespace Slic3r::branchingtree
diff --git a/src/libslic3r/BranchingTree/BranchingTree.hpp b/src/libslic3r/BranchingTree/BranchingTree.hpp
new file mode 100644
index 000000000..5313e30d9
--- /dev/null
+++ b/src/libslic3r/BranchingTree/BranchingTree.hpp
@@ -0,0 +1,157 @@
+#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  = 20.;
+
+    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;
+
+    // Tracking the weight of each junction, which is essentially the sum of
+    // the lenghts of all branches emanating from this junction.
+    float weight;
+
+    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;
+};
+
+// 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.
+bool build_tree(const indexed_triangle_set & its,
+                const std::vector<Node> &support_leafs,
+                Builder &                    builder,
+                const Properties &           properties = {});
+
+inline bool build_tree(const indexed_triangle_set & its,
+                       const std::vector<Node> &support_leafs,
+                       Builder &&                   builder,
+                       const Properties &           properties = {})
+{
+    return build_tree(its, support_leafs, builder, properties);
+}
+
+//class PointCloud;
+
+//bool build_tree(PointCloud &pcloud, Builder &builder);
+
+// Helper function to derive a bed polygon only from the model bounding box.
+inline 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
+
+#endif // SUPPORTTREEBRANCHING_HPP
diff --git a/src/libslic3r/BranchingTree/PointCloud.cpp b/src/libslic3r/BranchingTree/PointCloud.cpp
new file mode 100644
index 000000000..a4855ee53
--- /dev/null
+++ b/src/libslic3r/BranchingTree/PointCloud.cpp
@@ -0,0 +1,172 @@
+#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        max_slope)
+{
+    Vec3f Da = (B - A).normalized(), Db = -Da;
+    auto [polar_da, azim_da] = Geometry::dir_to_spheric(Da);
+    auto [polar_db, azim_db] = Geometry::dir_to_spheric(Db);
+    polar_da = std::max(polar_da, float(PI) - max_slope);
+    polar_db = std::max(polar_db, float(PI) - max_slope);
+
+    Da = Geometry::spheric_to_dir<float>(polar_da, azim_da);
+    Db = Geometry::spheric_to_dir<float>(polar_db, azim_db);
+
+    double t1 =
+        (A.z() * Db.x() + Db.z() * B.x() - B.z() * Db.x() - Db.z() * A.x()) /
+        (Da.x() * Db.z() - Da.z() * Db.x());
+
+    double t2 = (A.x() + Da.x() * t1 - B.x()) / Db.x();
+
+    return t1 > 0. && t2 > 0. ? A + t1 * Da : std::optional<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++) {
+        Vec3i face = its.indices[FI(i)];
+
+        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)
+    : m_leafs{std::move(support_leafs)}
+    , m_meshpoints{sample_mesh(M, props.sampling_radius())}
+    , m_bedpoints{sample_bed(props.bed_shape(),
+                             props.ground_level(),
+                             props.sampling_radius())}
+    , m_props{props}
+    , cos2bridge_slope{std::cos(props.max_slope()) *
+                       std::abs(std::cos(props.max_slope()))}
+    , MESHPTS_BEGIN{m_bedpoints.size()}
+    , SUPP_BEGIN{MESHPTS_BEGIN + m_meshpoints.size()}
+    , JUNCTIONS_BEGIN{SUPP_BEGIN + m_leafs.size()}
+    , m_searchable_indices(JUNCTIONS_BEGIN, true)
+    , m_queue_indices(JUNCTIONS_BEGIN, UNQUEUED)
+    , m_reachable_cnt{JUNCTIONS_BEGIN}
+    , m_ktree{CoordFn{this}, SUPP_BEGIN} // Only for bed and mesh points
+{
+    for (size_t i = 0; i < m_bedpoints.size(); ++i)
+        m_bedpoints[i].id = int(i);
+
+    for (size_t i = 0; i < m_meshpoints.size(); ++i)
+        m_meshpoints[i].id = int(MESHPTS_BEGIN + i);
+
+    for (size_t i = 0; i < m_leafs.size(); ++i)
+        m_leafs[i].id = int(SUPP_BEGIN + i);
+}
+
+float PointCloud::get_distance(const Vec3f &p, size_t node)
+{
+    auto t = get_type(node);
+    auto ret = std::numeric_limits<float>::infinity();
+
+    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, get(node).pos))
+            ret = std::numeric_limits<float>::infinity();
+        else
+            ret  = (get(node).pos - p).norm();
+
+        break;
+    }
+    case SUPP:
+    case JUNCTION:{
+        auto mergept = find_merge_pt(p, get(node).pos, m_props.max_slope());
+        if (!mergept || mergept->z() < (m_props.ground_level() + 2 * get(node).Rmin))
+            ret = std::numeric_limits<float>::infinity();
+        else
+            ret = (p - *mergept).norm();
+
+        break;
+    }
+    case NONE:
+        ;
+    }
+
+       // Setting the ret val to infinity will effectively discard this
+       // connection of nodes. max_branch_length property if 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..6f583759e
--- /dev/null
+++ b/src/libslic3r/BranchingTree/PointCloud.hpp
@@ -0,0 +1,214 @@
+#ifndef POINTCLOUD_HPP
+#define POINTCLOUD_HPP
+
+#include "BranchingTree.hpp"
+
+#include "libslic3r/KDTreeIndirect.hpp"
+#include "libslic3r/MutablePriorityQueue.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 { SUPP, MESH, BED, JUNCTION, NONE };
+
+// 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, SUPP_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));
+        }
+    };
+
+    KDTreeIndirect<3, float, CoordFn> m_ktree;
+
+    bool is_outside_support_cone(const Vec3f &supp, const Vec3f &pt)
+    {
+        Vec3d D = (pt - supp).cast<double>();
+        double dot_sq = -D.z() * std::abs(-D.z());
+
+        return dot_sq < D.squaredNorm() * cos2bridge_slope;
+    }
+
+    static constexpr auto UNQUEUED = size_t(-1);
+
+    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 SUPP: ret = &pc.m_leafs [id - pc.SUPP_BEGIN]; break;
+        case JUNCTION: ret = &pc.m_junctions[id - pc.JUNCTIONS_BEGIN]; break;
+        case NONE: assert(false);
+        }
+
+        return ret;
+    }
+
+public:
+
+    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);
+
+    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 < SUPP_BEGIN && !m_meshpoints.empty()) ret = MESH;
+        else if (node_id < JUNCTIONS_BEGIN && !m_leafs.empty()) ret = SUPP;
+        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 >= SUPP_BEGIN && node_id < JUNCTIONS_BEGIN ?
+                   node_id - SUPP_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);
+
+    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_searchable_indices.emplace_back(true);
+        m_queue_indices.emplace_back(UNQUEUED);
+        ++m_reachable_cnt;
+
+        return new_id;
+    }
+
+    void remove_node(size_t node_id)
+    {
+        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)
+    {
+        auto closest_anchors =
+            find_closest_points<5>(m_ktree, pos, [this, &pos](size_t id) {
+                return m_searchable_indices[id] &&
+                       !is_outside_support_cone(pos, get(id).pos);
+            });
+
+        for (size_t anchor : closest_anchors)
+            if (anchor != m_ktree.npos)
+                visitor(anchor, get_distance(pos, anchor));
+
+        for (size_t i = SUPP_BEGIN; i < m_searchable_indices.size(); ++i)
+            if (m_searchable_indices[i])
+                visitor(i, get_distance(pos, i));
+    }
+
+    auto start_queue()
+    {
+        auto ptsqueue = make_mutable_priority_queue<size_t, false>(
+            [this](size_t el, size_t idx) { m_queue_indices[el] = idx; },
+            ZCompareFn{this});
+
+        ptsqueue.reserve(m_leafs.size());
+        size_t iend = SUPP_BEGIN + m_leafs.size();
+        for (size_t i = SUPP_BEGIN; i < iend; ++i)
+            ptsqueue.push(i);
+
+        return ptsqueue;
+    }
+
+    const Properties & properties() const { return m_props; }
+};
+
+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);
+    }
+}
+
+bool 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/SLA/BranchingTreeSLA.cpp b/src/libslic3r/SLA/BranchingTreeSLA.cpp
new file mode 100644
index 000000000..751544fb5
--- /dev/null
+++ b/src/libslic3r/SLA/BranchingTreeSLA.cpp
@@ -0,0 +1,226 @@
+#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.2;
+
+    double get_radius(const branchingtree::Node &j)
+    {
+        double w = WIDENING_SCALE * m_sm.cfg.pillar_widening_factor * j.weight;
+
+        return std::max(double(j.Rmin), std::min(m_sm.cfg.base_radius_mm, 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);
+                m_builder.add_diffbridge(from.pos.cast<double>(),
+                                         to.pos.cast<double>(),
+                                         get_radius(from),
+                                         get_radius(to));
+            }
+        });
+    }
+
+    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 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 = calculate_anchor_placement(ex_tbb, m_sm,
+                                             fromj,
+                                             to.pos.cast<double>());
+
+    if (anchor) {
+        m_builder.add_diffbridge(fromj.pos, anchor->junction_point(), fromj.r,
+                                 anchor->r_back_mm);
+
+        m_builder.add_anchor(*anchor);
+
+        build_subtree(from.id);
+    }
+
+    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](size_t i) {
+            heads[i] = calculate_pinhead_placement(ex_seq, sm, i);
+        },
+        execution::max_concurrency(ex_tbb)
+    );
+
+    for (auto &h : heads)
+        if (h && h->is_valid()) {
+            leafs.emplace_back(h->junction_point().cast<float>(), h->r_back_mm);
+            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);
+
+    branchingtree::PointCloud nodes{its, 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/tests/sla_print/sla_test_utils.cpp b/tests/sla_print/sla_test_utils.cpp
index 07f19300c..a33a23c88 100644
--- a/tests/sla_print/sla_test_utils.cpp
+++ b/tests/sla_print/sla_test_utils.cpp
@@ -2,6 +2,7 @@
 #include "libslic3r/TriangleMeshSlicer.hpp"
 #include "libslic3r/SLA/AGGRaster.hpp"
 #include "libslic3r/SLA/DefaultSupportTree.hpp"
+#include "libslic3r/SLA/BranchingTreeSLA.hpp"
 
 #include <iomanip>
 
@@ -160,6 +161,11 @@ 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:;
     }