Merge remote-tracking branch 'upstream/master'

Former-commit-id: f36db27018db03e30c01a567709a730c06261bed

README.md

@@ -34,6 +34,9 @@ We use the [Eigen](http://eigen.tuxfamily.org) library heavily in our code. Our
 group prototypes a lot in MATLAB, and we have a useful [MATLAB to libigl+Eigen
 conversion table](matlab-to-eigen.html).
 
+We regularly test compiling our library on Mac OS X with clang, Linux with gcc
+and Windows with Visual Studio 2015 Community Edition.
+
 ## Tutorial
 
 As of version 1.0, libigl includes an introductory
@@ -191,10 +194,12 @@ Eurographics/ACM Symposium on Geometry Processing software award. Here are a
 few labs/companies/institutions using libigl:
 
  - [Adobe Research](http://www.adobe.com/technology/)  
+ - [Mesh](http://meshconsultants.ca/), consultants, Canada
  - [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
+ - [Czech Technical University in Prague](http://dcgi.felk.cvut.cz/), Czech
  - EPF Lausanne, [Computer Graphics and Geometry Laboratory](http://lgg.epfl.ch/people.php), Switzerland
  - ETH Zurich, [Interactive Geometry Lab](http://igl.ethz.ch/) and [Advanced Technologies Lab](http://ait.inf.ethz.ch/), Swizterland
  - George Mason University, [CraGL](http://cs.gmu.edu/~ygingold/), USA
@@ -204,6 +209,7 @@ few labs/companies/institutions using libigl:
  - NYUPoly, [Game Innovation Lab](http://game.engineering.nyu.edu/), USA
  - [TU Berlin](https://www.cg.tu-berlin.de), Germany
  - [TU Delft](http://www.tudelft.nl/en/), Netherlands
+ - [TU Wien](https://www.tuwien.ac.at/en/tuwien_home/), Austria
  - [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
  - [University College London](http://vecg.cs.ucl.ac.uk/), England
@@ -214,6 +220,7 @@ few labs/companies/institutions using libigl:
  - [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
+ - [Zhejiang University](http://www.math.zju.edu.cn/cagd/), China
 
 
 ## Contact

include/igl/AABB.h

@@ -284,14 +284,6 @@ private:
         Scalar & sqr_d,
         int & i,
         RowVectorDIMS & c) const;
-public:
-      static
-      inline void barycentric_coordinates(
-        const RowVectorDIMS & p, 
-        const RowVectorDIMS & a, 
-        const RowVectorDIMS & b, 
-        const RowVectorDIMS & c,
-        Eigen::Matrix<Scalar,1,3> & bary);
 public:
       EIGEN_MAKE_ALIGNED_OPERATOR_NEW
     };
@@ -300,10 +292,12 @@ public:
 // Implementation
 #include "EPS.h"
 #include "barycenter.h"
+#include "barycentric_coordinates.h"
 #include "colon.h"
 #include "colon.h"
 #include "doublearea.h"
 #include "matlab_format.h"
+#include "point_simplex_squared_distance.h"
 #include "project_to_line_segment.h"
 #include "sort.h"
 #include "volume.h"
@@ -398,6 +392,10 @@ inline void igl::AABB<DerivedV,DIM>::init(
 {
   using namespace Eigen;
   using namespace std;
+  if(V.size() == 0 || Ele.size() == 0 || I.size() == 0)
+  {
+    return;
+  }
   assert(DIM == V.cols() && "V.cols() should matched declared dimension");
   //const Scalar inf = numeric_limits<Scalar>::infinity();
   m_box = AlignedBox<Scalar,DIM>();
@@ -1026,101 +1024,11 @@ inline void igl::AABB<DerivedV,DIM>::leaf_squared_distance(
 {
   using namespace Eigen;
   using namespace std;
-
-  // Simplex size
-  const size_t ss = Ele.cols();
-  // Only one element per node
-  // plane unit normal
-  bool inside_triangle = false;
-  Scalar d_j = std::numeric_limits<Scalar>::infinity();
-  RowVectorDIMS pp;
-  // Only consider triangles, and non-degenerate triangles at that
-  if(ss == 3 && 
-      Ele(m_primitive,0) != Ele(m_primitive,1) && 
-      Ele(m_primitive,1) != Ele(m_primitive,2) && 
-      Ele(m_primitive,2) != Ele(m_primitive,0))
-  {
-    assert(DIM == 3 && "Only implemented for 3D triangles");
-    typedef Eigen::Matrix<Scalar,1,3> RowVector3S;
-    // can't be const because of annoying DIM template
-    RowVector3S v10(0,0,0);
-    v10.head(DIM) = (V.row(Ele(m_primitive,1))- V.row(Ele(m_primitive,0)));
-    RowVector3S v20(0,0,0);
-    v20.head(DIM) = (V.row(Ele(m_primitive,2))- V.row(Ele(m_primitive,0)));
-    const RowVectorDIMS n = (v10.cross(v20)).head(DIM);
-    Scalar n_norm = n.norm();
-    if(n_norm > 0)
-    {
-      const RowVectorDIMS un = n/n.norm();
-      // vector to plane
-      const RowVectorDIMS bc = 
-        1./3.*
-        ( V.row(Ele(m_primitive,0))+
-          V.row(Ele(m_primitive,1))+
-          V.row(Ele(m_primitive,2)));
-      const auto & v = p-bc;
-      // projected point on plane
-      d_j = v.dot(un);
-      pp = p - d_j*un;
-      // determine if pp is inside triangle
-      Eigen::Matrix<Scalar,1,3> b;
-      barycentric_coordinates(
-            pp,
-            V.row(Ele(m_primitive,0)),
-            V.row(Ele(m_primitive,1)),
-            V.row(Ele(m_primitive,2)),
-            b);
-      inside_triangle = fabs(fabs(b(0)) + fabs(b(1)) + fabs(b(2)) - 1.) <= 1e-10;
-    }
-  }
-  const auto & point_point_squared_distance = [&](const RowVectorDIMS & s)
-  {
-    const Scalar sqr_d_s = (p-s).squaredNorm();
-    set_min(p,sqr_d_s,m_primitive,s,sqr_d,i,c);
-  };
-  if(inside_triangle)
-  {
-    // point-triangle squared distance
-    const Scalar sqr_d_j = d_j*d_j;
-    //cout<<"point-triangle..."<<endl;
-    set_min(p,sqr_d_j,m_primitive,pp,sqr_d,i,c);
-  }else
-  {
-    if(ss >= 2)
-    {
-      // point-segment distance
-      // number of edges
-      size_t ne = ss==3?3:1;
-      for(size_t x = 0;x<ne;x++)
-      {
-        const size_t e1 = Ele(m_primitive,(x+1)%ss);
-        const size_t e2 = Ele(m_primitive,(x+2)%ss);
-        const RowVectorDIMS & s = V.row(e1);
-        const RowVectorDIMS & d = V.row(e2);
-        // Degenerate edge
-        if(e1 == e2 || (s-d).squaredNorm()==0)
-        {
-          // only consider once
-          if(e1 < e2)
-          {
-            point_point_squared_distance(s);
-          }
-          continue;
-        }
-        Matrix<Scalar,1,1> sqr_d_j_x(1,1);
-        Matrix<Scalar,1,1> t(1,1);
-        project_to_line_segment(p,s,d,t,sqr_d_j_x);
-        const RowVectorDIMS q = s+t(0)*(d-s);
-        set_min(p,sqr_d_j_x(0),m_primitive,q,sqr_d,i,c);
-      }
-    }else
-    {
-      // So then Ele is just a list of points...
-      assert(ss == 1);
-      const RowVectorDIMS & s = V.row(Ele(m_primitive,0));
-      point_point_squared_distance(s);
-    }
-  }
+  RowVectorDIMS c_candidate;
+  Scalar sqr_d_candidate;
+  igl::point_simplex_squared_distance<DIM>(
+    p,V,Ele,m_primitive,sqr_d_candidate,c_candidate);
+  set_min(p,sqr_d_candidate,m_primitive,c_candidate,sqr_d,i,c);
 }
 
 
@@ -1153,30 +1061,6 @@ inline void igl::AABB<DerivedV,DIM>::set_min(
 }
 
 
-template <typename DerivedV, int DIM>
-inline void
-igl::AABB<DerivedV,DIM>::barycentric_coordinates(
-  const RowVectorDIMS & p, 
-  const RowVectorDIMS & a, 
-  const RowVectorDIMS & b, 
-  const RowVectorDIMS & c,
-  Eigen::Matrix<Scalar,1,3> & bary)
-{
-  // http://gamedev.stackexchange.com/a/23745
-  const RowVectorDIMS v0 = b - a;
-  const RowVectorDIMS v1 = c - a;
-  const RowVectorDIMS v2 = p - a;
-  Scalar d00 = v0.dot(v0);
-  Scalar d01 = v0.dot(v1);
-  Scalar d11 = v1.dot(v1);
-  Scalar d20 = v2.dot(v0);
-  Scalar d21 = v2.dot(v1);
-  Scalar denom = d00 * d11 - d01 * d01;
-  bary(1) = (d11 * d20 - d01 * d21) / denom;
-  bary(2) = (d00 * d21 - d01 * d20) / denom;
-  bary(0) = 1.0f - bary(1) - bary(2);
-}
-
 template <typename DerivedV, int DIM>
 inline bool 
 igl::AABB<DerivedV,DIM>::intersect_ray(
@@ -1201,7 +1085,7 @@ igl::AABB<DerivedV,DIM>::intersect_ray(
     // 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);
+    return ray_mesh_intersect(origin,dir,V,Ele.row(m_primitive),hits);
   }
   std::vector<igl::Hit> left_hits;
   std::vector<igl::Hit> right_hits;
@@ -1248,7 +1132,7 @@ igl::AABB<DerivedV,DIM>::intersect_ray(
       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)&&
+        ray_mesh_intersect(origin,dir,V,Ele.row(tree->m_primitive),leaf_hit)&&
         leaf_hit.t < hit.t)
       {
         hit = leaf_hit;
@@ -1302,7 +1186,7 @@ igl::AABB<DerivedV,DIM>::intersect_ray(
     // 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);
+    return ray_mesh_intersect(origin,dir,V,Ele.row(m_primitive),hit);
   }
 
   // Doesn't seem like smartly choosing left before/after right makes a

include/igl/ambient_occlusion.cpp

@@ -32,9 +32,9 @@ IGL_INLINE void igl::ambient_occlusion(
   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
+#pragma omp parallel for
   for(int p = 0;p<n;p++)
   {
     const Vector3f origin = P.row(p).template cast<float>();

include/igl/barycentric_coordinates.cpp

@@ -55,33 +55,55 @@ template <
   typename DerivedC,
   typename DerivedL>
 IGL_INLINE void igl::barycentric_coordinates(
-  const Eigen::PlainObjectBase<DerivedP> & P,
-  const Eigen::PlainObjectBase<DerivedA> & A,
-  const Eigen::PlainObjectBase<DerivedB> & B,
-  const Eigen::PlainObjectBase<DerivedC> & C,
+  const Eigen::MatrixBase<DerivedP> & P,
+  const Eigen::MatrixBase<DerivedA> & A,
+  const Eigen::MatrixBase<DerivedB> & B,
+  const Eigen::MatrixBase<DerivedC> & C,
   Eigen::PlainObjectBase<DerivedL> & L)
 {
   using namespace Eigen;
-  assert(P.cols() == 2 && "query must be in 2d");
-  assert(A.cols() == 2 && "corners must be in 2d");
-  assert(B.cols() == 2 && "corners must be in 2d");
-  assert(C.cols() == 2 && "corners must be in 2d");
+#ifndef NDEBUG
+  const int DIM = P.cols();
+  assert(A.cols() == DIM && "corners must be in same dimension as query");
+  assert(B.cols() == DIM && "corners must be in same dimension as query");
+  assert(C.cols() == DIM && "corners must be in same dimension as query");
   assert(P.rows() == A.rows() && "Must have same number of queries as corners");
   assert(A.rows() == B.rows() && "Corners must be same size");
   assert(A.rows() == C.rows() && "Corners must be same size");
-  typedef Matrix<typename DerivedL::Scalar,DerivedL::RowsAtCompileTime,1> 
-    VectorXS;
-  // Total area
-  VectorXS dblA,LA,LB,LC;
-  doublearea(P,B,C,LA);
-  doublearea(A,P,C,LB);
-  doublearea(A,B,P,LC);
-  doublearea(A,B,C,dblA);
+#endif
+
+  // http://gamedev.stackexchange.com/a/23745
+  typedef 
+    Eigen::Array<
+      typename DerivedP::Scalar,
+               DerivedP::RowsAtCompileTime,
+               DerivedP::ColsAtCompileTime>
+    ArrayS;
+  typedef 
+    Eigen::Array<
+      typename DerivedP::Scalar,
+               DerivedP::RowsAtCompileTime,
+               1>
+    VectorS;
+
+  const ArrayS v0 = B.array() - A.array();
+  const ArrayS v1 = C.array() - A.array();
+  const ArrayS v2 = P.array() - A.array();
+  VectorS d00 = (v0*v0).rowwise().sum();
+  VectorS d01 = (v0*v1).rowwise().sum();
+  VectorS d11 = (v1*v1).rowwise().sum();
+  VectorS d20 = (v2*v0).rowwise().sum();
+  VectorS d21 = (v2*v1).rowwise().sum();
+  VectorS denom = d00 * d11 - d01 * d01;
   L.resize(P.rows(),3);
-  L<<LA,LB,LC;
-  L.array().colwise() /= dblA.array();
+  L.col(1) = (d11 * d20 - d01 * d21) / denom;
+  L.col(2) = (d00 * d21 - d01 * d20) / denom;
+  L.col(0) = 1.0f -(L.col(1) + L.col(2)).array();
 }
 
 #ifdef IGL_STATIC_LIBRARY
 template void igl::barycentric_coordinates<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::barycentric_coordinates<Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(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::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&);
+template void igl::barycentric_coordinates<Eigen::Matrix<double, 1, 2, 1, 1, 2>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&);
+template void igl::barycentric_coordinates<Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&);
 #endif

include/igl/barycentric_coordinates.h

@@ -39,15 +39,13 @@ namespace igl
   // Compute barycentric coordinates in a triangle
   //
   // Inputs:
-  //   P  #P by 2 Query points in 2d
-  //   A  #P by 2 Triangle corners in 2d
-  //   B  #P by 2 Triangle corners in 2d
-  //   C  #P by 2 Triangle corners in 2d
+  //   P  #P by dim Query points
+  //   A  #P by dim Triangle corners
+  //   B  #P by dim Triangle corners
+  //   C  #P by dim Triangle corners
   // Outputs:
-  //   L  #P by e list of barycentric coordinates
+  //   L  #P by 3 list of barycentric coordinates
   //   
-  // Known bugs: this code is not tested (and probably will not work) for
-  // triangles and queries in 3D even if the query lives in/on the triangle.
   template <
     typename DerivedP,
     typename DerivedA,
@@ -55,10 +53,10 @@ namespace igl
     typename DerivedC,
     typename DerivedL>
   IGL_INLINE void barycentric_coordinates(
-    const Eigen::PlainObjectBase<DerivedP> & P,
-    const Eigen::PlainObjectBase<DerivedA> & A,
-    const Eigen::PlainObjectBase<DerivedB> & B,
-    const Eigen::PlainObjectBase<DerivedC> & C,
+    const Eigen::MatrixBase<DerivedP> & P,
+    const Eigen::MatrixBase<DerivedA> & A,
+    const Eigen::MatrixBase<DerivedB> & B,
+    const Eigen::MatrixBase<DerivedC> & C,
     Eigen::PlainObjectBase<DerivedL> & L);
 
 }

include/igl/copyleft/cgal/piecewise_constant_winding_number.cpp

@@ -24,6 +24,6 @@ IGL_INLINE bool igl::copyleft::cgal::piecewise_constant_winding_number(
   // resolve intersections
   remesh_self_intersections(V,F,{},VV,FF,IF,J,IM);
   // _apply_ duplicate vertex mapping IM to FF
-  for_each(FF.data(),FF.data()+FF.size(),[&IM](int & a){a=IM(a);});
-  return piecewise_constant_winding_number(FF);
+  std::for_each(FF.data(),FF.data()+FF.size(),[&IM](int & a){a=IM(a);});
+  return igl::piecewise_constant_winding_number(FF);
 }

include/igl/copyleft/cgal/submesh_aabb_tree.cpp

@@ -50,5 +50,5 @@ IGL_INLINE void igl::copyleft::cgal::submesh_aabb_tree(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instanciation
-template void igl::copyleft::cgal::submesh_aabb_tree<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, CGAL::Epeck>(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, CGAL::AABB_tree<CGAL::AABB_traits<CGAL::Epeck, CGAL::AABB_triangle_primitive<CGAL::Epeck, std::__1::vector<CGAL::Epeck::Triangle_3, std::__1::allocator<CGAL::Epeck::Triangle_3> >::iterator, CGAL::Boolean_tag<false> > > >&, std::__1::vector<CGAL::Epeck::Triangle_3, std::__1::allocator<CGAL::Epeck::Triangle_3> >&, std::__1::vector<bool, std::__1::allocator<bool> >&);
+template void igl::copyleft::cgal::submesh_aabb_tree<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, CGAL::Epeck>(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, CGAL::AABB_tree<CGAL::AABB_traits<CGAL::Epeck, CGAL::AABB_triangle_primitive<CGAL::Epeck, std::vector<CGAL::Epeck::Triangle_3, std::allocator<CGAL::Epeck::Triangle_3> >::iterator, CGAL::Boolean_tag<false> > > >&, std::vector<CGAL::Epeck::Triangle_3, std::allocator<CGAL::Epeck::Triangle_3> >&, std::vector<bool, std::allocator<bool> >&);
 #endif

include/igl/embree/unproject_onto_mesh.cpp

@@ -7,8 +7,7 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "unproject_onto_mesh.h"
 #include "EmbreeIntersector.h"
-#include <igl/unproject.h>
-#include <igl/embree/unproject_in_mesh.h>
+#include "../unproject_onto_mesh.h"
 #include <vector>
 
 IGL_INLINE bool igl::embree::unproject_onto_mesh(
@@ -23,20 +22,14 @@ IGL_INLINE bool igl::embree::unproject_onto_mesh(
 {
   using namespace std;
   using namespace Eigen;
-  MatrixXd obj;
-  vector<igl::Hit> hits;
-
-  // This is lazy, it will find more than just the first hit
-  unproject_in_mesh(pos,model,proj,viewport,ei,obj,hits);
-
-  if (hits.size()> 0)
+  const auto & shoot_ray = [&ei](
+    const Eigen::Vector3f& s,
+    const Eigen::Vector3f& dir,
+    igl::Hit & hit)->bool
   {
-    bc << 1.0-hits[0].u-hits[0].v, hits[0].u, hits[0].v;
-    fid = hits[0].id;
-    return true;
-  }
-
-  return false;
+    return ei.intersectRay(s,dir,hit);
+  };
+  return igl::unproject_onto_mesh(pos,model,proj,viewport,shoot_ray,fid,bc);
 }
 
 IGL_INLINE bool igl::embree::unproject_onto_mesh(
@@ -50,14 +43,14 @@ IGL_INLINE bool igl::embree::unproject_onto_mesh(
   int& vid)
 {
   Eigen::Vector3f bc;
-  bool hit = unproject_onto_mesh(pos,F,model,proj,viewport,ei,fid,bc);
-  int i;
-  if (hit)
+  if(!igl::embree::unproject_onto_mesh(pos,F,model,proj,viewport,ei,fid,bc))
   {
-    bc.maxCoeff(&i);
-    vid = F(fid,i);
+    return false;
   }
-  return hit;
+  int i;
+  bc.maxCoeff(&i);
+  vid = F(fid,i);
+  return true;
 }
 
 

include/igl/extract_manifold_patches.h

@@ -6,13 +6,9 @@
 #include <vector>
 
 namespace igl {
-    // Extract a set of maximal manifold patches from a given mesh.
-    // A maximal manifold patch is a subset of the input faces that is
-    // connected and manifold, and it cannot be expanded further by
-    // including more faces.
-    //
-    // Assumes the input mesh have all self-intersection resolved.  See
-    // ``igl::cgal::remesh_self_intersection`` for more details.
+    // Extract a set of maximal patches from a given mesh.
+    // A maximal patch is a subset of the input faces that are connected via
+    // manifold edges; a patch is as large as possible.
     //
     // Inputs:
     //   F  #F by 3 list representing triangles.
@@ -25,15 +21,15 @@ namespace igl {
     // Returns:
     //   number of manifold patches.
     template <
-        typename DerivedF,
-        typename DerivedEMAP,
-        typename uE2EType,
-        typename DerivedP>
+      typename DerivedF,
+      typename DerivedEMAP,
+      typename uE2EType,
+      typename DerivedP>
     IGL_INLINE size_t extract_manifold_patches(
-            const Eigen::PlainObjectBase<DerivedF>& F,
-            const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
-            const std::vector<std::vector<uE2EType> >& uE2E,
-            Eigen::PlainObjectBase<DerivedP>& P);
+      const Eigen::PlainObjectBase<DerivedF>& F,
+      const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+      const std::vector<std::vector<uE2EType> >& uE2E,
+      Eigen::PlainObjectBase<DerivedP>& P);
 }
 #ifndef IGL_STATIC_LIBRARY
 #  include "extract_manifold_patches.cpp"

include/igl/ply.h.REMOVED.git-id

@@ -1 +1 @@
-a11a359439a307b36f2320e726048ce3dd344fda
+aa5b1239a610dbd02506563452ec69172ae4165a

include/igl/point_simplex_squared_distance.cpp

@@ -0,0 +1,148 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 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 "point_simplex_squared_distance.h"
+#include "project_to_line_segment.h"
+#include "barycentric_coordinates.h"
+#include <Eigen/Geometry>
+#include <limits>
+#include <cassert>
+
+template <
+  int DIM,
+  typename Derivedp,
+  typename DerivedV,
+  typename DerivedEle,
+  typename Derivedsqr_d,
+  typename Derivedc>
+IGL_INLINE void igl::point_simplex_squared_distance(
+  const Eigen::PlainObjectBase<Derivedp> & p,
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedEle> & Ele,
+  const typename DerivedEle::Index primitive,
+  Derivedsqr_d & sqr_d,
+  Eigen::PlainObjectBase<Derivedc> & c)
+{
+  assert(p.size() == DIM);
+  assert(V.cols() == DIM);
+  assert(Ele.cols() <= DIM+1);
+  assert(Ele.cols() <= 3 && "Only simplices up to triangles are considered");
+
+  typedef Derivedsqr_d Scalar;
+  typedef typename Eigen::Matrix<Scalar,1,DIM> RowVectorDIMS;
+  sqr_d = std::numeric_limits<Scalar>::infinity();
+  const auto & set_min = 
+    [&sqr_d,&c](const Scalar & sqr_d_candidate, const RowVectorDIMS & c_candidate)
+  {
+    if(sqr_d_candidate < sqr_d)
+    {
+      sqr_d = sqr_d_candidate;
+      c = c_candidate;
+    }
+  };
+
+  // Simplex size
+  const size_t ss = Ele.cols();
+  // Only one element per node
+  // plane unit normal
+  bool inside_triangle = false;
+  Scalar d_j = std::numeric_limits<Scalar>::infinity();
+  RowVectorDIMS pp;
+  // Only consider triangles, and non-degenerate triangles at that
+  if(ss == 3 && 
+      Ele(primitive,0) != Ele(primitive,1) && 
+      Ele(primitive,1) != Ele(primitive,2) && 
+      Ele(primitive,2) != Ele(primitive,0))
+  {
+    assert(DIM == 3 && "Only implemented for 3D triangles");
+    typedef Eigen::Matrix<Scalar,1,3> RowVector3S;
+    // can't be const because of annoying DIM template
+    RowVector3S v10(0,0,0);
+    v10.head(DIM) = (V.row(Ele(primitive,1))- V.row(Ele(primitive,0)));
+    RowVector3S v20(0,0,0);
+    v20.head(DIM) = (V.row(Ele(primitive,2))- V.row(Ele(primitive,0)));
+    const RowVectorDIMS n = (v10.cross(v20)).head(DIM);
+    Scalar n_norm = n.norm();
+    if(n_norm > 0)
+    {
+      const RowVectorDIMS un = n/n.norm();
+      // vector to plane
+      const RowVectorDIMS bc = 
+        1./3.*
+        ( V.row(Ele(primitive,0))+
+          V.row(Ele(primitive,1))+
+          V.row(Ele(primitive,2)));
+      const auto & v = p-bc;
+      // projected point on plane
+      d_j = v.dot(un);
+      pp = p - d_j*un;
+      // determine if pp is inside triangle
+      Eigen::Matrix<Scalar,1,3> b;
+      barycentric_coordinates(
+            pp,
+            V.row(Ele(primitive,0)),
+            V.row(Ele(primitive,1)),
+            V.row(Ele(primitive,2)),
+            b);
+      inside_triangle = fabs(fabs(b(0)) + fabs(b(1)) + fabs(b(2)) - 1.) <= 1e-10;
+    }
+  }
+  const auto & point_point_squared_distance = [&](const RowVectorDIMS & s)
+  {
+    const Scalar sqr_d_s = (p-s).squaredNorm();
+    set_min(sqr_d_s,s);
+  };
+  if(inside_triangle)
+  {
+    // point-triangle squared distance
+    const Scalar sqr_d_j = d_j*d_j;
+    //cout<<"point-triangle..."<<endl;
+    set_min(sqr_d_j,pp);
+  }else
+  {
+    if(ss >= 2)
+    {
+      // point-segment distance
+      // number of edges
+      size_t ne = ss==3?3:1;
+      for(size_t x = 0;x<ne;x++)
+      {
+        const size_t e1 = Ele(primitive,(x+1)%ss);
+        const size_t e2 = Ele(primitive,(x+2)%ss);
+        const RowVectorDIMS & s = V.row(e1);
+        const RowVectorDIMS & d = V.row(e2);
+        // Degenerate edge
+        if(e1 == e2 || (s-d).squaredNorm()==0)
+        {
+          // only consider once
+          if(e1 < e2)
+          {
+            point_point_squared_distance(s);
+          }
+          continue;
+        }
+        Eigen::Matrix<Scalar,1,1> sqr_d_j_x(1,1);
+        Eigen::Matrix<Scalar,1,1> t(1,1);
+        project_to_line_segment(p,s,d,t,sqr_d_j_x);
+        const RowVectorDIMS q = s+t(0)*(d-s);
+        set_min(sqr_d_j_x(0),q);
+      }
+    }else
+    {
+      // So then Ele is just a list of points...
+      assert(ss == 1);
+      const RowVectorDIMS & s = V.row(Ele(primitive,0));
+      point_point_squared_distance(s);
+    }
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instanciation
+template void igl::point_simplex_squared_distance<2, Eigen::Matrix<double, 1, 2, 1, 1, 2>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, Eigen::Matrix<double, 1, 2, 1, 1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, double&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> >&);
+template void igl::point_simplex_squared_distance<3, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, double&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&);
+#endif

include/igl/point_simplex_squared_distance.h

@@ -0,0 +1,43 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 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_POINT_SIMPLEX_SQUARED_DISTANCE_H
+#define IGL_POINT_SIMPLEX_SQUARED_DISTANCE_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Determine squared distance from a point to linear simplex
+  //
+  // Inputs:
+  //   p  d-long query point
+  //   V  #V by d list of vertices
+  //   Ele  #Ele by ss<=d+1 list of simplex indices into V
+  //   i  index into Ele of simplex
+  // Outputs:
+  //   sqr_d  squared distance of Ele(i) to p
+  //   c  closest point on Ele(i) 
+  //
+  template <
+    int DIM,
+    typename Derivedp,
+    typename DerivedV,
+    typename DerivedEle,
+    typename Derivedsqr_d,
+    typename Derivedc>
+  IGL_INLINE void point_simplex_squared_distance(
+    const Eigen::PlainObjectBase<Derivedp> & p,
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedEle> & Ele,
+    const typename DerivedEle::Index i,
+    Derivedsqr_d & sqr_d,
+    Eigen::PlainObjectBase<Derivedc> & c);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "point_simplex_squared_distance.cpp"
+#endif
+#endif

include/igl/pseudonormal_test.cpp

@@ -6,7 +6,8 @@
 // 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 "pseudonormal_test.h"
-#include "AABB.h"
+#include "barycentric_coordinates.h"
+#include "doublearea.h"
 #include "project_to_line_segment.h"
 #include <cassert>
 
@@ -49,7 +50,7 @@ IGL_INLINE void igl::pseudonormal_test(
   const double epsilon = 1e-12;
   if(area>MIN_DOUBLE_AREA)
   {
-    AABB<Eigen::MatrixXd,3>::barycentric_coordinates( c,A,B,C,b);
+    barycentric_coordinates( c,A,B,C,b);
     // Determine which normal to use
     const int type = (b.array()<=epsilon).cast<int>().sum();
     switch(type)

include/igl/ray_mesh_intersect.cpp

@@ -11,10 +11,10 @@ template <
   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,
+  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;
@@ -53,10 +53,10 @@ template <
   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,
+  const Eigen::MatrixBase<Derivedsource> & source,
+  const Eigen::MatrixBase<Deriveddir> & dir,
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
   igl::Hit & hit)
 {
   std::vector<igl::Hit> hits;
@@ -70,3 +70,8 @@ IGL_INLINE bool igl::ray_mesh_intersect(
     return false;
   }
 }
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instanciation
+template bool igl::ray_mesh_intersect<Eigen::Matrix<float, 3, 1, 0, 3, 1>, Eigen::Matrix<float, 3, 1, 0, 3, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<float, 3, 1, 0, 3, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 3, 1, 0, 3, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<igl::Hit, std::allocator<igl::Hit> >&);
+#endif

include/igl/ray_mesh_intersect.h

@@ -23,10 +23,10 @@ namespace igl
     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,
+    const Eigen::MatrixBase<Derivedsource> & source,
+    const Eigen::MatrixBase<Deriveddir> & dir,
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
     std::vector<igl::Hit> & hits);
   // Outputs:
   //   hit  first hit, set only if it exists
@@ -37,10 +37,10 @@ namespace igl
     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,
+    const Eigen::MatrixBase<Derivedsource> & source,
+    const Eigen::MatrixBase<Deriveddir> & dir,
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
     igl::Hit & hit);
 }
 #ifndef IGL_STATIC_LIBRARY

include/igl/readPLY.cpp

@@ -153,7 +153,6 @@ IGL_INLINE bool igl::readPLY(
     }
   }
   ply_close(in_ply);
-  fclose(fp);
   return true;
 }
 

include/igl/unproject_in_mesh.h

@@ -49,8 +49,9 @@ namespace igl
   //    model      model matrix
   //    proj       projection matrix
   //    viewport   vieweport vector
-  //    shoot_ray  function handle that outputs hits of a given ray against a
-  //      mesh (embedded in function handles as captured variable/data)
+  //    shoot_ray  function handle that outputs first hit of a given ray
+  //      against a mesh (embedded in function handles as captured
+  //      variable/data)
   // Outputs:
   //    obj        3d unprojected mouse point in mesh
   //    hits       vector of hits

include/igl/unproject_onto_mesh.cpp

@@ -0,0 +1,77 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Alec Jacobson
+//
+// 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 "unproject_onto_mesh.h"
+#include "unproject.h"
+#include "unproject_ray.h"
+#include "ray_mesh_intersect.h"
+#include <vector>
+
+template < typename DerivedV, typename DerivedF, typename Derivedbc>
+IGL_INLINE bool igl::unproject_onto_mesh(
+  const Eigen::Vector2f& pos,
+  const Eigen::Matrix4f& model,
+  const Eigen::Matrix4f& proj,
+  const Eigen::Vector4f& viewport,
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  int & fid,
+  Eigen::PlainObjectBase<Derivedbc> & bc)
+{
+  using namespace std;
+  using namespace Eigen;
+  const auto & shoot_ray = [&V,&F](
+    const Eigen::Vector3f& s,
+    const Eigen::Vector3f& dir,
+    igl::Hit & hit)->bool
+  {
+    std::vector<igl::Hit> hits;
+    if(!ray_mesh_intersect(s,dir,V,F,hits))
+    {
+      return false;
+    }
+    hit = hits[0];
+    return true;
+  };
+  return unproject_onto_mesh(pos,model,proj,viewport,shoot_ray,fid,bc);
+}
+
+template <typename Derivedbc>
+IGL_INLINE bool igl::unproject_onto_mesh(
+  const Eigen::Vector2f& pos,
+  const Eigen::Matrix4f& model,
+  const Eigen::Matrix4f& proj,
+  const Eigen::Vector4f& viewport,
+  const std::function<
+    bool(
+      const Eigen::Vector3f&,
+      const Eigen::Vector3f&,
+      igl::Hit &)
+      > & shoot_ray,
+  int & fid,
+  Eigen::PlainObjectBase<Derivedbc> & bc)
+{
+  using namespace std;
+  using namespace Eigen;
+  Vector3f s,dir;
+  unproject_ray(pos,model,proj,viewport,s,dir);
+  Hit hit;
+  if(!shoot_ray(s,dir,hit))
+  {
+    return false;
+  }
+  bc.resize(3);
+  bc << 1.0-hit.u-hit.v, hit.u, hit.v;
+  fid = hit.id;
+  return true;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+template bool igl::unproject_onto_mesh<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<float, 3, 1, 0, 3, 1> >(Eigen::Matrix<float, 2, 1, 0, 2, 1> const&, Eigen::Matrix<float, 4, 4, 0, 4, 4> const&, Eigen::Matrix<float, 4, 4, 0, 4, 4> const&, Eigen::Matrix<float, 4, 1, 0, 4, 1> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int&, Eigen::PlainObjectBase<Eigen::Matrix<float, 3, 1, 0, 3, 1> >&);
+#endif
+

include/igl/unproject_onto_mesh.h

@@ -0,0 +1,74 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 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_UNPROJECT_ONTO_MESH
+#define IGL_UNPROJECT_ONTO_MESH
+#include "igl_inline.h"
+#include "Hit.h"
+#include <Eigen/Core>
+#include <functional>
+
+namespace igl
+{
+  // Unproject a screen location (using current opengl viewport, projection, and
+  // model view) to a 3D position _onto_ a given mesh, if the ray through the
+  // given screen location (x,y) _hits_ the mesh.
+  //
+  // Inputs:
+  //    pos        screen space coordinates
+  //    model      model matrix
+  //    proj       projection matrix
+  //    viewport   vieweport vector
+  //    V   #V by 3 list of mesh vertex positions
+  //    F   #F by 3 list of mesh triangle indices into V
+  // Outputs:
+  //    fid  id of the first face hit
+  //    bc  barycentric coordinates of hit
+  // Returns true if there's a hit
+  template < typename DerivedV, typename DerivedF, typename Derivedbc>
+  IGL_INLINE bool unproject_onto_mesh(
+    const Eigen::Vector2f& pos,
+    const Eigen::Matrix4f& model,
+    const Eigen::Matrix4f& proj,
+    const Eigen::Vector4f& viewport,
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    int & fid,
+    Eigen::PlainObjectBase<Derivedbc> & bc);
+  //
+  // Inputs:
+  //    pos        screen space coordinates
+  //    model      model matrix
+  //    proj       projection matrix
+  //    viewport   vieweport vector
+  //    shoot_ray  function handle that outputs hits of a given ray against a
+  //      mesh (embedded in function handles as captured variable/data)
+  // Outputs:
+  //    fid  id of the first face hit
+  //    bc  barycentric coordinates of hit
+  // Returns true if there's a hit
+  template <typename Derivedbc>
+  IGL_INLINE bool unproject_onto_mesh(
+    const Eigen::Vector2f& pos,
+    const Eigen::Matrix4f& model,
+    const Eigen::Matrix4f& proj,
+    const Eigen::Vector4f& viewport,
+    const std::function<
+      bool(
+        const Eigen::Vector3f&,
+        const Eigen::Vector3f&,
+        igl::Hit  &)
+        > & shoot_ray,
+    int & fid,
+    Eigen::PlainObjectBase<Derivedbc> & bc);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "unproject_onto_mesh.cpp"
+#endif
+#endif
+
+

include/igl/viewer/Viewer.cpp

@@ -306,6 +306,22 @@ namespace viewer
     callback_mouse_scroll_data  = nullptr;
     callback_key_down_data      = nullptr;
     callback_key_up_data        = nullptr;
+
+#ifndef IGL_VIEWER_VIEWER_QUIET
+    const std::string usage(R"(igl::viewer::Viewer usage:
+  [drag]  Rotate scene
+  A,a     Toggle animation (tight draw loop)
+  I,i     Toggle invert normals
+  L,l     Toggle wireframe
+  O,o     Toggle orthographic/perspective projection
+  T,t     Toggle filled faces
+  Z       Snap to canonical view
+  [,]     Toggle between rotation control types (e.g. trackball, two-axis
+          valuator with fixed up)
+
+)");
+    std::cout<<usage;
+#endif
   }
 
   IGL_INLINE void Viewer::init_plugins()
@@ -456,18 +472,13 @@ namespace viewer
         return true;
 
     for (unsigned int i = 0; i<plugins.size(); ++i)
+    {
       if (plugins[i]->key_pressed(unicode_key, modifiers))
+      {
         return true;
-    const std::string usage(R"(igl::viewer::Viewer:
-A,a  Toggle animation (tight draw loop)
-I,i  Toggle invert normals
-L,l  Toggle wireframe
-O,o  Toggle orthographic/perspective projection
-T,t  Toggle filled faces
-Z    Snap to canonical view
-[,]  Toggle between rotation control types (e.g. trackball, two-axis valuator
-     with fixed up)
-)");
+      }
+    }
+
     switch(unicode_key)
     {
       case 'A':

include/igl/writePLY.cpp

@@ -132,7 +132,6 @@ IGL_INLINE bool igl::writePLY(
   }
 
   ply_close(ply);
-  fclose(fp);
   for(size_t i = 0;i<(size_t)F.rows();i++)
   {
     delete[] flist[i].verts;

python/tutorial/305_QuadraticProgramming.py

@@ -80,7 +80,7 @@ ux = igl.eigen.MatrixXd.Ones(V.rows(),1)
 
 # Equality constraint constrain solution to sum to 1
 Beq = igl.eigen.MatrixXd([[0.08]])
-Aeq = M.diagonal().transpose().sparseView()
+Aeq = M.diagonal().sparseView().transpose()
 
 # (Empty inequality constraints)
 solve(viewer)

scripts/pre-commit.sh

@@ -0,0 +1,33 @@
+#!/bin/sh
+#
+# To enable this hook, issue:
+#
+#     ln -s ../../scripts/pre-commit.sh .git/hooks/pre-commit
+#
+RED='\033[0;31m'
+NC='\033[0m'
+
+## Throw error if any files contain "std::__1"
+
+# list of changed files
+CHANGED=$(git diff --cached --name-only --diff-filter=ACM)
+# Ignore this file!
+CHANGED=$(echo "$CHANGED" | grep -v "scripts/pre-commit.sh")
+# Changed files containing the namespace "std::__1"
+STD1=$(grep -Il "std::__1" $CHANGED)
+if  [ $? -eq 0 ]; then
+  STD1=$(echo "$STD1" | sed -e "s/^/  /")
+  >&2 echo "[pre-commit hook] Error: Commit prohibited.
+
+The following files contain the offensive \"std::__1\" namespace: 
+
+${RED}$STD1${NC}
+
+Consider issueing:
+
+    sed -i '' -e \"s/std::__1/std/g\"" $STD1
+
+  exit 1
+else
+  exit 0
+fi

tutorial/305_QuadraticProgramming/main.cpp

@@ -86,7 +86,7 @@ int main(int argc, char *argv[])
   // Equality constraint constrain solution to sum to 1
   Beq.resize(1,1);
   Beq(0) = 0.08;
-  Aeq = M.diagonal().transpose().sparseView();
+  Aeq = M.diagonal().sparseView().transpose();
   // (Empty inequality constraints)
   solve(viewer);
   cout<<

tutorial/603_MEX/readOBJ_mex.cpp

@@ -36,9 +36,9 @@ void mexFunction(
   switch(nlhs)
   {
     case 2:
-      igl::prepare_lhs_index(F,plhs+1);
+      igl::matlab::prepare_lhs_index(F,plhs+1);
     case 1:
-      igl::prepare_lhs_double(V,plhs);
+      igl::matlab::prepare_lhs_double(V,plhs);
     default: break;
   }
 

tutorial/607_Picking/main.cpp

@@ -1,74 +1,50 @@
+#include "tutorial_shared_path.h"
 #include <igl/readOFF.h>
+#include <igl/unproject_onto_mesh.h>
 #include <igl/viewer/Viewer.h>
-#include <igl/embree/EmbreeIntersector.h>
-#include <igl/embree/unproject_onto_mesh.h>
-
-#include <algorithm>
-
-#include "tutorial_shared_path.h"
-
-// Mesh
-Eigen::MatrixXd V;
-Eigen::MatrixXi F;
-
-// Per-vertex color
-Eigen::MatrixXd C;
-
-igl::embree::EmbreeIntersector* ei;
-
-std::vector<int> picked_vertices;
-
-bool mouse_down(igl::viewer::Viewer& viewer, int button, int modifier)
-{
-
-  using namespace Eigen;
-  using namespace std;
-
-  int vid, fid;
-
-  // Cast a ray in the view direction starting from the mouse position
-  double x = viewer.current_mouse_x;
-  double y = viewer.core.viewport(3) - viewer.current_mouse_y;
-  bool hit = igl::embree::unproject_onto_mesh(Vector2f(x,y),
-                                F,
-                                viewer.core.view * viewer.core.model,
-                                viewer.core.proj,
-                                viewer.core.viewport,
-                                *ei,
-                                fid,
-                                vid);
-
-  // If the ray hits a face, assign a red color to the closest vertex
-  if (hit)
-  {
-    cerr << "Picked face(vertex): " << fid << " (" << vid << ")" << endl;
-    C.row(vid) << 1,0,0;
-    viewer.data.set_colors(C);
-    return true;
-  }
-
-  return false;
-}
+#include <iostream>
 
 int main(int argc, char *argv[])
 {
-
-  using namespace Eigen;
-  using namespace std;
+  // Mesh with per-face color
+  Eigen::MatrixXd V, C;
+  Eigen::MatrixXi F;
   // Load a mesh in OFF format
   igl::readOFF(TUTORIAL_SHARED_PATH "/fertility.off", V, F);
-
-  // Create a BVH for raycasting
-  ei = new igl::embree::EmbreeIntersector();
-  ei->init(V.cast<float>(),F);
-
-  // Initialize the colors to white for all vertices
-  C = MatrixXd::Constant(V.rows(),3,1);
-
-  // Show mesh
+  // Initialize white
+  C = Eigen::MatrixXd::Constant(F.rows(),3,1);
   igl::viewer::Viewer viewer;
+  viewer.callback_mouse_down = 
+    [&V,&F,&C](igl::viewer::Viewer& viewer, int, int)->bool
+  {
+    int fid;
+    Eigen::Vector3f bc;
+    // Cast a ray in the view direction starting from the mouse position
+    double x = viewer.current_mouse_x;
+    double y = viewer.core.viewport(3) - viewer.current_mouse_y;
+    if(igl::unproject_onto_mesh(
+      Eigen::Vector2f(x,y),
+      viewer.core.view * viewer.core.model,
+      viewer.core.proj,
+      viewer.core.viewport,
+      V,
+      F,
+      fid,
+      bc))
+    {
+      // paint hit red
+      C.row(fid)<<1,0,0;
+      viewer.data.set_colors(C);
+      return true;
+    }
+    return false;
+  };
+  std::cout<<R"(Usage:
+  [click]  Pick face on shape
+
+)";
+  // Show mesh
   viewer.data.set_mesh(V, F);
-  viewer.callback_mouse_down = &mouse_down;
   viewer.data.set_colors(C);
   viewer.core.show_lines = false;
   viewer.launch();