Measurement: optimize feature extraction:

- do not call ransac when not needed
- decrease number of ransac iterations
- precalculate inverse of transformation matrix
- pass the min error from ransac algorithm, do not calculate it again
This commit is contained in:
Lukas Matena 2022-12-05 14:50:25 +01:00
parent bed64cb7c7
commit a37f09edaa
5 changed files with 40 additions and 55 deletions

View File

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

View File

@ -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); Circled circle_taubin_newton(const Vec2ds& input, size_t cycles = 20);
// Find circle using RANSAC randomized algorithm. // 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. // Randomized algorithm by Emo Welzl, working with squared radii for efficiency. The returned circle radius is inflated by epsilon.
template<typename Vector, typename Points> template<typename Vector, typename Points>

View File

@ -14,7 +14,7 @@ namespace Measure {
constexpr double feature_hover_limit = 0.5; // how close to a feature the mouse must be to highlight it 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; Vec2ds out;
double z = 0.; double z = 0.;
@ -24,18 +24,17 @@ static std::pair<Vec3d, double> get_center_and_radius(const std::vector<Vec3d>&
out.emplace_back(pt_transformed.x(), pt_transformed.y()); 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) static std::array<Vec3d, 3> orthonormal_basis(const Vec3d& v)
{ {
@ -66,7 +65,6 @@ public:
float area; float area;
}; };
std::vector<SurfaceFeature> get_all_features() const;
std::optional<SurfaceFeature> get_feature(size_t face_idx, const Vec3d& point) const; std::optional<SurfaceFeature> get_feature(size_t face_idx, const Vec3d& point) const;
int get_num_of_planes() const; int get_num_of_planes() const;
const std::vector<int>& get_plane_triangle_indices(int idx) const; const std::vector<int>& get_plane_triangle_indices(int idx) const;
@ -268,6 +266,7 @@ void MeasuringImpl::extract_features()
q.setFromTwoVectors(plane.normal, Vec3d::UnitZ()); q.setFromTwoVectors(plane.normal, Vec3d::UnitZ());
Transform3d trafo = Transform3d::Identity(); Transform3d trafo = Transform3d::Identity();
trafo.rotate(q); trafo.rotate(q);
const Transform3d trafo_inv = trafo.inverse();
for (const std::vector<Vec3d>& border : plane.borders) { for (const std::vector<Vec3d>& border : plane.borders) {
if (border.size() <= 1) if (border.size() <= 1)
@ -275,28 +274,31 @@ void MeasuringImpl::extract_features()
bool done = false; 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);
// The whole border is one circle. Just add it into the list of features
// and we are done.
bool is_polygon = border.size()>4 && border.size()<=8; if (err < 0.05) {
bool lengths_match = std::all_of(border.begin()+2, border.end(), [is_polygon](const Vec3d& pt) { // The whole border is one circle. Just add it into the list of features
return Slic3r::is_approx((pt - *((&pt)-1)).squaredNorm(), (*((&pt)-1) - *((&pt)-2)).squaredNorm(), is_polygon ? 0.01 : 0.01); // and we are done.
});
if (lengths_match && (is_polygon || border.size() > 8)) { bool is_polygon = border.size()>4 && border.size()<=8;
if (is_polygon) { bool lengths_match = std::all_of(border.begin()+2, border.end(), [is_polygon](const Vec3d& pt) {
// This is a polygon, add the separate edges with the center. return Slic3r::is_approx((pt - *((&pt)-1)).squaredNorm(), (*((&pt)-1) - *((&pt)-2)).squaredNorm(), is_polygon ? 0.01 : 0.01);
for (int j=0; j<int(border.size()); ++j) });
plane.surface_features.emplace_back(SurfaceFeature(SurfaceFeatureType::Edge,
border[j==0 ? border.size()-1 : j-1], border[j], if (lengths_match && (is_polygon || border.size() > 8)) {
std::make_optional(center))); if (is_polygon) {
} else { // This is a polygon, add the separate edges with the center.
// The fit went well and it has more than 8 points - let's consider this a circle. for (int j=0; j<int(border.size()); ++j)
plane.surface_features.emplace_back(SurfaceFeature(SurfaceFeatureType::Circle, center, plane.normal, std::nullopt, radius)); plane.surface_features.emplace_back(SurfaceFeature(SurfaceFeatureType::Edge,
border[j==0 ? border.size()-1 : j-1], border[j],
std::make_optional(center)));
} else {
// The fit went well and it has more than 8 points - let's consider this a circle.
plane.surface_features.emplace_back(SurfaceFeature(SurfaceFeatureType::Circle, center, plane.normal, std::nullopt, radius));
}
done = true;
} }
done = true;
} }
} }
@ -391,11 +393,11 @@ void MeasuringImpl::extract_features()
} }
if (accept_circle) { 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 // Check that the fit went well. The tolerance is high, only to
// reject complete failures. // 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. // If the segment subtends less than 90 degrees, throw it away.
accept_circle &= single_circle_length / radius > 0.9*M_PI/2.; accept_circle &= single_circle_length / radius > 0.9*M_PI/2.;
@ -485,16 +487,6 @@ void MeasuringImpl::extract_features()
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;
}
@ -596,11 +588,6 @@ Measuring::Measuring(const indexed_triangle_set& its)
Measuring::~Measuring() {} 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 std::optional<SurfaceFeature> Measuring::get_feature(size_t face_idx, const Vec3d& point) const
{ {

View File

@ -94,9 +94,6 @@ public:
explicit Measuring(const indexed_triangle_set& its); explicit Measuring(const indexed_triangle_set& its);
~Measuring(); ~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 // Given a face_idx where the mouse cursor points, return a feature that
// should be highlighted (if any). // should be highlighted (if any).

View File

@ -654,7 +654,6 @@ void GLGizmoMeasure::on_render()
GLModel::Geometry init_data = init_plane_data(its, plane_triangles); GLModel::Geometry init_data = init_plane_data(its, plane_triangles);
m_plane.reset(); 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.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) }); m_raycasters.insert({ PLANE_ID, m_parent.add_raycaster_for_picking(SceneRaycaster::EType::Gizmo, PLANE_ID, *m_plane.mesh_raycaster) });