From 99514ba42b4b18f36236a1f1546b783b75c6d2ec Mon Sep 17 00:00:00 2001 From: Vojtech Bubnik Date: Thu, 21 May 2020 17:27:06 +0200 Subject: [PATCH] WIP: AABBTreeIndirect - optimized ray_box_intersect_invdir() test, sandbox for comparing the AABBTreeIndirect with libigl::AABB --- sandboxes/CMakeLists.txt | 1 + sandboxes/aabb-evaluation/CMakeLists.txt | 2 + sandboxes/aabb-evaluation/aabb-evaluation.cpp | 225 ++++++++++++++++++ src/libslic3r/AABBTreeIndirect.hpp | 65 ++++- 4 files changed, 284 insertions(+), 9 deletions(-) create mode 100644 sandboxes/aabb-evaluation/CMakeLists.txt create mode 100644 sandboxes/aabb-evaluation/aabb-evaluation.cpp diff --git a/sandboxes/CMakeLists.txt b/sandboxes/CMakeLists.txt index a2bd13bb0..23c15f089 100644 --- a/sandboxes/CMakeLists.txt +++ b/sandboxes/CMakeLists.txt @@ -2,3 +2,4 @@ #add_subdirectory(openvdb) add_subdirectory(meshboolean) add_subdirectory(opencsg) +#add_subdirectory(aabb-evaluation) \ No newline at end of file diff --git a/sandboxes/aabb-evaluation/CMakeLists.txt b/sandboxes/aabb-evaluation/CMakeLists.txt new file mode 100644 index 000000000..20011e345 --- /dev/null +++ b/sandboxes/aabb-evaluation/CMakeLists.txt @@ -0,0 +1,2 @@ +add_executable(aabb-evaluation aabb-evaluation.cpp) +target_link_libraries(aabb-evaluation libslic3r ${Boost_LIBRARIES} ${TBB_LIBRARIES} ${Boost_LIBRARIES} ${CMAKE_DL_LIBS}) diff --git a/sandboxes/aabb-evaluation/aabb-evaluation.cpp b/sandboxes/aabb-evaluation/aabb-evaluation.cpp new file mode 100644 index 000000000..b81fbc96c --- /dev/null +++ b/sandboxes/aabb-evaluation/aabb-evaluation.cpp @@ -0,0 +1,225 @@ +#include +#include +#include + +#include +#include +#include + +#include + +#ifdef _MSC_VER +#pragma warning(push) +#pragma warning(disable: 4244) +#pragma warning(disable: 4267) +#endif +#include +#include +#include +#include +#include +#include +#ifdef _MSC_VER +#pragma warning(pop) +#endif + +const std::string USAGE_STR = { + "Usage: aabb-evaluation stlfilename.stl" +}; + +using namespace Slic3r; + +void profile(const TriangleMesh &mesh) +{ + Eigen::MatrixXd V; + Eigen::MatrixXi F; + Eigen::MatrixXd vertex_normals; + sla::to_eigen_mesh(mesh, V, F); + igl::per_vertex_normals(V, F, vertex_normals); + + static constexpr int num_samples = 100; + const int num_vertices = std::min(10000, int(mesh.its.vertices.size())); + const Eigen::MatrixXd dirs = igl::random_dir_stratified(num_samples).cast(); + + Eigen::MatrixXd occlusion_output0; + { + AABBTreeIndirect::Tree3f tree; + { + PROFILE_BLOCK(AABBIndirect_Init); + tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(mesh.its.vertices, mesh.its.indices); + } + { + PROFILE_BLOCK(EigenMesh3D_AABBIndirectF_AmbientOcclusion); + occlusion_output0.resize(num_vertices, 1); + for (int ivertex = 0; ivertex < num_vertices; ++ ivertex) { + const Eigen::Vector3d origin = mesh.its.vertices[ivertex].template cast(); + const Eigen::Vector3d normal = vertex_normals.row(ivertex).template cast(); + int num_hits = 0; + for (int s = 0; s < num_samples; s++) { + Eigen::Vector3d d = dirs.row(s); + if(d.dot(normal) < 0) { + // reverse ray + d *= -1; + } + igl::Hit hit; + if (AABBTreeIndirect::intersect_ray_first_hit(mesh.its.vertices, mesh.its.indices, tree, (origin + 1e-4 * d).eval(), d, hit)) + ++ num_hits; + } + occlusion_output0(ivertex) = (double)num_hits/(double)num_samples; + } + } + + { + PROFILE_BLOCK(EigenMesh3D_AABBIndirectFF_AmbientOcclusion); + occlusion_output0.resize(num_vertices, 1); + for (int ivertex = 0; ivertex < num_vertices; ++ ivertex) { + const Eigen::Vector3d origin = mesh.its.vertices[ivertex].template cast(); + const Eigen::Vector3d normal = vertex_normals.row(ivertex).template cast(); + int num_hits = 0; + for (int s = 0; s < num_samples; s++) { + Eigen::Vector3d d = dirs.row(s); + if(d.dot(normal) < 0) { + // reverse ray + d *= -1; + } + igl::Hit hit; + if (AABBTreeIndirect::intersect_ray_first_hit(mesh.its.vertices, mesh.its.indices, tree, + Eigen::Vector3f((origin + 1e-4 * d).template cast()), + Eigen::Vector3f(d.template cast()), hit)) + ++ num_hits; + } + occlusion_output0(ivertex) = (double)num_hits/(double)num_samples; + } + } + } + + Eigen::MatrixXd occlusion_output1; + { + std::vector vertices; + std::vector triangles; + for (int i = 0; i < V.rows(); ++ i) + vertices.emplace_back(V.row(i).transpose()); + for (int i = 0; i < F.rows(); ++ i) + triangles.emplace_back(F.row(i).transpose()); + AABBTreeIndirect::Tree3d tree; + { + PROFILE_BLOCK(AABBIndirectD_Init); + tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(vertices, triangles); + } + + { + PROFILE_BLOCK(EigenMesh3D_AABBIndirectD_AmbientOcclusion); + occlusion_output1.resize(num_vertices, 1); + for (int ivertex = 0; ivertex < num_vertices; ++ ivertex) { + const Eigen::Vector3d origin = V.row(ivertex).template cast(); + const Eigen::Vector3d normal = vertex_normals.row(ivertex).template cast(); + int num_hits = 0; + for (int s = 0; s < num_samples; s++) { + Eigen::Vector3d d = dirs.row(s); + if(d.dot(normal) < 0) { + // reverse ray + d *= -1; + } + igl::Hit hit; + if (AABBTreeIndirect::intersect_ray_first_hit(vertices, triangles, tree, Eigen::Vector3d(origin + 1e-4 * d), d, hit)) + ++ num_hits; + } + occlusion_output1(ivertex) = (double)num_hits/(double)num_samples; + } + } + } + + // Build the AABB accelaration tree + + Eigen::MatrixXd occlusion_output2; + { + igl::AABB AABB; + { + PROFILE_BLOCK(EigenMesh3D_AABB_Init); + AABB.init(V, F); + } + { + PROFILE_BLOCK(EigenMesh3D_AABB_AmbientOcclusion); + occlusion_output2.resize(num_vertices, 1); + for (int ivertex = 0; ivertex < num_vertices; ++ ivertex) { + const Eigen::Vector3d origin = V.row(ivertex).template cast(); + const Eigen::Vector3d normal = vertex_normals.row(ivertex).template cast(); + int num_hits = 0; + for (int s = 0; s < num_samples; s++) { + Eigen::Vector3d d = dirs.row(s); + if(d.dot(normal) < 0) { + // reverse ray + d *= -1; + } + igl::Hit hit; + if (AABB.intersect_ray(V, F, origin + 1e-4 * d, d, hit)) + ++ num_hits; + } + occlusion_output2(ivertex) = (double)num_hits/(double)num_samples; + } + } + } + + Eigen::MatrixXd occlusion_output3; + { + typedef Eigen::Map> MapMatrixXfUnaligned; + typedef Eigen::Map> MapMatrixXiUnaligned; + igl::AABB AABB; + auto vertices = MapMatrixXfUnaligned(mesh.its.vertices.front().data(), mesh.its.vertices.size(), 3); + auto faces = MapMatrixXiUnaligned(mesh.its.indices.front().data(), mesh.its.indices.size(), 3); + { + PROFILE_BLOCK(EigenMesh3D_AABBf_Init); + AABB.init( + vertices, + faces); + } + + { + PROFILE_BLOCK(EigenMesh3D_AABBf_AmbientOcclusion); + occlusion_output3.resize(num_vertices, 1); + for (int ivertex = 0; ivertex < num_vertices; ++ ivertex) { + const Eigen::Vector3d origin = mesh.its.vertices[ivertex].template cast(); + const Eigen::Vector3d normal = vertex_normals.row(ivertex).template cast(); + int num_hits = 0; + for (int s = 0; s < num_samples; s++) { + Eigen::Vector3d d = dirs.row(s); + if(d.dot(normal) < 0) { + // reverse ray + d *= -1; + } + igl::Hit hit; + if (AABB.intersect_ray(vertices, faces, (origin + 1e-4 * d).eval().template cast(), d.template cast(), hit)) + ++ num_hits; + } + occlusion_output3(ivertex) = (double)num_hits/(double)num_samples; + } + } + } + + PROFILE_UPDATE(); + PROFILE_OUTPUT(nullptr); +} + +int main(const int argc, const char *argv[]) +{ + if(argc < 2) { + std::cout << USAGE_STR << std::endl; + return EXIT_SUCCESS; + } + + TriangleMesh mesh; + if (! mesh.ReadSTLFile(argv[1])) { + std::cerr << "Error loading " << argv[1] << std::endl; + return -1; + } + + mesh.repair(); + if (mesh.facets_count() == 0) { + std::cerr << "Error loading " << argv[1] << " . It is empty." << std::endl; + return -1; + } + + profile(mesh); + + return EXIT_SUCCESS; +} diff --git a/src/libslic3r/AABBTreeIndirect.hpp b/src/libslic3r/AABBTreeIndirect.hpp index 134b15b71..8d9dfbdef 100644 --- a/src/libslic3r/AABBTreeIndirect.hpp +++ b/src/libslic3r/AABBTreeIndirect.hpp @@ -19,8 +19,6 @@ extern "C" } // Definition of the ray intersection hit structure. #include -// Find intersection parameters of a ray with axis aligned bounding box. -#include namespace Slic3r { namespace AABBTreeIndirect { @@ -210,6 +208,11 @@ private: std::vector m_nodes; }; +using Tree2f = Tree<2, float>; +using Tree3f = Tree<3, float>; +using Tree2d = Tree<2, double>; +using Tree3d = Tree<3, double>; + namespace detail { template struct RayIntersector { @@ -224,6 +227,7 @@ namespace detail { const VectorType origin; const VectorType dir; + const VectorType invdir; }; template @@ -231,6 +235,48 @@ namespace detail { std::vector hits; }; + template + inline bool ray_box_intersect_invdir( + const Eigen::MatrixBase &origin, + const Eigen::MatrixBase &inv_dir, + Eigen::AlignedBox box, + const Scalar &t0, + const Scalar &t1, + Scalar &tmin, + Scalar &tmax) { + // http://people.csail.mit.edu/amy/papers/box-jgt.pdf + // "An Efficient and Robust Ray–Box Intersection Algorithm" + if (inv_dir.x() < 0) + std::swap(box.min().x(), box.max().x()); + if (inv_dir.y() < 0) + std::swap(box.min().y(), box.max().y()); + tmin = (box.min().x() - origin.x()) * inv_dir.x(); + Scalar tymax = (box.max().y() - origin.y()) * inv_dir.y(); + if (tmin > tymax) + return false; + tmax = (box.max().x() - origin.x()) * inv_dir.x(); + Scalar tymin = (box.min().y() - origin.y()) * inv_dir.y(); + if (tymin > tmax) + return false; + if (tymin > tmin) + tmin = tymin; + if (tymax < tmax) + tmax = tymax; + if (inv_dir.z() < 0) + std::swap(box.min().z(), box.max().z()); + Scalar tzmin = (box.min().z() - origin.z()) * inv_dir.z(); + if (tzmin > tmax) + return false; + Scalar tzmax = (box.max().z() - origin.z()) * inv_dir.z(); + if (tmin > tzmax) + return false; + if (tzmin > tmin) + tmin = tzmin; + if (tzmax < tmax) + tmax = tzmax; + return tmin < t1 && tmax > t0; + } + template static inline bool intersect_ray_recursive_first_hit( RayIntersectorType &ray_intersector, @@ -243,7 +289,7 @@ namespace detail { { Scalar t_start, t_end; - if (! igl::ray_box_intersect(ray_intersector.origin, ray_intersector.dir, node.bbox.template cast(), Scalar(0), min_t, t_start, t_end)) + if (! ray_box_intersect_invdir(ray_intersector.origin, ray_intersector.invdir, node.bbox.template cast(), Scalar(0), min_t, t_start, t_end)) return false; } @@ -293,7 +339,7 @@ namespace detail { { Scalar t_start, t_end; - if (! igl::ray_box_intersect(ray_intersector.origin, ray_intersector.dir, node.bbox.template cast(), + if (! ray_box_intersect_invdir(ray_intersector.origin, ray_intersector.invdir, node.bbox.template cast(), Scalar(0), std::numeric_limits::infinity(), t_start, t_end)) return; } @@ -483,13 +529,14 @@ inline Tree<3, typename VertexType::Scalar> build_aabb_tree_over_indexed_triangl // Indexed triangle set - 3D vertices. const std::vector &vertices, // Indexed triangle set - triangular faces, references to vertices. - const std::vector &faces) + const std::vector &faces, + //FIXME do we want to apply an epsilon? + const typename VertexType::Scalar eps = 0) { using TreeType = Tree<3, typename VertexType::Scalar>; using CoordType = typename TreeType::CoordType; using VectorType = typename TreeType::VectorType; using BoundingBox = typename TreeType::BoundingBox; - static constexpr CoordType eps = CoordType(1e-4); struct InputType { size_t idx() const { return m_idx; } @@ -503,7 +550,7 @@ inline Tree<3, typename VertexType::Scalar> build_aabb_tree_over_indexed_triangl std::vector input; input.reserve(faces.size()); - VectorType veps(eps, eps, eps); + const VectorType veps(eps, eps, eps); for (size_t i = 0; i < faces.size(); ++ i) { const IndexedFaceType &face = faces[i]; const VertexType &v1 = vertices[face(0)]; @@ -546,7 +593,7 @@ inline bool intersect_ray_first_hit( using Scalar = typename VectorType::Scalar; auto ray_intersector = detail::RayIntersector { vertices, faces, tree, - origin, dir + origin, dir, VectorType(dir.cwiseInverse()) }; return ! tree.empty() && detail::intersect_ray_recursive_first_hit( ray_intersector, size_t(0), std::numeric_limits::infinity(), hit); @@ -574,7 +621,7 @@ inline bool intersect_ray_all_hits( { auto ray_intersector = detail::RayIntersectorHits { vertices, faces, tree, - origin, dir + origin, dir, VectorType(dir.cwiseInverse()) }; if (! tree.empty()) { ray_intersector.hits.reserve(8);