Measure: Initial porting of Measure Gizmo

src/libslic3r/CMakeLists.txt

@@ -208,6 +208,9 @@ set(lisbslic3r_sources
     ModelArrange.cpp
     MultiMaterialSegmentation.cpp
     MultiMaterialSegmentation.hpp
+    Measure.hpp
+    Measure.cpp
+    MeasureUtils.hpp
     CustomGCode.cpp
     CustomGCode.hpp
     Arrange.hpp
@@ -307,6 +310,7 @@ set(lisbslic3r_sources
     Surface.hpp
     SurfaceCollection.cpp
     SurfaceCollection.hpp
+    SurfaceMesh.hpp
     SVG.cpp
     SVG.hpp
     Technologies.hpp

src/libslic3r/Geometry.hpp

@@ -1,3 +1,18 @@
+///|/ Copyright (c) Prusa Research 2016 - 2023 Vojtěch Bubník @bubnikv, Enrico Turri @enricoturri1966, Tomáš Mészáros @tamasmeszaros, Lukáš Matěna @lukasmatena, Filip Sykala @Jony01, Lukáš Hejl @hejllukas
+///|/ Copyright (c) 2017 Eyal Soha @eyal0
+///|/ Copyright (c) Slic3r 2013 - 2016 Alessandro Ranellucci @alranel
+///|/
+///|/ ported from lib/Slic3r/Geometry.pm:
+///|/ Copyright (c) Prusa Research 2017 - 2022 Vojtěch Bubník @bubnikv
+///|/ Copyright (c) Slic3r 2011 - 2015 Alessandro Ranellucci @alranel
+///|/ Copyright (c) 2013 Jose Luis Perez Diez
+///|/ Copyright (c) 2013 Anders Sundman
+///|/ Copyright (c) 2013 Jesse Vincent
+///|/ Copyright (c) 2012 Mike Sheldrake @mesheldrake
+///|/ Copyright (c) 2012 Mark Hindess
+///|/
+///|/ PrusaSlicer is released under the terms of the AGPLv3 or higher
+///|/
 #ifndef slic3r_Geometry_hpp_
 #define slic3r_Geometry_hpp_
 

src/libslic3r/Geometry/Circle.cpp

@@ -1,3 +1,7 @@
+///|/ Copyright (c) Prusa Research 2021 - 2022 Lukáš Matěna @lukasmatena, Filip Sykala @Jony01, Vojtěch Bubník @bubnikv
+///|/
+///|/ PrusaSlicer is released under the terms of the AGPLv3 or higher
+///|/
 #include "Circle.hpp"
 
 #include "../Polygon.hpp"
@@ -108,7 +112,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 +136,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

@@ -1,3 +1,7 @@
+///|/ Copyright (c) Prusa Research 2021 - 2022 Lukáš Matěna @lukasmatena, Filip Sykala @Jony01, Vojtěch Bubník @bubnikv, Enrico Turri @enricoturri1966
+///|/
+///|/ PrusaSlicer is released under the terms of the AGPLv3 or higher
+///|/
 #ifndef slic3r_Geometry_Circle_hpp_
 #define slic3r_Geometry_Circle_hpp_
 
@@ -102,7 +106,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.hpp

@@ -0,0 +1,200 @@
+///|/ Copyright (c) Prusa Research 2022 - 2023 Lukáš Matěna @lukasmatena, Enrico Turri @enricoturri1966, Vojtěch Bubník @bubnikv
+///|/
+///|/ PrusaSlicer is released under the terms of the AGPLv3 or higher
+///|/
+#ifndef Slic3r_Measure_hpp_
+#define Slic3r_Measure_hpp_
+
+#include <optional>
+#include <memory>
+
+#include "Point.hpp"
+
+
+struct indexed_triangle_set;
+
+
+
+namespace Slic3r {
+
+class TriangleMesh;
+
+namespace Measure {
+
+
+enum class SurfaceFeatureType : int {
+    Undef  = 0,
+    Point  = 1 << 0,
+    Edge   = 1 << 1,
+    Circle = 1 << 2,
+    Plane  = 1 << 3
+};
+
+class SurfaceFeature {
+public:
+    SurfaceFeature(SurfaceFeatureType type, const Vec3d& pt1, const Vec3d& pt2, std::optional<Vec3d> pt3 = std::nullopt, double value = 0.0)
+        : m_type(type), m_pt1(pt1), m_pt2(pt2), m_pt3(pt3), m_value(value) {}
+
+    explicit SurfaceFeature(const Vec3d& pt)
+    : m_type{SurfaceFeatureType::Point}, m_pt1{pt} {}
+
+    // Get type of this feature.
+    SurfaceFeatureType get_type() const { return m_type; }
+
+    // For points, return the point.
+    Vec3d get_point() const { assert(m_type == SurfaceFeatureType::Point); return m_pt1; }
+
+    // For edges, return start and end.
+    std::pair<Vec3d, Vec3d> get_edge() const { assert(m_type == SurfaceFeatureType::Edge); return std::make_pair(m_pt1, m_pt2); }    
+
+    // For circles, return center, radius and normal.
+    std::tuple<Vec3d, double, Vec3d> get_circle() const { assert(m_type == SurfaceFeatureType::Circle); return std::make_tuple(m_pt1, m_value, m_pt2); }
+
+    // For planes, return index into vector provided by Measuring::get_plane_triangle_indices, normal and point.
+    std::tuple<int, Vec3d, Vec3d> get_plane() const { assert(m_type == SurfaceFeatureType::Plane); return std::make_tuple(int(m_value), m_pt1, m_pt2); }
+
+    // For anything, return an extra point that should also be considered a part of this.
+    std::optional<Vec3d> get_extra_point() const { assert(m_type != SurfaceFeatureType::Undef); return m_pt3; }
+
+    bool operator == (const SurfaceFeature& other) const {
+        if (this->m_type != other.m_type) return false;
+        switch (this->m_type)
+        {
+        case SurfaceFeatureType::Undef: { break; }
+        case SurfaceFeatureType::Point: { return (this->m_pt1.isApprox(other.m_pt1)); }
+        case SurfaceFeatureType::Edge: {
+            return (this->m_pt1.isApprox(other.m_pt1) && this->m_pt2.isApprox(other.m_pt2)) ||
+                   (this->m_pt1.isApprox(other.m_pt2) && this->m_pt2.isApprox(other.m_pt1));
+        }
+        case SurfaceFeatureType::Plane:
+        case SurfaceFeatureType::Circle: {
+            return (this->m_pt1.isApprox(other.m_pt1) && this->m_pt2.isApprox(other.m_pt2) && std::abs(this->m_value - other.m_value) < EPSILON);
+        }
+        }
+
+        return false;
+    }
+
+    bool operator != (const SurfaceFeature& other) const {
+        return !operator == (other);
+    }
+
+private:
+    SurfaceFeatureType m_type{ SurfaceFeatureType::Undef };
+    Vec3d m_pt1{ Vec3d::Zero() };
+    Vec3d m_pt2{ Vec3d::Zero() };
+    std::optional<Vec3d> m_pt3;
+    double m_value{ 0.0 };
+};
+
+
+
+class MeasuringImpl;
+
+
+class Measuring {
+public:
+    // Construct the measurement object on a given its.
+    explicit Measuring(const indexed_triangle_set& its);
+    ~Measuring();
+
+
+    // 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;
+
+    // Return total number of planes.
+    int get_num_of_planes() 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;
+
+    // Returns the mesh used for measuring
+    const indexed_triangle_set& get_its() const;
+
+private: 
+    std::unique_ptr<MeasuringImpl> priv;
+};
+
+
+struct DistAndPoints {
+    DistAndPoints(double dist_, Vec3d from_, Vec3d to_) : dist(dist_), from(from_), to(to_) {}
+    double dist;
+    Vec3d from;
+    Vec3d to;
+};
+
+struct AngleAndEdges {
+    AngleAndEdges(double angle_, const Vec3d& center_, const std::pair<Vec3d, Vec3d>& e1_, const std::pair<Vec3d, Vec3d>& e2_, double radius_, bool coplanar_)
+        : angle(angle_), center(center_), e1(e1_), e2(e2_), radius(radius_), coplanar(coplanar_) {}
+    double angle;
+    Vec3d center;
+    std::pair<Vec3d, Vec3d> e1;
+    std::pair<Vec3d, Vec3d> e2;
+    double radius;
+    bool coplanar;
+
+    static const AngleAndEdges Dummy;
+};
+
+struct MeasurementResult {
+    std::optional<AngleAndEdges> angle;
+    std::optional<DistAndPoints> distance_infinite;
+    std::optional<DistAndPoints> distance_strict;
+    std::optional<Vec3d>  distance_xyz;
+
+    bool has_distance_data() const {
+        return distance_infinite.has_value() || distance_strict.has_value();
+    }
+
+    bool has_any_data() const {
+        return angle.has_value() || distance_infinite.has_value() || distance_strict.has_value() || distance_xyz.has_value();
+    }
+};
+
+// Returns distance/angle between two SurfaceFeatures.
+MeasurementResult get_measurement(const SurfaceFeature& a, const SurfaceFeature& b, const Measuring* measuring = nullptr);
+
+inline Vec3d edge_direction(const Vec3d& from, const Vec3d& to) { return (to - from).normalized(); }
+inline Vec3d edge_direction(const std::pair<Vec3d, Vec3d>& e) { return edge_direction(e.first, e.second); }
+inline Vec3d edge_direction(const SurfaceFeature& edge) {
+    assert(edge.get_type() == SurfaceFeatureType::Edge);
+    return edge_direction(edge.get_edge());
+}
+
+inline Vec3d plane_normal(const SurfaceFeature& plane) {
+    assert(plane.get_type() == SurfaceFeatureType::Plane);
+    return std::get<1>(plane.get_plane());
+}
+
+inline bool are_parallel(const Vec3d& v1, const Vec3d& v2) { return std::abs(std::abs(v1.dot(v2)) - 1.0) < EPSILON; }
+inline bool are_perpendicular(const Vec3d& v1, const Vec3d& v2) { return std::abs(v1.dot(v2)) < EPSILON; }
+
+inline bool are_parallel(const std::pair<Vec3d, Vec3d>& e1, const std::pair<Vec3d, Vec3d>& e2) {
+    return are_parallel(e1.second - e1.first, e2.second - e2.first);
+}
+inline bool are_parallel(const SurfaceFeature& f1, const SurfaceFeature& f2) {
+    if (f1.get_type() == SurfaceFeatureType::Edge && f2.get_type() == SurfaceFeatureType::Edge)
+        return are_parallel(edge_direction(f1), edge_direction(f2));
+    else if (f1.get_type() == SurfaceFeatureType::Edge && f2.get_type() == SurfaceFeatureType::Plane)
+        return are_perpendicular(edge_direction(f1), plane_normal(f2));
+    else
+        return false;
+}
+
+inline bool are_perpendicular(const SurfaceFeature& f1, const SurfaceFeature& f2) {
+    if (f1.get_type() == SurfaceFeatureType::Edge && f2.get_type() == SurfaceFeatureType::Edge)
+        return are_perpendicular(edge_direction(f1), edge_direction(f2));
+    else if (f1.get_type() == SurfaceFeatureType::Edge && f2.get_type() == SurfaceFeatureType::Plane)
+        return are_parallel(edge_direction(f1), plane_normal(f2));
+    else
+        return false;
+}
+
+} // namespace Measure
+} // namespace Slic3r
+
+#endif // Slic3r_Measure_hpp_

src/libslic3r/MeasureUtils.hpp

@@ -0,0 +1,390 @@
+///|/ Copyright (c) Prusa Research 2022 Enrico Turri @enricoturri1966
+///|/
+///|/ PrusaSlicer is released under the terms of the AGPLv3 or higher
+///|/
+#ifndef Slic3r_MeasureUtils_hpp_
+#define Slic3r_MeasureUtils_hpp_
+
+#include <initializer_list>
+
+namespace Slic3r {
+namespace Measure {
+
+// Utility class used to calculate distance circle-circle
+// Adaptation of code found in: 
+// https://github.com/davideberly/GeometricTools/blob/master/GTE/Mathematics/Polynomial1.h
+
+class Polynomial1
+{
+public:
+    Polynomial1(std::initializer_list<double> values)
+    {
+        // C++ 11 will call the default constructor for
+        // Polynomial1<Real> p{}, so it is guaranteed that
+        // values.size() > 0.
+        m_coefficient.resize(values.size());
+        std::copy(values.begin(), values.end(), m_coefficient.begin());
+        EliminateLeadingZeros();
+    }
+
+    // Construction and destruction.  The first constructor creates a
+    // polynomial of the specified degree but sets all coefficients to
+    // zero (to ensure initialization).  You are responsible for setting
+    // the coefficients, presumably with the degree-term set to a nonzero
+    // number.  In the second constructor, the degree is the number of
+    // initializers plus 1, but then adjusted so that coefficient[degree]
+    // is not zero (unless all initializer values are zero).
+    explicit Polynomial1(uint32_t degree)
+        : m_coefficient(static_cast<size_t>(degree) + 1, 0.0)
+    {}
+
+    // Eliminate any leading zeros in the polynomial, except in the case
+    // the degree is 0 and the coefficient is 0.  The elimination is
+    // necessary when arithmetic operations cause a decrease in the degree
+    // of the result.  For example, (1 + x + x^2) + (1 + 2*x - x^2) =
+    // (2 + 3*x).  The inputs both have degree 2, so the result is created
+    // with degree 2.  After the addition we find that the degree is in
+    // fact 1 and resize the array of coefficients.  This function is
+    // called internally by the arithmetic operators, but it is exposed in
+    // the public interface in case you need it for your own purposes.
+    void EliminateLeadingZeros()
+    {
+        const size_t size = m_coefficient.size();
+        if (size > 1) {
+            const double zero = 0.0;
+            int32_t leading;
+            for (leading = static_cast<int32_t>(size) - 1; leading > 0; --leading) {
+                if (m_coefficient[leading] != zero)
+                    break;
+            }
+
+            m_coefficient.resize(++leading);
+        }
+    }
+
+    // Set all coefficients to the specified value.
+    void SetCoefficients(double value)
+    {
+        std::fill(m_coefficient.begin(), m_coefficient.end(), value);
+    }
+
+    inline uint32_t GetDegree() const
+    {
+        // By design, m_coefficient.size() > 0.
+        return static_cast<uint32_t>(m_coefficient.size() - 1);
+    }
+
+    inline const double& operator[](uint32_t i) const { return m_coefficient[i]; }
+    inline double& operator[](uint32_t i) { return m_coefficient[i]; }
+
+    // Evaluate the polynomial.  If the polynomial is invalid, the
+    // function returns zero.
+    double operator()(double t) const
+    {
+        int32_t i = static_cast<int32_t>(m_coefficient.size());
+        double result = m_coefficient[--i];
+        for (--i; i >= 0; --i) {
+            result *= t;
+            result += m_coefficient[i];
+        }
+        return result;
+    }
+
+protected:
+    // The class is designed so that m_coefficient.size() >= 1.
+    std::vector<double> m_coefficient;
+};
+
+inline Polynomial1 operator * (const Polynomial1& p0, const Polynomial1& p1)
+{
+    const uint32_t p0Degree = p0.GetDegree();
+    const uint32_t p1Degree = p1.GetDegree();
+    Polynomial1 result(p0Degree + p1Degree);
+    result.SetCoefficients(0.0);
+    for (uint32_t i0 = 0; i0 <= p0Degree; ++i0) {
+        for (uint32_t i1 = 0; i1 <= p1Degree; ++i1) {
+            result[i0 + i1] += p0[i0] * p1[i1];
+        }
+    }
+    return result;
+}
+
+inline Polynomial1 operator + (const Polynomial1& p0, const Polynomial1& p1)
+{
+    const uint32_t p0Degree = p0.GetDegree();
+    const uint32_t p1Degree = p1.GetDegree();
+    uint32_t i;
+    if (p0Degree >= p1Degree) {
+        Polynomial1 result(p0Degree);
+        for (i = 0; i <= p1Degree; ++i) {
+            result[i] = p0[i] + p1[i];
+        }
+        for (/**/; i <= p0Degree; ++i) {
+            result[i] = p0[i];
+        }
+        result.EliminateLeadingZeros();
+        return result;
+    }
+    else {
+        Polynomial1 result(p1Degree);
+        for (i = 0; i <= p0Degree; ++i) {
+            result[i] = p0[i] + p1[i];
+        }
+        for (/**/; i <= p1Degree; ++i) {
+            result[i] = p1[i];
+        }
+        result.EliminateLeadingZeros();
+        return result;
+    }
+}
+
+inline Polynomial1 operator - (const Polynomial1& p0, const Polynomial1& p1)
+{
+    const uint32_t p0Degree = p0.GetDegree();
+    const uint32_t p1Degree = p1.GetDegree();
+    uint32_t i;
+    if (p0Degree >= p1Degree) {
+        Polynomial1 result(p0Degree);
+        for (i = 0; i <= p1Degree; ++i) {
+            result[i] = p0[i] - p1[i];
+        }
+        for (/**/; i <= p0Degree; ++i) {
+            result[i] = p0[i];
+        }
+        result.EliminateLeadingZeros();
+        return result;
+    }
+    else {
+        Polynomial1 result(p1Degree);
+        for (i = 0; i <= p0Degree; ++i) {
+            result[i] = p0[i] - p1[i];
+        }
+        for (/**/; i <= p1Degree; ++i) {
+            result[i] = -p1[i];
+        }
+        result.EliminateLeadingZeros();
+        return result;
+    }
+}
+
+inline Polynomial1 operator * (double scalar, const Polynomial1& p)
+{
+    const uint32_t degree = p.GetDegree();
+    Polynomial1 result(degree);
+    for (uint32_t i = 0; i <= degree; ++i) {
+        result[i] = scalar * p[i];
+    }
+    return result;
+}
+
+// Utility class used to calculate distance circle-circle
+// Adaptation of code found in: 
+// https://github.com/davideberly/GeometricTools/blob/master/GTE/Mathematics/RootsPolynomial.h
+
+class RootsPolynomial
+{
+public:
+    // General equations: sum_{i=0}^{d} c(i)*t^i = 0.  The input array 'c'
+    // must have at least d+1 elements and the output array 'root' must
+    // have at least d elements.
+
+    // Find the roots on (-infinity,+infinity).
+    static int32_t Find(int32_t degree, const double* c, uint32_t maxIterations, double* roots)
+    {
+        if (degree >= 0 && c != nullptr) {
+            const double zero = 0.0;
+            while (degree >= 0 && c[degree] == zero) {
+                --degree;
+            }
+
+            if (degree > 0) {
+                // Compute the Cauchy bound.
+                const double one = 1.0;
+                const double invLeading = one / c[degree];
+                double maxValue = zero;
+                for (int32_t i = 0; i < degree; ++i) {
+                    const double value = std::fabs(c[i] * invLeading);
+                    if (value > maxValue)
+                        maxValue = value;
+                }
+                const double bound = one + maxValue;
+
+                return FindRecursive(degree, c, -bound, bound, maxIterations, roots);
+            }
+            else if (degree == 0)
+                // The polynomial is a nonzero constant.
+                return 0;
+            else {
+                // The polynomial is identically zero.
+                roots[0] = zero;
+                return 1;
+            }
+        }
+        else
+            // Invalid degree or c.
+            return 0;
+    }
+
+    // If you know that p(tmin) * p(tmax) <= 0, then there must be at
+    // least one root in [tmin, tmax].  Compute it using bisection.
+    static bool Find(int32_t degree, const double* c, double tmin, double tmax, uint32_t maxIterations, double& root)
+    {
+        const double zero = 0.0;
+        double pmin = Evaluate(degree, c, tmin);
+        if (pmin == zero) {
+            root = tmin;
+            return true;
+        }
+        double pmax = Evaluate(degree, c, tmax);
+        if (pmax == zero) {
+            root = tmax;
+            return true;
+        }
+
+        if (pmin * pmax > zero)
+            // It is not known whether the interval bounds a root.
+            return false;
+
+        if (tmin >= tmax)
+            // Invalid ordering of interval endpoitns. 
+            return false;
+
+        for (uint32_t i = 1; i <= maxIterations; ++i) {
+            root = 0.5 * (tmin + tmax);
+
+            // This test is designed for 'float' or 'double' when tmin
+            // and tmax are consecutive floating-point numbers.
+            if (root == tmin || root == tmax)
+                break;
+
+            const double p = Evaluate(degree, c, root);
+            const double product = p * pmin;
+            if (product < zero) {
+                tmax = root;
+                pmax = p;
+            }
+            else if (product > zero) {
+                tmin = root;
+                pmin = p;
+            }
+            else
+                break;
+        }
+
+        return true;
+    }
+
+    // Support for the Find functions.
+    static int32_t FindRecursive(int32_t degree, double const* c, double tmin, double tmax, uint32_t maxIterations, double* roots)
+    {
+        // The base of the recursion.
+        const double zero = 0.0;
+        double root = zero;
+        if (degree == 1) {
+            int32_t numRoots;
+            if (c[1] != zero) {
+                root = -c[0] / c[1];
+                numRoots = 1;
+            }
+            else if (c[0] == zero) {
+                root = zero;
+                numRoots = 1;
+            }
+            else
+                numRoots = 0;
+
+            if (numRoots > 0 && tmin <= root && root <= tmax) {
+                roots[0] = root;
+                return 1;
+            }
+            return 0;
+        }
+
+        // Find the roots of the derivative polynomial scaled by 1/degree.
+        // The scaling avoids the factorial growth in the coefficients;
+        // for example, without the scaling, the high-order term x^d
+        // becomes (d!)*x through multiple differentiations.  With the
+        // scaling we instead get x.  This leads to better numerical
+        // behavior of the root finder.
+        const int32_t derivDegree = degree - 1;
+        std::vector<double> derivCoeff(static_cast<size_t>(derivDegree) + 1);
+        std::vector<double> derivRoots(derivDegree);
+        for (int32_t i = 0, ip1 = 1; i <= derivDegree; ++i, ++ip1) {
+            derivCoeff[i] = c[ip1] * (double)(ip1) / (double)degree;
+        }
+        const int32_t numDerivRoots = FindRecursive(degree - 1, &derivCoeff[0], tmin, tmax, maxIterations, &derivRoots[0]);
+
+        int32_t numRoots = 0;
+        if (numDerivRoots > 0) {
+            // Find root on [tmin,derivRoots[0]].
+            if (Find(degree, c, tmin, derivRoots[0], maxIterations, root))
+                roots[numRoots++] = root;
+
+            // Find root on [derivRoots[i],derivRoots[i+1]].
+            for (int32_t i = 0, ip1 = 1; i <= numDerivRoots - 2; ++i, ++ip1) {
+                if (Find(degree, c, derivRoots[i], derivRoots[ip1], maxIterations, root))
+                    roots[numRoots++] = root;
+            }
+
+            // Find root on [derivRoots[numDerivRoots-1],tmax].
+            if (Find(degree, c, derivRoots[static_cast<size_t>(numDerivRoots) - 1], tmax, maxIterations, root))
+                roots[numRoots++] = root;
+        }
+        else {
+            // The polynomial is monotone on [tmin,tmax], so has at most one root.
+            if (Find(degree, c, tmin, tmax, maxIterations, root))
+                roots[numRoots++] = root;
+        }
+        return numRoots;
+    }
+
+    static double Evaluate(int32_t degree, const double* c, double t)
+    {
+        int32_t i = degree;
+        double result = c[i];
+        while (--i >= 0) {
+            result = t * result + c[i];
+        }
+        return result;
+    }
+};
+
+// Adaptation of code found in: 
+// https://github.com/davideberly/GeometricTools/blob/master/GTE/Mathematics/Vector.h
+
+// Construct a single vector orthogonal to the nonzero input vector.  If
+// the maximum absolute component occurs at index i, then the orthogonal
+// vector U has u[i] = v[i+1], u[i+1] = -v[i], and all other components
+// zero.  The index addition i+1 is computed modulo N.
+inline Vec3d get_orthogonal(const Vec3d& v, bool unitLength)
+{
+    double cmax = std::fabs(v[0]);
+    int32_t imax = 0;
+    for (int32_t i = 1; i < 3; ++i) {
+        double c = std::fabs(v[i]);
+        if (c > cmax) {
+            cmax = c;
+            imax = i;
+        }
+    }
+
+    Vec3d result = Vec3d::Zero();
+    int32_t inext = imax + 1;
+    if (inext == 3)
+        inext = 0;
+
+    result[imax] = v[inext];
+    result[inext] = -v[imax];
+    if (unitLength) {
+        const double sqrDistance = result[imax] * result[imax] + result[inext] * result[inext];
+        const double invLength = 1.0 / std::sqrt(sqrDistance);
+        result[imax] *= invLength;
+        result[inext] *= invLength;
+    }
+    return result;
+}
+
+} // namespace Slic3r
+} // namespace Measure
+
+#endif // Slic3r_MeasureUtils_hpp_

src/libslic3r/SurfaceMesh.hpp

@@ -0,0 +1,167 @@
+///|/ Copyright (c) Prusa Research 2022 Lukáš Matěna @lukasmatena
+///|/
+///|/ PrusaSlicer is released under the terms of the AGPLv3 or higher
+///|/
+#ifndef slic3r_SurfaceMesh_hpp_
+#define slic3r_SurfaceMesh_hpp_
+
+#include <admesh/stl.h>
+#include <libslic3r/TriangleMesh.hpp>
+
+#include "boost/container/small_vector.hpp"
+
+namespace Slic3r {
+
+class TriangleMesh;
+
+
+
+enum Face_index : int;
+
+class Halfedge_index {
+    friend class SurfaceMesh;
+
+public:
+    Halfedge_index() : m_face(Face_index(-1)), m_side(0) {}
+    Face_index face() const { return m_face; }
+    unsigned char side() const { return m_side; }
+    bool is_invalid() const { return int(m_face) < 0; }
+    bool operator!=(const Halfedge_index& rhs) const { return ! ((*this) == rhs); }
+    bool operator==(const Halfedge_index& rhs) const { return m_face == rhs.m_face && m_side == rhs.m_side; }
+
+private:
+    Halfedge_index(int face_idx, unsigned char side_idx) : m_face(Face_index(face_idx)), m_side(side_idx) {}
+
+    Face_index m_face;
+    unsigned char m_side;
+};
+
+
+
+class Vertex_index {
+    friend class SurfaceMesh;
+
+public:
+    Vertex_index() : m_face(Face_index(-1)), m_vertex_idx(0) {}
+    bool is_invalid() const { return int(m_face) < 0; }
+    bool operator==(const Vertex_index& rhs) const = delete; // Use SurfaceMesh::is_same_vertex.
+
+private:
+    Vertex_index(int face_idx, unsigned char vertex_idx) : m_face(Face_index(face_idx)), m_vertex_idx(vertex_idx) {}
+
+    Face_index m_face;
+    unsigned char m_vertex_idx;
+};
+
+
+
+class SurfaceMesh {
+public:
+    explicit SurfaceMesh(const indexed_triangle_set& its)
+    : m_its(its),
+      m_face_neighbors(its_face_neighbors_par(its))
+    {}
+    SurfaceMesh(const SurfaceMesh&)            = delete;
+    SurfaceMesh& operator=(const SurfaceMesh&) = delete;
+
+    Vertex_index source(Halfedge_index h) const { assert(! h.is_invalid()); return Vertex_index(h.m_face, h.m_side); }
+    Vertex_index target(Halfedge_index h) const { assert(! h.is_invalid()); return Vertex_index(h.m_face, h.m_side == 2 ? 0 : h.m_side + 1); }
+    Face_index face(Halfedge_index h) const { assert(! h.is_invalid()); return h.m_face; }
+
+    Halfedge_index next(Halfedge_index h)     const { assert(! h.is_invalid()); h.m_side = (h.m_side + 1) % 3; return h; }
+    Halfedge_index prev(Halfedge_index h)     const { assert(! h.is_invalid()); h.m_side = (h.m_side == 0 ? 2 : h.m_side - 1); return h; }
+    Halfedge_index halfedge(Vertex_index v)   const { return Halfedge_index(v.m_face, (v.m_vertex_idx == 0 ? 2 : v.m_vertex_idx - 1)); }  
+    Halfedge_index halfedge(Face_index f)     const { return Halfedge_index(f, 0); }  
+    Halfedge_index opposite(Halfedge_index h) const {
+        if (h.is_invalid())
+            return h;
+
+        int face_idx = m_face_neighbors[h.m_face][h.m_side];
+        Halfedge_index h_candidate = halfedge(Face_index(face_idx));
+
+        if (h_candidate.is_invalid())
+            return Halfedge_index(); // invalid
+
+        for (int i=0; i<3; ++i) {
+            if (is_same_vertex(source(h_candidate), target(h))) {
+                // Meshes in PrusaSlicer should be fixed enough for the following not to happen.
+                assert(is_same_vertex(target(h_candidate), source(h)));
+                return h_candidate;
+            }
+            h_candidate = next(h_candidate);
+        }
+        return Halfedge_index(); // invalid
+    }
+
+    Halfedge_index next_around_target(Halfedge_index h) const { return opposite(next(h)); }
+    Halfedge_index prev_around_target(Halfedge_index h) const { Halfedge_index op = opposite(h); return (op.is_invalid() ? Halfedge_index() : prev(op)); }
+    Halfedge_index next_around_source(Halfedge_index h) const { Halfedge_index op = opposite(h); return (op.is_invalid() ? Halfedge_index() : next(op)); }
+    Halfedge_index prev_around_source(Halfedge_index h) const { return opposite(prev(h)); }
+    Halfedge_index halfedge(Vertex_index source, Vertex_index target) const
+    {
+        Halfedge_index hi(source.m_face, source.m_vertex_idx);
+        assert(! hi.is_invalid());
+
+        const Vertex_index orig_target = this->target(hi);
+        Vertex_index current_target = orig_target;
+
+        while (! is_same_vertex(current_target, target)) {
+            hi = next_around_source(hi);
+            if (hi.is_invalid())
+                break;
+            current_target = this->target(hi);
+            if (is_same_vertex(current_target, orig_target))
+                return Halfedge_index(); // invalid
+        }
+
+        return hi;
+    }
+
+    const stl_vertex& point(Vertex_index v) const { return m_its.vertices[m_its.indices[v.m_face][v.m_vertex_idx]]; }
+
+    size_t degree(Vertex_index v) const
+    {
+        // In case the mesh is broken badly, the loop might end up to be infinite,
+        // never getting back to the first halfedge. Remember list of all half-edges
+        // and trip if any is encountered for the second time.
+        Halfedge_index h_first = halfedge(v);
+        boost::container::small_vector<Halfedge_index, 10> he_visited;
+        Halfedge_index h = next_around_target(h_first);
+        size_t degree = 2;
+        while (! h.is_invalid() && h != h_first) {
+            he_visited.emplace_back(h);
+            h = next_around_target(h);
+            if (std::find(he_visited.begin(), he_visited.end(), h) == he_visited.end())
+                return 0;
+            ++degree;
+        }
+        return h.is_invalid() ? 0 : degree - 1;
+    }
+
+    size_t degree(Face_index f) const {
+        size_t total = 0;
+        for (unsigned char i=0; i<3; ++i) {
+            size_t d = degree(Vertex_index(f, i));
+            if (d == 0)
+                return 0;
+            total += d;
+        }
+        assert(total - 6 >= 0);
+        return total - 6; // we counted 3 halfedges from f, and one more for each neighbor
+    }
+
+    bool is_border(Halfedge_index h) const { return m_face_neighbors[h.m_face][h.m_side] == -1; }
+
+    bool is_same_vertex(const Vertex_index& a, const Vertex_index& b) const { return m_its.indices[a.m_face][a.m_vertex_idx] == m_its.indices[b.m_face][b.m_vertex_idx]; }
+    Vec3i get_face_neighbors(Face_index face_id) const { assert(int(face_id) < int(m_face_neighbors.size())); return m_face_neighbors[face_id]; }
+
+
+
+private:
+    const std::vector<Vec3i> m_face_neighbors;
+    const indexed_triangle_set& m_its;
+};
+
+} //namespace Slic3r
+
+#endif // slic3r_SurfaceMesh_hpp_