Add test for issue/2270 and a fix (#2271)

* Add test for issue/2270 and a (dirty) fix

* restore func; simplify test (failing)

* more robust fix

* Implement solution from T.Ize in "Robust BVH Ray Traversal" 4.2

* Implement solution from T.Ize in "Robust BVH Ray Traversal" 4.1

* Move nextafter in its own file

* Add std:: to other nextafter calls

* Add test for nextafter

* Use unaryExpr for nextafter on Eigen elements

* Fix and rename into increment_ulp

* missing templates

* correct overload name; missing templates

---------

Co-authored-by: Alec Jacobson <[email protected]>

include/igl/AABB.cpp

@@ -10,6 +10,7 @@
 #include "barycenter.h"
 #include "colon.h"
 #include "doublearea.h"
+#include "increment_ulp.h"
 #include "point_simplex_squared_distance.h"
 #include "project_to_line_segment.h"
 #include "sort.h"
@@ -24,6 +25,7 @@
 #include <queue>
 #include <stack>
 
+
 template <typename DerivedV, int DIM>
 template <typename DerivedEle, typename Derivedbb_mins, typename Derivedbb_maxs, typename Derivedelements>
 IGL_INLINE void igl::AABB<DerivedV,DIM>::init(
@@ -831,37 +833,10 @@ igl::AABB<DerivedV,DIM>::intersect_ray(
   const RowVectorDIMS & dir,
   std::vector<igl::Hit> & hits) const
 {
-  hits.clear();
-  const Scalar t0 = 0;
-  const Scalar t1 = std::numeric_limits<Scalar>::infinity();
-  {
-    Scalar _1,_2;
-    if(!ray_box_intersect(origin,dir,m_box,t0,t1,_1,_2))
-    {
-      return false;
-    }
-  }
-  if(this->is_leaf())
-  {
-    // Actually process elements
-    assert((Ele.size() == 0 || Ele.cols() == 3) && "Elements should be triangles");
-    // Cheesecake way of hitting element
-    bool ret = ray_mesh_intersect(origin,dir,V,Ele.row(m_primitive),hits);
-    // Since we only gave ray_mesh_intersect a single face, it will have set
-    // any hits to id=0. Set these to this primitive's id
-    for(auto & hit : hits)
-    {
-      hit.id = m_primitive;
-    }
-    return ret;
-  }
-  std::vector<igl::Hit> left_hits;
-  std::vector<igl::Hit> right_hits;
-  const bool left_ret = m_left->intersect_ray(V,Ele,origin,dir,left_hits);
-  const bool right_ret = m_right->intersect_ray(V,Ele,origin,dir,right_hits);
-  hits.insert(hits.end(),left_hits.begin(),left_hits.end());
-  hits.insert(hits.end(),right_hits.begin(),right_hits.end());
-  return left_ret || right_ret;
+  RowVectorDIMS inv_dir = dir.cwiseInverse();
+  RowVectorDIMS inv_dir_pad = inv_dir;
+  igl::increment_ulp(inv_dir_pad, 2);
+  return intersect_ray_opt(V, Ele, origin, dir, inv_dir, inv_dir_pad, hits);
 }
 
 template <typename DerivedV, int DIM>
@@ -933,22 +908,76 @@ igl::AABB<DerivedV,DIM>::intersect_ray(
   const Scalar _min_t,
   igl::Hit & hit) const
 {
-  //// Naive, slow
-  //std::vector<igl::Hit> hits;
-  //intersect_ray(V,Ele,origin,dir,hits);
-  //if(hits.size() > 0)
-  //{
-  //  hit = hits.front();
-  //  return true;
-  //}else
-  //{
-  //  return false;
-  //}
+  RowVectorDIMS inv_dir = dir.cwiseInverse();
+  RowVectorDIMS inv_dir_pad = inv_dir;
+  igl::increment_ulp(inv_dir_pad, 2);
+  return intersect_ray_opt(V, Ele, origin, dir, inv_dir, inv_dir_pad, _min_t, hit);
+}
+
+template <typename DerivedV, int DIM>
+template <typename DerivedEle>
+IGL_INLINE bool
+igl::AABB<DerivedV,DIM>::intersect_ray_opt(
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedEle> & Ele,
+  const RowVectorDIMS & origin,
+  const RowVectorDIMS & dir,
+  const RowVectorDIMS & inv_dir,
+  const RowVectorDIMS & inv_dir_pad,
+  std::vector<igl::Hit> & hits) const
+{
+  hits.clear();
+  const Scalar t0 = 0;
+  const Scalar t1 = std::numeric_limits<Scalar>::infinity();
+  {
+    Scalar _1,_2;
+
+    if(!ray_box_intersect(origin,inv_dir,inv_dir_pad,m_box,t0,t1,_1,_2))
+    {
+      return false;
+    }
+  }
+  if(this->is_leaf())
+  {
+    // Actually process elements
+    assert((Ele.size() == 0 || Ele.cols() == 3) && "Elements should be triangles");
+    // Cheesecake way of hitting element
+    bool ret = ray_mesh_intersect(origin,dir,V,Ele.row(m_primitive),hits);
+    // Since we only gave ray_mesh_intersect a single face, it will have set
+    // any hits to id=0. Set these to this primitive's id
+    for(auto & hit: hits)
+    {
+      hit.id = m_primitive;
+    }
+    return ret;
+  }
+  std::vector<igl::Hit> left_hits;
+  std::vector<igl::Hit> right_hits;
+  const bool left_ret = m_left->intersect_ray_opt(V,Ele,origin,dir,inv_dir,inv_dir_pad,left_hits);
+  const bool right_ret = m_right->intersect_ray_opt(V,Ele,origin,dir,inv_dir,inv_dir_pad,right_hits);
+  hits.insert(hits.end(),left_hits.begin(),left_hits.end());
+  hits.insert(hits.end(),right_hits.begin(),right_hits.end());
+  return left_ret || right_ret;
+}
+
+template <typename DerivedV, int DIM>
+template <typename DerivedEle>
+IGL_INLINE bool
+igl::AABB<DerivedV,DIM>::intersect_ray_opt(
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedEle> & Ele,
+  const RowVectorDIMS & origin,
+  const RowVectorDIMS & dir,
+  const RowVectorDIMS & inv_dir,
+  const RowVectorDIMS & inv_dir_pad,
+  const Scalar _min_t,
+  igl::Hit & hit) const
+{
   Scalar min_t = _min_t;
   const Scalar t0 = 0;
   {
     Scalar _1,_2;
-    if(!ray_box_intersect(origin,dir,m_box,t0,min_t,_1,_2))
+    if(!ray_box_intersect(origin,inv_dir,inv_dir_pad,m_box,t0,min_t,_1,_2))
     {
       return false;
     }
@@ -967,7 +996,7 @@ igl::AABB<DerivedV,DIM>::intersect_ray(
   // differnce
   igl::Hit left_hit;
   igl::Hit right_hit;
-  bool left_ret = m_left->intersect_ray(V,Ele,origin,dir,min_t,left_hit);
+  bool left_ret = m_left->intersect_ray_opt(V,Ele,origin,dir,inv_dir,inv_dir_pad,min_t,left_hit);
   if(left_ret && left_hit.t<min_t)
   {
     // It's scary that this line doesn't seem to matter....
@@ -978,7 +1007,7 @@ igl::AABB<DerivedV,DIM>::intersect_ray(
   {
     left_ret = false;
   }
-  bool right_ret = m_right->intersect_ray(V,Ele,origin,dir,min_t,right_hit);
+  bool right_ret = m_right->intersect_ray_opt(V,Ele,origin,dir,inv_dir,inv_dir_pad,min_t,right_hit);
   if(right_ret && right_hit.t<min_t)
   {
     min_t = right_hit.t;

include/igl/AABB.h

@@ -443,6 +443,43 @@ private:
         Scalar & sqr_d,
         int & i,
         Eigen::PlainObjectBase<RowVectorDIMS> & c) const;
+
+      /// Intersect a ray with the mesh return all hits
+      ///
+      /// @param[in]  V  #V by dim list of vertex positions
+      /// @param[in]  Ele  #Ele by dim list of simplex indices
+      /// @param[in]  origin  dim-long ray origin
+      /// @param[in]  dir  dim-long ray direction
+      /// @param[out]  hits  list of hits
+      /// @return  true if any hits
+      template <typename DerivedEle>
+      IGL_INLINE bool intersect_ray_opt(
+        const Eigen::MatrixBase<DerivedV> & V,
+        const Eigen::MatrixBase<DerivedEle> & Ele,
+        const RowVectorDIMS & origin,
+        const RowVectorDIMS & dir,
+        const RowVectorDIMS & inv_dir,
+        const RowVectorDIMS & inv_dir_pad,
+        std::vector<igl::Hit> & hits) const;
+      /// Intersect a ray with the mesh return first hit farther than `min_t`
+      ///
+      /// @param[in]  V  #V by dim list of vertex positions
+      /// @param[in]  Ele  #Ele by dim list of simplex indices
+      /// @param[in]  origin  dim-long ray origin
+      /// @param[in]  dir  dim-long ray direction
+      /// @param[in]  min_t  minimum t value to consider
+      /// @param[out]  hit  first hit
+      /// @return  true if any hit
+      template <typename DerivedEle>
+      IGL_INLINE bool intersect_ray_opt(
+        const Eigen::MatrixBase<DerivedV> & V,
+        const Eigen::MatrixBase<DerivedEle> & Ele,
+        const RowVectorDIMS & origin,
+        const RowVectorDIMS & dir,
+        const RowVectorDIMS & inv_dir,
+        const RowVectorDIMS & inv_dir_pad,
+        const Scalar min_t,
+        igl::Hit & hit) const;
 public:
       EIGEN_MAKE_ALIGNED_OPERATOR_NEW
     };

include/igl/copyleft/cgal/assign_scalar.cpp

@@ -7,6 +7,8 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "assign_scalar.h"
 
+#include <cmath>
+
 template <>
 IGL_INLINE void igl::copyleft::cgal::assign_scalar(
   const CGAL::Epeck::FT & rhs,
@@ -75,7 +77,7 @@ IGL_INLINE void igl::copyleft::cgal::assign_scalar(
   const auto interval = CGAL::to_interval(cgal);
   d = interval.first;
   do {
-      const double next = nextafter(d, interval.second);
+      const double next = std::nextafter(d, interval.second);
       if (CGAL::abs(cgal-d) < CGAL::abs(cgal-next)) break;
       d = next;
   } while (d < interval.second);
@@ -90,7 +92,7 @@ IGL_INLINE void igl::copyleft::cgal::assign_scalar(
   const auto interval = CGAL::to_interval(cgal);
   d = interval.first;
   do {
-      const float next = nextafter(d, float(interval.second));
+      const float next = std::nextafter(d, float(interval.second));
       if (CGAL::abs(cgal-d) < CGAL::abs(cgal-next)) break;
       d = next;
   } while (d < float(interval.second));
@@ -131,7 +133,7 @@ IGL_INLINE void igl::copyleft::cgal::assign_scalar(
   const auto interval = CGAL::to_interval(cgal);
   d = interval.first;
   do {
-      const double next = nextafter(d, interval.second);
+      const double next = std::nextafter(d, interval.second);
       if (CGAL::abs(cgal-d) < CGAL::abs(cgal-next)) break;
       d = next;
   } while (d < interval.second);
@@ -144,7 +146,7 @@ IGL_INLINE void igl::copyleft::cgal::assign_scalar(
   const auto interval = CGAL::to_interval(cgal);
   d = interval.first;
   do {
-      const float next = nextafter(d, float(interval.second));
+      const float next = std::nextafter(d, float(interval.second));
       if (CGAL::abs(cgal-d) < CGAL::abs(cgal-next)) break;
       d = next;
   } while (d < float(interval.second));
@@ -166,7 +168,7 @@ IGL_INLINE void igl::copyleft::cgal::assign_scalar(
   const auto interval = CGAL::to_interval(cgal);
   d = interval.first;
   do {
-      const double next = nextafter(d, interval.second);
+      const double next = std::nextafter(d, interval.second);
       if (CGAL::abs(cgal-d) < CGAL::abs(cgal-next)) break;
       d = next;
   } while (d < interval.second);
@@ -179,7 +181,7 @@ IGL_INLINE void igl::copyleft::cgal::assign_scalar(
   const auto interval = CGAL::to_interval(cgal);
   d = interval.first;
   do {
-      const float next = nextafter(d, float(interval.second));
+      const float next = std::nextafter(d, float(interval.second));
       if (CGAL::abs(cgal-d) < CGAL::abs(cgal-next)) break;
       d = next;
   } while (d < float(interval.second));

include/igl/increment_ulp.cpp

@@ -0,0 +1,39 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2013 Alec Jacobson <[email protected]>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "increment_ulp.h"
+#include <cmath>
+#include <limits>
+
+
+// see "Robust BVH Ray Traversal" by Thiago Ize, section 3:
+// for why we need this
+template <typename Derived>
+IGL_INLINE void igl::increment_ulp(
+    Eigen::MatrixBase<Derived>& inout,
+    int it
+    )
+{
+    typedef typename Derived::Scalar Scalar;
+
+    inout = inout.unaryExpr([&it](Scalar v){
+              for (int k = 0; k < it; ++k) {
+                v = std::nextafter(v, std::signbit(v) ? -std::numeric_limits<Scalar>::infinity(): std::numeric_limits<Scalar>::infinity());
+              }
+              return v;
+            });
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::increment_ulp<Eigen::Matrix<float, 1, -1, 1, 1, -1>>(Eigen::MatrixBase<Eigen::Matrix<float, 1, -1, 1, 1, -1>>&, int);
+// generated by autoexplicit.sh
+template void igl::increment_ulp<Eigen::Matrix<float, 1, 3, 1, 1, 3>>(Eigen::MatrixBase<Eigen::Matrix<float, 1, 3, 1, 1, 3>>&, int);
+// generated by autoexplicit.sh
+template void igl::increment_ulp<Eigen::Matrix<double, 1, 3, 1, 1, 3>>(Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3>>&, int);
+#endif

include/igl/increment_ulp.h

@@ -0,0 +1,30 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Alec Jacobson <[email protected]>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_INCREMENT_ULP_H
+#define IGL_INCREMENT_ULP_H
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+/// Increment Unit in Last Place of a matrix
+///
+/// @param[in,out]  inout  input matrix
+/// @param[in]      it  number of increments
+template <typename Derived>
+IGL_INLINE void increment_ulp(
+    Eigen::MatrixBase<Derived>& inout,
+    int it);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "increment_ulp.cpp"
+#endif
+
+#endif

include/igl/ray_box_intersect.cpp

@@ -8,6 +8,7 @@
 #include "ray_box_intersect.h"
 #include <array>
 #include <igl/matlab_format.h>
+#include <igl/increment_ulp.h>
 
 template <
   typename Derivedsource,
@@ -15,102 +16,26 @@ template <
   typename Scalar>
 IGL_INLINE bool igl::ray_box_intersect(
   const Eigen::MatrixBase<Derivedsource> & origin,
-  const Eigen::MatrixBase<Deriveddir> & dir,
+  const Eigen::MatrixBase<Deriveddir> & inv_dir,
+  const Eigen::MatrixBase<Deriveddir> & inv_dir_pad,
   const Eigen::AlignedBox<Scalar,3> & box,
   const Scalar & t0,
   const Scalar & t1,
   Scalar & tmin,
   Scalar & tmax)
 {
-#ifdef false
-  // https://github.com/RMonica/basic_next_best_view/blob/master/src/RayTracer.cpp
-  const auto & intersectRayBox = [](
-    const Eigen::Vector3f& rayo,
-    const Eigen::Vector3f& rayd,
-    const Eigen::Vector3f& bmin,
-    const Eigen::Vector3f& bmax,
-    float & tnear,
-    float & tfar
-    )->bool
-  {
-    Eigen::Vector3f bnear;
-    Eigen::Vector3f bfar;
-    // Checks for intersection testing on each direction coordinate
-    // Computes
-    float t1, t2;
-    tnear = -1e+6f, tfar = 1e+6f; //, tCube;
-    bool intersectFlag = true;
-    for (int i = 0; i < 3; ++i) {
-  //    std::cout << "coordinate " << i << ": bmin " << bmin(i) << ", bmax " << bmax(i) << std::endl;
-      assert(bmin(i) <= bmax(i));
-      if (::fabs(rayd(i)) < 1e-6) {   // Ray parallel to axis i-th
-        if (rayo(i) < bmin(i) || rayo(i) > bmax(i)) {
-          intersectFlag = false;
-        }
-      }
-      else {
-        // Finds the nearest and the farthest vertices of the box from the ray origin
-        if (::fabs(bmin(i) - rayo(i)) < ::fabs(bmax(i) - rayo(i))) {
-          bnear(i) = bmin(i);
-          bfar(i) = bmax(i);
-        }
-        else {
-          bnear(i) = bmax(i);
-          bfar(i) = bmin(i);
-        }
-  //      std::cout << "  bnear " << bnear(i) << ", bfar " << bfar(i) << std::endl;
-        // Finds the distance parameters t1 and t2 of the two ray-box intersections:
-        // t1 must be the closest to the ray origin rayo.
-        t1 = (bnear(i) - rayo(i)) / rayd(i);
-        t2 = (bfar(i) - rayo(i)) / rayd(i);
-        if (t1 > t2) {
-          std::swap(t1,t2);
-        }
-        // The two intersection values are used to saturate tnear and tfar
-        if (t1 > tnear) {
-          tnear = t1;
-        }
-        if (t2 < tfar) {
-          tfar = t2;
-        }
-  //      std::cout << "  t1 " << t1 << ", t2 " << t2 << ", tnear " << tnear << ", tfar " << tfar
-  //        << "  tnear > tfar? " << (tnear > tfar) << ", tfar < 0? " << (tfar < 0) << std::endl;
-        if(tnear > tfar) {
-          intersectFlag = false;
-        }
-        if(tfar < 0) {
-        intersectFlag = false;
-        }
-      }
-    }
-    // Checks whether intersection occurs or not
-    return intersectFlag;
-  };
-  float tmin_f, tmax_f;
-  bool ret = intersectRayBox(
-      origin.   template cast<float>(),
-      dir.      template cast<float>(),
-      box.min().template cast<float>(),
-      box.max().template cast<float>(),
-      tmin_f,
-      tmax_f);
-  tmin = tmin_f;
-  tmax = tmax_f;
-  return ret;
-#else
   using namespace Eigen;
-  // This should be precomputed and provided as input
   typedef Matrix<Scalar,1,3>  RowVector3S;
-  const RowVector3S inv_dir( 1./dir(0),1./dir(1),1./dir(2));
   const std::array<bool, 3> sign = { inv_dir(0)<0, inv_dir(1)<0, inv_dir(2)<0};
   // http://people.csail.mit.edu/amy/papers/box-jgt.pdf
   // "An Efficient and Robust Ray–Box Intersection Algorithm"
+  // corrected in "Robust BVH Ray Traversal" by Thiago Ize, section 3:
   Scalar tymin, tymax, tzmin, tzmax;
   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);
+  tmax = ( bounds[1-sign[0]](0) - origin(0)) * inv_dir_pad(0);
   tymin = (bounds[sign[1]](1)   - origin(1)) * inv_dir(1);
-  tymax = (bounds[1-sign[1]](1) - origin(1)) * inv_dir(1);
+  tymax = (bounds[1-sign[1]](1) - origin(1)) * inv_dir_pad(1);
   // NaN-safe min and max
   const auto berger_perrin_min = [&](
       const Scalar a,
@@ -124,32 +49,51 @@ IGL_INLINE bool igl::ray_box_intersect(
   {
     return (a > b) ? a : b;
   };
-
   if ( (tmin > tymax) || (tymin > tmax) )
   {
     return false;
   }
-  tmin = berger_perrin_max(tmin,tymin);
-  tmax = berger_perrin_min(tmax,tymax);
+  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;
-  }
-  tmin = berger_perrin_max(tmin,tzmin);
-  tmax = berger_perrin_min(tmax,tzmax);
-  if(!( (tmin < t1) && (tmax > t0) ))
+  tzmax = (bounds[1-sign[2]](2) - origin(2)) * inv_dir_pad(2);
+  if ((tmin > tzmax) || (tzmin > tmax))
   {
     return false;
   }
-  return true;
-#endif
+  tmin = berger_perrin_max(tmin, tzmin);
+  tmax = berger_perrin_min(tmax, tzmax);
+  return ((tmin < t1) && (tmax > t0));
+}
+
+template <
+  typename Derivedsource,
+  typename Deriveddir,
+  typename Scalar>
+IGL_INLINE bool igl::ray_box_intersect(
+  const Eigen::MatrixBase<Derivedsource> & origin,
+  const Eigen::MatrixBase<Deriveddir> & dir,
+  const Eigen::AlignedBox<Scalar,3> & box,
+  const Scalar & t0,
+  const Scalar & t1,
+  Scalar & tmin,
+  Scalar & tmax)
+{
+    // precompute the inv_dir
+    Eigen::Matrix<Scalar, 1, 3> inv_dir = dir.cwiseInverse();
+    // see "Robust BVH Ray Traversal" by Thiago Ize, section 3:
+    // for why we need this
+    Eigen::Matrix<Scalar, 1, 3> inv_dir_pad = inv_dir;
+    igl::increment_ulp(inv_dir_pad, 2);
+    return igl::ray_box_intersect(origin, inv_dir, inv_dir_pad, box, t0, t1, tmin, tmax);
 }
 
 #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::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::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::AlignedBox<double, 3> const&, double const&, double const&, double&, double&);
+template bool igl::ray_box_intersect<Eigen::Matrix<float, 3, 1, 0, 3, 1>, Eigen::Matrix<float, 3, 1, 0, 3, 1>, float>(Eigen::MatrixBase<Eigen::Matrix<float, 3, 1, 0, 3, 1>> const&, Eigen::MatrixBase<Eigen::Matrix<float, 3, 1, 0, 3, 1>> const&, Eigen::AlignedBox<float, 3> const&, float const&, float const&, float&, float&);
 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_box_intersect.h

@@ -35,6 +35,21 @@ namespace igl
     const Scalar & t1,
     Scalar & tmin,
     Scalar & tmax);
+  /// \overload
+  /// \brief same with direction inverse precomputed
+  template <
+    typename Derivedsource,
+    typename Deriveddir,
+    typename Scalar>
+  IGL_INLINE bool ray_box_intersect(
+    const Eigen::MatrixBase<Derivedsource> & source,
+    const Eigen::MatrixBase<Deriveddir> & inv_dir,
+    const Eigen::MatrixBase<Deriveddir> & inv_dir_pad,
+    const Eigen::AlignedBox<Scalar,3> & box,
+    const Scalar & t0,
+    const Scalar & t1,
+    Scalar & tmin,
+    Scalar & tmax);
 }
 #ifndef IGL_STATIC_LIBRARY
 #  include "ray_box_intersect.cpp"

tests/include/igl/increment_ulp.cpp

@@ -0,0 +1,15 @@
+#include <test_common.h>
+#include <igl/increment_ulp.h>
+#include <limits>
+
+TEST_CASE("increment_ulp: dense", "[igl]" )
+{
+  Eigen::RowVectorXf v(5), ref(5);
+  v << std::numeric_limits<float>::lowest(), 0.f, std::numeric_limits<float>::max(), std::numeric_limits<float>::infinity(), std::numeric_limits<float>::quiet_NaN();
+  ref << -std::numeric_limits<float>::infinity(), std::numeric_limits<float>::denorm_min(), std::numeric_limits<float>::infinity(), std::numeric_limits<float>::infinity(), std::numeric_limits<float>::quiet_NaN();
+
+  igl::increment_ulp(v, 1);
+
+  REQUIRE(v.head(4) == ref.head(4));
+  REQUIRE((std::isnan(v(4)) && std::isnan(ref(4))));
+}

tests/include/igl/ray_mesh_intersect.cpp

@@ -1,5 +1,6 @@
 #include <test_common.h>
 #include <igl/ray_mesh_intersect.h>
+#include <igl/ray_box_intersect.h>
 #include <igl/AABB.h>
 
 
@@ -29,7 +30,6 @@ TEST_CASE("ray_mesh_intersect: one_triangle", "[igl]")
   IGL_POP_FPE;
 }
 
-
 TEST_CASE("ray_mesh_intersect: corner-case", "[igl]")
 {
   IGL_PUSH_FPE;
@@ -79,3 +79,36 @@ TEST_CASE("ray_mesh_intersect: corner-case", "[igl]")
   }
   IGL_POP_FPE;
 }
+
+TEST_CASE("ray_mesh_intersect: corner-case2", "[igl]")
+{
+  IGL_PUSH_FPE;
+
+  Eigen::MatrixXf vertices(3, 3);
+  vertices <<
+      -2.891303300857544, 0.7025225162506104, 1.157850384712219,
+      -2.870383024215698, 0.7444183230400085, 1.18663215637207,
+      -2.890183448791504, 0.7462523579597473, 1.157822966575623;
+
+
+  Eigen::MatrixXi faces(1, 3);
+  faces <<
+      0, 2, 1;
+
+  Eigen::Vector3f origin;
+  Eigen::Vector3f direction;
+
+  origin << -5.411622047424316, -0.02165498770773411, 0.7916983366012573;
+  direction << 0.9475222229957581, 0.2885690927505493, 0.1375846415758133;
+
+  std::vector<igl::Hit> hits, hits_bvh;
+  bool is_hit = igl::ray_mesh_intersect(origin, direction, vertices, faces, hits);
+  Eigen::AlignedBox3f box;
+  box.extend(vertices.row(0).transpose());
+  box.extend(vertices.row(1).transpose());
+  box.extend(vertices.row(2).transpose());
+
+  float tmin, tmax;
+  bool is_hit_box = igl::ray_box_intersect(origin, direction, box, 0.0f, std::numeric_limits<float>::max(), tmin, tmax);
+  REQUIRE (is_hit == is_hit_box);
+}