Merge remote-tracking branch 'upstream/master' into develop

Former-commit-id: b62250263d007a5e3f34bbeb56571ab470c43397

README.md

@@ -190,6 +190,8 @@ Libigl is used by many research groups around the world. In 2015, it won the
 Eurographics/ACM Symposium on Geometry Processing software award. Here are a
 few labs/companies/institutions using libigl:
 
+ - [Adobe Research](http://www.adobe.com/technology/)  
+ - [Pixar Research](http://graphics.pixar.com/research/)
  - [Spine by Esoteric Software](http://esotericsoftware.com/) is an animation tool dedicated to 2D characters.
  - Columbia University, [Columbia Computer Graphics Group](http://www.cs.columbia.edu/cg/), USA
  - [Cornell University](http://www.graphics.cornell.edu/), USA
@@ -200,13 +202,18 @@ few labs/companies/institutions using libigl:
  - [National Institute of Informatics](http://www.nii.ac.jp/en/), Japan
  - New York University, [Media Research Lab](http://mrl.nyu.edu/), USA
  - NYUPoly, [Game Innovation Lab](http://game.engineering.nyu.edu/), USA
- - [Telecom ParisTech](http://www.telecom-paristech.fr/en/formation-et-innovation-dans-le-numerique.html), Paris, France
+ - [TU Berlin](https://www.cg.tu-berlin.de), Germany
  - [TU Delft](http://www.tudelft.nl/en/), Netherlands
+ - [Telecom ParisTech](http://www.telecom-paristech.fr/en/formation-et-innovation-dans-le-numerique.html), Paris, France
  - [Universidade Federal de Santa Catarina](http://mtm.ufsc.br/~leo/), Brazil
- - [Università della Svizzera Italiana](http://www.usi.ch/en), Switzerland
  - [University College London](http://vecg.cs.ucl.ac.uk/), England
+ - [University of California Berkeley](http://vis.berkeley.edu/), USA
  - [University of Cambridge](http://www.cam.ac.uk/), England
  - [University of Pennsylvania](http://cg.cis.upenn.edu/), USA
+ - [University of Texas at Austin](http://www.cs.utexas.edu/users/evouga/), USA
+ - [University of Victoria](https://www.csc.uvic.ca/Research/graphics/), Canada
+ - [Università della Svizzera Italiana](http://www.usi.ch/en), Switzerland
+ - [Université Toulouse III Paul Sabatier](http://www.univ-tlse3.fr/), France
 
 
 ## Contact

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,10 +307,14 @@ 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>
 #include <list>
+#include <queue>
+#include <stack>
 
 template <typename DerivedV, int DIM>
   template <typename Derivedbb_mins, typename Derivedbb_maxs>
@@ -1149,4 +1179,160 @@ 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();
+  {
+    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
+    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
+{
+#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>
+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;
+  {
+    Scalar _1,_2;
+    if(!ray_box_intersect(origin,dir,m_box,t0,min_t,_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
+    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;
+  }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,135 @@
+// 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);
+  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;
+    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);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+// generated by autoexplicit.sh
+template void igl::ambient_occlusion<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(std::__1::function<bool (Eigen::Matrix<float, 3, 1, 0, 3, 1> const&, Eigen::Matrix<float, 3, 1, 0, 3, 1> const&)> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+// generated by autoexplicit.sh
+template void igl::ambient_occlusion<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(std::__1::function<bool (Eigen::Matrix<float, 3, 1, 0, 3, 1> const&, Eigen::Matrix<float, 3, 1, 0, 3, 1> const&)> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+// generated by autoexplicit.sh
+template void igl::ambient_occlusion<Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(std::__1::function<bool (Eigen::Matrix<float, 3, 1, 0, 3, 1> const&, Eigen::Matrix<float, 3, 1, 0, 3, 1> const&)> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+#endif

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,138 @@
+#include "ray_box_intersect.h"
+#include <vector>
+
+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,
+  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::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 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;
+  }
+  if(!( (tmin < t1) && (tmax > t0) ))
+  {
+    return false;
+  }
+  return true;
+#endif
+}

include/igl/ray_box_intersect.h

@@ -0,0 +1,44 @@
+// 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
+  // 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,
+    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,
+    Scalar & tmin,
+    Scalar & tmax);
+}
+#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);
 }

include/igl/xml/XMLSerializable.h

@@ -1,6 +1,7 @@
 #ifndef IGL_XML_XMLSERIALIZABLE_H
 #define IGL_XML_XMLSERIALIZABLE_H
 
+#include "serialize_xml.h"
 #include "../igl_inline.h"
 #include "../serialize.h"
 

include/igl/xml/serialization_test.cpp

@@ -8,7 +8,8 @@
 //#define IGL_SERIALIZATION_TEST_H
 
 //#include <igl/Timer.h>
-#include <igl/xml/serialize_xml.h>
+#include "serialize_xml.h"
+#include "XMLSerializable.h"
 
 namespace igl
 {

include/igl/xml/serialize_xml.cpp

@@ -9,6 +9,7 @@
 #include "serialize_xml.h"
 #include "../STR.h"
 #include "../serialize.h"
+#include "XMLSerializable.h"
 
 #include <iterator>
 #include <limits>

include/igl/xml/serialize_xml.h

@@ -17,7 +17,6 @@
 
 #include "../igl_inline.h"
 
-#include "XMLSerializable.h"
 
 #include <Eigen/Dense>
 #include <Eigen/Sparse>
@@ -38,6 +37,7 @@ namespace igl
 {
   namespace xml
   {
+    struct XMLSerializableBase;
     // serializes the given object either to a xml file or to the provided doc data
     //
     // Templates:

tutorial/601_Serialization/main.cpp

@@ -8,22 +8,22 @@
 Eigen::MatrixXd V;
 Eigen::MatrixXi F;
 
-//// derive from igl::Serializable to serialize your own type
-//struct State : public igl::Serializable
-//{
-//  Eigen::MatrixXd V;
-//  Eigen::MatrixXi F;
-//  std::vector<int> ids;
-//
-//  // You have to define this function to
-//  // register the fields you want to serialize
-//  void InitSerialization()
-//  {
-//    this->Add(V  , "V");
-//    this->Add(F  , "F");
-//    this->Add(ids, "ids");
-//  }
-//};
+// derive from igl::Serializable to serialize your own type
+struct State : public igl::Serializable
+{
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  std::vector<int> ids;
+
+  // You have to define this function to
+  // register the fields you want to serialize
+  void InitSerialization()
+  {
+    this->Add(V  , "V");
+    this->Add(F  , "F");
+    this->Add(ids, "ids");
+  }
+};
 
 //// alternatively you can do it like the following to have
 //// a non-intrusive serialization:

tutorial/cmake/FindMATLAB.cmake

@@ -203,6 +203,11 @@ SET(MATLAB_LIBRARIES
 
 IF(MATLAB_INCLUDE_DIR AND MATLAB_LIBRARIES)
   SET(MATLAB_FOUND 1)
+  message(STATUS "Found MATLAB: ${MATLAB_INCLUDE_DIR}")
+else()
+  if (NOT MATLAB_FIND_QUIETLY)
+    message(FATAL_ERROR "could NOT find MATLAB")
+  endif (NOT MATLAB_FIND_QUIETLY)
 ENDIF(MATLAB_INCLUDE_DIR AND MATLAB_LIBRARIES)
 
 MARK_AS_ADVANCED(