Slic3r::Measure Namespace Reference

Classes
struct	AngleAndEdges
struct	DistAndPoints
struct	MeasurementResult
class	Measuring
class	MeasuringImpl
class	Polynomial1
class	RootsPolynomial
class	SurfaceFeature

Enumerations
enum class	SurfaceFeatureType : int {   Undef = 0 , Point = 1 << 0 , Edge = 1 << 1 , Circle = 1 << 2 ,   Plane = 1 << 3 }

Functions
static std::tuple< Vec3d, double, double >	get_center_and_radius (const std::vector< Vec3d > &points, const Transform3d &trafo, const Transform3d &trafo_inv)
static std::array< Vec3d, 3 >	orthonormal_basis (const Vec3d &v)
static AngleAndEdges	angle_edge_edge (const std::pair< Vec3d, Vec3d > &e1, const std::pair< Vec3d, Vec3d > &e2)
static AngleAndEdges	angle_edge_plane (const std::pair< Vec3d, Vec3d > &e, const std::tuple< int, Vec3d, Vec3d > &p)
static AngleAndEdges	angle_plane_plane (const std::tuple< int, Vec3d, Vec3d > &p1, const std::tuple< int, Vec3d, Vec3d > &p2)
MeasurementResult	get_measurement (const SurfaceFeature &a, const SurfaceFeature &b, const Measuring *measuring)
Vec3d	edge_direction (const Vec3d &from, const Vec3d &to)
Vec3d	edge_direction (const std::pair< Vec3d, Vec3d > &e)
Vec3d	edge_direction (const SurfaceFeature &edge)
Vec3d	plane_normal (const SurfaceFeature &plane)
bool	are_parallel (const Vec3d &v1, const Vec3d &v2)
bool	are_perpendicular (const Vec3d &v1, const Vec3d &v2)
bool	are_parallel (const std::pair< Vec3d, Vec3d > &e1, const std::pair< Vec3d, Vec3d > &e2)
bool	are_parallel (const SurfaceFeature &f1, const SurfaceFeature &f2)
bool	are_perpendicular (const SurfaceFeature &f1, const SurfaceFeature &f2)
Polynomial1	operator* (const Polynomial1 &p0, const Polynomial1 &p1)
Polynomial1	operator+ (const Polynomial1 &p0, const Polynomial1 &p1)
Polynomial1	operator- (const Polynomial1 &p0, const Polynomial1 &p1)
Polynomial1	operator* (double scalar, const Polynomial1 &p)
Vec3d	get_orthogonal (const Vec3d &v, bool unitLength)

Variables
constexpr double	feature_hover_limit = 0.5

Enumeration Type Documentation

SurfaceFeatureType

enum class Slic3r::Measure::SurfaceFeatureType : int

strong

Enumerator
Undef
Point
Edge
Circle
Plane

                              : int {
    Undef  = 0,
    Point  = 1 << 0,
    Edge   = 1 << 1,
    Circle = 1 << 2,
    Plane  = 1 << 3
};

Function Documentation

angle_edge_edge()

static AngleAndEdges Slic3r::Measure::angle_edge_edge ( const std::pair< Vec3d, Vec3d > &  e1, const std::pair< Vec3d, Vec3d > &  e2  )

static

{
    if (are_parallel(e1, e2))
        return AngleAndEdges::Dummy;
 
    Vec3d e1_unit = edge_direction(e1.first, e1.second);
    Vec3d e2_unit = edge_direction(e2.first, e2.second);
 
    // project edges on the plane defined by them
    Vec3d normal = e1_unit.cross(e2_unit).normalized();
    const Eigen::Hyperplane<double, 3> plane(normal, e1.first);
    Vec3d e11_proj = plane.projection(e1.first);
    Vec3d e12_proj = plane.projection(e1.second);
    Vec3d e21_proj = plane.projection(e2.first);
    Vec3d e22_proj = plane.projection(e2.second);
 
    const bool coplanar = (e2.first - e21_proj).norm() < EPSILON && (e2.second - e22_proj).norm() < EPSILON;
 
    // rotate the plane to become the XY plane
    auto qp = Eigen::Quaternion<double>::FromTwoVectors(normal, Vec3d::UnitZ());
    auto qp_inverse = qp.inverse();
    const Vec3d e11_rot = qp * e11_proj;
    const Vec3d e12_rot = qp * e12_proj;
    const Vec3d e21_rot = qp * e21_proj;
    const Vec3d e22_rot = qp * e22_proj;
 
    // discard Z
    const Vec2d e11_rot_2d = Vec2d(e11_rot.x(), e11_rot.y());
    const Vec2d e12_rot_2d = Vec2d(e12_rot.x(), e12_rot.y());
    const Vec2d e21_rot_2d = Vec2d(e21_rot.x(), e21_rot.y());
    const Vec2d e22_rot_2d = Vec2d(e22_rot.x(), e22_rot.y());
 
    // find intersection (arc center) of edges in XY plane
    const Eigen::Hyperplane<double, 2> e1_rot_2d_line = Eigen::Hyperplane<double, 2>::Through(e11_rot_2d, e12_rot_2d);
    const Eigen::Hyperplane<double, 2> e2_rot_2d_line = Eigen::Hyperplane<double, 2>::Through(e21_rot_2d, e22_rot_2d);
    const Vec2d center_rot_2d = e1_rot_2d_line.intersection(e2_rot_2d_line);
 
    // arc center in original coordinate
    const Vec3d center = qp_inverse * Vec3d(center_rot_2d.x(), center_rot_2d.y(), e11_rot.z());
 
    // ensure the edges are pointing away from the center
    std::pair<Vec3d, Vec3d> out_e1 = e1;
    std::pair<Vec3d, Vec3d> out_e2 = e2;
    if ((center_rot_2d - e11_rot_2d).squaredNorm() > (center_rot_2d - e12_rot_2d).squaredNorm()) {
        std::swap(e11_proj, e12_proj);
        std::swap(out_e1.first, out_e1.second);
        e1_unit = -e1_unit;
    }
    if ((center_rot_2d - e21_rot_2d).squaredNorm() > (center_rot_2d - e22_rot_2d).squaredNorm()) {
        std::swap(e21_proj, e22_proj);
        std::swap(out_e2.first, out_e2.second);
        e2_unit = -e2_unit;
    }
 
    // arc angle
    const double angle = std::acos(std::clamp(e1_unit.dot(e2_unit), -1.0, 1.0));
    // arc radius
    const Vec3d e1_proj_mid = 0.5 * (e11_proj + e12_proj);
    const Vec3d e2_proj_mid = 0.5 * (e21_proj + e22_proj);
    const double radius = std::min((center - e1_proj_mid).norm(), (center - e2_proj_mid).norm());
 
    return { angle, center, out_e1, out_e2, radius, coplanar };
}

References Slic3r::angle(), are_parallel(), Slic3r::Measure::AngleAndEdges::Dummy, edge_direction(), EPSILON, Eigen::Quaternion< _Scalar, _Options >::FromTwoVectors(), Eigen::Hyperplane< _Scalar, _AmbientDim, _Options >::intersection(), Eigen::Hyperplane< _Scalar, _AmbientDim, _Options >::projection(), and Eigen::Hyperplane< _Scalar, _AmbientDim, _Options >::Through().

Referenced by angle_edge_plane(), angle_plane_plane(), and get_measurement().

Here is the call graph for this function:

Here is the caller graph for this function:

angle_edge_plane()

static AngleAndEdges Slic3r::Measure::angle_edge_plane ( const std::pair< Vec3d, Vec3d > &  e, const std::tuple< int, Vec3d, Vec3d > &  p  )

static

{
    const auto& [idx, normal, origin] = p;
    Vec3d e1e2_unit = edge_direction(e);
    if (are_perpendicular(e1e2_unit, normal))
        return AngleAndEdges::Dummy;
 
    // ensure the edge is pointing away from the intersection
    // 1st calculate instersection between edge and plane
    const Eigen::Hyperplane<double, 3> plane(normal, origin);
    const Eigen::ParametrizedLine<double, 3> line = Eigen::ParametrizedLine<double, 3>::Through(e.first, e.second);
    const Vec3d inters = line.intersectionPoint(plane);
 
    // then verify edge direction and revert it, if needed
    Vec3d e1 = e.first;
    Vec3d e2 = e.second;
    if ((e1 - inters).squaredNorm() > (e2 - inters).squaredNorm()) {
        std::swap(e1, e2);
        e1e2_unit = -e1e2_unit;
    }
 
    if (are_parallel(e1e2_unit, normal)) {
        const std::array<Vec3d, 3> basis = orthonormal_basis(e1e2_unit);
        const double radius = (0.5 * (e1 + e2) - inters).norm();
        const Vec3d edge_on_plane_dir = (basis[1].dot(origin - inters) >= 0.0) ? basis[1] : -basis[1];
        std::pair<Vec3d, Vec3d> edge_on_plane = std::make_pair(inters, inters + radius * edge_on_plane_dir);
        if (!inters.isApprox(e1)) {
            edge_on_plane.first  += radius * edge_on_plane_dir;
            edge_on_plane.second += radius * edge_on_plane_dir;
        }
        return AngleAndEdges(0.5 * double(PI), inters, std::make_pair(e1, e2), edge_on_plane, radius, inters.isApprox(e1));
    }
 
    const Vec3d e1e2 = e2 - e1;
    const double e1e2_len = e1e2.norm();
 
    // calculate 2nd edge (on the plane)
    const Vec3d temp = normal.cross(e1e2);
    const Vec3d edge_on_plane_unit = normal.cross(temp).normalized();
    std::pair<Vec3d, Vec3d> edge_on_plane = { origin, origin + e1e2_len * edge_on_plane_unit };
 
    // ensure the 2nd edge is pointing in the correct direction
    const Vec3d test_edge = (edge_on_plane.second - edge_on_plane.first).cross(e1e2);
    if (test_edge.dot(temp) < 0.0)
        edge_on_plane = { origin, origin - e1e2_len * edge_on_plane_unit };
 
    AngleAndEdges ret = angle_edge_edge({ e1, e2 }, edge_on_plane);
    ret.radius = (inters - 0.5 * (e1 + e2)).norm();
    return ret;
}

References angle_edge_edge(), are_parallel(), are_perpendicular(), Slic3r::cross(), Slic3r::Measure::AngleAndEdges::Dummy, edge_direction(), Eigen::ParametrizedLine< _Scalar, _AmbientDim, _Options >::intersectionPoint(), orthonormal_basis(), PI, Slic3r::Measure::AngleAndEdges::radius, and Eigen::ParametrizedLine< _Scalar, _AmbientDim, _Options >::Through().

Referenced by get_measurement().

Here is the call graph for this function:

Here is the caller graph for this function:

angle_plane_plane()

static AngleAndEdges Slic3r::Measure::angle_plane_plane ( const std::tuple< int, Vec3d, Vec3d > &  p1, const std::tuple< int, Vec3d, Vec3d > &  p2  )

static

{
    const auto& [idx1, normal1, origin1] = p1;
    const auto& [idx2, normal2, origin2] = p2;
 
    // are planes parallel ?
    if (are_parallel(normal1, normal2))
        return AngleAndEdges::Dummy;
 
    auto intersection_plane_plane = [](const Vec3d& n1, const Vec3d& o1, const Vec3d& n2, const Vec3d& o2) {
        Eigen::MatrixXd m(2, 3);
        m << n1.x(), n1.y(), n1.z(), n2.x(), n2.y(), n2.z();
        Eigen::VectorXd b(2);
        b << o1.dot(n1), o2.dot(n2);
        Eigen::VectorXd x = m.colPivHouseholderQr().solve(b);
        return std::make_pair(n1.cross(n2).normalized(), Vec3d(x(0), x(1), x(2)));
    };
 
    // Calculate intersection line between planes
    const auto [intersection_line_direction, intersection_line_origin] = intersection_plane_plane(normal1, origin1, normal2, origin2);
 
    // Project planes' origin on intersection line
    const Eigen::ParametrizedLine<double, 3> intersection_line = Eigen::ParametrizedLine<double, 3>(intersection_line_origin, intersection_line_direction);
    const Vec3d origin1_proj = intersection_line.projection(origin1);
    const Vec3d origin2_proj = intersection_line.projection(origin2);
 
    // Calculate edges on planes
    const Vec3d edge_on_plane1_unit = (origin1 - origin1_proj).normalized();
    const Vec3d edge_on_plane2_unit = (origin2 - origin2_proj).normalized();
    const double radius = std::max(10.0, std::max((origin1 - origin1_proj).norm(), (origin2 - origin2_proj).norm()));
    const std::pair<Vec3d, Vec3d> edge_on_plane1 = { origin1_proj + radius * edge_on_plane1_unit, origin1_proj + 2.0 * radius * edge_on_plane1_unit };
    const std::pair<Vec3d, Vec3d> edge_on_plane2 = { origin2_proj + radius * edge_on_plane2_unit, origin2_proj + 2.0 * radius * edge_on_plane2_unit };
 
    AngleAndEdges ret = angle_edge_edge(edge_on_plane1, edge_on_plane2);
    ret.radius = radius;
    return ret;
}

References angle_edge_edge(), are_parallel(), Slic3r::Measure::AngleAndEdges::Dummy, Eigen::ParametrizedLine< _Scalar, _AmbientDim, _Options >::projection(), and Slic3r::Measure::AngleAndEdges::radius.

Referenced by get_measurement().

Here is the call graph for this function:

Here is the caller graph for this function:

are_parallel() [1/3]

bool Slic3r::Measure::are_parallel ( const std::pair< Vec3d, Vec3d > &  e1, const std::pair< Vec3d, Vec3d > &  e2  )

inline

                                                                                           {
    return are_parallel(e1.second - e1.first, e2.second - e2.first);
}

References are_parallel().

Here is the call graph for this function:

are_parallel() [2/3]

bool Slic3r::Measure::are_parallel ( const SurfaceFeature &  f1, const SurfaceFeature &  f2  )

inline

                                                                             {
    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;
}

References are_parallel(), are_perpendicular(), Edge, edge_direction(), Slic3r::Measure::SurfaceFeature::get_type(), Plane, and plane_normal().

Here is the call graph for this function:

are_parallel() [3/3]

bool Slic3r::Measure::are_parallel ( const Vec3d &  v1, const Vec3d &  v2  )

inline

{ return std::abs(std::abs(v1.dot(v2)) - 1.0) < EPSILON; }

References EPSILON.

Referenced by angle_edge_edge(), angle_edge_plane(), angle_plane_plane(), are_parallel(), are_parallel(), are_perpendicular(), and get_measurement().

Here is the caller graph for this function:

are_perpendicular() [1/2]

bool Slic3r::Measure::are_perpendicular ( const SurfaceFeature &  f1, const SurfaceFeature &  f2  )

inline

                                                                                  {
    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;
}

References are_parallel(), are_perpendicular(), Edge, edge_direction(), Slic3r::Measure::SurfaceFeature::get_type(), Plane, and plane_normal().

Here is the call graph for this function:

are_perpendicular() [2/2]

bool Slic3r::Measure::are_perpendicular ( const Vec3d &  v1, const Vec3d &  v2  )

inline

{ return std::abs(v1.dot(v2)) < EPSILON; }

References EPSILON.

Referenced by angle_edge_plane(), are_parallel(), are_perpendicular(), and get_measurement().

Here is the caller graph for this function:

edge_direction() [1/3]

Vec3d Slic3r::Measure::edge_direction ( const std::pair< Vec3d, Vec3d > &  e )

inline

{ return edge_direction(e.first, e.second); }

References edge_direction().

Here is the call graph for this function:

edge_direction() [2/3]

Vec3d Slic3r::Measure::edge_direction ( const SurfaceFeature &  edge )

inline

                                                        {
    assert(edge.get_type() == SurfaceFeatureType::Edge);
    return edge_direction(edge.get_edge());
}

References Edge, edge_direction(), Slic3r::Measure::SurfaceFeature::get_edge(), and Slic3r::Measure::SurfaceFeature::get_type().

Here is the call graph for this function:

edge_direction() [3/3]

Vec3d Slic3r::Measure::edge_direction ( const Vec3d &  from, const Vec3d &  to  )

inline

{ return (to - from).normalized(); }

Referenced by angle_edge_edge(), angle_edge_plane(), are_parallel(), are_perpendicular(), edge_direction(), edge_direction(), and Slic3r::GUI::GLGizmoMeasure::render_dimensioning().

Here is the caller graph for this function:

get_center_and_radius()

static std::tuple< Vec3d, double, double > Slic3r::Measure::get_center_and_radius ( const std::vector< Vec3d > &  points, const Transform3d &  trafo, const Transform3d &  trafo_inv  )

static

{
    Vec2ds out;
    double z = 0.;
    for (const Vec3d& pt : points) {
        Vec3d pt_transformed = trafo * pt;
        z = pt_transformed.z();
        out.emplace_back(pt_transformed.x(), pt_transformed.y());
    }
 
    const int iter = points.size() < 10  ? 2 :
                     points.size() < 100 ? 4 :
                     6;
 
    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);
}

References Slic3r::Geometry::circle_ransac(), error, and Eigen::Transform< _Scalar, _Dim, _Mode, _Options >::inverse().

Referenced by Slic3r::Measure::MeasuringImpl::extract_features().

Here is the call graph for this function:

Here is the caller graph for this function:

get_measurement()

MeasurementResult Slic3r::Measure::get_measurement	(	const SurfaceFeature &	a,
		const SurfaceFeature &	b,
		const Measuring *	measuring
	)

{
    assert(a.get_type() != SurfaceFeatureType::Undef && b.get_type() != SurfaceFeatureType::Undef);
 
    const bool swap = int(a.get_type()) > int(b.get_type());
    const SurfaceFeature& f1 = swap ? b : a;
    const SurfaceFeature& f2 = swap ? a : b;
 
    MeasurementResult result;
 
    if (f1.get_type() == SurfaceFeatureType::Point) {
        if (f2.get_type() == SurfaceFeatureType::Point) {
            Vec3d diff = (f2.get_point() - f1.get_point());
            result.distance_strict = std::make_optional(DistAndPoints{diff.norm(), f1.get_point(), f2.get_point()});
            result.distance_xyz = diff.cwiseAbs();
 
        } else if (f2.get_type() == SurfaceFeatureType::Edge) {
            const auto [s,e] = f2.get_edge();
            const Eigen::ParametrizedLine<double, 3> line(s, (e-s).normalized());
            const double dist_inf = line.distance(f1.get_point());
            const Vec3d proj = line.projection(f1.get_point());
            const double len_sq = (e-s).squaredNorm();
            const double dist_start_sq = (proj-s).squaredNorm();
            const double dist_end_sq = (proj-e).squaredNorm();
            if (dist_start_sq < len_sq && dist_end_sq < len_sq) {
                // projection falls on the line - the strict distance is the same as infinite
                result.distance_strict = std::make_optional(DistAndPoints{dist_inf, f1.get_point(), proj});
            } else { // the result is the closer of the endpoints
                const bool s_is_closer = dist_start_sq < dist_end_sq;
                result.distance_strict = std::make_optional(DistAndPoints{std::sqrt(std::min(dist_start_sq, dist_end_sq) + sqr(dist_inf)), f1.get_point(), s_is_closer ? s : e});
            }
            result.distance_infinite = std::make_optional(DistAndPoints{dist_inf, f1.get_point(), proj});
        } else if (f2.get_type() == SurfaceFeatureType::Circle) {
            // Find a plane containing normal, center and the point.
            const auto [c, radius, n] = f2.get_circle();
            const Eigen::Hyperplane<double, 3> circle_plane(n, c);
            const Vec3d proj = circle_plane.projection(f1.get_point());
            if (proj.isApprox(c)) {
                const Vec3d p_on_circle = c + radius * get_orthogonal(n, true);
                result.distance_strict = std::make_optional(DistAndPoints{ radius, c, p_on_circle });
            }
            else {
                const Eigen::Hyperplane<double, 3> circle_plane(n, c);
                const Vec3d proj = circle_plane.projection(f1.get_point());
                const double dist = std::sqrt(std::pow((proj - c).norm() - radius, 2.) +
                    (f1.get_point() - proj).squaredNorm());
 
                const Vec3d p_on_circle = c + radius * (proj - c).normalized();
                result.distance_strict = std::make_optional(DistAndPoints{ dist, f1.get_point(), p_on_circle }); // TODO
            }
        } else if (f2.get_type() == SurfaceFeatureType::Plane) {
            const auto [idx, normal, pt] = f2.get_plane();
            Eigen::Hyperplane<double, 3> plane(normal, pt);
            result.distance_infinite = std::make_optional(DistAndPoints{plane.absDistance(f1.get_point()), f1.get_point(), plane.projection(f1.get_point())}); // TODO
            // TODO: result.distance_strict =
        }
    }
    else if (f1.get_type() == SurfaceFeatureType::Edge) {
        if (f2.get_type() == SurfaceFeatureType::Edge) {
            std::vector<DistAndPoints> distances;
 
            auto add_point_edge_distance = [&distances](const Vec3d& v, const std::pair<Vec3d, Vec3d>& e) {
                const MeasurementResult res = get_measurement(SurfaceFeature(v), SurfaceFeature(SurfaceFeatureType::Edge, e.first, e.second));
                double distance = res.distance_strict->dist;
                Vec3d v2 = res.distance_strict->to;
 
                const Vec3d e1e2 = e.second - e.first;
                const Vec3d e1v2 = v2 - e.first;
                if (e1v2.dot(e1e2) >= 0.0 && e1v2.norm() < e1e2.norm())
                    distances.emplace_back(distance, v, v2);
            };
 
            std::pair<Vec3d, Vec3d> e1 = f1.get_edge();
            std::pair<Vec3d, Vec3d> e2 = f2.get_edge();
 
            distances.emplace_back((e2.first - e1.first).norm(), e1.first, e2.first);
            distances.emplace_back((e2.second - e1.first).norm(), e1.first, e2.second);
            distances.emplace_back((e2.first - e1.second).norm(), e1.second, e2.first);
            distances.emplace_back((e2.second - e1.second).norm(), e1.second, e2.second);
            add_point_edge_distance(e1.first, e2);
            add_point_edge_distance(e1.second, e2);
            add_point_edge_distance(e2.first, e1);
            add_point_edge_distance(e2.second, e1);
            auto it = std::min_element(distances.begin(), distances.end(),
                [](const DistAndPoints& item1, const DistAndPoints& item2) {
                    return item1.dist < item2.dist;
                });
            result.distance_infinite = std::make_optional(*it);
 
            result.angle = angle_edge_edge(f1.get_edge(), f2.get_edge());
        } else if (f2.get_type() == SurfaceFeatureType::Circle) {
            const std::pair<Vec3d, Vec3d> e = f1.get_edge();
            const auto& [center, radius, normal] = f2.get_circle();
            const Vec3d e1e2 = (e.second - e.first);
            const Vec3d e1e2_unit = e1e2.normalized();
 
            std::vector<DistAndPoints> distances;
            distances.emplace_back(*get_measurement(SurfaceFeature(e.first), f2).distance_strict);
            distances.emplace_back(*get_measurement(SurfaceFeature(e.second), f2).distance_strict);
 
            const Eigen::Hyperplane<double, 3> plane(e1e2_unit, center);
            const Eigen::ParametrizedLine<double, 3> line = Eigen::ParametrizedLine<double, 3>::Through(e.first, e.second);
            const Vec3d inter = line.intersectionPoint(plane);
            const Vec3d e1inter = inter - e.first;
            if (e1inter.dot(e1e2) >= 0.0 && e1inter.norm() < e1e2.norm())
                distances.emplace_back(*get_measurement(SurfaceFeature(inter), f2).distance_strict);
 
            auto it = std::min_element(distances.begin(), distances.end(),
                [](const DistAndPoints& item1, const DistAndPoints& item2) {
                    return item1.dist < item2.dist;
                });
            result.distance_infinite = std::make_optional(DistAndPoints{it->dist, it->from, it->to});
        } else if (f2.get_type() == SurfaceFeatureType::Plane) {
            assert(measuring != nullptr);
 
            const auto [from, to] = f1.get_edge();
            const auto [idx, normal, origin] = f2.get_plane();
 
            const Vec3d edge_unit = (to - from).normalized();
            if (are_perpendicular(edge_unit, normal)) {
                std::vector<DistAndPoints> distances;
                const Eigen::Hyperplane<double, 3> plane(normal, origin);
                distances.push_back(DistAndPoints{ plane.absDistance(from), from, plane.projection(from) });
                distances.push_back(DistAndPoints{ plane.absDistance(to), to, plane.projection(to) });
                auto it = std::min_element(distances.begin(), distances.end(),
                    [](const DistAndPoints& item1, const DistAndPoints& item2) {
                        return item1.dist < item2.dist;
                    });
                result.distance_infinite = std::make_optional(DistAndPoints{ it->dist, it->from, it->to });
            }
            else {
                const std::vector<SurfaceFeature>& plane_features = measuring->get_plane_features(idx);
                std::vector<DistAndPoints> distances;
                for (const SurfaceFeature& sf : plane_features) {
                    if (sf.get_type() == SurfaceFeatureType::Edge) {
                        const auto m = get_measurement(sf, f1);
                        if (!m.distance_infinite.has_value()) {
                            distances.clear();
                            break;
                        }
                        else
                            distances.push_back(*m.distance_infinite);
                    }
                }
                if (!distances.empty()) {
                    auto it = std::min_element(distances.begin(), distances.end(),
                        [](const DistAndPoints& item1, const DistAndPoints& item2) {
                            return item1.dist < item2.dist;
                        });
                    result.distance_infinite = std::make_optional(DistAndPoints{ it->dist, it->from, it->to });
                }
            }
            result.angle = angle_edge_plane(f1.get_edge(), f2.get_plane());
        }
    } else if (f1.get_type() == SurfaceFeatureType::Circle) {
        if (f2.get_type() == SurfaceFeatureType::Circle) {
            const auto [c0, r0, n0] = f1.get_circle();
            const auto [c1, r1, n1] = f2.get_circle();
 
            // The following code is an adaptation of the algorithm found in: 
            // https://github.com/davideberly/GeometricTools/blob/master/GTE/Mathematics/DistCircle3Circle3.h
            // and described in:
            // https://www.geometrictools.com/Documentation/DistanceToCircle3.pdf
 
            struct ClosestInfo
            {
                double sqrDistance{ 0.0 };
                Vec3d circle0Closest{ Vec3d::Zero() };
                Vec3d circle1Closest{ Vec3d::Zero() };
 
                inline bool operator < (const ClosestInfo& other) const { return sqrDistance < other.sqrDistance; }
            };
            std::array<ClosestInfo, 16> candidates{};
 
            const double zero = 0.0;
 
            const Vec3d D = c1 - c0;
 
            if (!are_parallel(n0, n1)) {
                // Get parameters for constructing the degree-8 polynomial phi.
                const double one = 1.0;
                const double two = 2.0;
                const double  r0sqr = sqr(r0);
                const double  r1sqr = sqr(r1);
 
                // Compute U1 and V1 for the plane of circle1.
                const std::array<Vec3d, 3> basis = orthonormal_basis(n1);
                const Vec3d U1 = basis[0];
                const Vec3d V1 = basis[1];
 
                // Construct the polynomial phi(cos(theta)).
                const Vec3d N0xD = n0.cross(D);
                const Vec3d N0xU1 = n0.cross(U1);
                const Vec3d N0xV1 = n0.cross(V1);
                const double a0 = r1 * D.dot(U1);
                const double a1 = r1 * D.dot(V1);
                const double a2 = N0xD.dot(N0xD);
                const double a3 = r1 * N0xD.dot(N0xU1);
                const double a4 = r1 * N0xD.dot(N0xV1);
                const double a5 = r1sqr * N0xU1.dot(N0xU1);
                const double a6 = r1sqr * N0xU1.dot(N0xV1);
                const double a7 = r1sqr * N0xV1.dot(N0xV1);
                Polynomial1 p0{ a2 + a7, two * a3, a5 - a7 };
                Polynomial1 p1{ two * a4, two * a6 };
                Polynomial1 p2{ zero, a1 };
                Polynomial1 p3{ -a0 };
                Polynomial1 p4{ -a6, a4, two * a6 };
                Polynomial1 p5{ -a3, a7 - a5 };
                Polynomial1 tmp0{ one, zero, -one };
                Polynomial1 tmp1 = p2 * p2 + tmp0 * p3 * p3;
                Polynomial1 tmp2 = two * p2 * p3;
                Polynomial1 tmp3 = p4 * p4 + tmp0 * p5 * p5;
                Polynomial1 tmp4 = two * p4 * p5;
                Polynomial1 p6 = p0 * tmp1 + tmp0 * p1 * tmp2 - r0sqr * tmp3;
                Polynomial1 p7 = p0 * tmp2 + p1 * tmp1 - r0sqr * tmp4;
 
                // Parameters for polynomial root finding. The roots[] array
                // stores the roots. We need only the unique ones, which is
                // the responsibility of the set uniqueRoots. The pairs[]
                // array stores the (cosine,sine) information mentioned in the
                // PDF. TODO: Choose the maximum number of iterations for root
                // finding based on specific polynomial data?
                const uint32_t maxIterations = 128;
                int32_t degree = 0;
                size_t numRoots = 0;
                std::array<double, 8> roots{};
                std::set<double> uniqueRoots{};
                size_t numPairs = 0;
                std::array<std::pair<double, double>, 16> pairs{};
                double temp = zero;
                double sn = zero;
 
                if (p7.GetDegree() > 0 || p7[0] != zero) {
                    // H(cs,sn) = p6(cs) + sn * p7(cs)
                    Polynomial1 phi = p6 * p6 - tmp0 * p7 * p7;
                    degree = static_cast<int32_t>(phi.GetDegree());
                    assert(degree > 0);
                    numRoots = RootsPolynomial::Find(degree, &phi[0], maxIterations, roots.data());
                    for (size_t i = 0; i < numRoots; ++i) {
                        uniqueRoots.insert(roots[i]);
                    }
 
                    for (auto const& cs : uniqueRoots) {
                        if (std::fabs(cs) <= one) {
                            temp = p7(cs);
                            if (temp != zero) {
                                sn = -p6(cs) / temp;
                                pairs[numPairs++] = std::make_pair(cs, sn);
                            }
                            else {
                                temp = std::max(one - sqr(cs), zero);
                                sn = std::sqrt(temp);
                                pairs[numPairs++] = std::make_pair(cs, sn);
                                if (sn != zero)
                                    pairs[numPairs++] = std::make_pair(cs, -sn);
                            }
                        }
                    }
                }
                else {
                    // H(cs,sn) = p6(cs)
                    degree = static_cast<int32_t>(p6.GetDegree());
                    assert(degree > 0);
                    numRoots = RootsPolynomial::Find(degree, &p6[0], maxIterations, roots.data());
                    for (size_t i = 0; i < numRoots; ++i) {
                        uniqueRoots.insert(roots[i]);
                    }
 
                    for (auto const& cs : uniqueRoots) {
                        if (std::fabs(cs) <= one) {
                            temp = std::max(one - sqr(cs), zero);
                            sn = std::sqrt(temp);
                            pairs[numPairs++] = std::make_pair(cs, sn);
                            if (sn != zero)
                                pairs[numPairs++] = std::make_pair(cs, -sn);
                        }
                    }
                }
 
                for (size_t i = 0; i < numPairs; ++i) {
                    ClosestInfo& info = candidates[i];
                    Vec3d delta = D + r1 * (pairs[i].first * U1 + pairs[i].second * V1);
                    info.circle1Closest = c0 + delta;
                    const double N0dDelta = n0.dot(delta);
                    const double lenN0xDelta = n0.cross(delta).norm();
                    if (lenN0xDelta > 0.0) {
                        const double diff = lenN0xDelta - r0;
                        info.sqrDistance = sqr(N0dDelta) + sqr(diff);
                        delta -= N0dDelta * n0;
                        delta.normalize();
                        info.circle0Closest = c0 + r0 * delta;
                    }
                    else {
                        const Vec3d r0U0 = r0 * get_orthogonal(n0, true);
                        const Vec3d diff = delta - r0U0;
                        info.sqrDistance = diff.dot(diff);
                        info.circle0Closest = c0 + r0U0;
                    }
                }
 
                std::sort(candidates.begin(), candidates.begin() + numPairs);
            }
            else {
                ClosestInfo& info = candidates[0];
            
                const double N0dD = n0.dot(D);
                const Vec3d normProj = N0dD * n0;
                const Vec3d compProj = D - normProj;
                Vec3d U = compProj;
                const double d = U.norm();
                U.normalize();
 
                // The configuration is determined by the relative location of the
                // intervals of projection of the circles on to the D-line.
                // Circle0 projects to [-r0,r0] and circle1 projects to
                // [d-r1,d+r1].
                const double dmr1 = d - r1;
                double distance;
                if (dmr1 >= r0) {
                    // d >= r0 + r1
                    // The circles are separated (d > r0 + r1) or tangent with one
                    // outside the other (d = r0 + r1).
                    distance = dmr1 - r0;
                    info.circle0Closest = c0 + r0 * U;
                    info.circle1Closest = c1 - r1 * U;
                }
                else {
                    // d < r0 + r1
                    // The cases implicitly use the knowledge that d >= 0.
                    const double dpr1 = d + r1;
                    if (dpr1 <= r0) {
                        // Circle1 is inside circle0.
                        distance = r0 - dpr1;
                        if (d > 0.0) {
                            info.circle0Closest = c0 + r0 * U;
                            info.circle1Closest = c1 + r1 * U;
                        }
                        else {
                            // The circles are concentric, so U = (0,0,0).
                            // Construct a vector perpendicular to N0 to use for
                            // closest points.
                            U = get_orthogonal(n0, true);
                            info.circle0Closest = c0 + r0 * U;
                            info.circle1Closest = c1 + r1 * U;
                        }
                    }
                    else if (dmr1 <= -r0) {
                        // Circle0 is inside circle1.
                        distance = -r0 - dmr1;
                        if (d > 0.0) {
                            info.circle0Closest = c0 - r0 * U;
                            info.circle1Closest = c1 - r1 * U;
                        }
                        else {
                            // The circles are concentric, so U = (0,0,0).
                            // Construct a vector perpendicular to N0 to use for
                            // closest points.
                            U = get_orthogonal(n0, true);
                            info.circle0Closest = c0 + r0 * U;
                            info.circle1Closest = c1 + r1 * U;
                        }
                    }
                    else {
                        distance = (c1 - c0).norm();
                        info.circle0Closest = c0;
                        info.circle1Closest = c1;
                    }
                }
 
                info.sqrDistance = distance * distance + N0dD * N0dD;
            }
 
            result.distance_infinite = std::make_optional(DistAndPoints{ std::sqrt(candidates[0].sqrDistance), candidates[0].circle0Closest, candidates[0].circle1Closest }); // TODO
        } else if (f2.get_type() == SurfaceFeatureType::Plane) {
            assert(measuring != nullptr);
 
            const auto [center, radius, normal1] = f1.get_circle();
            const auto [idx2, normal2, origin2] = f2.get_plane();
 
            const bool coplanar = are_parallel(normal1, normal2) && Eigen::Hyperplane<double, 3>(normal1, center).absDistance(origin2) < EPSILON;
            if (!coplanar) {
                const std::vector<SurfaceFeature>& plane_features = measuring->get_plane_features(idx2);
                std::vector<DistAndPoints> distances;
                for (const SurfaceFeature& sf : plane_features) {
                    if (sf.get_type() == SurfaceFeatureType::Edge) {
                        const auto m = get_measurement(sf, f1);
                        if (!m.distance_infinite.has_value()) {
                            distances.clear();
                            break;
                        }
                        else
                            distances.push_back(*m.distance_infinite);
                    }
                }
                if (!distances.empty()) {
                    auto it = std::min_element(distances.begin(), distances.end(),
                        [](const DistAndPoints& item1, const DistAndPoints& item2) {
                            return item1.dist < item2.dist;
                        });
                    result.distance_infinite = std::make_optional(DistAndPoints{ it->dist, it->from, it->to });
                }
            }
        }
    } else if (f1.get_type() == SurfaceFeatureType::Plane) {
        const auto [idx1, normal1, pt1] = f1.get_plane();
        const auto [idx2, normal2, pt2] = f2.get_plane();
 
        if (are_parallel(normal1, normal2)) {
            // The planes are parallel, calculate distance.
            const Eigen::Hyperplane<double, 3> plane(normal1, pt1);
            result.distance_infinite = std::make_optional(DistAndPoints{ plane.absDistance(pt2), pt2, plane.projection(pt2) }); // TODO
        }
        else
            result.angle = angle_plane_plane(f1.get_plane(), f2.get_plane());
    }
    
    return result;
}

References Eigen::Hyperplane< _Scalar, _AmbientDim, _Options >::absDistance(), Slic3r::Measure::MeasurementResult::angle, angle_edge_edge(), angle_edge_plane(), angle_plane_plane(), are_parallel(), are_perpendicular(), Circle, Slic3r::diff(), Slic3r::Measure::DistAndPoints::dist, Eigen::ParametrizedLine< _Scalar, _AmbientDim, _Options >::distance(), Slic3r::Measure::MeasurementResult::distance_infinite, Slic3r::Measure::MeasurementResult::distance_strict, Slic3r::Measure::MeasurementResult::distance_xyz, Edge, EPSILON, Slic3r::Measure::RootsPolynomial::Find(), Slic3r::Measure::SurfaceFeature::get_circle(), Slic3r::Measure::SurfaceFeature::get_edge(), get_measurement(), get_orthogonal(), Slic3r::Measure::SurfaceFeature::get_plane(), Slic3r::Measure::Measuring::get_plane_features(), Slic3r::Measure::SurfaceFeature::get_point(), Slic3r::Measure::SurfaceFeature::get_type(), Slic3r::Measure::Polynomial1::GetDegree(), Eigen::ParametrizedLine< _Scalar, _AmbientDim, _Options >::intersectionPoint(), orthonormal_basis(), Plane, Point, Eigen::Hyperplane< _Scalar, _AmbientDim, _Options >::projection(), Eigen::ParametrizedLine< _Scalar, _AmbientDim, _Options >::projection(), Slic3r::sqr(), Eigen::ParametrizedLine< _Scalar, _AmbientDim, _Options >::Through(), and Undef.

Referenced by Slic3r::Measure::MeasuringImpl::get_feature(), get_measurement(), and Slic3r::GUI::GLGizmoMeasure::update_measurement_result().

Here is the call graph for this function:

Here is the caller graph for this function:

get_orthogonal()

Vec3d Slic3r::Measure::get_orthogonal ( const Vec3d &  v, bool  unitLength  )

inline

{
    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;
}

Referenced by get_measurement(), and Slic3r::GUI::GLGizmoMeasure::on_render_input_window().

Here is the caller graph for this function:

operator*() [1/2]

Polynomial1 Slic3r::Measure::operator* ( const Polynomial1 &  p0, const Polynomial1 &  p1  )

inline

{
    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;
}

References Slic3r::Measure::Polynomial1::GetDegree(), and Slic3r::Measure::Polynomial1::SetCoefficients().

Here is the call graph for this function:

operator*() [2/2]

Polynomial1 Slic3r::Measure::operator* ( double  scalar, const Polynomial1 &  p  )

inline

{
    const uint32_t degree = p.GetDegree();
    Polynomial1 result(degree);
    for (uint32_t i = 0; i <= degree; ++i) {
        result[i] = scalar * p[i];
    }
    return result;
}

References Slic3r::Measure::Polynomial1::GetDegree().

Here is the call graph for this function:

operator+()

Polynomial1 Slic3r::Measure::operator+ ( const Polynomial1 &  p0, const Polynomial1 &  p1  )

inline

{
    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;
    }
}

References Slic3r::Measure::Polynomial1::EliminateLeadingZeros(), and Slic3r::Measure::Polynomial1::GetDegree().

Here is the call graph for this function:

operator-()

Polynomial1 Slic3r::Measure::operator- ( const Polynomial1 &  p0, const Polynomial1 &  p1  )

inline

{
    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;
    }
}

References Slic3r::Measure::Polynomial1::EliminateLeadingZeros(), and Slic3r::Measure::Polynomial1::GetDegree().

Here is the call graph for this function:

orthonormal_basis()

static std::array< Vec3d, 3 > Slic3r::Measure::orthonormal_basis ( const Vec3d &  v )

static

{
    std::array<Vec3d, 3> ret;
    ret[2] = v.normalized();
    int index;
    ret[2].cwiseAbs().maxCoeff(&index);
    switch (index)
    {
    case 0: { ret[0] = Vec3d(ret[2].y(), -ret[2].x(), 0.0).normalized(); break; }
    case 1: { ret[0] = Vec3d(0.0, ret[2].z(), -ret[2].y()).normalized(); break; }
    case 2: { ret[0] = Vec3d(-ret[2].z(), 0.0, ret[2].x()).normalized(); break; }
    }
    ret[1] = ret[2].cross(ret[0]).normalized();
    return ret;
}

Referenced by angle_edge_plane(), and get_measurement().

Here is the caller graph for this function:

plane_normal()

Vec3d Slic3r::Measure::plane_normal ( const SurfaceFeature &  plane )

inline

                                                       {
    assert(plane.get_type() == SurfaceFeatureType::Plane);
    return std::get<1>(plane.get_plane());
}

References Slic3r::Measure::SurfaceFeature::get_plane(), Slic3r::Measure::SurfaceFeature::get_type(), and Plane.

Referenced by are_parallel(), and are_perpendicular().

Here is the call graph for this function:

Here is the caller graph for this function:

Variable Documentation

feature_hover_limit

constexpr double Slic3r::Measure::feature_hover_limit = 0.5

constexpr

Referenced by Slic3r::Measure::MeasuringImpl::get_feature().