expose ambient occlusion without embree

Former-commit-id: 2a42010680d7ad7b4f7a5860c81b256ae5d8f5e9

include/igl/AABB.h

@@ -8,6 +8,7 @@
 #ifndef IGL_AABB_H
 #define IGL_AABB_H
 
+#include "Hit.h"
 #include <Eigen/Core>
 #include <Eigen/Geometry>
 #include <vector>
@@ -134,7 +135,8 @@ public:
           const Eigen::VectorXi & I);
       // Return whether at leaf node
       inline bool is_leaf() const;
-      // Find the indices of elements containing given point.
+      // Find the indices of elements containing given point: this makes sense
+      // when Ele is a co-dimension 0 simplex (tets in 3D, triangles in 2D).
       //
       // Inputs:
       //   V  #V by dim list of mesh vertex positions. **Should be same as used to
@@ -185,19 +187,43 @@ public:
       // Known bugs: currently assumes Elements are triangles regardless of
       // dimension.
       inline Scalar squared_distance(
-          const Eigen::PlainObjectBase<DerivedV> & V,
-          const Eigen::MatrixXi & Ele, 
-          const RowVectorDIMS & p,
-          int & i,
-          RowVectorDIMS & c) const;
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::MatrixXi & Ele, 
+        const RowVectorDIMS & p,
+        int & i,
+        RowVectorDIMS & c) const;
 //private:
       inline Scalar squared_distance(
-          const Eigen::PlainObjectBase<DerivedV> & V,
-          const Eigen::MatrixXi & Ele, 
-          const RowVectorDIMS & p,
-          const Scalar min_sqr_d,
-          int & i,
-          RowVectorDIMS & c) const;
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::MatrixXi & Ele, 
+        const RowVectorDIMS & p,
+        const Scalar min_sqr_d,
+        int & i,
+        RowVectorDIMS & c) const;
+      // All hits
+      inline bool intersect_ray(
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::MatrixXi & Ele, 
+        const RowVectorDIMS & origin,
+        const RowVectorDIMS & dir,
+        std::vector<igl::Hit> & hits) const;
+      // First hit
+      inline bool intersect_ray(
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::MatrixXi & Ele, 
+        const RowVectorDIMS & origin,
+        const RowVectorDIMS & dir,
+        igl::Hit & hit) const;
+//private:
+      inline bool intersect_ray(
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::MatrixXi & Ele, 
+        const RowVectorDIMS & origin,
+        const RowVectorDIMS & dir,
+        const Scalar min_t,
+        igl::Hit & hit) const;
+
+
 public:
       template <
         typename DerivedP, 
@@ -281,6 +307,8 @@ public:
 #include "project_to_line_segment.h"
 #include "sort.h"
 #include "volume.h"
+#include "ray_box_intersect.h"
+#include "ray_mesh_intersect.h"
 #include <iostream>
 #include <iomanip>
 #include <limits>
@@ -1149,4 +1177,108 @@ igl::AABB<DerivedV,DIM>::barycentric_coordinates(
   bary(0) = 1.0f - bary(1) - bary(2);
 }
 
+template <typename DerivedV, int DIM>
+inline bool 
+igl::AABB<DerivedV,DIM>::intersect_ray(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::MatrixXi & Ele, 
+  const RowVectorDIMS & origin,
+  const RowVectorDIMS & dir,
+  std::vector<igl::Hit> & hits) const
+{
+  hits.clear();
+  const Scalar t0 = 0;
+  const Scalar t1 = std::numeric_limits<Scalar>::infinity();
+  if(!ray_box_intersect(origin,dir,m_box,t0,t1))
+  {
+    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
+    return ray_mesh_intersect(origin,dir,V,Ele.row(m_primitive).eval(),hits);
+  }
+  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;
+}
+
+template <typename DerivedV, int DIM>
+inline bool 
+igl::AABB<DerivedV,DIM>::intersect_ray(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::MatrixXi & Ele, 
+  const RowVectorDIMS & origin,
+  const RowVectorDIMS & dir,
+  igl::Hit & hit) const
+{
+  return intersect_ray(
+    V,Ele,origin,dir,std::numeric_limits<Scalar>::infinity(),hit);
+}
+
+template <typename DerivedV, int DIM>
+inline bool 
+igl::AABB<DerivedV,DIM>::intersect_ray(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::MatrixXi & Ele, 
+  const RowVectorDIMS & origin,
+  const RowVectorDIMS & dir,
+  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;
+  //}
+  Scalar min_t = _min_t;
+  const Scalar t0 = 0;
+  if(!ray_box_intersect(origin,dir,m_box,t0,min_t))
+  {
+    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
+    return ray_mesh_intersect(origin,dir,V,Ele.row(m_primitive).eval(),hit);
+  }
+  igl::Hit left_hit;
+  igl::Hit right_hit;
+  bool left_ret = m_left->intersect_ray(V,Ele,origin,dir,min_t,left_hit);
+  if(left_ret && left_hit.t<min_t)
+  {
+    min_t = left_hit.t;
+    hit = left_hit;
+    left_ret = true;
+  }else
+  {
+    left_ret = false;
+  }
+  bool right_ret = m_right->intersect_ray(V,Ele,origin,dir,min_t,right_hit);
+  if(right_ret && right_hit.t<min_t)
+  {
+    min_t = right_hit.t;
+    hit = right_hit;
+    right_ret = true;
+  }else
+  {
+    right_ret = false;
+  }
+  return left_ret || right_ret;
+}
+
 #endif

include/igl/ambient_occlusion.cpp

@@ -0,0 +1,124 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// 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 "ambient_occlusion.h"
+#include "random_dir.h"
+#include "ray_mesh_intersect.h"
+#include "EPS.h"
+#include "Hit.h"
+
+template <
+  typename DerivedP,
+  typename DerivedN,
+  typename DerivedS >
+IGL_INLINE void igl::ambient_occlusion(
+  const std::function<
+    bool(
+      const Eigen::Vector3f&,
+      const Eigen::Vector3f&)
+      > & shoot_ray,
+  const Eigen::PlainObjectBase<DerivedP> & P,
+  const Eigen::PlainObjectBase<DerivedN> & N,
+  const int num_samples,
+  Eigen::PlainObjectBase<DerivedS> & S)
+{
+  using namespace Eigen;
+  using namespace igl;
+  const int n = P.rows();
+  // Resize output
+  S.resize(n,1);
+  // Embree seems to be parallel when constructing but not when tracing rays
+#pragma omp parallel for
+  // loop over mesh vertices
+  for(int p = 0;p<n;p++)
+  {
+    const Vector3f origin = P.row(p).template cast<float>();
+    const Vector3f normal = N.row(p).template cast<float>();
+    int num_hits = 0;
+    MatrixXf D = random_dir_stratified(num_samples).cast<float>();
+    for(int s = 0;s<num_samples;s++)
+    {
+      //Vector3d d = random_dir();
+      Vector3f d = D.row(s);
+      if(d.dot(normal) < 0)
+      {
+        // reverse ray
+        d *= -1;
+      }
+      if(shoot_ray(origin,d))
+      {
+        num_hits++;
+      }
+    }
+    S(p) = (double)num_hits/(double)num_samples;
+  }
+}
+
+template <
+  typename DerivedV,
+  int DIM,
+  typename DerivedF,
+  typename DerivedP,
+  typename DerivedN,
+  typename DerivedS >
+IGL_INLINE void igl::ambient_occlusion(
+  const igl::AABB<DerivedV,DIM> & aabb,
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  const Eigen::PlainObjectBase<DerivedP> & P,
+  const Eigen::PlainObjectBase<DerivedN> & N,
+  const int num_samples,
+  Eigen::PlainObjectBase<DerivedS> & S)
+{
+  const auto & shoot_ray = [&aabb,&V,&F](
+    const Eigen::Vector3f& _s,
+    const Eigen::Vector3f& dir)->bool
+  {
+    Eigen::Vector3f s = _s+1e-4*dir;
+    igl::Hit hit;
+    return aabb.intersect_ray(
+      V,
+      F,
+      s  .cast<typename DerivedV::Scalar>().eval(),
+      dir.cast<typename DerivedV::Scalar>().eval(),
+      hit);
+  };
+  return ambient_occlusion(shoot_ray,P,N,num_samples,S);
+
+}
+
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedP,
+  typename DerivedN,
+  typename DerivedS >
+IGL_INLINE void igl::ambient_occlusion(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  const Eigen::PlainObjectBase<DerivedP> & P,
+  const Eigen::PlainObjectBase<DerivedN> & N,
+  const int num_samples,
+  Eigen::PlainObjectBase<DerivedS> & S)
+{
+  if(F.rows() < 100)
+  {
+    // Super naive
+    const auto & shoot_ray = [&V,&F](
+      const Eigen::Vector3f& _s,
+      const Eigen::Vector3f& dir)->bool
+    {
+      Eigen::Vector3f s = _s+1e-4*dir;
+      igl::Hit hit;
+      return ray_mesh_intersect(s,dir,V,F,hit);
+    };
+    return ambient_occlusion(shoot_ray,P,N,num_samples,S);
+  }
+  AABB<DerivedV,3> aabb;
+  aabb.init(V,F);
+  return ambient_occlusion(aabb,V,F,P,N,num_samples,S);
+}

include/igl/ambient_occlusion.h

@@ -0,0 +1,80 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// 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_AMBIENT_OCCLUSION_H
+#define IGL_AMBIENT_OCCLUSION_H
+#include "igl_inline.h"
+#include "AABB.h"
+#include <Eigen/Core>
+#include <functional>
+namespace igl
+{
+  // Compute ambient occlusion per given point
+  //
+  // Inputs:
+  //    shoot_ray  function handle that outputs hits of a given ray against a
+  //      mesh (embedded in function handles as captured variable/data)
+  //    P  #P by 3 list of origin points
+  //    N  #P by 3 list of origin normals
+  // Outputs:
+  //    S  #P list of ambient occlusion values between 1 (fully occluded) and
+  //      0 (not occluded)
+  //
+  template <
+    typename DerivedP,
+    typename DerivedN,
+    typename DerivedS >
+  IGL_INLINE void ambient_occlusion(
+    const std::function<
+      bool(
+        const Eigen::Vector3f&,
+        const Eigen::Vector3f&)
+        > & shoot_ray,
+    const Eigen::PlainObjectBase<DerivedP> & P,
+    const Eigen::PlainObjectBase<DerivedN> & N,
+    const int num_samples,
+    Eigen::PlainObjectBase<DerivedS> & S);
+  // Inputs:
+  //   AABB  axis-aligned bounding box hierarchy around (V,F)
+  template <
+    typename DerivedV,
+    int DIM,
+    typename DerivedF,
+    typename DerivedP,
+    typename DerivedN,
+    typename DerivedS >
+  IGL_INLINE void ambient_occlusion(
+    const igl::AABB<DerivedV,DIM> & aabb,
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    const Eigen::PlainObjectBase<DerivedP> & P,
+    const Eigen::PlainObjectBase<DerivedN> & N,
+    const int num_samples,
+    Eigen::PlainObjectBase<DerivedS> & S);
+  // Inputs:
+  //    V  #V by 3 list of mesh vertex positions
+  //    F  #F by 3 list of mesh face indices into V
+  template <
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedP,
+    typename DerivedN,
+    typename DerivedS >
+  IGL_INLINE void ambient_occlusion(
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    const Eigen::PlainObjectBase<DerivedP> & P,
+    const Eigen::PlainObjectBase<DerivedN> & N,
+    const int num_samples,
+    Eigen::PlainObjectBase<DerivedS> & S);
+
+};
+#ifndef IGL_STATIC_LIBRARY
+#  include "ambient_occlusion.cpp"
+#endif
+
+#endif

include/igl/embree/ambient_occlusion.cpp

@@ -6,10 +6,9 @@
 // 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 "ambient_occlusion.h"
+#include "../ambient_occlusion.h"
 #include "EmbreeIntersector.h"
 #include "../Hit.h"
-#include "../random_dir.h"
-#include "../EPS.h"
 
 template <
   typename DerivedP,
@@ -22,38 +21,15 @@ IGL_INLINE void igl::embree::ambient_occlusion(
   const int num_samples,
   Eigen::PlainObjectBase<DerivedS> & S)
 {
-  using namespace Eigen;
-  using namespace igl;
-  const int n = P.rows();
-  // Resize output
-  S.resize(n,1);
-  // Embree seems to be parallel when constructing but not when tracing rays
-#pragma omp parallel for
-  // loop over mesh vertices
-  for(int p = 0;p<n;p++)
+  const auto & shoot_ray = [&ei](
+    const Eigen::Vector3f& s,
+    const Eigen::Vector3f& dir)->bool
   {
-    const Vector3f origin = P.row(p).template cast<float>();
-    const Vector3f normal = N.row(p).template cast<float>();
-    int num_hits = 0;
-    MatrixXf D = random_dir_stratified(num_samples).cast<float>();
-    for(int s = 0;s<num_samples;s++)
-    {
-      //Vector3d d = random_dir();
-      Vector3f d = D.row(s);
-      if(d.dot(normal) < 0)
-      {
-        // reverse ray
-        d *= -1;
-      }
-      igl::Hit hit;
-      const float tnear = 1e-4f;
-      if(ei.intersectRay(origin,d,hit,tnear))
-      {
-        num_hits++;
-      }
-    }
-    S(p) = (double)num_hits/(double)num_samples;
-  }
+    igl::Hit hit;
+    const float tnear = 1e-4f;
+    return ei.intersectRay(s,dir,hit,tnear);
+  };
+  return igl::ambient_occlusion(shoot_ray,P,N,num_samples,S);
 }
 
 template <

include/igl/ray_box_intersect.cpp

@@ -0,0 +1,53 @@
+#include "ray_box_intersect.h"
+template <
+  typename Derivedsource,
+  typename Deriveddir,
+  typename Scalar>
+IGL_INLINE bool igl::ray_box_intersect(
+  const Eigen::PlainObjectBase<Derivedsource> & origin,
+  const Eigen::PlainObjectBase<Deriveddir> & dir,
+  const Eigen::AlignedBox<Scalar,3> & box,
+  const Scalar & t0,
+  const Scalar & t1)
+{
+  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::vector<bool> 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"
+  Scalar tmin, tmax, tymin, tymax, tzmin, tzmax;
+  std::vector<RowVector3S> bounds = {box.min(),box.max()};
+  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) )
+  {
+    return false;
+  }
+  if (tymin > tmin)
+  {
+    tmin = tymin;
+  }
+  if (tymax < tmax)
+  {
+    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;
+  }
+  return ( (tmin < t1) && (tmax > t0) );
+}

include/igl/ray_box_intersect.h

@@ -0,0 +1,39 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// 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_RAY_BOX_INTERSECT_H
+#define IGL_RAY_BOX_INTERSECT_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <Eigen/Geometry>
+namespace igl
+{
+  // Determine whether a ray origin+t*dir and box intersect within the ray's parameterized
+  // range (t0,t1)
+  //
+  // Inputs:
+  //   source  3-vector origin of ray
+  //   dir  3-vector direction of ray
+  //   box  axis aligned box
+  //   t0  hit only if hit.t less than t0
+  //   t1  hit only if hit.t greater than t1
+  // Returns true if hit
+  template <
+    typename Derivedsource,
+    typename Deriveddir,
+    typename Scalar>
+  IGL_INLINE bool ray_box_intersect(
+    const Eigen::PlainObjectBase<Derivedsource> & source,
+    const Eigen::PlainObjectBase<Deriveddir> & dir,
+    const Eigen::AlignedBox<Scalar,3> & box,
+    const Scalar & t0,
+    const Scalar & t1);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "ray_box_intersect.cpp"
+#endif
+#endif

include/igl/ray_mesh_intersect.cpp

@@ -0,0 +1,73 @@
+#include "ray_mesh_intersect.h"
+
+extern "C"
+{
+#include "raytri.c"
+}
+
+template <
+  typename Derivedsource,
+  typename Deriveddir,
+  typename DerivedV, 
+  typename DerivedF> 
+IGL_INLINE bool igl::ray_mesh_intersect(
+  const Eigen::PlainObjectBase<Derivedsource> & s,
+  const Eigen::PlainObjectBase<Deriveddir> & dir,
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  std::vector<igl::Hit> & hits)
+{
+  using namespace Eigen;
+  using namespace igl;
+  using namespace std;
+  // Should be but can't be const 
+  Vector3d s_d = s.template cast<double>();
+  Vector3d dir_d = dir.template cast<double>();
+  hits.clear();
+  // 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)
+    {
+      hits.push_back({(int)f,(int)-1,(float)u,(float)v,(float)t});
+    }
+  }
+  // Sort hits based on distance
+  std::sort(
+    hits.begin(),
+    hits.end(),
+    [](const Hit & a, const Hit & b)->bool{ return a.t < b.t;});
+  return hits.size() > 0;
+}
+
+template <
+  typename Derivedsource,
+  typename Deriveddir,
+  typename DerivedV, 
+  typename DerivedF> 
+IGL_INLINE bool igl::ray_mesh_intersect(
+  const Eigen::PlainObjectBase<Derivedsource> & source,
+  const Eigen::PlainObjectBase<Deriveddir> & dir,
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  igl::Hit & hit)
+{
+  std::vector<igl::Hit> hits;
+  ray_mesh_intersect(source,dir,V,F,hits);
+  if(hits.size() > 0)
+  {
+    hit = hits.front();
+    return true;
+  }else
+  {
+    return false;
+  }
+}

include/igl/ray_mesh_intersect.h

@@ -0,0 +1,49 @@
+#ifndef IGL_RAY_MESH_INTERSECT_H
+#define IGL_RAY_MESH_INTERSECT_H
+#include "igl_inline.h"
+#include "Hit.h"
+#include <Eigen/Core>
+#include <vector>
+namespace igl
+{
+  // Shoot a ray against a mesh (V,F) and collect all hits.
+  //
+  // Inputs:
+  //   source  3-vector origin of ray
+  //   dir  3-vector direction of ray
+  //   V  #V by 3 list of mesh vertex positions
+  //   F  #F by 3 list of mesh face indices into V
+  // Outputs:
+  //    hits  **sorted** list of hits
+  // Returns true if there were any hits (hits.size() > 0)
+  //
+  template <
+    typename Derivedsource,
+    typename Deriveddir,
+    typename DerivedV, 
+    typename DerivedF> 
+  IGL_INLINE bool ray_mesh_intersect(
+    const Eigen::PlainObjectBase<Derivedsource> & source,
+    const Eigen::PlainObjectBase<Deriveddir> & dir,
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    std::vector<igl::Hit> & hits);
+  // Outputs:
+  //   hit  first hit, 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_mesh_intersect(
+    const Eigen::PlainObjectBase<Derivedsource> & source,
+    const Eigen::PlainObjectBase<Deriveddir> & dir,
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    igl::Hit & hit);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "ray_mesh_intersect.cpp"
+#endif
+#endif

include/igl/unproject_in_mesh.cpp

@@ -7,6 +7,7 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "unproject_in_mesh.h"
 #include "unproject_ray.h"
+#include "ray_mesh_intersect.h"
 
 template < typename Derivedobj>
   IGL_INLINE int igl::unproject_in_mesh(
@@ -72,30 +73,7 @@ template < typename DerivedV, typename DerivedF, typename Derivedobj>
     const Eigen::Vector3f& dir,
     std::vector<igl::Hit> & hits)
   {
-    // Should be but can't be const 
-    Vector3d s_d = s.template cast<double>();
-    Vector3d dir_d = dir.template cast<double>();
-    hits.clear();
-    // 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))
-      {
-        hits.push_back({(int)f,(int)-1,(float)u,(float)v,(float)t});
-      }
-    }
-    // Sort hits based on distance
-    std::sort(
-      hits.begin(),
-      hits.end(),
-      [](const Hit & a, const Hit & b)->bool{ return a.t < b.t;});
+    ray_mesh_intersect(s,dir,V,F,hits);
   };
   return unproject_in_mesh(pos,model,proj,viewport,shoot_ray,obj,hits);
 }