Make Tree-Support Deterministic (#15565)

* Make tree support deterministic without giving up its parallelism

* Break equal-distance ties in the tree support MST by coordinates

* test: cover the determinism this PR fixes

The MST unit tests here cover the tie-break, but the drop_nodes rework
has no test.

Adds two cases to the tree support suite. The thread-scheduling one
slices five configs twice each and compares the support point sequence,
which is what the node ordering moves. The MST tie one pins the branch
diameter and line width that carry Prim's equal-distance ties into the
toolpaths.

slice_with_tree_support takes an optional config list so the second case
can add the tree parameters it needs, and the double-slice comparison is
shared rather than written twice.

Both fail on main without this PR. The first passes from 60d1ceb580, the
second from e148865dd6.

---------

Co-authored-by: raistlin7447 <kris.austin@gmail.com>

src/libslic3r/MinimumSpanningTree.cpp

@@ -60,10 +60,15 @@ auto MinimumSpanningTree::prim(std::vector<Point> vertices) const -> AdjacencyGr
         //This search is O(V) right now, which can be made down to O(log(V)). This reduces the overall time complexity from O(V*V) to O(V*log(E)).
         //However that requires an implementation of a heap that supports the decreaseKey operation, which is not in the std library.
         //TODO: Implement this?
+        // Break equal-distance ties on coordinates: the map is keyed by address, so its
+        // iteration order (and therefore the first minimum) would otherwise depend on where
+        // the vertices were allocated.
         using MapValue = std::pair<const Point*, coordf_t>;
         const auto closest = std::min_element(smallest_distance.begin(), smallest_distance.end(),
                                               [](const MapValue& a, const MapValue& b) {
-                                                  return a.second < b.second;
+                                                  if (a.second != b.second)
+                                                      return a.second < b.second;
+                                                  return *a.first < *b.first;
                                               });
 
         //Add this point to the graph and remove it from the candidates.

src/libslic3r/Support/TreeSupport.cpp

@@ -2846,7 +2846,9 @@ void TreeSupport::drop_nodes()
             const MinimumSpanningTree& mst = spanning_trees[group_index];
             //In the first pass, merge all nodes that are close together.
             std::vector<std::pair<const Point, SupportNode*>> nodes_vec(nodes_this_part.begin(), nodes_this_part.end());
-            tbb::parallel_for_each(nodes_vec.begin(), nodes_vec.end(), [&](const std::pair<const Point, SupportNode*>& entry) {
+            // Sequential: nodes merge into and invalidate each other in place, so parallel execution
+            // makes the merge order (and thus the result) depend on thread scheduling.
+            std::for_each(nodes_vec.begin(), nodes_vec.end(), [&](const std::pair<const Point, SupportNode*>& entry) {
                 SupportNode* p_node = entry.second;
                 SupportNode& node = *p_node;
                 if (!p_node->valid)
@@ -2934,7 +2936,32 @@ void TreeSupport::drop_nodes()
             );
 
             //In the second pass, move all middle nodes.
-            tbb::parallel_for_each(nodes_vec.begin(), nodes_vec.end(), [&](const std::pair<const Point, SupportNode*>& entry) {
+            // Still parallel: this pass only reads other nodes. Side effects (invalidation, new
+            // nodes, contact_nodes/unsupported_branch_leaves updates) are recorded per node and
+            // applied afterwards in node order. Node creation must be deferred too, since
+            // SupportNode's constructor writes `parent->child = this` on other nodes.
+            struct PendingNode {
+                Point        position;
+                int          distance_to_top           = 0;
+                int          support_roof_layers_below = 0;
+                bool         to_buildplate             = false;
+                SupportNode *parent                    = nullptr;
+                bool         zero_max_move             = false;
+                bool         has_overhang              = false;
+                ExPolygon    overhang;
+                bool         clamp_radius              = false;
+                coordf_t     parent_radius             = 0;
+                double       dist_to_outer             = 0;
+            };
+            struct PassTwoResult {
+                bool                     invalidate       = false;
+                bool                     unsupported_leaf = false;
+                std::vector<PendingNode> pending;
+            };
+            std::vector<PassTwoResult> pass2_results(nodes_vec.size());
+            auto pass2_body = [&](size_t node_idx) {
+                const std::pair<const Point, SupportNode*>& entry = nodes_vec[node_idx];
+                PassTwoResult& pass2_out = pass2_results[node_idx];
 
                 SupportNode* p_node = entry.second;
                 const SupportNode& node = *p_node;
@@ -2949,14 +2976,16 @@ void TreeSupport::drop_nodes()
                     ExPolygons overhangs_next = diff_clipped({ node.overhang }, get_collision(0, obj_layer_nr_next));
                     for(auto& overhang:overhangs_next) {
                         Point        next_pt     = overhang.contour.centroid();
-                        SupportNode *next_node   = m_ts_data->create_node(next_pt, p_node->distance_to_top + 1, obj_layer_nr_next,
+                        PendingNode pending;
-                                                                          p_node->support_roof_layers_below - (p_node->distance_to_top >= 0 ? 1 : 0),
+                        pending.position                  = next_pt;
-                                                                          to_buildplate, p_node, print_z_next, height_next);
+                        pending.distance_to_top           = p_node->distance_to_top + 1;
-                        next_node->max_move_dist = 0;
+                        pending.support_roof_layers_below = p_node->support_roof_layers_below - (p_node->distance_to_top >= 0 ? 1 : 0);
-                        next_node->overhang = std::move(overhang);
+                        pending.to_buildplate             = to_buildplate;
-                        m_ts_data->m_mutex.lock();
+                        pending.parent                    = p_node;
-                        contact_nodes[layer_nr_next].emplace_back(next_node);
+                        pending.zero_max_move             = true;
-                        m_ts_data->m_mutex.unlock();
+                        pending.has_overhang              = true;
+                        pending.overhang                  = std::move(overhang);
+                        pass2_out.pending.emplace_back(std::move(pending));
 
                     }
                     return;
@@ -2973,17 +3002,17 @@ void TreeSupport::drop_nodes()
                     {
                         if (support_on_buildplate_only)
                         {
-                            unsupported_branch_leaves.push_front({ layer_nr, p_node });
+                            pass2_out.unsupported_leaf = true;
                         }
                         else {
-                            p_node->valid = false;
+                            pass2_out.invalidate = true;
                         }
                         return;
                     }
                     // if the link between parent and current is cut by contours, mark current as bottom contact node
                     if (p_node->parent && intersection_ln({p_node->position, p_node->parent->position}, layer_contours).empty()==false)
                     {
-                        p_node->valid = false;
+                        pass2_out.invalidate = true;
                         return;
                     }
                 }
@@ -3096,20 +3125,47 @@ void TreeSupport::drop_nodes()
                 }
                 auto              next_collision = get_collision(0, obj_layer_nr_next);
                 const bool   to_buildplate  = !is_inside_ex(m_ts_data->m_layer_outlines[obj_layer_nr_next], next_layer_vertex);
-                SupportNode *     next_node     = m_ts_data->create_node(next_layer_vertex, node.distance_to_top + 1, obj_layer_nr_next,
-                    node.support_roof_layers_below - (node.distance_to_top >= 0 ? 1 : 0),
-                    to_buildplate, p_node, print_z_next, height_next);
                 // don't increase radius if next node will collide partially with the object (STUDIO-7883)
-                to_outside             = projection_onto(next_collision, next_node->position);
+                to_outside             = projection_onto(next_collision, next_layer_vertex);
                 direction_to_outer     = to_outside - node.position;
                 double dist_to_outer   = unscale_(direction_to_outer.cast<double>().norm());
-                next_node->radius      = std::max(node.radius, std::min(next_node->radius, dist_to_outer));
+                PendingNode pending;
-                get_max_move_dist(next_node);
+                pending.position                  = next_layer_vertex;
-                m_ts_data->m_mutex.lock();
+                pending.distance_to_top           = node.distance_to_top + 1;
-                contact_nodes[layer_nr_next].push_back(next_node);
+                pending.support_roof_layers_below = node.support_roof_layers_below - (node.distance_to_top >= 0 ? 1 : 0);
-                m_ts_data->m_mutex.unlock();
+                pending.to_buildplate             = to_buildplate;
+                pending.parent                    = p_node;
+                pending.clamp_radius              = true;
+                pending.parent_radius             = node.radius;
+                pending.dist_to_outer             = dist_to_outer;
+                pass2_out.pending.emplace_back(std::move(pending));
+            };
+            tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes_vec.size()),
+                [&pass2_body](const tbb::blocked_range<size_t>& node_range) {
+                    for (size_t node_idx = node_range.begin(); node_idx < node_range.end(); ++ node_idx)
+                        pass2_body(node_idx);
+                });
+            // Apply the recorded side effects in node order.
+            for (size_t node_idx = 0; node_idx < nodes_vec.size(); ++ node_idx) {
+                PassTwoResult& pass2_out = pass2_results[node_idx];
+                for (PendingNode& pending : pass2_out.pending) {
+                    SupportNode* next_node = m_ts_data->create_node(pending.position, pending.distance_to_top, obj_layer_nr_next,
+                        pending.support_roof_layers_below, pending.to_buildplate, pending.parent, print_z_next, height_next);
+                    if (pending.zero_max_move)
+                        next_node->max_move_dist = 0;
+                    if (pending.has_overhang)
+                        next_node->overhang = std::move(pending.overhang);
+                    if (pending.clamp_radius) {
+                        next_node->radius = std::max(pending.parent_radius, std::min(next_node->radius, pending.dist_to_outer));
+                        get_max_move_dist(next_node);
+                    }
+                    contact_nodes[layer_nr_next].push_back(next_node);
+                }
+                if (pass2_out.unsupported_leaf)
+                    unsupported_branch_leaves.push_front({ layer_nr, nodes_vec[node_idx].second });
+                if (pass2_out.invalidate)
+                    nodes_vec[node_idx].second->valid = false;
             }
-            );
         }
 
 #ifdef SUPPORT_TREE_DEBUG_TO_SVG

src/libslic3r/Support/TreeSupport3D.cpp

@@ -2382,13 +2382,10 @@ static void merge_influence_areas(
     size_t num_buckets_initial;
     {
         // How many buckets per first merge iteration?
-        const size_t num_threads     = tbb::this_task_arena::max_concurrency();
+        // Fixed at 4: merging is not associative, so sizing buckets off max_concurrency() made
-        // 4 buckets per thread if possible,
+        // results depend on the core count of the slicing machine.
-        const size_t num_buckets_min = (input_size + 2) / 4;
+        const size_t bucket_size     = 4;
-        // 2 buckets per thread otherwise.
+        num_buckets_initial          = (input_size + 2) / 4;
-        const size_t num_buckets_max = input_size / 2;
-        num_buckets_initial          = num_buckets_min >= num_threads ? num_buckets_min : num_buckets_max;
-        const size_t bucket_size     = num_buckets_min >= num_threads ? 4 : 2;
         // Fill in the buckets.
         SupportElementMerging *it = influence_areas.data();
         // Reserve one more bucket to keep a single influence area which will not be merged in the first iteration.

tests/fff_print/test_tree_support.cpp

@@ -1,5 +1,7 @@
 #include <catch2/catch_all.hpp>
 
+#include <algorithm>
+
 #include "libslic3r/Layer.hpp"
 #include "libslic3r/TriangleMesh.hpp"
 
@@ -33,10 +35,13 @@ TriangleMesh scaled(TestMesh id, float scale)
     return mesh;
 }
 
+// `extra` is applied last, so a caller can add or override any key.
 void slice_with_tree_support(const TriangleMesh &mesh, Slic3r::Print &print, const char *style,
-                             int threshold_angle = 30, int build_plate_only = 0, int raft_layers = 0)
+                             int threshold_angle = 30, int build_plate_only = 0, int raft_layers = 0,
+                             std::initializer_list<Slic3r::ConfigBase::SetDeserializeItem> extra = {})
 {
-    Slic3r::Test::init_and_process_print({ mesh }, print, {
+    DynamicPrintConfig config = DynamicPrintConfig::full_print_config();
+    config.set_deserialize_strict({
         { "enable_support",              1 },
         { "support_type",                "tree(auto)" },
         { "support_style",               style },
@@ -45,6 +50,8 @@ void slice_with_tree_support(const TriangleMesh &mesh, Slic3r::Print &print, con
         { "raft_layers",                 raft_layers },
         { "layer_height",                0.2 },
     });
+    config.set_deserialize_strict(extra);
+    Slic3r::Test::init_and_process_print({ mesh }, print, config);
 }
 
 Points support_points(const Slic3r::Print &print)
@@ -63,6 +70,32 @@ size_t support_point_count(const TriangleMesh &mesh, const char *style, int thre
     return support_points(print).size();
 }
 
+// Index of the first differing point, or the common length when they match. An index keeps a
+// failure readable; comparing the vectors themselves dumps thousands of points.
+size_t first_difference(const Points &a, const Points &b)
+{
+    const size_t common = std::min(a.size(), b.size());
+    for (size_t i = 0; i < common; ++i)
+        if (a[i] != b[i])
+            return i;
+    return common;
+}
+
+// Slice `mesh` twice and require an identical support point sequence. Point counts and total
+// length are order insensitive, so the sequence is what a reordering shows up in.
+void sliced_twice_matches(const TriangleMesh &mesh, int build_plate_only, const char *style = "tree_slim",
+                          std::initializer_list<Slic3r::ConfigBase::SetDeserializeItem> extra = {})
+{
+    Slic3r::Print first_print, second_print;
+    slice_with_tree_support(mesh, first_print, style, 30, build_plate_only, 0, extra);
+    slice_with_tree_support(mesh, second_print, style, 30, build_plate_only, 0, extra);
+    const Points first  = support_points(first_print);
+    const Points second = support_points(second_print);
+    REQUIRE(first.size() > 1000); // without support the comparison below passes vacuously
+    REQUIRE(second.size() == first.size());
+    REQUIRE(first_difference(first, second) == first.size());
+}
+
 } // namespace
 
 TEST_CASE("Tree support is generated for an overhang and not for a plain cube", "[TreeSupport]")
@@ -123,3 +156,35 @@ TEST_CASE("A raft is still generated under tree support", "[TreeSupport]")
     // The raft goes under the object.
     REQUIRE(rafted_object->layers().front()->print_z > unrafted_object->layers().front()->print_z);
 }
+
+// drop_nodes() decides the node merges and spawns the next layer's nodes in parallel. Every one of
+// those decisions has to be applied in a fixed order, or the same model gives different branches on
+// each slice.
+TEST_CASE("Tree support toolpaths do not depend on thread scheduling", "[TreeSupport][Regression]")
+{
+    // Scaled up so that a layer holds enough nodes for the parallel range to be split. At stock
+    // size it stays in one chunk and the order never varies.
+    SECTION("overhang")            { sliced_twice_matches(scaled(TestMesh::overhang, 2.f), 0); }
+    SECTION("bridge with hole")    { sliced_twice_matches(scaled(TestMesh::bridge_with_hole, 3.f), 0); }
+    // Dropping every branch that cannot reach the bed leaves the survivors dense enough that the
+    // neighbour merge fires in bulk.
+    SECTION("on the build plate")  { sliced_twice_matches(scaled(TestMesh::overhang, 4.f), 1); }
+    // Branches resting on the model are what put nodes in a part group other than 0, which is the
+    // only way to reach the prune in the second pass. tree_hybrid additionally builds polygon
+    // nodes, so it is the only style that exercises the overhang merge.
+    SECTION("resting on the model") { sliced_twice_matches(two_tier_mesh(), 0); }
+    SECTION("hybrid on the model")  { sliced_twice_matches(two_tier_mesh(), 0, "tree_hybrid"); }
+}
+
+// Prim breaks equal-distance ties by heap address. A 1 mm branch diameter puts neighbours close
+// enough to tie, and an explicit line width pins max_move_dist, so the moved tie winner reaches
+// the support toolpaths.
+TEST_CASE("Tree support toolpaths do not depend on the MST tie order", "[TreeSupport][Regression]")
+{
+    sliced_twice_matches(two_tier_mesh(), 0, "tree_hybrid", {
+        { "tree_support_branch_diameter", 1.0 },
+        { "tree_support_branch_distance", 5.0 },
+        { "tree_support_branch_angle",    40 },
+        { "support_line_width",           0.4 },
+    });
+}

tests/libslic3r/CMakeLists.txt

@@ -29,6 +29,7 @@ add_executable(${_TEST_NAME}_tests
     test_polygon.cpp
     test_mutable_polygon.cpp
     test_mutable_priority_queue.cpp
+    test_minimum_spanning_tree.cpp
     test_nozzle_volume_type.cpp
     test_step.cpp
     test_stl.cpp

tests/libslic3r/test_minimum_spanning_tree.cpp

@@ -0,0 +1,66 @@
+#include <catch2/catch_all.hpp>
+
+#include <algorithm>
+
+#include "libslic3r/MinimumSpanningTree.hpp"
+#include "libslic3r/Point.hpp"
+
+using namespace Slic3r;
+
+// A 5x5 lattice: at every step of Prim's algorithm several candidates sit at the same
+// distance from the tree, so the tie-break decides the tree's shape.
+static std::vector<Point> lattice()
+{
+    std::vector<Point> vertices;
+    for (int y = 0; y < 5; ++y)
+        for (int x = 0; x < 5; ++x)
+            vertices.emplace_back(Point::new_scale(x, y));
+    return vertices;
+}
+
+static std::vector<Point> sorted_neighbours(const MinimumSpanningTree &mst, const Point &vertex)
+{
+    std::vector<Point> neighbours = mst.adjacent_nodes(vertex);
+    std::sort(neighbours.begin(), neighbours.end());
+    return neighbours;
+}
+
+TEST_CASE("Minimum spanning tree connects every vertex", "[MinimumSpanningTree]")
+{
+    const std::vector<Point> vertices = lattice();
+    const MinimumSpanningTree mst(vertices);
+
+    REQUIRE(mst.vertices().size() == vertices.size());
+    size_t adjacency_entries = 0;
+    for (const Point &vertex : vertices) {
+        const std::vector<Point> neighbours = mst.adjacent_nodes(vertex);
+        REQUIRE(! neighbours.empty());
+        adjacency_entries += neighbours.size();
+    }
+    // A tree on n vertices has n - 1 edges, each listed from both ends.
+    REQUIRE(adjacency_entries == 2 * (vertices.size() - 1));
+}
+
+TEST_CASE("Minimum spanning tree does not depend on the order of the non-root vertices", "[MinimumSpanningTree][Regression]")
+{
+    const std::vector<Point> vertices = lattice();
+    const MinimumSpanningTree reference(vertices);
+
+    // The root stays first: Prim's tree legitimately depends on where it starts.
+    // Every other order of the remaining vertices must give the same tree.
+    std::vector<std::vector<Point>> orders;
+    orders.emplace_back(vertices);
+    std::reverse(orders.back().begin() + 1, orders.back().end());
+    for (size_t shift = 1; shift + 1 < vertices.size(); ++shift) {
+        orders.emplace_back(vertices);
+        std::rotate(orders.back().begin() + 1, orders.back().begin() + 1 + shift, orders.back().end());
+    }
+
+    for (const std::vector<Point> &order : orders) {
+        const MinimumSpanningTree mst(order);
+        for (const Point &vertex : vertices) {
+            INFO("vertex " << vertex.x() << "," << vertex.y());
+            REQUIRE(sorted_neighbours(mst, vertex) == sorted_neighbours(reference, vertex));
+        }
+    }
+}