Reduce dynamic allocations for AABB queries (#2001) [ci skip]

* Use static num of cols if available.

* Don't allocate a new vector for only one triangle.

* bug in test

---------

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

include/igl/point_simplex_squared_distance.cpp

@@ -149,8 +149,8 @@ IGL_INLINE void igl::point_simplex_squared_distance(
   Derivedsqr_d & sqr_d,
   Eigen::MatrixBase<Derivedc> & c)
 {
-  // Use Dynamic because we don't know Ele.cols() at compile time.
-  Eigen::Matrix<typename Derivedc::Scalar,1,Eigen::Dynamic> b;
+  // DerivedEle::ColsAtCompileTime will fallback to Eigen::Dynamic if unknown.
+  Eigen::Matrix<typename Derivedc::Scalar,1,DerivedEle::ColsAtCompileTime> b;
   point_simplex_squared_distance<DIM>( p, V, Ele, primitive, sqr_d, c, b );
 }
 

include/igl/ray_mesh_intersect.cpp

@@ -12,40 +12,68 @@ extern "C"
 #include "raytri.c"
 }
 
+namespace igl {
 template <
   typename Derivedsource,
   typename Deriveddir,
   typename DerivedV,
   typename DerivedF>
-IGL_INLINE bool igl::ray_mesh_intersect(
+IGL_INLINE bool ray_triangle_intersect(
   const Eigen::MatrixBase<Derivedsource> & s,
   const Eigen::MatrixBase<Deriveddir> & dir,
   const Eigen::MatrixBase<DerivedV> & V,
   const Eigen::MatrixBase<DerivedF> & F,
-  std::vector<igl::Hit> & hits)
+  const int f,
+  igl::Hit& hit)
 {
   using namespace Eigen;
-  using namespace std;
-  // Should be but can't be const
+
+  // intersect_triangle1 needs non-const inputs.
   Vector3d s_d = s.template cast<double>();
   Vector3d dir_d = dir.template cast<double>();
+  RowVector3d v0 = V.row(F(f,0)).template cast<double>();
+  RowVector3d v1 = V.row(F(f,1)).template cast<double>();
+  RowVector3d v2 = V.row(F(f,2)).template cast<double>();
+
+  // shoot ray, record hit
+  double t,u,v;
+  if(intersect_triangle1(
+    s_d.data(), dir_d.data(), v0.data(), v1.data(), v2.data(), &t, &u, &v) &&
+    t>0)
+  {
+    hit = {static_cast<int>(f),static_cast<int>(-1),
+           static_cast<float>(u),static_cast<float>(v),
+           static_cast<float>(t)};
+    return true;
+  }
+
+  return false;
+}
+
+template <
+  typename Derivedsource,
+  typename Deriveddir,
+  typename DerivedV,
+  typename DerivedF>
+IGL_INLINE bool ray_mesh_intersect(
+  const Eigen::MatrixBase<Derivedsource> & s,
+  const Eigen::MatrixBase<Deriveddir> & dir,
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
+  std::vector<igl::Hit> & hits)
+{
+  using namespace Eigen;
+
   hits.clear();
   hits.reserve(F.rows());
 
   // loop over all triangles
   for(int f = 0;f<F.rows();f++)
   {
-    // Should be but can't be const
-    RowVector3d v0 = V.row(F(f,0)).template cast<double>();
-    RowVector3d v1 = V.row(F(f,1)).template cast<double>();
-    RowVector3d v2 = V.row(F(f,2)).template cast<double>();
-    // shoot ray, record hit
-    double t,u,v;
-    if(intersect_triangle1(
-      s_d.data(), dir_d.data(), v0.data(), v1.data(), v2.data(), &t, &u, &v) &&
-      t>0)
+    igl::Hit hit;
+    if(ray_triangle_intersect(s, dir, V, F, f, hit))
     {
-      hits.push_back({(int)f,(int)-1,(float)u,(float)v,(float)t});
+      hits.push_back(hit);
     }
   }
   // Sort hits based on distance
@@ -61,24 +89,32 @@ template <
   typename Deriveddir,
   typename DerivedV,
   typename DerivedF>
-IGL_INLINE bool igl::ray_mesh_intersect(
-  const Eigen::MatrixBase<Derivedsource> & source,
+IGL_INLINE bool ray_mesh_intersect(
+  const Eigen::MatrixBase<Derivedsource> & s,
   const Eigen::MatrixBase<Deriveddir> & dir,
   const Eigen::MatrixBase<DerivedV> & V,
   const Eigen::MatrixBase<DerivedF> & F,
   igl::Hit & hit)
 {
-  std::vector<igl::Hit> hits;
-  ray_mesh_intersect(source,dir,V,F,hits);
-  if(hits.size() > 0)
+  if (F.rows() == 1)
   {
-    hit = hits.front();
-    return true;
+    // Shortcut for AABB based queries.
+    return ray_triangle_intersect(s, dir, V, F, 0, hit);
   }else
   {
-    return false;
+    std::vector<igl::Hit> hits;
+    ray_mesh_intersect(s,dir,V,F,hits);
+    if(hits.size() > 0)
+    {
+      hit = hits.front();
+      return true;
+    }else
+    {
+      return false;
+    }
   }
 }
+}
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation

include/igl/ray_mesh_intersect.h

@@ -48,6 +48,24 @@ namespace igl
     const Eigen::MatrixBase<DerivedV> & V,
     const Eigen::MatrixBase<DerivedF> & F,
     igl::Hit & hit);
+
+  // Explicit function to check one triangle (given by index f) of the mesh.
+  // This function is used by ray_mesh_intersect to check each triangle.
+  // Outputs:
+  //   hit  hit with the given triangle, set only if it exists
+  // Returns true if there was a hit
+  template <
+    typename Derivedsource,
+    typename Deriveddir,
+    typename DerivedV,
+    typename DerivedF>
+  IGL_INLINE bool ray_triangle_intersect(
+    const Eigen::MatrixBase<Derivedsource> & source,
+    const Eigen::MatrixBase<Deriveddir> & dir,
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const int f,
+    igl::Hit& hit);
 }
 #ifndef IGL_STATIC_LIBRARY
 #  include "ray_mesh_intersect.cpp"

tests/include/igl/ray_mesh_intersect.cpp

@@ -0,0 +1,25 @@
+#include <test_common.h>
+#include <igl/ray_mesh_intersect.h>
+
+TEST_CASE("ray_mesh_intersect: one_triangle", "[igl]")
+{
+  Eigen::MatrixXd V(3,3);
+  V.row(0) << 0.0, 0.0, 0.0;
+  V.row(1) << 1.0, 0.0, 0.0;
+  V.row(2) << 0.5, 1.0, 0.0;
+  
+  Eigen::MatrixXi F(1,3);
+  F.row(0) << 0,1,2;
+  
+  Eigen::Vector3f source{0.5, 0.5, -1.0};
+  Eigen::Vector3f direction{0.0, 0.0, 1.0};
+  
+  igl::Hit hit;
+  REQUIRE(igl::ray_mesh_intersect(source, direction, V, F, hit) == true);
+  REQUIRE(hit.t == Approx(1.0));
+  
+  std::vector<igl::Hit> hits;
+  REQUIRE(igl::ray_mesh_intersect(source, direction, V, F, hits) == true);
+  REQUIRE(hits.size() == 1);
+  REQUIRE(hits.front().t == Approx(1.0));
+}