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

AABB: Some further polishing and a reference to an SSE implementation

Browse Source 
of the 3D Box vs. ray intersection implementation.
This commit is contained in:
Vojtech Bubnik 2020-05-22 11:35:49 +02:00
parent 925bf1af70
commit ac1f24e5c9
2 changed files with 54 additions and 61 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
sandboxes/aabb-evaluation/aabb-evaluation.cpp
View File

@ -16,7 +16,6 @@
#include <igl/ray_mesh_intersect.h> #include <igl/ray_mesh_intersect.h>
#include <igl/point_mesh_squared_distance.h> #include <igl/point_mesh_squared_distance.h>
#include <igl/remove_duplicate_vertices.h> #include <igl/remove_duplicate_vertices.h>
#include <igl/collapse_small_triangles.h>
#include <igl/signed_distance.h> #include <igl/signed_distance.h>
#include <igl/random_dir.h> #include <igl/random_dir.h>
#ifdef _MSC_VER #ifdef _MSC_VER

58
src/libslic3r/AABBTreeIndirect.hpp
View File

@ -236,26 +236,29 @@ namespace detail {
std::vector<igl::Hit> hits; std::vector<igl::Hit> hits;
}; };
//FIXME implement SSE for float AABB trees with float ray queries.
// SSE/SSE2 is supported by any Intel/AMD x64 processor.
// SSE support requires 16 byte alignment of the AABB nodes, representing the bounding boxes with 4+4 floats,
// storing the node index as the 4th element of the bounding box min value etc.
// https://www.flipcode.com/archives/SSE_RayBox_Intersection_Test.shtml
template <typename Derivedsource, typename Deriveddir, typename Scalar> template <typename Derivedsource, typename Deriveddir, typename Scalar>
inline bool ray_box_intersect_invdir( inline bool ray_box_intersect_invdir(
const Eigen::MatrixBase<Derivedsource> &origin, const Eigen::MatrixBase<Derivedsource> &origin,
const Eigen::MatrixBase<Deriveddir> &inv_dir, const Eigen::MatrixBase<Deriveddir> &inv_dir,
Eigen::AlignedBox<Scalar,3> box, Eigen::AlignedBox<Scalar,3> box,
const Scalar &t0, const Scalar &t0,
const Scalar &t1, const Scalar &t1) {
Scalar &tmin,
Scalar &tmax) {
// http://people.csail.mit.edu/amy/papers/box-jgt.pdf // http://people.csail.mit.edu/amy/papers/box-jgt.pdf
// "An Efficient and Robust Ray–Box Intersection Algorithm" // "An Efficient and Robust Ray–Box Intersection Algorithm"
if (inv_dir.x() < 0) if (inv_dir.x() < 0)
std::swap(box.min().x(), box.max().x()); std::swap(box.min().x(), box.max().x());
if (inv_dir.y() < 0) if (inv_dir.y() < 0)
std::swap(box.min().y(), box.max().y()); std::swap(box.min().y(), box.max().y());
tmin = (box.min().x() - origin.x()) * inv_dir.x(); Scalar tmin = (box.min().x() - origin.x()) * inv_dir.x();
Scalar tymax = (box.max().y() - origin.y()) * inv_dir.y(); Scalar tymax = (box.max().y() - origin.y()) * inv_dir.y();
if (tmin > tymax) if (tmin > tymax)
return false; return false;
tmax = (box.max().x() - origin.x()) * inv_dir.x(); Scalar tmax = (box.max().x() - origin.x()) * inv_dir.x();
Scalar tymin = (box.min().y() - origin.y()) * inv_dir.y(); Scalar tymin = (box.min().y() - origin.y()) * inv_dir.y();
if (tymin > tmax) if (tymin > tmax)
return false; return false;
@ -328,11 +331,8 @@ namespace detail {
const auto &node = ray_intersector.tree.node(node_idx); const auto &node = ray_intersector.tree.node(node_idx);
assert(node.is_valid()); assert(node.is_valid());
{ if (! ray_box_intersect_invdir(ray_intersector.origin, ray_intersector.invdir, node.bbox.template cast<Scalar>(), Scalar(0), min_t))
Scalar 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; return false;
}
if (node.is_leaf()) { if (node.is_leaf()) {
// shoot ray, record hit // shoot ray, record hit
@ -345,10 +345,9 @@ namespace detail {
&& t > 0.) { && t > 0.) {
hit = igl::Hit { int(node.idx), -1, float(u), float(v), float(t) }; hit = igl::Hit { int(node.idx), -1, float(u), float(v), float(t) };
return true; return true;
} } else
return false; return false;
} } else {
// Left / right child node index. // Left / right child node index.
size_t left = node_idx * 2 + 1; size_t left = node_idx * 2 + 1;
size_t right = left + 1; size_t right = left + 1;
@ -367,6 +366,7 @@ namespace detail {
right_ret = false; right_ret = false;
return left_ret || right_ret; return left_ret || right_ret;
} }
}
template<typename RayIntersectorType> template<typename RayIntersectorType>
static inline void intersect_ray_recursive_all_hits(RayIntersectorType &ray_intersector, size_t node_idx) static inline void intersect_ray_recursive_all_hits(RayIntersectorType &ray_intersector, size_t node_idx)
@ -376,34 +376,28 @@ namespace detail {
const auto &node = ray_intersector.tree.node(node_idx); const auto &node = ray_intersector.tree.node(node_idx);
assert(node.is_valid()); assert(node.is_valid());
{
Scalar t_start, t_end;
if (! ray_box_intersect_invdir(ray_intersector.origin, ray_intersector.invdir, 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)) Scalar(0), std::numeric_limits<Scalar>::infinity()))
return; return;
}
if (node.is_leaf()) { if (node.is_leaf()) {
using Vector = Eigen::Matrix<double, 3, 1>;
Vector origin_d = ray_intersector.origin.template cast<double>();
Vector dir_d = ray_intersector.dir .template cast<double>();
auto face = ray_intersector.faces[node.idx]; auto face = ray_intersector.faces[node.idx];
Vector v0 = ray_intersector.vertices[face(0)].template cast<double>();
Vector v1 = ray_intersector.vertices[face(1)].template cast<double>();
Vector v2 = ray_intersector.vertices[face(2)].template cast<double>();
// shoot ray, record hit
double t, u, v; double t, u, v;
if (intersect_triangle1(origin_d.data(), dir_d.data(), v0.data(), v1.data(), v2.data(), &t, &u, &v) && t > 0.) 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 > 0.) {
ray_intersector.hits.emplace_back(igl::Hit{ int(node.idx), -1, float(u), float(v), float(t) }); ray_intersector.hits.emplace_back(igl::Hit{ int(node.idx), -1, float(u), float(v), float(t) });
return;
} }
} else {
// Left / right child node index. // Left / right child node index.
size_t left = node_idx * 2 + 1; size_t left = node_idx * 2 + 1;
size_t right = left + 1; size_t right = left + 1;
intersect_ray_recursive_all_hits(ray_intersector, left); intersect_ray_recursive_all_hits(ray_intersector, left);
intersect_ray_recursive_all_hits(ray_intersector, right); intersect_ray_recursive_all_hits(ray_intersector, right);
} }
}
// Nothing to do with COVID-19 social distancing. // Nothing to do with COVID-19 social distancing.
template<typename AVertexType, typename AIndexedFaceType, typename ATreeType, typename AVectorType> template<typename AVertexType, typename AIndexedFaceType, typename ATreeType, typename AVectorType>
@ -431,17 +425,17 @@ namespace detail {
Vector ap = p - a; Vector ap = p - a;
Scalar d1 = ab.dot(ap); Scalar d1 = ab.dot(ap);
Scalar d2 = ac.dot(ap); Scalar d2 = ac.dot(ap);
if (d1 <= Scalar(0) && d2 <= Scalar(0)) if (d1 <= 0 && d2 <= 0)
return a; return a;
// Check if P in vertex region outside B // Check if P in vertex region outside B
Vector bp = p - b; Vector bp = p - b;
Scalar d3 = ab.dot(bp); Scalar d3 = ab.dot(bp);
Scalar d4 = ac.dot(bp); Scalar d4 = ac.dot(bp);
if (d3 >= Scalar(0) && d4 <= d3) if (d3 >= 0 && d4 <= d3)
return b; return b;
// Check if P in edge region of AB, if so return projection of P onto AB // Check if P in edge region of AB, if so return projection of P onto AB
Scalar vc = d1*d4 - d3*d2; Scalar vc = d1*d4 - d3*d2;
if (a != b && vc <= Scalar(0) && d1 >= Scalar(0) && d3 <= Scalar(0)) { if (a != b && vc <= 0 && d1 >= 0 && d3 <= 0) {
Scalar v = d1 / (d1 - d3); Scalar v = d1 / (d1 - d3);
return a + v * ab; return a + v * ab;
} }
@ -449,17 +443,17 @@ namespace detail {
Vector cp = p - c; Vector cp = p - c;
Scalar d5 = ab.dot(cp); Scalar d5 = ab.dot(cp);
Scalar d6 = ac.dot(cp); Scalar d6 = ac.dot(cp);
if (d6 >= Scalar(0) && d5 <= d6) if (d6 >= 0 && d5 <= d6)
return c; return c;
// Check if P in edge region of AC, if so return projection of P onto AC // Check if P in edge region of AC, if so return projection of P onto AC
Scalar vb = d5*d2 - d1*d6; Scalar vb = d5*d2 - d1*d6;
if (vb <= Scalar(0) && d2 >= Scalar(0) && d6 <= Scalar(0)) { if (vb <= 0 && d2 >= 0 && d6 <= 0) {
Scalar w = d2 / (d2 - d6); Scalar w = d2 / (d2 - d6);
return a + w * ac; return a + w * ac;
} }
// Check if P in edge region of BC, if so return projection of P onto BC // Check if P in edge region of BC, if so return projection of P onto BC
Scalar va = d3*d6 - d5*d4; Scalar va = d3*d6 - d5*d4;
if (va <= Scalar(0) && (d4 - d3) >= Scalar(0) && (d5 - d6) >= Scalar(0)) { if (va <= 0 && (d4 - d3) >= 0 && (d5 - d6) >= 0) {
Scalar w = (d4 - d3) / ((d4 - d3) + (d5 - d6)); Scalar w = (d4 - d3) / ((d4 - d3) + (d5 - d6));
return b + w * (c - b); return b + w * (c - b);
} }