3DPrinting/PrusaSlicer-NonPlainar
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

WIP: AABBTreeIndirect - optimized ray_box_intersect_invdir() test,

Browse source 
sandbox for comparing the AABBTreeIndirect with libigl::AABB
This commit is contained in:
Vojtech Bubnik 2020-05-21 17:27:06 +02:00
parent 2b8f655020
commit 99514ba42b
4 changed files with 284 additions and 9 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

1
sandboxes/CMakeLists.txt
View file

@ -2,3 +2,4 @@
#add_subdirectory(openvdb)
add_subdirectory(meshboolean)
add_subdirectory(opencsg)
#add_subdirectory(aabb-evaluation)

2
sandboxes/aabb-evaluation/CMakeLists.txt Normal file
View file

@ -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})

225
sandboxes/aabb-evaluation/aabb-evaluation.cpp Normal file
View file

@ -0,0 +1,225 @@
#include <iostream>
#include <fstream>
#include <string>
#include <libslic3r/TriangleMesh.hpp>
#include <libslic3r/AABBTreeIndirect.hpp>
#include <libslic3r/SLA/EigenMesh3D.hpp>
#include <Shiny/Shiny.h>
#ifdef _MSC_VER
#pragma warning(push)
#pragma warning(disable: 4244)
#pragma warning(disable: 4267)
#endif
#include <igl/ray_mesh_intersect.h>
#include <igl/point_mesh_squared_distance.h>
#include <igl/remove_duplicate_vertices.h>
#include <igl/collapse_small_triangles.h>
#include <igl/signed_distance.h>
#include <igl/random_dir.h>
#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<double>();
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<double>();
const Eigen::Vector3d normal = vertex_normals.row(ivertex).template cast<double>();
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<double>();
const Eigen::Vector3d normal = vertex_normals.row(ivertex).template cast<double>();
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<float>()),
Eigen::Vector3f(d.template cast<float>()), hit))
++ num_hits;
}
occlusion_output0(ivertex) = (double)num_hits/(double)num_samples;
}
}
}
Eigen::MatrixXd occlusion_output1;
{
std::vector<Vec3d> vertices;
std::vector<Vec3i> 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<double>();
const Eigen::Vector3d normal = vertex_normals.row(ivertex).template cast<double>();
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<Eigen::MatrixXd, 3> 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<double>();
const Eigen::Vector3d normal = vertex_normals.row(ivertex).template cast<double>();
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<const Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor | Eigen::DontAlign>> MapMatrixXfUnaligned;
typedef Eigen::Map<const Eigen::Matrix<int, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor | Eigen::DontAlign>> MapMatrixXiUnaligned;
igl::AABB<MapMatrixXfUnaligned, 3> 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<double>();
const Eigen::Vector3d normal = vertex_normals.row(ivertex).template cast<double>();
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<float>(), d.template cast<float>(), 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;
}

65
src/libslic3r/AABBTreeIndirect.hpp
View file

@ -19,8 +19,6 @@ extern "C"
}
// Definition of the ray intersection hit structure.
#include <igl/Hit.h>
// Find intersection parameters of a ray with axis aligned bounding box.
#include <igl/ray_box_intersect.h>
namespace Slic3r {
namespace AABBTreeIndirect {
@ -210,6 +208,11 @@ private:
std::vector<Node> 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<typename AVertexType, typename AIndexedFaceType, typename ATreeType, typename AVectorType>
struct RayIntersector {
@ -224,6 +227,7 @@ namespace detail {
const VectorType origin;
const VectorType dir;
const VectorType invdir;
};
template<typename VertexType, typename IndexedFaceType, typename TreeType, typename VectorType>
@ -231,6 +235,48 @@ namespace detail {
std::vector<igl::Hit> hits;
};
template <typename Derivedsource, typename Deriveddir, typename Scalar>
inline bool ray_box_intersect_invdir(
const Eigen::MatrixBase<Derivedsource> &origin,
const Eigen::MatrixBase<Deriveddir> &inv_dir,
Eigen::AlignedBox<Scalar,3> 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<typename RayIntersectorType, typename Scalar>
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>(), Scalar(0), min_t, t_start, t_end))
if (! ray_box_intersect_invdir(ray_intersector.origin, ray_intersector.invdir, node.bbox.template cast<Scalar>(), 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<Scalar>(),
if (! ray_box_intersect_invdir(ray_intersector.origin, ray_intersector.invdir, node.bbox.template cast<Scalar>(),
Scalar(0), std::numeric_limits<Scalar>::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<VertexType> &vertices,
// Indexed triangle set - triangular faces, references to vertices.
const std::vector<IndexedFaceType> &faces)
const std::vector<IndexedFaceType> &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<InputType> 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<VertexType, IndexedFaceType, TreeType, VectorType> {
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<Scalar>::infinity(), hit);
@ -574,7 +621,7 @@ inline bool intersect_ray_all_hits(
{
auto ray_intersector = detail::RayIntersectorHits<VertexType, IndexedFaceType, TreeType, VectorType> {
vertices, faces, tree,
origin, dir
origin, dir, VectorType(dir.cwiseInverse())
};
if (! tree.empty()) {
ray_intersector.hits.reserve(8);