Merge branch 'lm_optimize_measurement'

src/libslic3r/Geometry/Circle.cpp

@@ -108,7 +108,7 @@ Circled circle_taubin_newton(const Vec2ds& input, size_t cycles)
     return out;
 }
 
-Circled circle_ransac(const Vec2ds& input, size_t iterations)
+Circled circle_ransac(const Vec2ds& input, size_t iterations, double* min_error)
 {
     if (input.size() < 3)
         return Circled::make_invalid();
@@ -132,6 +132,8 @@ Circled circle_ransac(const Vec2ds& input, size_t iterations)
             circle_best = c;
         }
     }
+    if (min_error)
+        *min_error = err_min;
     return circle_best;
 }
 

src/libslic3r/Geometry/Circle.hpp

@@ -102,7 +102,7 @@ inline Vec2d circle_center_taubin_newton(const Vec2ds& input, size_t cycles = 20
 Circled circle_taubin_newton(const Vec2ds& input, size_t cycles = 20);
 
 // Find circle using RANSAC randomized algorithm.
-Circled circle_ransac(const Vec2ds& input, size_t iterations = 20);
+Circled circle_ransac(const Vec2ds& input, size_t iterations = 20, double* min_error = nullptr);
 
 // Randomized algorithm by Emo Welzl, working with squared radii for efficiency. The returned circle radius is inflated by epsilon.
 template<typename Vector, typename Points>

src/libslic3r/Measure.cpp

@@ -8,13 +8,15 @@
 
 #include <numeric>
 
+#define DEBUG_EXTRACT_ALL_FEATURES_AT_ONCE 0
+
 namespace Slic3r {
 namespace Measure {
 
 
 constexpr double feature_hover_limit = 0.5; // how close to a feature the mouse must be to highlight it
 
-static std::pair<Vec3d, double> get_center_and_radius(const std::vector<Vec3d>& points, const Transform3d& trafo)
+static std::tuple<Vec3d, double, double> get_center_and_radius(const std::vector<Vec3d>& points, const Transform3d& trafo, const Transform3d& trafo_inv)
 {
     Vec2ds out;
     double z = 0.;
@@ -24,18 +26,17 @@ static std::pair<Vec3d, double> get_center_and_radius(const std::vector<Vec3d>&
         out.emplace_back(pt_transformed.x(), pt_transformed.y());
     }
 
-    auto circle = Geometry::circle_ransac(out, 20); // FIXME: iterations?
+    const int iter = points.size() < 10  ? 2 :
+                     points.size() < 100 ? 4 :
+                     6;
 
-    return std::make_pair(trafo.inverse() * Vec3d(circle.center.x(), circle.center.y(), z), circle.radius);
+    double error = std::numeric_limits<double>::max();
+    auto circle = Geometry::circle_ransac(out, iter, &error);
+    
+    return std::make_tuple(trafo.inverse() * Vec3d(circle.center.x(), circle.center.y(), z), circle.radius, error);
 }
 
-static bool circle_fit_is_ok(const std::vector<Vec3d>& pts, const Vec3d& center, double radius)
+
-{
-    for (const Vec3d& pt : pts)
-        if (std::abs((pt - center).norm() - radius) > 0.05)
-            return false;
-    return true;
-}
 
 static std::array<Vec3d, 3> orthonormal_basis(const Vec3d& v)
 {
@@ -64,17 +65,18 @@ public:
         std::vector<SurfaceFeature> surface_features;
         Vec3d normal;
         float area;
+        bool features_extracted = false;
     };
 
-    std::vector<SurfaceFeature> get_all_features() const;
+    std::optional<SurfaceFeature> get_feature(size_t face_idx, const Vec3d& point);
-    std::optional<SurfaceFeature> get_feature(size_t face_idx, const Vec3d& point) const;
+    int get_num_of_planes() const;
-    std::vector<std::vector<int>> get_planes_triangle_indices() const;
+    const std::vector<int>& get_plane_triangle_indices(int idx) const;
-    const std::vector<SurfaceFeature>& get_plane_features(unsigned int plane_id) const;
+    const std::vector<SurfaceFeature>& get_plane_features(unsigned int plane_id);
     const TriangleMesh& get_mesh() const;
 
 private:
     void update_planes();
-    void extract_features();
+    void extract_features(int plane_idx);
     
     std::vector<PlaneData> m_planes;
     std::vector<size_t>    m_face_to_plane;
@@ -90,7 +92,13 @@ MeasuringImpl::MeasuringImpl(const indexed_triangle_set& its)
 : m_mesh(its)
 {
     update_planes();
-    extract_features();
+
+    // Extracting features will be done as needed.
+    // To extract all planes at once, run the following:
+#if DEBUG_EXTRACT_ALL_FEATURES_AT_ONCE
+    for (int i=0; i<int(m_planes.size()); ++i)
+        extract_features(i);
+#endif
 }
 
 
@@ -252,14 +260,11 @@ void MeasuringImpl::update_planes()
 
 
 
-void MeasuringImpl::extract_features()
+void MeasuringImpl::extract_features(int plane_idx)
 {
-    std::vector<double> angles; // placed in outer scope to prevent reallocations
+    assert(! m_planes[plane_idx].features_extracted);
-    std::vector<double> lengths;
 
-
+    PlaneData& plane = m_planes[plane_idx];
-    for (int i=0; i<(int)m_planes.size(); ++i) {
-        PlaneData& plane = m_planes[i];
     plane.surface_features.clear();
     const Vec3d& normal = plane.normal;
 
@@ -267,6 +272,10 @@ void MeasuringImpl::extract_features()
     q.setFromTwoVectors(plane.normal, Vec3d::UnitZ());
     Transform3d trafo = Transform3d::Identity();
     trafo.rotate(q);
+    const Transform3d trafo_inv = trafo.inverse();
+
+    std::vector<double> angles; // placed in outer scope to prevent reallocations
+    std::vector<double> lengths;
 
     for (const std::vector<Vec3d>& border : plane.borders) {
         if (border.size() <= 1)
@@ -274,8 +283,10 @@ void MeasuringImpl::extract_features()
 
         bool done = false;
 
-            if (const auto& [center, radius] = get_center_and_radius(border, trafo);
+        if (border.size() > 4) {
-                (border.size()>4) && circle_fit_is_ok(border, center, radius)) {
+            const auto& [center, radius, err] = get_center_and_radius(border, trafo, trafo_inv);
+
+            if (err < 0.05) {
                 // The whole border is one circle. Just add it into the list of features
                 // and we are done.
 
@@ -298,6 +309,7 @@ void MeasuringImpl::extract_features()
                     done = true;
                 }
             }
+        }
 
         if (! done) {
             // In this case, the border is not a circle and may contain circular
@@ -390,11 +402,11 @@ void MeasuringImpl::extract_features()
                         }
 
                         if (accept_circle) {
-                                const auto& [center, radius] = get_center_and_radius(single_circle, trafo);
+                            const auto& [center, radius, err] = get_center_and_radius(single_circle, trafo, trafo_inv);
 
                             // Check that the fit went well. The tolerance is high, only to
                             // reject complete failures.
-                                accept_circle &= circle_fit_is_ok(single_circle, center, radius);
+                            accept_circle &= err < 0.05;
 
                             // If the segment subtends less than 90 degrees, throw it away.
                             accept_circle &= single_circle_length / radius > 0.9*M_PI/2.;
@@ -475,23 +487,12 @@ void MeasuringImpl::extract_features()
     }
     cog /= double(counter);
     plane.surface_features.emplace_back(SurfaceFeature(SurfaceFeatureType::Plane,
-            plane.normal, cog, std::optional<Vec3d>(), i + 0.0001));
+        plane.normal, cog, std::optional<Vec3d>(), plane_idx + 0.0001));
 
     plane.borders.clear();
     plane.borders.shrink_to_fit();
-    }
-}
 
-
+    plane.features_extracted = true;
-
-std::vector<SurfaceFeature> MeasuringImpl::get_all_features() const
-{
-    std::vector<SurfaceFeature> features;
-    //PlaneData& plane = m_planes[0];
-    for (const PlaneData& plane : m_planes)    
-        for (const SurfaceFeature& feature : plane.surface_features)
-            features.emplace_back(feature);
-    return features;        
 }
 
 
@@ -499,13 +500,18 @@ std::vector<SurfaceFeature> MeasuringImpl::get_all_features() const
 
 
 
-std::optional<SurfaceFeature> MeasuringImpl::get_feature(size_t face_idx, const Vec3d& point) const
+
+
+std::optional<SurfaceFeature> MeasuringImpl::get_feature(size_t face_idx, const Vec3d& point)
 {
     if (face_idx >= m_face_to_plane.size())
         return std::optional<SurfaceFeature>();
 
     const PlaneData& plane = m_planes[m_face_to_plane[face_idx]];
 
+    if (! plane.features_extracted)
+        extract_features(m_face_to_plane[face_idx]);
+    
     size_t closest_feature_idx = size_t(-1);
     double min_dist = std::numeric_limits<double>::max();
 
@@ -554,17 +560,24 @@ std::optional<SurfaceFeature> MeasuringImpl::get_feature(size_t face_idx, const
 
 
 
-std::vector<std::vector<int>> MeasuringImpl::get_planes_triangle_indices() const
+int MeasuringImpl::get_num_of_planes() const
 {
-    std::vector<std::vector<int>> out;
+    return (m_planes.size());
-    for (const PlaneData& plane : m_planes)
-        out.emplace_back(plane.facets);        
-    return out;
 }
 
-const std::vector<SurfaceFeature>& MeasuringImpl::get_plane_features(unsigned int plane_id) const
+
+
+const std::vector<int>& MeasuringImpl::get_plane_triangle_indices(int idx) const
+{
+    assert(idx >= 0 && idx < int(m_planes.size()));
+    return m_planes[idx].facets;
+}
+
+const std::vector<SurfaceFeature>& MeasuringImpl::get_plane_features(unsigned int plane_id)
 {
     assert(plane_id < m_planes.size());
+    if (! m_planes[plane_id].features_extracted)
+        extract_features(plane_id);
     return m_planes[plane_id].surface_features;
 }
 
@@ -590,11 +603,6 @@ Measuring::Measuring(const indexed_triangle_set& its)
 Measuring::~Measuring() {}
 
 
-std::vector<SurfaceFeature> Measuring::get_all_features() const
-{
-    return priv->get_all_features();
-}
-
 
 std::optional<SurfaceFeature> Measuring::get_feature(size_t face_idx, const Vec3d& point) const
 {
@@ -602,10 +610,15 @@ std::optional<SurfaceFeature> Measuring::get_feature(size_t face_idx, const Vec3
 }
 
 
-
+int Measuring::get_num_of_planes() const
-std::vector<std::vector<int>> Measuring::get_planes_triangle_indices() const
 {
-    return priv->get_planes_triangle_indices();
+    return priv->get_num_of_planes();
+}
+
+
+const std::vector<int>& Measuring::get_plane_triangle_indices(int idx) const
+{
+    return priv->get_plane_triangle_indices(idx);
 }
 
 const std::vector<SurfaceFeature>& Measuring::get_plane_features(unsigned int plane_id) const

src/libslic3r/Measure.hpp

@@ -94,18 +94,16 @@ public:
     explicit Measuring(const indexed_triangle_set& its);
     ~Measuring();
 
-    // Return a reference to a list of all features identified on the its.
-    // Use only for debugging. Expensive, do not call often.
-    std::vector<SurfaceFeature> get_all_features() const;
 
     // Given a face_idx where the mouse cursor points, return a feature that
     // should be highlighted (if any).
     std::optional<SurfaceFeature> get_feature(size_t face_idx, const Vec3d& point) const;
 
-    // Returns a list of triangle indices for each identified plane. Each
+    // Return total number of planes.
-    // Plane object contains an index into this vector. Expensive, do not
+    int get_num_of_planes() const;
-    // call too often.
+
-    std::vector<std::vector<int>> get_planes_triangle_indices() const;
+    // Returns a list of triangle indices for given plane.
+    const std::vector<int>& get_plane_triangle_indices(int idx) const;
 
     // Returns the surface features of the plane with the given index
     const std::vector<SurfaceFeature>& get_plane_features(unsigned int plane_id) const;

src/slic3r/GUI/Gizmos/GLGizmoMeasure.cpp

@@ -93,11 +93,8 @@ static std::string center_on_feature_type_as_string(Measure::SurfaceFeatureType
     return ret;
 }
 
-static GLModel::Geometry init_plane_data(const indexed_triangle_set& its, const std::vector<std::vector<int>>& planes_triangles, int idx)
+static GLModel::Geometry init_plane_data(const indexed_triangle_set& its, const std::vector<int>& triangle_indices)
 {
-    assert(0 <= idx && idx < (int)planes_triangles.size());
-    const std::vector<int>& triangle_indices = planes_triangles[idx];
-
     GLModel::Geometry init_data;
     init_data.format = { GUI::GLModel::Geometry::EPrimitiveType::Triangles, GLModel::Geometry::EVertexLayout::P3N3 };
     unsigned int i = 0;
@@ -653,11 +650,10 @@ void GLGizmoMeasure::on_render()
                     if (m_last_plane_idx != idx) {
                         m_last_plane_idx = idx;
                         const indexed_triangle_set& its = m_measuring->get_mesh().its;
-                        const std::vector<std::vector<int>> planes_triangles = m_measuring->get_planes_triangle_indices();
+                        const std::vector<int>& plane_triangles = m_measuring->get_plane_triangle_indices(idx);
-                        GLModel::Geometry init_data = init_plane_data(its, planes_triangles, idx);
+                        GLModel::Geometry init_data = init_plane_data(its, plane_triangles);
                         m_plane.reset();
                         m_plane.mesh_raycaster = std::make_unique<MeshRaycaster>(std::make_shared<const TriangleMesh>(init_data.get_as_indexed_triangle_set()));
-                        m_plane.model.init_from(std::move(init_data));
                     }
 
                     m_raycasters.insert({ PLANE_ID, m_parent.add_raycaster_for_picking(SceneRaycaster::EType::Gizmo, PLANE_ID, *m_plane.mesh_raycaster) });
@@ -1041,10 +1037,9 @@ void GLGizmoMeasure::update_if_needed()
 {
     auto update_plane_models_cache = [this](const indexed_triangle_set& its) {
         m_plane_models_cache.clear();
-        const std::vector<std::vector<int>> planes_triangles = m_measuring->get_planes_triangle_indices();
+        for (int idx = 0; idx < m_measuring->get_num_of_planes(); ++idx) {
-        for (int idx = 0; idx < (int)planes_triangles.size(); ++idx) {
             m_plane_models_cache.emplace_back(GLModel());
-            GLModel::Geometry init_data = init_plane_data(its, planes_triangles, idx);
+            GLModel::Geometry init_data = init_plane_data(its, m_measuring->get_plane_triangle_indices(idx));
             m_plane_models_cache.back().init_from(std::move(init_data));
         }
     };