refactoring raycaster, gathering only hitpoints and their normals.

src/libslic3r/GCode.cpp

@@ -2586,7 +2586,7 @@ std::string GCode::extrude_loop(ExtrusionLoop loop, std::string description, dou
         loop.split_at(last_pos, false);
     }
     else
-        m_seam_placer.place_seam(loop, this->last_pos(), m_config.external_perimeters_first,
+        m_seam_placer.place_seam(m_layer->object(), loop, m_layer->slice_z, this->last_pos(), m_config.external_perimeters_first,
                                  EXTRUDER_CONFIG(nozzle_diameter), lower_layer_edge_grid ? lower_layer_edge_grid->get() : nullptr);
 
     // clip the path to avoid the extruder to get exactly on the first point of the loop;
@@ -2715,10 +2715,10 @@ std::string GCode::extrude_perimeters(const Print &print, const std::vector<Obje
             if (m_layer->lower_layer && ! lower_layer_edge_grid)
                 lower_layer_edge_grid = calculate_layer_edge_grid(*m_layer->lower_layer);
 
-            m_seam_placer.plan_perimeters(std::vector<const ExtrusionEntity*>(region.perimeters.begin(), region.perimeters.end()),
+//            m_seam_placer.plan_perimeters(std::vector<const ExtrusionEntity*>(region.perimeters.begin(), region.perimeters.end()),
-                *m_layer, m_config.seam_position, this->last_pos(), EXTRUDER_CONFIG(nozzle_diameter),
+//                *m_layer, m_config.seam_position, this->last_pos(), EXTRUDER_CONFIG(nozzle_diameter),
-                (m_layer == NULL ? nullptr : m_layer->object()),
+//                (m_layer == NULL ? nullptr : m_layer->object()),
-                (lower_layer_edge_grid ? lower_layer_edge_grid.get() : nullptr));
+//                (lower_layer_edge_grid ? lower_layer_edge_grid.get() : nullptr));
 
             for (const ExtrusionEntity* ee : region.perimeters)
                 gcode += this->extrude_entity(*ee, "perimeter", -1., &lower_layer_edge_grid);

src/libslic3r/GCode/SeamPlacerNG.cpp

@@ -6,7 +6,6 @@
 #include <boost/log/trivial.hpp>
 #include <random>
 #include <algorithm>
-#include <atomic>
 
 #include "libslic3r/ExtrusionEntity.hpp"
 #include "libslic3r/Print.hpp"
@@ -16,7 +15,7 @@
 #include "libslic3r/SVG.hpp"
 #include "libslic3r/Layer.hpp"
 
-//TODO remove
+// TODO remove
 #include <boost/nowide/cstdio.hpp>
 
 namespace Slic3r {
@@ -24,60 +23,58 @@ namespace Slic3r {
 namespace SeamPlacerImpl {
 
 /// Coordinate frame
-class Frame
+class Frame {
-{
 public:
-    Frame()
+    Frame() {
-    {
         mX = Vec3d(1, 0, 0);
         mY = Vec3d(0, 1, 0);
         mZ = Vec3d(0, 0, 1);
     }
 
-    Frame(const Vec3d &x, const Vec3d &y, const Vec3d &z)
+    Frame(const Vec3d &x, const Vec3d &y, const Vec3d &z) :
-        : mX(x), mY(y), mZ(z)
+            mX(x), mY(y), mZ(z) {
-    {}
+    }
 
-    void set_from_z(const Vec3d &z)
+    void set_from_z(const Vec3d &z) {
-    {
         mZ = z.normalized();
         Vec3d tmpZ = mZ;
-        Vec3d tmpX = (std::abs(tmpZ.x()) > 0.99f) ? Vec3d(0, 1, 0) :
+        Vec3d tmpX = (std::abs(tmpZ.x()) > 0.99f) ? Vec3d(0, 1, 0) : Vec3d(1, 0, 0);
-                                                    Vec3d(1, 0, 0);
         mY = (tmpZ.cross(tmpX)).normalized();
         mX = mY.cross(tmpZ);
     }
 
-    Vec3d to_world(const Vec3d &a) const
+    Vec3d to_world(const Vec3d &a) const {
-    {
         return a.x() * mX + a.y() * mY + a.z() * mZ;
     }
 
-    Vec3d to_local(const Vec3d &a) const
+    Vec3d to_local(const Vec3d &a) const {
-    {
         return Vec3d(mX.dot(a), mY.dot(a), mZ.dot(a));
     }
 
-    const Vec3d &binormal() const { return mX; }
+    const Vec3d& binormal() const {
+        return mX;
+    }
 
-    const Vec3d &tangent() const { return mY; }
+    const Vec3d& tangent() const {
+        return mY;
+    }
 
-    const Vec3d &normal() const { return mZ; }
+    const Vec3d& normal() const {
+        return mZ;
+    }
 
 private:
     Vec3d mX, mY, mZ;
 };
 
-Vec3d sample_sphere_uniform(const Vec2f &samples)
+Vec3d sample_sphere_uniform(const Vec2f &samples) {
-{
     float term_one = 2.0f * M_PIf32 * samples.x();
     float term_two = 2.0f * sqrt(samples.y() - samples.y() * samples.y());
     return {cos(term_one) * term_two, sin(term_one) * term_two,
         1.0f - 2.0f * samples.y()};
 }
 
-Vec3d sample_power_cosine_hemisphere(const Vec2f &samples, float power)
+Vec3d sample_power_cosine_hemisphere(const Vec2f &samples, float power) {
-{
     const float term1 = 2.f * M_PIf32 * samples.x();
     const float term2 = pow(samples.y(), 1.f / (power + 1.f));
     const float term3 = sqrt(1.f - term2 * term2);
@@ -85,55 +82,60 @@ Vec3d sample_power_cosine_hemisphere(const Vec2f &samples, float power)
     return Vec3d(cos(term1) * term3, sin(term1) * term3, term2);
 }
 
-void raycast_visibility(
+void resolve_geometry_hit(const igl::Hit &hitpoint,
-    size_t                                  ray_count,
-    const AABBTreeIndirect::Tree<3, float> &raycasting_tree,
         const indexed_triangle_set &triangles,
-    const KDTreeIndirect<3, coordf_t, KDTreeCoordinateFunctor>
+        const KDTreeIndirect<3, coordf_t, KDTreeCoordinateFunctor> &perimeter_points_tree,
-                               &perimeter_points_tree,
+        const std::vector<SeamCandidate> &perimeter_points,
-    std::vector<SeamCandidate> &perimeter_points)
+        std::vector<std::atomic<size_t>> &visibility_counters) {
-{
+    auto face = triangles.indices[hitpoint.id];
+    auto edge1 = triangles.vertices[face[1]] - triangles.vertices[face[0]];
+    auto edge2 = triangles.vertices[face[2]] - triangles.vertices[face[0]];
+
+    Vec3d hit_pos = (triangles.vertices[face[0]] + edge1 * hitpoint.u + edge2 * hitpoint.v).cast<double>();
+    auto perimeter_point_index = find_closest_point(perimeter_points_tree, hit_pos);
+    auto perimeter_point_pos = perimeter_points[perimeter_point_index].position;
+    auto dist_squared = (perimeter_point_pos - hit_pos).squaredNorm();
+
+    visibility_counters[perimeter_point_index].fetch_add(1.0 / dist_squared, std::memory_order_relaxed);
+}
+
+std::vector<HitInfo> raycast_visibility(size_t ray_count,
+        const AABBTreeIndirect::Tree<3, float> &raycasting_tree,
+        const indexed_triangle_set &triangles) {
     auto bbox = raycasting_tree.node(0).bbox;
     Vec3d vision_sphere_center = bbox.center().cast<double>();
-    float vision_sphere_raidus = (bbox.sizes().maxCoeff() *
+    float vision_sphere_raidus = (bbox.sizes().maxCoeff() * 0.55); // 0.5 (half) covers whole object,
-                                  0.55); // 0.5 (half) covers whole object,
                                                                    // 0.05 added to avoid corner cases
 
     // Prepare random samples per ray
     std::random_device rnd_device;
-    std::mt19937                          mersenne_engine{rnd_device()};
+    std::mt19937 mersenne_engine { rnd_device() };
-    std::uniform_real_distribution<float> dist{0, 1};
+    std::uniform_real_distribution<float> dist { 0, 1 };
 
     auto gen = [&dist, &mersenne_engine]() {
         return Vec2f(dist(mersenne_engine), dist(mersenne_engine));
     };
 
-    BOOST_LOG_TRIVIAL(debug) << "PM: generate random samples: start";
+    BOOST_LOG_TRIVIAL(debug)
+    << "PM: generate random samples: start";
     std::vector<Vec2f> global_dir_random_samples(ray_count);
-    generate(begin(global_dir_random_samples), end(global_dir_random_samples),
+    generate(begin(global_dir_random_samples), end(global_dir_random_samples), gen);
-             gen);
     std::vector<Vec2f> local_dir_random_samples(ray_count);
-    generate(begin(local_dir_random_samples), end(local_dir_random_samples),
+    generate(begin(local_dir_random_samples), end(local_dir_random_samples), gen);
-             gen);
 
-    BOOST_LOG_TRIVIAL(debug) << "PM: generate random samples: end";
+    BOOST_LOG_TRIVIAL(debug)
-
+    << "PM: generate random samples: end";
-    std::vector<std::atomic<size_t>> visibility_counters(
-        perimeter_points.size());
 
     BOOST_LOG_TRIVIAL(debug)
     << "PM: raycast visibility for " << ray_count << " rays: start";
     // raycast visibility
-    tbb::parallel_for(
+    std::vector<HitInfo> hit_points = tbb::parallel_reduce(tbb::blocked_range<size_t>(0, ray_count),
-        tbb::blocked_range<size_t>(0, ray_count),
+            std::vector<HitInfo> { },
-        [&](tbb::blocked_range<size_t> r) {
+            [&](tbb::blocked_range<size_t> r, std::vector<HitInfo> init) {
                 for (size_t index = r.begin(); index < r.end(); ++index) {
-                Vec3d global_ray_dir = sample_sphere_uniform(
+                    Vec3d global_ray_dir = sample_sphere_uniform(global_dir_random_samples[index]);
-                    global_dir_random_samples[index]);
+                    Vec3d ray_origin = (vision_sphere_center - global_ray_dir * vision_sphere_raidus);
-                Vec3d ray_origin = (vision_sphere_center -
+                    Vec3d local_dir = sample_power_cosine_hemisphere(local_dir_random_samples[index], 1.5);
-                                    global_ray_dir * vision_sphere_raidus);
-                Vec3d local_dir  = sample_power_cosine_hemisphere(
-                     local_dir_random_samples[index], 1.0);
 
                     Frame f;
                     f.set_from_z(global_ray_dir);
@@ -142,71 +144,57 @@ void raycast_visibility(
                     igl::Hit hitpoint;
                     // FIXME: This query will not compile for float ray origin and
                     // direction for some reason
-                auto hit = AABBTreeIndirect::intersect_ray_first_hit(
+                    auto hit = AABBTreeIndirect::intersect_ray_first_hit(triangles.vertices,
-                    triangles.vertices, triangles.indices, raycasting_tree,
+                            triangles.indices, raycasting_tree, ray_origin, final_ray_dir, hitpoint);
-                    ray_origin, final_ray_dir, hitpoint);
 
                     if (hit) {
                         auto face = triangles.indices[hitpoint.id];
-                    auto edge1 = triangles.vertices[face[1]] -
+                        auto edge1 = triangles.vertices[face[1]] - triangles.vertices[face[0]];
-                                 triangles.vertices[face[0]];
+                        auto edge2 = triangles.vertices[face[2]] - triangles.vertices[face[0]];
-                    auto edge2 = triangles.vertices[face[2]] -
-                                 triangles.vertices[face[0]];
 
-                    Vec3d hit_pos = (triangles.vertices[face[0]] +
+                        Vec3d hit_pos = (triangles.vertices[face[0]] + edge1 * hitpoint.u + edge2 * hitpoint.v).cast<double>();
-                                     edge1 * hitpoint.u + edge2 * hitpoint.v)
+                        Vec3d surface_normal = edge1.cross(edge2).cast<double>();
-                                        .cast<double>();
 
-                    auto perimeter_point_index =
+                        init.push_back(HitInfo { hit_pos, surface_normal });
-                        find_closest_point(perimeter_points_tree, hit_pos);
-
-                    visibility_counters[perimeter_point_index]
-                        .fetch_add(1, std::memory_order_relaxed);
                     }
                 }
-        });
+                return init;
+            },
+            [](std::vector<HitInfo> left, std::vector<HitInfo> right) {
+                left.insert(left.end(), right.begin(), right.end());
+                return left;
+            }
+            );
 
     BOOST_LOG_TRIVIAL(debug)
     << "PM: raycast visibility for " << ray_count << " rays: end";
 
-    BOOST_LOG_TRIVIAL(debug)
+//TODO disable, only debug code
-        << "PM: assign visibility to perimenter points : start";
+//#ifdef 0
-
-    tbb::parallel_for(tbb::blocked_range<size_t>(0, perimeter_points.size()),
-                      [&](tbb::blocked_range<size_t> r) {
-                          for (size_t index = r.begin(); index < r.end();
-                               ++index) {
-                              perimeter_points[index].visibility =
-                                  visibility_counters[index];
-                          }
-                      });
-    BOOST_LOG_TRIVIAL(debug)
-        << "PM: assign visibility to perimenter points : end";
-
     its_write_obj(triangles, "triangles.obj");
 
     Slic3r::CNumericLocalesSetter locales_setter;
     FILE *fp = boost::nowide::fopen("perimeter.obj", "w");
     if (fp == nullptr) {
-        BOOST_LOG_TRIVIAL(error) << "Couldn't open "
+        BOOST_LOG_TRIVIAL(error)
-                                 << "perimeter.obj"
+        << "Couldn't open " << "perimeter.obj" << " for writing";
-                                 << " for writing";
     }
 
-    for (size_t i = 0; i < perimeter_points.size(); ++i)
+    for (size_t i = 0; i < hit_points.size(); ++i)
-        fprintf(fp, "v %f %f %f %zu\n", perimeter_points[i].position[0],
+        fprintf(fp, "v %f %f %f \n", hit_points[i].m_position[0], hit_points[i].m_position[1],
-                perimeter_points[i].position[1],
+                hit_points[i].m_position[2]);
-                perimeter_points[i].position[2],
-                perimeter_points[i].visibility);
     fclose(fp);
+//#endif
+
+    return hit_points;
 }
 
 } // namespace SeamPlacerImpl
 
-void SeamPlacer::init(const Print &print)
+void SeamPlacer::init(const Print &print) {
-{
     using namespace SeamPlacerImpl;
-    seam_candidates_trees.clear();
+    m_perimeter_points_trees_per_object.clear();
+    m_perimeter_points_per_object.clear();
 
     for (const PrintObject *po : print.objects()) {
         BOOST_LOG_TRIVIAL(debug)
@@ -216,51 +204,98 @@ void SeamPlacer::init(const Print &print)
         auto triangle_set = po->model_object()->raw_indexed_triangle_set();
         its_transform(triangle_set, obj_transform);
 
-        auto raycasting_tree =
+        auto raycasting_tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(triangle_set.vertices,
-            AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(
+                triangle_set.indices);
-                triangle_set.vertices, triangle_set.indices);
 
-        BOOST_LOG_TRIVIAL(debug) << "PM: build AABB tree for raycasting: end";
+        BOOST_LOG_TRIVIAL(debug)
+        << "PM: build AABB tree for raycasting: end";
+
+        std::vector<HitInfo> hit_points = raycast_visibility(ray_count_per_object, raycasting_tree, triangle_set);
+        HitInfoCoordinateFunctor hit_points_functor { &hit_points };
+        KDTreeIndirect<3, coordf_t, HitInfoCoordinateFunctor> hit_points_tree { hit_points_functor, hit_points.size() };
 
         BOOST_LOG_TRIVIAL(debug)
         << "PM: gather and build KD tree with seam candidates: start";
         // gather seam candidates (perimeter points)
-        auto seam_candidates = tbb::parallel_reduce(
+        m_perimeter_points_per_object[po] = tbb::parallel_reduce(tbb::blocked_range<size_t>(0, po->layers().size()),
-            tbb::blocked_range<size_t>(0, po->layers().size()),
+                std::vector<SeamCandidate> { }, [&](tbb::blocked_range<size_t> r, std::vector<SeamCandidate> init) {
-            std::vector<SeamCandidate>{},
-            [&](tbb::blocked_range<size_t> r,
-                std::vector<SeamCandidate> init) {
                     for (size_t index = r.begin(); index < r.end(); ++index) {
                         const auto layer = po->layers()[index];
                         auto unscaled_z = layer->slice_z;
-                    for (Points points : layer->lslices) {
+                        for (const ExPolygon &expoly : layer->lslices) {
-                        for (Point point : points) {
+                            // Contour - insert first point marked as Polygon
-                            auto unscaled_p        = unscale(point);
+                            // start, then insert rest sequentially.
-                            auto unscaled_position = Vec3d{unscaled_p.x(),
+                            {
-                                                           unscaled_p.y(),
+                                auto unscaled_p = unscale(expoly.contour[0]);
-                                                           unscaled_z};
+                                init.emplace_back(Vec3d { unscaled_p.x(), unscaled_p.y(), unscaled_z }, true);
-                            init.emplace_back(unscaled_position);
+                            }
+                            for (size_t index = 1; index < expoly.contour.size(); ++index) {
+                                auto unscaled_p = unscale(expoly.contour[index]);
+                                init.emplace_back(Vec3d { unscaled_p.x(), unscaled_p.y(), unscaled_z });
+                            }
+
+                            for (const Polygon &hole : expoly.holes) {
+                                // Perform the same for each hole
+                                {
+                                    auto unscaled_p = unscale(hole[0]);
+                                    init.emplace_back(Vec3d { unscaled_p.x(), unscaled_p.y(), unscaled_z }, true);
+                                }
+                                for (size_t index = 1; index < hole.size(); ++index) {
+                                    auto unscaled_p = unscale(hole[index]);
+                                    init.emplace_back(Vec3d { unscaled_p.x(), unscaled_p.y(), unscaled_z });
+                                }
                             }
                         }
                     }
                     return init;
                 },
-            [](std::vector<SeamCandidate> prev,
+                [](std::vector<SeamCandidate> prev, std::vector<SeamCandidate> next) {
-               std::vector<SeamCandidate> next) {
                     prev.insert(prev.end(), next.begin(), next.end());
                     return prev;
                 });
+        auto &perimeter_points = m_perimeter_points_per_object[po];
         // Build KD tree with seam candidates
-        auto functor = KDTreeCoordinateFunctor{&seam_candidates};
+        auto functor = KDTreeCoordinateFunctor { &perimeter_points };
-        auto perimeter_points_tree = KDTreeIndirect<
+        m_perimeter_points_trees_per_object.emplace(std::piecewise_construct, std::forward_as_tuple(po),
-            3, coordf_t, KDTreeCoordinateFunctor>{functor,
+                std::forward_as_tuple(functor, m_perimeter_points_per_object[po].size()));
-                                                  seam_candidates.size()};
+        SeamPlacer::PointTree &perimeter_points_tree = m_perimeter_points_trees_per_object.find(po)->second;
 
         BOOST_LOG_TRIVIAL(debug)
         << "PM: gather and build KD tree with seam candidates: end";
-
-        raycast_visibility(100000, raycasting_tree, triangle_set,
-                           perimeter_points_tree, seam_candidates);
     }
 }
+
+void SeamPlacer::place_seam(const PrintObject *po, ExtrusionLoop &loop, coordf_t unscaled_z, const Point &last_pos,
+        bool external_first,
+        double nozzle_diameter, const EdgeGrid::Grid *lower_layer_edge_grid) {
+    assert(m_perimeter_points_trees_per_object.find(po) != nullptr);
+    assert(m_perimeter_points_per_object.find(po) != nullptr);
+    const auto &perimeter_points_tree = m_perimeter_points_trees_per_object.find(po)->second;
+    const auto &perimeter_points = m_perimeter_points_per_object.find(po)->second;
+
+    Points loop_points { };
+    loop.collect_points(loop_points);
+
+    // vector of pairs: first-> index into perimeter points, second-> index into loop points
+    auto closest_perimeter_point_indices = std::vector<std::pair<size_t, size_t>>(loop_points.size());
+    for (size_t p_index = 0; p_index < loop_points.size(); ++p_index) {
+        auto unscaled_p = unscale(loop_points[p_index]);
+        closest_perimeter_point_indices.emplace_back(find_closest_point(perimeter_points_tree,
+                Vec3d { unscaled_p.x(), unscaled_p.y(), unscaled_z }), p_index);
+    }
+
+    std::sort(closest_perimeter_point_indices.begin(), closest_perimeter_point_indices.end(),
+            [&](const std::pair<size_t, size_t> left, const std::pair<size_t, size_t> right) {
+                return perimeter_points[left.first].visibility < perimeter_points[right.first].visibility;
+            });
+
+    loop.split_at_vertex(loop_points[closest_perimeter_point_indices[0].second]);
+//
+//    Point seam = last_pos;
+//    if (!loop.split_at_vertex(seam))
+//        // The point is not in the original loop.
+//        // Insert it.
+//        loop.split_at(seam, true);
+}
+
 } // namespace Slic3r

src/libslic3r/GCode/SeamPlacerNG.hpp

@@ -24,57 +24,59 @@ class Grid;
 
 namespace SeamPlacerImpl {
 
-struct SeamCandidate
+struct SeamCandidate {
-{
+    explicit SeamCandidate(const Vec3d &pos) :
-    SeamCandidate(const Vec3d &pos) : position(pos), visibility(0) {}
+            position(pos), visibility(0.0), polygon_start(false) {
-    // TODO is there some equivalent for Vec with coordf_t type?
+    }
+    SeamCandidate(const Vec3d &pos, bool polygon_start) :
+            position(pos), visibility(0.0), polygon_start(polygon_start) {
+    }
     Vec3d position;
-    size_t visibility;
+    float visibility;
+    bool polygon_start;
 };
 
-struct KDTreeCoordinateFunctor
+struct HitInfo {
-{
+    Vec3d m_position;
-    KDTreeCoordinateFunctor(std::vector<SeamCandidate> *seam_candidates)
+    Vec3d m_surface_normal;
-        : seam_candidates(seam_candidates)
+};
-    {}
+
+
+struct KDTreeCoordinateFunctor {
+    KDTreeCoordinateFunctor(std::vector<SeamCandidate> *seam_candidates) :
+            seam_candidates(seam_candidates) {
+    }
     std::vector<SeamCandidate> *seam_candidates;
-    float                       operator()(size_t index, size_t dim) const
+    float operator()(size_t index, size_t dim) const {
-    {
         return seam_candidates->operator[](index).position[dim];
     }
 };
+
+struct HitInfoCoordinateFunctor {
+    HitInfoCoordinateFunctor(std::vector<HitInfo> *hit_points) :
+        m_hit_points(hit_points) {
+    }
+    std::vector<HitInfo> *m_hit_points;
+    float operator()(size_t index, size_t dim) const {
+        return m_hit_points->operator[](index).m_position[dim];
+    }
+};
 } // namespace SeamPlacerImpl
 
-class SeamPlacer
+class SeamPlacer {
-{
+    using PointTree =
+    KDTreeIndirect<3, coordf_t, SeamPlacerImpl::KDTreeCoordinateFunctor>;
+    const size_t ray_count_per_object = 100000;
+
 public:
-    std::vector<
+    std::unordered_map<const PrintObject*, PointTree> m_perimeter_points_trees_per_object;
-        KDTreeIndirect<3, coordf_t, SeamPlacerImpl::KDTreeCoordinateFunctor>>
+    std::unordered_map<const PrintObject*, std::vector<SeamPlacerImpl::SeamCandidate>> m_perimeter_points_per_object;
-        seam_candidates_trees;
 
     void init(const Print &print);
 
-    void plan_perimeters(const std::vector<const ExtrusionEntity *> perimeters,
+    void place_seam(const PrintObject *po, ExtrusionLoop &loop, coordf_t unscaled_z, const Point &last_pos,
-                         const Layer                               &layer,
-                         SeamPosition          seam_position,
-                         Point                 last_pos,
-                         coordf_t              nozzle_dmr,
-                         const PrintObject    *po,
-                         const EdgeGrid::Grid *lower_layer_edge_grid)
-    {}
-
-    void place_seam(ExtrusionLoop        &loop,
-                    const Point          &last_pos,
             bool external_first,
-                    double                nozzle_diameter,
+            double nozzle_diameter, const EdgeGrid::Grid *lower_layer_edge_grid);
-                    const EdgeGrid::Grid *lower_layer_edge_grid)
-    {
-        Point seam = last_pos;
-        if (!loop.split_at_vertex(seam))
-            // The point is not in the original loop.
-            // Insert it.
-            loop.split_at(seam, true);
-    }
 };
 
 } // namespace Slic3r

src/libslic3r/KDTreeIndirect.hpp

@@ -231,6 +231,53 @@ size_t find_closest_point(const KDTreeIndirectType& kdtree, const PointType& poi
 	return find_closest_point(kdtree, point, [](size_t) { return true; });
 }
 
+// Find a nearby points (spherical neighbourhood) using Euclidian metrics.
+template<typename KDTreeIndirectType, typename PointType, typename FilterFn>
+std::vector<size_t> find_nearby_points(const KDTreeIndirectType &kdtree, const PointType &center,
+        typename KDTreeIndirectType::CoordType& max_distance, FilterFn filter)
+        {
+    using CoordType = typename KDTreeIndirectType::CoordType;
+
+    struct Visitor {
+        const KDTreeIndirectType &kdtree;
+        const PointType &center;
+        const CoordType &max_distance_squared;
+        const FilterFn filter;
+        std::vector<size_t> result;
+
+        Visitor(const KDTreeIndirectType &kdtree, const PointType &center, const CoordType &max_distance,
+                FilterFn filter) :
+                kdtree(kdtree), center(center), max_distance_squared(max_distance*max_distance), filter(filter) {
+        }
+        unsigned int operator()(size_t idx, size_t dimension) {
+            if (this->filter(idx)) {
+                auto dist = CoordType(0);
+                for (size_t i = 0; i < KDTreeIndirectType::NumDimensions; ++i) {
+                    CoordType d = center[i] - kdtree.coordinate(idx, i);
+                    dist += d * d;
+                }
+                if (dist < max_distance_squared) {
+                    result.push_back(idx);
+                }
+            }
+            return kdtree.descent_mask(center[dimension], max_distance_squared, idx, dimension);
+        }
+    } visitor(kdtree, center, max_distance, filter);
+
+    kdtree.visit(visitor);
+    return visitor.result;
+}
+
+template<typename KDTreeIndirectType, typename PointType>
+std::vector<size_t> find_nearby_points(const KDTreeIndirectType &kdtree, const PointType &center,
+        typename KDTreeIndirectType::CoordType& max_distance)
+        {
+    return find_nearby_points(kdtree, center, max_distance, [](size_t) {
+        return true;
+    });
+}
+
+
 } // namespace Slic3r
 
 #endif /* slic3r_KDTreeIndirect_hpp_ */