robust ray box intersect (#2244) [ci skip]

* nan-proof min/max

* templates

include/igl/AABB.cpp

@@ -1025,6 +1025,10 @@ namespace igl
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template bool igl::AABB<Eigen::Matrix<float, -1, -1, 0, -1, -1>, 3>::intersect_ray<Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<float, 1, 3, 1, 1, 3> const&, Eigen::Matrix<float, 1, 3, 1, 1, 3> const&, std::vector<igl::Hit, std::allocator<igl::Hit> >&) const;
+// generated by autoexplicit.sh
+template void igl::AABB<Eigen::Matrix<float, -1, -1, 0, -1, -1>, 3>::init<Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&);
+// generated by autoexplicit.sh
 template void igl::AABB<Eigen::Matrix<double, -1, 3, 1, -1, 3>, 3>::squared_distance<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 1, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&) const;
 // generated by autoexplicit.sh
 template void igl::AABB<Eigen::Matrix<double, -1, 3, 1, -1, 3>, 3>::init<Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&);

include/igl/barycenter.cpp

@@ -34,6 +34,8 @@ IGL_INLINE void igl::barycenter(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::barycenter<Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> >&);
+// generated by autoexplicit.sh
 template void igl::barycenter<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
 template void igl::barycenter<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
 template void igl::barycenter<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);

include/igl/ray_box_intersect.cpp

@@ -7,6 +7,7 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "ray_box_intersect.h"
 #include <array>
+#include <igl/matlab_format.h>
 
 template <
   typename Derivedsource,
@@ -105,37 +106,39 @@ IGL_INLINE bool igl::ray_box_intersect(
   // http://people.csail.mit.edu/amy/papers/box-jgt.pdf
   // "An Efficient and Robust Ray–Box Intersection Algorithm"
   Scalar tymin, tymax, tzmin, tzmax;
-  std::array<RowVector3S, 2> bounds = {box.min(),box.max()};
+  std::array<RowVector3S, 2> bounds = {box.min().array(),box.max().array()};
   tmin = ( bounds[sign[0]](0)   - origin(0)) * inv_dir(0);
   tmax = ( bounds[1-sign[0]](0) - origin(0)) * inv_dir(0);
   tymin = (bounds[sign[1]](1)   - origin(1)) * inv_dir(1);
   tymax = (bounds[1-sign[1]](1) - origin(1)) * inv_dir(1);
-  if ( (tmin > tymax) || (tymin > tmax) )
+  // NaN-safe min and max
+  const auto berger_perrin_min = [&](
+      const Scalar a,
+      const Scalar b) -> Scalar
   {
-    return false;
-  }
-  if (tymin > tmin)
+    return (a < b) ? a : b;
+  };
+  const auto berger_perrin_max = [&](
+      const Scalar a,
+      const Scalar b) -> Scalar
   {
-    tmin = tymin;
-  }
-  if (tymax < tmax)
+    return (a > b) ? a : b;
+  };
+
+  if ( (tmin > tymax) || (tymin > tmax) )
   {
-    tmax = tymax;
+    return false;
   }
+  tmin = berger_perrin_max(tmin,tymin);
+  tmax = berger_perrin_min(tmax,tymax);
   tzmin = (bounds[sign[2]](2) - origin(2))   * inv_dir(2);
   tzmax = (bounds[1-sign[2]](2) - origin(2)) * inv_dir(2);
   if ( (tmin > tzmax) || (tzmin > tmax) )
   {
     return false;
   }
-  if (tzmin > tmin)
-  {
-    tmin = tzmin;
-  }
-  if (tzmax < tmax)
-  {
-    tmax = tzmax;
-  }
+  tmin = berger_perrin_max(tmin,tzmin);
+  tmax = berger_perrin_min(tmax,tzmax);
   if(!( (tmin < t1) && (tmax > t0) ))
   {
     return false;
@@ -146,5 +149,7 @@ IGL_INLINE bool igl::ray_box_intersect(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+// generated by autoexplicit.sh
+template bool igl::ray_box_intersect<Eigen::Matrix<float, 1, 3, 1, 1, 3>, Eigen::Matrix<float, 1, 3, 1, 1, 3>, float>(Eigen::MatrixBase<Eigen::Matrix<float, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 1, 3, 1, 1, 3> > const&, Eigen::AlignedBox<float, 3> const&, float const&, float const&, float&, float&);
 template bool igl::ray_box_intersect<Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, double>(Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::AlignedBox<double, 3> const&, double const&, double const&, double&, double&);
 #endif

include/igl/ray_mesh_intersect.cpp

@@ -119,6 +119,10 @@ IGL_INLINE bool ray_mesh_intersect(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template bool igl::ray_mesh_intersect<Eigen::Matrix<float, 1, 3, 1, 1, 3>, Eigen::Matrix<float, 1, 3, 1, 1, 3>, Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Block<Eigen::Matrix<int, -1, -1, 0, -1, -1> const, 1, -1, false> >(Eigen::MatrixBase<Eigen::Matrix<float, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<int, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, std::vector<igl::Hit, std::allocator<igl::Hit> >&);
+// generated by autoexplicit.sh
+template bool igl::ray_mesh_intersect<Eigen::Matrix<float, 3, 1, 0, 3, 1>, Eigen::Matrix<float, 3, 1, 0, 3, 1>, Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<float, 3, 1, 0, 3, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 3, 1, 0, 3, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<igl::Hit, std::allocator<igl::Hit> >&);
+// generated by autoexplicit.sh
 template bool igl::ray_mesh_intersect<Eigen::Matrix<float, 3, 1, 0, 3, 1>, Eigen::Matrix<float, 3, 1, 0, 3, 1>, Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<float, 3, 1, 0, 3, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 3, 1, 0, 3, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, std::vector<igl::Hit, std::allocator<igl::Hit> >&);
 template bool igl::ray_mesh_intersect<Eigen::Matrix<float, 3, 1, 0, 3, 1>, Eigen::Matrix<float, 3, 1, 0, 3, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<float, 3, 1, 0, 3, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 3, 1, 0, 3, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<igl::Hit, std::allocator<igl::Hit> >&);
 template bool igl::ray_mesh_intersect<Eigen::Matrix<float, 3, 1, 0, 3, 1>, Eigen::Matrix<float, 3, 1, 0, 3, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<float, 3, 1, 0, 3, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 3, 1, 0, 3, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::Hit&);

tests/include/igl/ray_mesh_intersect.cpp

@@ -1,5 +1,6 @@
 #include <test_common.h>
 #include <igl/ray_mesh_intersect.h>
+#include <igl/AABB.h>
 
 TEST_CASE("ray_mesh_intersect: one_triangle", "[igl]")
 {
@@ -23,3 +24,52 @@ TEST_CASE("ray_mesh_intersect: one_triangle", "[igl]")
   REQUIRE(hits.size() == 1);
   REQUIRE(hits.front().t == Approx(1.0));
 }
+
+
+TEST_CASE("ray_mesh_intersect: corner-case", "[igl]")
+{
+  //       .      //
+  //      /|\     //
+  //     / | \    //
+  //    /  |  \   //
+  //   /   |   \  //
+  //  .----x----. //      
+  //   \   |   /  //               
+  //    \  |  /   //               
+  //     \ | /    //               
+  //      \|/     //                
+  //       .      //
+
+  Eigen::MatrixXf vertices(5, 3);
+  vertices <<
+      -1., 0., 0.,
+      0., 1., 0.,
+      1., 0., 0.,
+      0., -1., 0.,
+      0., 0., 0.;
+
+  Eigen::MatrixXi faces(4, 3);
+  faces <<
+      0, 1, 4,
+      1, 2, 4,
+      2, 3, 4,
+      3, 0, 4;
+
+  igl::AABB<Eigen::MatrixXf, 3> mesh_bvh;
+
+  mesh_bvh.init(vertices, faces);
+
+  for (float eps: {1e-5f, 0.f})
+  {
+    Eigen::Vector3f origin(eps, eps, 1.f + eps);
+    Eigen::Vector3f direction(0.f, 0.f, -1.f);
+
+    std::vector<igl::Hit> hits, hits_bvh;
+    bool is_hit = igl::ray_mesh_intersect(origin, direction, vertices, faces, hits);
+    bool is_hit_bvh = mesh_bvh.intersect_ray(vertices, faces, origin, direction, hits_bvh);
+
+    REQUIRE (is_hit);
+    REQUIRE (is_hit == is_hit_bvh);
+    REQUIRE (hits.size() == hits_bvh.size());
+  }
+}