From 1c76df89eaf9dcbad50601a8a8aabd9b6a7e2f84 Mon Sep 17 00:00:00 2001 From: Vojtech Bubnik Date: Fri, 27 Aug 2021 21:04:11 +0200 Subject: [PATCH] Fix of paint on supports don't work for object that has been scaled up #6718 The triangle-ray intersection function used a hard coded epsilon, which did not work for triangle meshes, that were either too small or too large. Newly the epsilon may be provided to the AABBTreeIndirect search functions externally and IndexedMesh calculates a suitable epsilon on demand from an average triangle mesh edge length. --- src/libslic3r/AABBTreeIndirect.hpp | 125 ++++++++++++++++++++--------- src/libslic3r/SLA/IndexedMesh.cpp | 26 +++--- src/libslic3r/SLA/IndexedMesh.hpp | 8 +- src/libslic3r/TriangleMesh.cpp | 15 ++++ src/libslic3r/TriangleMesh.hpp | 1 + src/slic3r/GUI/MeshUtils.hpp | 2 +- 6 files changed, 127 insertions(+), 50 deletions(-) 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)