Fix horrible complexity of custom seam lookup (#5067)

- polygons are offset individually
- custom areas are kept separately for each PrintObject
- AABB tree is used to get logN lookup complexity

src/libslic3r/AABBTreeIndirect.hpp

@@ -726,6 +726,36 @@ inline bool is_any_triangle_in_radius(
     return hit_point.allFinite();
 }
 
+
+// Traverse the tree and return the index of an entity whose bounding box
+// contains a given point. Returns size_t(-1) when the point is outside.
+template<typename TreeType, typename VectorType>
+size_t get_candidate_idx(const TreeType& tree, const VectorType& v)
+{
+    if (tree.empty() || ! tree.node(0).bbox.contains(v))
+        return size_t(-1);
+
+    size_t node_idx = 0;
+    while (true) {
+        decltype(tree.node(node_idx)) node = tree.node(node_idx);
+        static_assert(std::is_reference<decltype(node)>::value,
+                      "Nodes shall be addressed by reference.");
+        assert(node.is_valid());
+        assert(node.bbox.contains(v));
+
+        if (! node.is_leaf()) {
+            if (tree.left_child(node_idx).bbox.contains(v))
+                node_idx = tree.left_child_idx(node_idx);
+            else if (tree.right_child(node_idx).bbox.contains(v))
+                node_idx = tree.right_child_idx(node_idx);
+            else
+                return size_t(-1);
+        } else
+            return node.idx;
+    }
+}
+
+
 } // namespace AABBTreeIndirect
 } // namespace Slic3r
 

src/libslic3r/GCode/SeamPlacer.cpp

@@ -191,19 +191,93 @@ void SeamPlacer::init(const Print& print)
 {
     m_enforcers.clear();
     m_blockers.clear();
-    //m_last_seam_position.clear();
     m_seam_history.clear();
+    m_po_list.clear();
 
-   for (const PrintObject* po : print.objects()) {
-       po->project_and_append_custom_facets(true, EnforcerBlockerType::ENFORCER, m_enforcers);
-       po->project_and_append_custom_facets(true, EnforcerBlockerType::BLOCKER, m_blockers);
-   }
     const std::vector<double>& nozzle_dmrs = print.config().nozzle_diameter.values;
     float max_nozzle_dmr = *std::max_element(nozzle_dmrs.begin(), nozzle_dmrs.end());
-   for (ExPolygons& explgs : m_enforcers)
+
-       explgs = Slic3r::offset_ex(explgs, scale_(max_nozzle_dmr));
+
-   for (ExPolygons& explgs : m_blockers)
+    std::vector<ExPolygons> temp_enf;
+    std::vector<ExPolygons> temp_blk;
+
+    for (const PrintObject* po : print.objects()) {
+        temp_enf.clear();
+        temp_blk.clear();
+        po->project_and_append_custom_facets(true, EnforcerBlockerType::ENFORCER, temp_enf);
+        po->project_and_append_custom_facets(true, EnforcerBlockerType::BLOCKER, temp_blk);
+
+        // Offset the triangles out slightly.
+        for (auto* custom_per_object : {&temp_enf, &temp_blk})
+            for (ExPolygons& explgs : *custom_per_object)
                 explgs = Slic3r::offset_ex(explgs, scale_(max_nozzle_dmr));
+
+//     FIXME: Offsetting should be done somehow cheaper, but following does not work
+//        for (auto* custom_per_object : {&temp_enf, &temp_blk}) {
+//            for (ExPolygons& plgs : *custom_per_object) {
+//                for (ExPolygon& plg : plgs) {
+//                    auto out = Slic3r::offset_ex(plg, scale_(max_nozzle_dmr));
+//                    plg = out.empty() ? ExPolygon() : out.front();
+//                    assert(out.empty() || out.size() == 1);
+//                }
+//            }
+//        }
+
+
+
+        // Remember this PrintObject and initialize a store of enforcers and blockers for it.
+        m_po_list.push_back(po);
+        size_t po_idx = m_po_list.size() - 1;
+        m_enforcers.emplace_back(std::vector<CustomTrianglesPerLayer>(temp_enf.size()));
+        m_blockers.emplace_back(std::vector<CustomTrianglesPerLayer>(temp_blk.size()));
+
+        // A helper class to store data to build the AABB tree from.
+        class CustomTriangleRef {
+        public:
+            CustomTriangleRef(size_t idx,
+                              Point&& centroid,
+                              BoundingBox&& bb)
+                : m_idx{idx}, m_centroid{centroid},
+                  m_bbox{AlignedBoxType(bb.min, bb.max)}
+            {}
+            size_t idx() const              { return m_idx;      }
+            const Point& centroid() const   { return m_centroid; }
+            const TreeType::BoundingBox& bbox() const { return m_bbox; }
+
+        private:
+            size_t m_idx;
+            Point m_centroid;
+            AlignedBoxType m_bbox;
+        };
+
+        // A lambda to extract the ExPolygons and save them into the member AABB tree.
+        // Will be called for enforcers and blockers separately.
+        auto add_custom = [](std::vector<ExPolygons>& src, std::vector<CustomTrianglesPerLayer>& dest) {
+            // Go layer by layer, and append all the ExPolygons into the AABB tree.
+            size_t layer_idx = 0;
+            for (ExPolygons& expolys_on_layer : src) {
+                CustomTrianglesPerLayer& layer_data = dest[layer_idx];
+                std::vector<CustomTriangleRef> triangles_data;
+                layer_data.polys.reserve(expolys_on_layer.size());
+                triangles_data.reserve(expolys_on_layer.size());
+
+                for (ExPolygon& expoly : expolys_on_layer) {
+                    if (expoly.empty())
+                        continue;
+                    layer_data.polys.emplace_back(std::move(expoly));
+                    triangles_data.emplace_back(layer_data.polys.size() - 1,
+                                                layer_data.polys.back().centroid(),
+                                                layer_data.polys.back().bounding_box());
+                }
+                // All polygons are saved, build the AABB tree for them.
+                layer_data.tree.build(std::move(triangles_data));
+                ++layer_idx;
+            }
+        };
+
+        add_custom(temp_enf, m_enforcers.at(po_idx));
+        add_custom(temp_blk, m_blockers.at(po_idx));
+    }
 }
 
 
@@ -216,7 +290,9 @@ Point SeamPlacer::get_seam(const size_t layer_idx, const SeamPosition seam_posit
     BoundingBox polygon_bb = polygon.bounding_box();
     const coord_t  nozzle_r   = coord_t(scale_(0.5 * nozzle_dmr) + 0.5);
 
-    if (this->is_custom_seam_on_layer(layer_idx)) {
+    size_t po_idx = std::find(m_po_list.begin(), m_po_list.end(), po) - m_po_list.begin();
+
+    if (this->is_custom_seam_on_layer(layer_idx, po_idx)) {
         // Seam enf/blockers can begin and end in between the original vertices.
         // Let add extra points in between and update the leghths.
         polygon.densify(MINIMAL_POLYGON_SIDE);
@@ -229,11 +305,10 @@ Point SeamPlacer::get_seam(const size_t layer_idx, const SeamPosition seam_posit
         if (seam_position == spAligned) {
             // Seam is aligned to the seam at the preceding layer.
             if (po != nullptr) {
-                std::optional<Point> pos = m_seam_history.get_last_seam(po, layer_idx, polygon_bb);
+                std::optional<Point> pos = m_seam_history.get_last_seam(m_po_list[po_idx], layer_idx, polygon_bb);
                 if (pos.has_value()) {
-                    //last_pos = m_last_seam_position[po];
                     last_pos = *pos;
-                    last_pos_weight = is_custom_enforcer_on_layer(layer_idx) ? 0.f : 1.f;
+                    last_pos_weight = is_custom_enforcer_on_layer(layer_idx, po_idx) ? 0.f : 1.f;
                 }
             }
         }
@@ -313,12 +388,12 @@ Point SeamPlacer::get_seam(const size_t layer_idx, const SeamPosition seam_posit
 
         // Custom seam. Huge (negative) constant penalty is applied inside
         // blockers (enforcers) to rule out points that should not win.
-        this->apply_custom_seam(polygon, penalties, lengths, layer_idx, seam_position);
+        this->apply_custom_seam(polygon, po_idx, penalties, lengths, layer_idx, seam_position);
 
         // Find a point with a minimum penalty.
         size_t idx_min = std::min_element(penalties.begin(), penalties.end()) - penalties.begin();
 
-        if (seam_position != spAligned || ! is_custom_enforcer_on_layer(layer_idx)) {
+        if (seam_position != spAligned || ! is_custom_enforcer_on_layer(layer_idx, po_idx)) {
             // Very likely the weight of idx_min is very close to the weight of last_pos_proj_idx.
             // In that case use last_pos_proj_idx instead.
             float penalty_aligned  = penalties[last_pos_proj_idx];
@@ -368,15 +443,15 @@ Point SeamPlacer::get_seam(const size_t layer_idx, const SeamPosition seam_posit
             // The other loops will get a seam close to the random point chosen
             // on the innermost contour.
             //FIXME This works correctly for inner contours first only.
-            last_pos = this->get_random_seam(layer_idx, polygon);
+            last_pos = this->get_random_seam(layer_idx, polygon, po_idx);
         }
-        if (loop.role() == erExternalPerimeter && is_custom_seam_on_layer(layer_idx)) {
+        if (loop.role() == erExternalPerimeter && is_custom_seam_on_layer(layer_idx, po_idx)) {
             // There is a possibility that the loop will be influenced by custom
             // seam enforcer/blocker. In this case do not inherit the seam
             // from internal loops (which may conflict with the custom selection
             // and generate another random one.
             bool saw_custom = false;
-            Point candidate = this->get_random_seam(layer_idx, polygon, &saw_custom);
+            Point candidate = this->get_random_seam(layer_idx, polygon, po_idx, &saw_custom);
             if (saw_custom)
                 last_pos = candidate;
         }
@@ -385,7 +460,7 @@ Point SeamPlacer::get_seam(const size_t layer_idx, const SeamPosition seam_posit
 }
 
 
-Point SeamPlacer::get_random_seam(size_t layer_idx, const Polygon& polygon,
+Point SeamPlacer::get_random_seam(size_t layer_idx, const Polygon& polygon, size_t po_idx,
                                   bool* saw_custom) const
 {
     // Parametrize the polygon by its length.
@@ -394,7 +469,7 @@ Point SeamPlacer::get_random_seam(size_t layer_idx, const Polygon& polygon,
     // Which of the points are inside enforcers/blockers?
     std::vector<size_t> enforcers_idxs;
     std::vector<size_t> blockers_idxs;
-    this->get_enforcers_and_blockers(layer_idx, polygon, enforcers_idxs, blockers_idxs);
+    this->get_enforcers_and_blockers(layer_idx, polygon, po_idx, enforcers_idxs, blockers_idxs);
 
     bool has_enforcers = ! enforcers_idxs.empty();
     bool has_blockers = ! blockers_idxs.empty();
@@ -444,32 +519,44 @@ Point SeamPlacer::get_random_seam(size_t layer_idx, const Polygon& polygon,
 
 void SeamPlacer::get_enforcers_and_blockers(size_t layer_id,
                              const Polygon& polygon,
+                             size_t po_idx,
                              std::vector<size_t>& enforcers_idxs,
                              std::vector<size_t>& blockers_idxs) const
 {
     enforcers_idxs.clear();
     blockers_idxs.clear();
 
-    // FIXME: This is quadratic and it should be improved, maybe by building
+    auto is_inside = [](const Point& pt,
-    // an AABB tree (or at least utilize bounding boxes).
+                        const CustomTrianglesPerLayer& custom_data) -> bool {
+        assert(! custom_data.polys.empty());
+        // Now ask the AABB tree which polygon we should check and check it.
+        size_t candidate = AABBTreeIndirect::get_candidate_idx(custom_data.tree, pt);
+        if (candidate != size_t(-1)
+         && custom_data.polys[candidate].contains(pt))
+            return true;
+        return false;
+    };
+
+    if (! m_enforcers[po_idx].empty()) {
+        const CustomTrianglesPerLayer& enforcers = m_enforcers[po_idx][layer_id];
+        if (! enforcers.polys.empty()) {
             for (size_t i=0; i<polygon.points.size(); ++i) {
-
+                if (is_inside(polygon.points[i], enforcers))
-        if (! m_enforcers.empty()) {
+                    enforcers_idxs.emplace_back(i);
-            assert(layer_id < m_enforcers.size());
+            }
-            for (const ExPolygon& explg : m_enforcers[layer_id]) {
-                if (explg.contains(polygon.points[i]))
-                    enforcers_idxs.push_back(i);
         }
     }
 
-        if (! m_blockers.empty()) {
+    if (! m_blockers[po_idx].empty()) {
-            assert(layer_id < m_blockers.size());
+        const CustomTrianglesPerLayer& blockers = m_blockers[po_idx][layer_id];
-            for (const ExPolygon& explg : m_blockers[layer_id]) {
+        if (! blockers.polys.empty()) {
-                if (explg.contains(polygon.points[i]))
+            for (size_t i=0; i<polygon.points.size(); ++i) {
-                    blockers_idxs.push_back(i);
+                if (is_inside(polygon.points[i], blockers))
+                    blockers_idxs.emplace_back(i);
             }
         }
     }
+
 }
 
 
@@ -543,17 +630,17 @@ static std::vector<size_t> find_enforcer_centers(const Polygon& polygon,
 
 
 
-void SeamPlacer::apply_custom_seam(const Polygon& polygon,
+void SeamPlacer::apply_custom_seam(const Polygon& polygon, size_t po_idx,
                                    std::vector<float>& penalties,
                                    const std::vector<float>& lengths,
                                    int layer_id, SeamPosition seam_position) const
 {
-    if (! is_custom_seam_on_layer(layer_id))
+    if (! is_custom_seam_on_layer(layer_id, po_idx))
         return;
 
     std::vector<size_t> enforcers_idxs;
     std::vector<size_t> blockers_idxs;
-    this->get_enforcers_and_blockers(layer_id, polygon, enforcers_idxs, blockers_idxs);
+    this->get_enforcers_and_blockers(layer_id, polygon, po_idx, enforcers_idxs, blockers_idxs);
 
     for (size_t i : enforcers_idxs) {
         assert(i < penalties.size());

src/libslic3r/GCode/SeamPlacer.hpp

@@ -3,9 +3,10 @@
 
 #include <optional>
 
-#include "libslic3r/ExPolygon.hpp"
+#include "libslic3r/Polygon.hpp"
 #include "libslic3r/PrintConfig.hpp"
 #include "libslic3r/BoundingBox.hpp"
+#include "libslic3r/AABBTreeIndirect.hpp"
 
 namespace Slic3r {
 
@@ -44,9 +45,20 @@ public:
                    coordf_t nozzle_diameter, const PrintObject* po,
                    bool was_clockwise, const EdgeGrid::Grid* lower_layer_edge_grid);
 
+    using TreeType = AABBTreeIndirect::Tree<2, coord_t>;
+    using AlignedBoxType = Eigen::AlignedBox<TreeType::CoordType, TreeType::NumDimensions>;
+
 private:
-    std::vector<ExPolygons> m_enforcers;
+
-    std::vector<ExPolygons> m_blockers;
+    struct CustomTrianglesPerLayer {
+        Polygons polys;
+        TreeType tree;
+    };
+
+
+    std::vector<std::vector<CustomTrianglesPerLayer>> m_enforcers;
+    std::vector<std::vector<CustomTrianglesPerLayer>> m_blockers;
+    std::vector<const PrintObject*> m_po_list;
 
     //std::map<const PrintObject*, Point>  m_last_seam_position;
     SeamHistory  m_seam_history;
@@ -54,32 +66,33 @@ private:
     // Get indices of points inside enforcers and blockers.
     void get_enforcers_and_blockers(size_t layer_id,
                                     const Polygon& polygon,
+                                    size_t po_id,
                                     std::vector<size_t>& enforcers_idxs,
                                     std::vector<size_t>& blockers_idxs) const;
 
     // Apply penalties to points inside enforcers/blockers.
-    void apply_custom_seam(const Polygon& polygon,
+    void apply_custom_seam(const Polygon& polygon, size_t po_id,
                            std::vector<float>& penalties,
                            const std::vector<float>& lengths,
                            int layer_id, SeamPosition seam_position) const;
 
     // Return random point of a polygon. The distribution will be uniform
     // along the contour and account for enforcers and blockers.
-    Point get_random_seam(size_t layer_idx, const Polygon& polygon,
+    Point get_random_seam(size_t layer_idx, const Polygon& polygon, size_t po_id,
                           bool* saw_custom = nullptr) const;
 
     // Is there any enforcer/blocker on this layer?
-    bool is_custom_seam_on_layer(size_t layer_id) const {
+    bool is_custom_seam_on_layer(size_t layer_id, size_t po_idx) const {
-        return is_custom_enforcer_on_layer(layer_id)
+        return is_custom_enforcer_on_layer(layer_id, po_idx)
-            || is_custom_blocker_on_layer(layer_id);
+            || is_custom_blocker_on_layer(layer_id, po_idx);
     }
 
-    bool is_custom_enforcer_on_layer(size_t layer_id) const {
+    bool is_custom_enforcer_on_layer(size_t layer_id, size_t po_idx) const {
-        return (! m_enforcers.empty() && ! m_enforcers[layer_id].empty());
+        return (! m_enforcers.at(po_idx).empty() && ! m_enforcers.at(po_idx)[layer_id].polys.empty());
     }
 
-    bool is_custom_blocker_on_layer(size_t layer_id) const {
+    bool is_custom_blocker_on_layer(size_t layer_id, size_t po_idx) const {
-        return (! m_blockers.empty() && ! m_blockers[layer_id].empty());
+        return (! m_blockers.at(po_idx).empty() && ! m_blockers.at(po_idx)[layer_id].polys.empty());
     }
 };
 

src/libslic3r/PrintObject.cpp

@@ -11,7 +11,6 @@
 #include "Slicing.hpp"
 #include "Tesselate.hpp"
 #include "Utils.hpp"
-#include "AABBTreeIndirect.hpp"
 #include "Fill/FillAdaptive.hpp"
 #include "Format/STL.hpp"