diff --git a/src/libslic3r/AABBTreeIndirect.hpp b/src/libslic3r/AABBTreeIndirect.hpp index 76aa36194..217166f8c 100644 --- a/src/libslic3r/AABBTreeIndirect.hpp +++ b/src/libslic3r/AABBTreeIndirect.hpp @@ -15,11 +15,6 @@ #include "Utils.hpp" // for next_highest_power_of_2() -extern "C" -{ -// Ray-Triangle Intersection Test Routines by Tomas Moller, May 2000 -#include -} // Definition of the ray intersection hit structure. #include @@ -231,6 +226,9 @@ namespace detail { const VectorType origin; const VectorType dir; const VectorType invdir; + + // epsilon for ray-triangle intersection, see intersect_triangle1() + const double eps; }; template @@ -283,44 +281,91 @@ namespace detail { return tmin < t1 && tmax > t0; } + // The following intersect_triangle() is derived from raytri.c routine intersect_triangle1() + // Ray-Triangle Intersection Test Routines + // Different optimizations of my and Ben Trumbore's + // code from journals of graphics tools (JGT) + // http://www.acm.org/jgt/ + // by Tomas Moller, May 2000 template - std::enable_if_t::value && std::is_same::value, bool> - intersect_triangle(const V &origin, const V &dir, const W &v0, const W &v1, const W &v2, double &t, double &u, double &v) { - return intersect_triangle1(const_cast(origin.data()), const_cast(dir.data()), - const_cast(v0.data()), const_cast(v1.data()), const_cast(v2.data()), - &t, &u, &v); + std::enable_if_t::value&& std::is_same::value, bool> + intersect_triangle(const V &orig, const V &dir, const W &vert0, const W &vert1, const W &vert2, double &t, double &u, double &v, double eps) + { + // find vectors for two edges sharing vert0 + const V edge1 = vert1 - vert0; + const V edge2 = vert2 - vert0; + // begin calculating determinant - also used to calculate U parameter + const V pvec = dir.cross(edge2); + // if determinant is near zero, ray lies in plane of triangle + const double det = edge1.dot(pvec); + V qvec; + + if (det > eps) { + // calculate distance from vert0 to ray origin + V tvec = orig - vert0; + // calculate U parameter and test bounds + u = tvec.dot(pvec); + if (u < 0.0 || u > det) + return false; + // prepare to test V parameter + qvec = tvec.cross(edge1); + // calculate V parameter and test bounds + v = dir.dot(qvec); + if (v < 0.0 || u + v > det) + return false; + } else if (det < -eps) { + // calculate distance from vert0 to ray origin + V tvec = orig - vert0; + // calculate U parameter and test bounds + u = tvec.dot(pvec); + if (u > 0.0 || u < det) + return false; + // prepare to test V parameter + qvec = tvec.cross(edge1); + // calculate V parameter and test bounds + v = dir.dot(qvec); + if (v > 0.0 || u + v < det) + return false; + } else + // ray is parallel to the plane of the triangle + return false; + + double inv_det = 1.0 / det; + // calculate t, ray intersects triangle + t = edge2.dot(qvec) * inv_det; + u *= inv_det; + v *= inv_det; + return true; } template std::enable_if_t::value && !std::is_same::value, bool> - intersect_triangle(const V &origin, const V &dir, const W &v0, const W &v1, const W &v2, double &t, double &u, double &v) { - using Vector = Eigen::Matrix; - Vector w0 = v0.template cast(); - Vector w1 = v1.template cast(); - Vector w2 = v2.template cast(); - return intersect_triangle1(const_cast(origin.data()), const_cast(dir.data()), - w0.data(), w1.data(), w2.data(), &t, &u, &v); + intersect_triangle(const V &origin, const V &dir, const W &v0, const W &v1, const W &v2, double &t, double &u, double &v, double eps) { + return intersect_triangle(origin, dir, v0.template cast(), v1.template cast(), v2.template cast(), t, u, v, eps); } template std::enable_if_t::value && std::is_same::value, bool> - intersect_triangle(const V &origin, const V &dir, const W &v0, const W &v1, const W &v2, double &t, double &u, double &v) { - using Vector = Eigen::Matrix; - Vector o = origin.template cast(); - Vector d = dir.template cast(); - return intersect_triangle1(o.data(), d.data(), const_cast(v0.data()), const_cast(v1.data()), const_cast(v2.data()), &t, &u, &v); + intersect_triangle(const V &origin, const V &dir, const W &v0, const W &v1, const W &v2, double &t, double &u, double &v, double eps) { + return intersect_triangle(origin.template cast(), dir.template cast(), v0, v1, v2, t, u, v, eps); } template std::enable_if_t::value && ! std::is_same::value, bool> - intersect_triangle(const V &origin, const V &dir, const W &v0, const W &v1, const W &v2, double &t, double &u, double &v) { - using Vector = Eigen::Matrix; - Vector o = origin.template cast(); - Vector d = dir.template cast(); - Vector w0 = v0.template cast(); - Vector w1 = v1.template cast(); - Vector w2 = v2.template cast(); - return intersect_triangle1(o.data(), d.data(), w0.data(), w1.data(), w2.data(), &t, &u, &v); + intersect_triangle(const V &origin, const V &dir, const W &v0, const W &v1, const W &v2, double &t, double &u, double &v, double eps) { + return intersect_triangle(origin.template cast(), dir.template cast(), v0.template cast(), v1.template cast(), v2.template cast(), t, u, v, eps); + } + + template + double intersect_triangle_epsilon(const Tree &tree) { + double eps = 0.000001; + if (! tree.empty()) { + const typename Tree::BoundingBox &bbox = tree.nodes().front().bbox; + double l = (bbox.max() - bbox.min()).cwiseMax(); + if (l > 0) + eps /= (l * l); + } + return eps; } template @@ -343,7 +388,7 @@ namespace detail { if (intersect_triangle( ray_intersector.origin, ray_intersector.dir, ray_intersector.vertices[face(0)], ray_intersector.vertices[face(1)], ray_intersector.vertices[face(2)], - t, u, v) + t, u, v, ray_intersector.eps) && t > 0.) { hit = igl::Hit { int(node.idx), -1, float(u), float(v), float(t) }; return true; @@ -388,7 +433,7 @@ namespace detail { if (intersect_triangle( ray_intersector.origin, ray_intersector.dir, ray_intersector.vertices[face(0)], ray_intersector.vertices[face(1)], ray_intersector.vertices[face(2)], - t, u, v) + t, u, v, ray_intersector.eps) && t > 0.) { ray_intersector.hits.emplace_back(igl::Hit{ int(node.idx), -1, float(u), float(v), float(t) }); } @@ -623,12 +668,15 @@ inline bool intersect_ray_first_hit( // Direction of the ray. const VectorType &dir, // First intersection of the ray with the indexed triangle set. - igl::Hit &hit) + igl::Hit &hit, + // Epsilon for the ray-triangle intersection, it should be proportional to an average triangle edge length. + const double eps = 0.000001) { using Scalar = typename VectorType::Scalar; - auto ray_intersector = detail::RayIntersector { + auto ray_intersector = detail::RayIntersector { vertices, faces, tree, - origin, dir, VectorType(dir.cwiseInverse()) + origin, dir, VectorType(dir.cwiseInverse()), + eps }; return ! tree.empty() && detail::intersect_ray_recursive_first_hit( ray_intersector, size_t(0), std::numeric_limits::infinity(), hit); @@ -652,11 +700,14 @@ inline bool intersect_ray_all_hits( // Direction of the ray. const VectorType &dir, // All intersections of the ray with the indexed triangle set, sorted by parameter t. - std::vector &hits) + std::vector &hits, + // Epsilon for the ray-triangle intersection, it should be proportional to an average triangle edge length. + const double eps = 0.000001) { auto ray_intersector = detail::RayIntersectorHits { { vertices, faces, {tree}, - origin, dir, VectorType(dir.cwiseInverse()) } + origin, dir, VectorType(dir.cwiseInverse()), + eps } }; if (! tree.empty()) { ray_intersector.hits.reserve(8); diff --git a/src/libslic3r/SLA/IndexedMesh.cpp b/src/libslic3r/SLA/IndexedMesh.cpp index 887ef1555..07c4203ab 100644 --- a/src/libslic3r/SLA/IndexedMesh.cpp +++ b/src/libslic3r/SLA/IndexedMesh.cpp @@ -17,10 +17,18 @@ namespace sla { class IndexedMesh::AABBImpl { private: AABBTreeIndirect::Tree3f m_tree; + double m_triangle_ray_epsilon; public: - void init(const indexed_triangle_set &its) + void init(const indexed_triangle_set &its, bool calculate_epsilon) { + m_triangle_ray_epsilon = 0.000001; + if (calculate_epsilon) { + // Calculate epsilon from average triangle edge length. + double l = its_average_edge_length(its); + if (l > 0) + m_triangle_ray_epsilon = 0.000001 * l * l; + } m_tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set( its.vertices, its.indices); } @@ -31,7 +39,7 @@ public: igl::Hit & hit) { AABBTreeIndirect::intersect_ray_first_hit(its.vertices, its.indices, - m_tree, s, dir, hit); + m_tree, s, dir, hit, m_triangle_ray_epsilon); } void intersect_ray(const indexed_triangle_set &its, @@ -40,7 +48,7 @@ public: std::vector & hits) { AABBTreeIndirect::intersect_ray_all_hits(its.vertices, its.indices, - m_tree, s, dir, hits); + m_tree, s, dir, hits, m_triangle_ray_epsilon); } double squared_distance(const indexed_triangle_set & its, @@ -60,25 +68,25 @@ public: } }; -template void IndexedMesh::init(const M &mesh) +template void IndexedMesh::init(const M &mesh, bool calculate_epsilon) { BoundingBoxf3 bb = bounding_box(mesh); m_ground_level += bb.min(Z); // Build the AABB accelaration tree - m_aabb->init(*m_tm); + m_aabb->init(*m_tm, calculate_epsilon); } -IndexedMesh::IndexedMesh(const indexed_triangle_set& tmesh) +IndexedMesh::IndexedMesh(const indexed_triangle_set& tmesh, bool calculate_epsilon) : m_aabb(new AABBImpl()), m_tm(&tmesh) { - init(tmesh); + init(tmesh, calculate_epsilon); } -IndexedMesh::IndexedMesh(const TriangleMesh &mesh) +IndexedMesh::IndexedMesh(const TriangleMesh &mesh, bool calculate_epsilon) : m_aabb(new AABBImpl()), m_tm(&mesh.its) { - init(mesh); + init(mesh, calculate_epsilon); } IndexedMesh::~IndexedMesh() {} diff --git a/src/libslic3r/SLA/IndexedMesh.hpp b/src/libslic3r/SLA/IndexedMesh.hpp index 25ab75b83..9348a97c9 100644 --- a/src/libslic3r/SLA/IndexedMesh.hpp +++ b/src/libslic3r/SLA/IndexedMesh.hpp @@ -42,12 +42,14 @@ class IndexedMesh { std::vector m_holes; #endif - template void init(const M &mesh); + template void init(const M &mesh, bool calculate_epsilon); public: - explicit IndexedMesh(const indexed_triangle_set&); - explicit IndexedMesh(const TriangleMesh &mesh); + // calculate_epsilon ... calculate epsilon for triangle-ray intersection from an average triangle edge length. + // If set to false, a default epsilon is used, which works for "reasonable" meshes. + explicit IndexedMesh(const indexed_triangle_set &tmesh, bool calculate_epsilon = false); + explicit IndexedMesh(const TriangleMesh &mesh, bool calculate_epsilon = false); IndexedMesh(const IndexedMesh& other); IndexedMesh& operator=(const IndexedMesh&); diff --git a/src/libslic3r/TriangleMesh.cpp b/src/libslic3r/TriangleMesh.cpp index d4baabc97..fa8f8bce6 100644 --- a/src/libslic3r/TriangleMesh.cpp +++ b/src/libslic3r/TriangleMesh.cpp @@ -1275,6 +1275,21 @@ float its_volume(const indexed_triangle_set &its) return volume; } +float its_average_edge_length(const indexed_triangle_set &its) +{ + if (its.indices.empty()) + return 0.f; + + double edge_length = 0.f; + for (size_t i = 0; i < its.indices.size(); ++ i) { + const its_triangle v = its_triangle_vertices(its, i); + edge_length += (v[1] - v[0]).cast().norm() + + (v[2] - v[0]).cast().norm() + + (v[1] - v[2]).cast().norm(); + } + return float(edge_length / (3 * its.indices.size())); +} + std::vector its_split(const indexed_triangle_set &its) { return its_split<>(its); diff --git a/src/libslic3r/TriangleMesh.hpp b/src/libslic3r/TriangleMesh.hpp index b7a1bebb1..c463af5a2 100644 --- a/src/libslic3r/TriangleMesh.hpp +++ b/src/libslic3r/TriangleMesh.hpp @@ -199,6 +199,7 @@ inline stl_normal its_unnormalized_normal(const indexed_triangle_set &its, } float its_volume(const indexed_triangle_set &its); +float its_average_edge_length(const indexed_triangle_set &its); void its_merge(indexed_triangle_set &A, const indexed_triangle_set &B); void its_merge(indexed_triangle_set &A, const std::vector &triangles); diff --git a/src/slic3r/GUI/MeshUtils.hpp b/src/slic3r/GUI/MeshUtils.hpp index ec6c337c0..65c326116 100644 --- a/src/slic3r/GUI/MeshUtils.hpp +++ b/src/slic3r/GUI/MeshUtils.hpp @@ -112,7 +112,7 @@ public: // The class references extern TriangleMesh, which must stay alive // during MeshRaycaster existence. MeshRaycaster(const TriangleMesh& mesh) - : m_emesh(mesh) + : m_emesh(mesh, true) // calculate epsilon for triangle-ray intersection from an average edge length { m_normals.reserve(mesh.stl.facet_start.size()); for (const stl_facet& facet : mesh.stl.facet_start)