minor changes with no effect to AABB ray

Former-commit-id: 1b127a8c08b3db2321a95d4cdae637530ae2429e

include/igl/AABB.h

@@ -313,6 +313,8 @@ public:
 #include <iomanip>
 #include <limits>
 #include <list>
+#include <queue>
+#include <stack>
 
 template <typename DerivedV, int DIM>
   template <typename Derivedbb_mins, typename Derivedbb_maxs>
@@ -1189,9 +1191,12 @@ igl::AABB<DerivedV,DIM>::intersect_ray(
   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;
+    Scalar _1,_2;
+    if(!ray_box_intersect(origin,dir,m_box,t0,t1,_1,_2))
+    {
+      return false;
+    }
   }
   if(this->is_leaf())
   {
@@ -1218,8 +1223,50 @@ igl::AABB<DerivedV,DIM>::intersect_ray(
   const RowVectorDIMS & dir,
   igl::Hit & hit) const
 {
+#if false
+  // BFS
+  std::queue<const AABB *> Q;
+  // Or DFS
+  //std::stack<const AABB *> Q;
+  Q.push(this);
+  bool any_hit = false;
+  hit.t = std::numeric_limits<Scalar>::infinity();
+  while(!Q.empty())
+  {
+    const AABB * tree = Q.front();
+    //const AABB * tree = Q.top();
+    Q.pop();
+    {
+      Scalar _1,_2;
+      if(!ray_box_intersect(
+        origin,dir,tree->m_box,Scalar(0),Scalar(hit.t),_1,_2))
+      {
+        continue;
+      }
+    }
+    if(tree->is_leaf())
+    {
+      // Actually process elements
+      assert((Ele.size() == 0 || Ele.cols() == 3) && "Elements should be triangles");
+      igl::Hit leaf_hit;
+      if(
+        ray_mesh_intersect(origin,dir,V,Ele.row(tree->m_primitive).eval(),leaf_hit)&&
+        leaf_hit.t < hit.t)
+      {
+        hit = leaf_hit;
+      }
+      continue;
+    }
+    // Add children to queue
+    Q.push(tree->m_left);
+    Q.push(tree->m_right);
+  }
+  return any_hit;
+#else
+  // DFS
   return intersect_ray(
     V,Ele,origin,dir,std::numeric_limits<Scalar>::infinity(),hit);
+#endif
 }
 
 template <typename DerivedV, int DIM>
@@ -1245,9 +1292,12 @@ igl::AABB<DerivedV,DIM>::intersect_ray(
   //}
   Scalar min_t = _min_t;
   const Scalar t0 = 0;
-  if(!ray_box_intersect(origin,dir,m_box,t0,min_t))
   {
-    return false;
+    Scalar _1,_2;
+    if(!ray_box_intersect(origin,dir,m_box,t0,min_t,_1,_2))
+    {
+      return false;
+    }
   }
   if(this->is_leaf())
   {
@@ -1256,11 +1306,15 @@ igl::AABB<DerivedV,DIM>::intersect_ray(
     // Cheesecake way of hitting element
     return ray_mesh_intersect(origin,dir,V,Ele.row(m_primitive).eval(),hit);
   }
+
+  // Doesn't seem like smartly choosing left before/after right makes a
+  // differnce
   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)
   {
+    // It's scary that this line doesn't seem to matter....
     min_t = left_hit.t;
     hit = left_hit;
     left_ret = true;

include/igl/ambient_occlusion.cpp

@@ -31,18 +31,19 @@ IGL_INLINE void igl::ambient_occlusion(
   const int n = P.rows();
   // Resize output
   S.resize(n,1);
+  VectorXi hits = VectorXi::Zero(n,1);
   // Embree seems to be parallel when constructing but not when tracing rays
 #pragma omp parallel for
+  const MatrixXf D = random_dir_stratified(num_samples).cast<float>();
   // 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();
+//      //Vector3d d = random_dir();
       Vector3f d = D.row(s);
       if(d.dot(normal) < 0)
       {

include/igl/ray_box_intersect.cpp

@@ -8,8 +8,86 @@ IGL_INLINE bool igl::ray_box_intersect(
   const Eigen::PlainObjectBase<Deriveddir> & dir,
   const Eigen::AlignedBox<Scalar,3> & box,
   const Scalar & t0,
-  const Scalar & t1)
+  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;
@@ -17,7 +95,7 @@ IGL_INLINE bool igl::ray_box_intersect(
   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;
+  Scalar 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);
@@ -49,5 +127,10 @@ IGL_INLINE bool igl::ray_box_intersect(
   {
     tmax = tzmax;
   }
-  return ( (tmin < t1) && (tmax > t0) );
+  if(!( (tmin < t1) && (tmax > t0) ))
+  {
+    return false;
+  }
+  return true;
+#endif
 }

include/igl/ray_box_intersect.h

@@ -21,6 +21,9 @@ namespace igl
   //   box  axis aligned box
   //   t0  hit only if hit.t less than t0
   //   t1  hit only if hit.t greater than t1
+  // Outputs:
+  //   tmin  minimum of interval of overlap within [t0,t1]
+  //   tmax  maximum of interval of overlap within [t0,t1]
   // Returns true if hit
   template <
     typename Derivedsource,
@@ -31,7 +34,9 @@ namespace igl
     const Eigen::PlainObjectBase<Deriveddir> & dir,
     const Eigen::AlignedBox<Scalar,3> & box,
     const Scalar & t0,
-    const Scalar & t1);
+    const Scalar & t1,
+    Scalar & tmin,
+    Scalar & tmax);
 }
 #ifndef IGL_STATIC_LIBRARY
 #  include "ray_box_intersect.cpp"