Merge branch 'master' of https://github.com/libigl/libigl

Former-commit-id: 31355bf27bb5dfb91afd001704da204bca53dc68

.gitignore

@@ -93,3 +93,4 @@ tests/bin
 python/build3
 *.pyc
 python/build4
+python/scripts/generated

.travis.yml

@@ -11,13 +11,13 @@ matrix:
         - cd python
         - mkdir build
         - cd build
-        - cmake -DCMAKE_CXX_COMPILER=g++-4.8 -DCMAKE_C_COMPILER=gcc-4.8 -DLIBIGL_WITH_EMBREE=OFF ../
+        - cmake -DCMAKE_CXX_COMPILER=g++-4.8 -DCMAKE_C_COMPILER=gcc-4.8 ../
         - make -j 2
         - cd ../../
         - cd tutorial
         - mkdir build
         - cd build
-        - cmake -DLIBIGL_USE_STATIC_LIBRARY=ON -DCMAKE_CXX_COMPILER=g++-4.8 -DCMAKE_C_COMPILER=gcc-4.8 -DLIBIGL_WITH_EMBREE=OFF ../
+        - cmake -DLIBIGL_USE_STATIC_LIBRARY=ON -DCMAKE_CXX_COMPILER=g++-4.8 -DCMAKE_C_COMPILER=gcc-4.8 ../
         - make -j 2
       addons:
         apt:

README.md

@@ -244,7 +244,7 @@ page](https://github.com/libigl/libigl/issues).
 
 ## Copyright
 2016 Alec Jacobson, Daniele Panozzo, Christian Schüller, Olga Diamanti, Qingnan
-Zhou, Nico Pietroni, Stefan Brugger, Kenshi Takayama, Wenzel Jakob, Nikolas De
+Zhou, Sebastian Koch, Nico Pietroni, Stefan Brugger, Kenshi Takayama, Wenzel Jakob, Nikolas De
 Giorgis, Luigi Rocca, Leonardo Sacht, Kevin Walliman, Olga Sorkine-Hornung, and others.
 
 Please see individual files for appropriate copyright notices.

include/igl/AABB.cpp

@@ -0,0 +1,953 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2013 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 "AABB.h"
+#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"
+#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>
+inline void igl::AABB<DerivedV,DIM>::init(
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::MatrixXi & Ele, 
+    const Eigen::PlainObjectBase<Derivedbb_mins> & bb_mins,
+    const Eigen::PlainObjectBase<Derivedbb_maxs> & bb_maxs,
+    const Eigen::VectorXi & elements,
+    const int i)
+{
+  using namespace std;
+  using namespace Eigen;
+  deinit();
+  if(bb_mins.size() > 0)
+  {
+    assert(bb_mins.rows() == bb_maxs.rows() && "Serial tree arrays must match");
+    assert(bb_mins.cols() == V.cols() && "Serial tree array dim must match V");
+    assert(bb_mins.cols() == bb_maxs.cols() && "Serial tree arrays must match");
+    assert(bb_mins.rows() == elements.rows() &&
+        "Serial tree arrays must match");
+    // construct from serialization
+    m_box.extend(bb_mins.row(i).transpose());
+    m_box.extend(bb_maxs.row(i).transpose());
+    m_primitive = elements(i);
+    // Not leaf then recurse
+    if(m_primitive == -1)
+    {
+      m_left = new AABB();
+      m_left->init( V,Ele,bb_mins,bb_maxs,elements,2*i+1);
+      m_right = new AABB();
+      m_right->init( V,Ele,bb_mins,bb_maxs,elements,2*i+2);
+      //m_depth = std::max( m_left->m_depth, m_right->m_depth)+1;
+    }
+  }else
+  {
+    VectorXi allI = colon<int>(0,Ele.rows()-1);
+    MatrixXDIMS BC;
+    if(Ele.cols() == 1)
+    {
+      // points
+      BC = V;
+    }else
+    {
+      // Simplices
+      barycenter(V,Ele,BC);
+    }
+    MatrixXi SI(BC.rows(),BC.cols());
+    {
+      MatrixXDIMS _;
+      MatrixXi IS;
+      igl::sort(BC,1,true,_,IS);
+      // Need SI(i) to tell which place i would be sorted into
+      const int dim = IS.cols();
+      for(int i = 0;i<IS.rows();i++)
+      {
+        for(int d = 0;d<dim;d++)
+        {
+          SI(IS(i,d),d) = i;
+        }
+      }
+    }
+    init(V,Ele,SI,allI);
+  }
+}
+
+  template <typename DerivedV, int DIM>
+inline void igl::AABB<DerivedV,DIM>::init(
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::MatrixXi & Ele)
+{
+  using namespace Eigen;
+  // deinit will be immediately called...
+  return init(V,Ele,MatrixXDIMS(),MatrixXDIMS(),VectorXi(),0);
+}
+
+  template <typename DerivedV, int DIM>
+inline void igl::AABB<DerivedV,DIM>::init(
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::MatrixXi & Ele, 
+    const Eigen::MatrixXi & SI,
+    const Eigen::VectorXi & I)
+{
+  using namespace Eigen;
+  using namespace std;
+  deinit();
+  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>();
+  // Compute bounding box
+  for(int i = 0;i<I.rows();i++)
+  {
+    for(int c = 0;c<Ele.cols();c++)
+    {
+      m_box.extend(V.row(Ele(I(i),c)).transpose());
+      m_box.extend(V.row(Ele(I(i),c)).transpose());
+    }
+  }
+  switch(I.size())
+  {
+    case 0:
+      {
+        assert(false);
+      }
+    case 1:
+      {
+        m_primitive = I(0);
+        break;
+      }
+    default:
+      {
+        // Compute longest direction
+        int max_d = -1;
+        m_box.diagonal().maxCoeff(&max_d);
+        // Can't use median on BC directly because many may have same value,
+        // but can use median on sorted BC indices
+        VectorXi SIdI(I.rows());
+        for(int i = 0;i<I.rows();i++)
+        {
+          SIdI(i) = SI(I(i),max_d);
+        }
+        // Since later I use <= I think I don't need to worry about odd/even
+        // Pass by copy to avoid changing input
+        const auto median = [](VectorXi A)->Scalar
+        {
+          size_t n = A.size()/2;
+          nth_element(A.data(),A.data()+n,A.data()+A.size());
+          if(A.rows() % 2 == 1)
+          {
+            return A(n);
+          }else
+          {
+            nth_element(A.data(),A.data()+n-1,A.data()+A.size());
+            return 0.5*(A(n)+A(n-1));
+          }
+        };
+        const Scalar med = median(SIdI);
+        VectorXi LI((I.rows()+1)/2),RI(I.rows()/2);
+        assert(LI.rows()+RI.rows() == I.rows());
+        // Distribute left and right
+        {
+          int li = 0;
+          int ri = 0;
+          for(int i = 0;i<I.rows();i++)
+          {
+            if(SIdI(i)<=med)
+            {
+              LI(li++) = I(i);
+            }else
+            {
+              RI(ri++) = I(i);
+            }
+          }
+        }
+        //m_depth = 0;
+        if(LI.rows()>0)
+        {
+          m_left = new AABB();
+          m_left->init(V,Ele,SI,LI);
+          //m_depth = std::max(m_depth, m_left->m_depth+1);
+        }
+        if(RI.rows()>0)
+        {
+          m_right = new AABB();
+          m_right->init(V,Ele,SI,RI);
+          //m_depth = std::max(m_depth, m_right->m_depth+1);
+        }
+      }
+  }
+}
+
+template <typename DerivedV, int DIM>
+inline bool igl::AABB<DerivedV,DIM>::is_leaf() const
+{
+  return m_primitive != -1;
+}
+
+template <typename DerivedV, int DIM>
+template <typename Derivedq>
+inline std::vector<int> igl::AABB<DerivedV,DIM>::find(
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::MatrixXi & Ele, 
+    const Eigen::PlainObjectBase<Derivedq> & q,
+    const bool first) const
+{
+  using namespace std;
+  using namespace Eigen;
+  assert(q.size() == DIM && 
+      "Query dimension should match aabb dimension");
+  assert(Ele.cols() == V.cols()+1 && 
+      "AABB::find only makes sense for (d+1)-simplices");
+  const Scalar epsilon = igl::EPS<Scalar>();
+  // Check if outside bounding box
+  bool inside = m_box.contains(q.transpose());
+  if(!inside)
+  {
+    return std::vector<int>();
+  }
+  assert(m_primitive==-1 || (m_left == NULL && m_right == NULL));
+  if(is_leaf())
+  {
+    // Initialize to some value > -epsilon
+    Scalar a1=0,a2=0,a3=0,a4=0;
+    switch(DIM)
+    {
+      case 3:
+        {
+          // Barycentric coordinates
+          typedef Eigen::Matrix<Scalar,1,3> RowVector3S;
+          const RowVector3S V1 = V.row(Ele(m_primitive,0));
+          const RowVector3S V2 = V.row(Ele(m_primitive,1));
+          const RowVector3S V3 = V.row(Ele(m_primitive,2));
+          const RowVector3S V4 = V.row(Ele(m_primitive,3));
+          a1 = volume_single(V2,V4,V3,(RowVector3S)q);
+          a2 = volume_single(V1,V3,V4,(RowVector3S)q);
+          a3 = volume_single(V1,V4,V2,(RowVector3S)q);
+          a4 = volume_single(V1,V2,V3,(RowVector3S)q);
+          break;
+        }
+      case 2:
+        {
+          // Barycentric coordinates
+          typedef Eigen::Matrix<Scalar,2,1> Vector2S;
+          const Vector2S V1 = V.row(Ele(m_primitive,0));
+          const Vector2S V2 = V.row(Ele(m_primitive,1));
+          const Vector2S V3 = V.row(Ele(m_primitive,2));
+          // Hack for now to keep templates simple. If becomes bottleneck
+          // consider using std::enable_if_t 
+          const Vector2S q2 = q.head(2);
+          a1 = doublearea_single(V1,V2,q2);
+          a2 = doublearea_single(V2,V3,q2);
+          a3 = doublearea_single(V3,V1,q2);
+          break;
+        }
+      default:assert(false);
+    }
+    // Normalization is important for correcting sign
+    Scalar sum = a1+a2+a3+a4;
+    a1 /= sum;
+    a2 /= sum;
+    a3 /= sum;
+    a4 /= sum;
+    if(
+        a1>=-epsilon && 
+        a2>=-epsilon && 
+        a3>=-epsilon && 
+        a4>=-epsilon)
+    {
+      return std::vector<int>(1,m_primitive);
+    }else
+    {
+      return std::vector<int>();
+    }
+  }
+  std::vector<int> left = m_left->find(V,Ele,q,first);
+  if(first && !left.empty())
+  {
+    return left;
+  }
+  std::vector<int> right = m_right->find(V,Ele,q,first);
+  if(first)
+  {
+    return right;
+  }
+  left.insert(left.end(),right.begin(),right.end());
+  return left;
+}
+
+template <typename DerivedV, int DIM>
+inline int igl::AABB<DerivedV,DIM>::subtree_size() const
+{
+  // 1 for self
+  int n = 1;
+  int n_left = 0,n_right = 0;
+  if(m_left != NULL)
+  {
+    n_left = m_left->subtree_size();
+  }
+  if(m_right != NULL)
+  {
+    n_right = m_right->subtree_size();
+  }
+  n += 2*std::max(n_left,n_right);
+  return n;
+}
+
+
+template <typename DerivedV, int DIM>
+template <typename Derivedbb_mins, typename Derivedbb_maxs>
+inline void igl::AABB<DerivedV,DIM>::serialize(
+    Eigen::PlainObjectBase<Derivedbb_mins> & bb_mins,
+    Eigen::PlainObjectBase<Derivedbb_maxs> & bb_maxs,
+    Eigen::VectorXi & elements,
+    const int i) const
+{
+  using namespace std;
+  using namespace Eigen;
+  // Calling for root then resize output
+  if(i==0)
+  {
+    const int m = subtree_size();
+    //cout<<"m: "<<m<<endl;
+    bb_mins.resize(m,DIM);
+    bb_maxs.resize(m,DIM);
+    elements.resize(m,1);
+  }
+  //cout<<i<<" ";
+  bb_mins.row(i) = m_box.min();
+  bb_maxs.row(i) = m_box.max();
+  elements(i) = m_primitive;
+  if(m_left != NULL)
+  {
+    m_left->serialize(bb_mins,bb_maxs,elements,2*i+1);
+  }
+  if(m_right != NULL)
+  {
+    m_right->serialize(bb_mins,bb_maxs,elements,2*i+2);
+  }
+}
+
+template <typename DerivedV, int DIM>
+inline typename igl::AABB<DerivedV,DIM>::Scalar 
+igl::AABB<DerivedV,DIM>::squared_distance(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::MatrixXi & Ele, 
+  const RowVectorDIMS & p,
+  int & i,
+  RowVectorDIMS & c) const
+{
+  return squared_distance(V,Ele,p,std::numeric_limits<Scalar>::infinity(),i,c);
+}
+
+
+template <typename DerivedV, int DIM>
+inline typename igl::AABB<DerivedV,DIM>::Scalar 
+igl::AABB<DerivedV,DIM>::squared_distance(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::MatrixXi & Ele, 
+  const RowVectorDIMS & p,
+  Scalar min_sqr_d,
+  int & i,
+  RowVectorDIMS & c) const
+{
+  using namespace Eigen;
+  using namespace std;
+  Scalar sqr_d = min_sqr_d;
+  //assert(DIM == 3 && "Code has only been tested for DIM == 3");
+  assert((Ele.cols() == 3 || Ele.cols() == 2 || Ele.cols() == 1)
+    && "Code has only been tested for simplex sizes 3,2,1");
+
+  assert(m_primitive==-1 || (m_left == NULL && m_right == NULL));
+  if(is_leaf())
+  {
+    leaf_squared_distance(V,Ele,p,sqr_d,i,c);
+  }else
+  {
+    bool looked_left = false;
+    bool looked_right = false;
+    const auto & look_left = [&]()
+    {
+      int i_left;
+      RowVectorDIMS c_left = c;
+      Scalar sqr_d_left = m_left->squared_distance(V,Ele,p,sqr_d,i_left,c_left);
+      set_min(p,sqr_d_left,i_left,c_left,sqr_d,i,c);
+      looked_left = true;
+    };
+    const auto & look_right = [&]()
+    {
+      int i_right;
+      RowVectorDIMS c_right = c;
+      Scalar sqr_d_right = m_right->squared_distance(V,Ele,p,sqr_d,i_right,c_right);
+      set_min(p,sqr_d_right,i_right,c_right,sqr_d,i,c);
+      looked_right = true;
+    };
+
+    // must look left or right if in box
+    if(m_left->m_box.contains(p.transpose()))
+    {
+      look_left();
+    }
+    if(m_right->m_box.contains(p.transpose()))
+    {
+      look_right();
+    }
+    // if haven't looked left and could be less than current min, then look
+    Scalar  left_min_sqr_d = m_left->m_box.squaredExteriorDistance(p.transpose());
+    Scalar right_min_sqr_d = m_right->m_box.squaredExteriorDistance(p.transpose());
+    if(left_min_sqr_d < right_min_sqr_d)
+    {
+      if(!looked_left && left_min_sqr_d<sqr_d)
+      {
+        look_left();
+      }
+      if( !looked_right && right_min_sqr_d<sqr_d)
+      {
+        look_right();
+      }
+    }else
+    {
+      if( !looked_right && right_min_sqr_d<sqr_d)
+      {
+        look_right();
+      }
+      if(!looked_left && left_min_sqr_d<sqr_d)
+      {
+        look_left();
+      }
+    }
+  }
+  return sqr_d;
+}
+
+template <typename DerivedV, int DIM>
+template <
+  typename DerivedP, 
+  typename DerivedsqrD, 
+  typename DerivedI, 
+  typename DerivedC>
+inline void igl::AABB<DerivedV,DIM>::squared_distance(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::MatrixXi & Ele, 
+  const Eigen::PlainObjectBase<DerivedP> & P,
+  Eigen::PlainObjectBase<DerivedsqrD> & sqrD,
+  Eigen::PlainObjectBase<DerivedI> & I,
+  Eigen::PlainObjectBase<DerivedC> & C) const
+{
+  assert(P.cols() == V.cols() && "cols in P should match dim of cols in V");
+  sqrD.resize(P.rows(),1);
+  I.resize(P.rows(),1);
+  C.resize(P.rows(),P.cols());
+  for(int p = 0;p<P.rows();p++)
+  {
+    RowVectorDIMS Pp = P.row(p), c;
+    int Ip;
+    sqrD(p) = squared_distance(V,Ele,Pp,Ip,c);
+    I(p) = Ip;
+    C.row(p).head(DIM) = c;
+  }
+}
+
+template <typename DerivedV, int DIM>
+template < 
+  typename Derivedother_V,
+  typename DerivedsqrD, 
+  typename DerivedI, 
+  typename DerivedC>
+inline void igl::AABB<DerivedV,DIM>::squared_distance(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::MatrixXi & Ele, 
+  const AABB<Derivedother_V,DIM> & other,
+  const Eigen::PlainObjectBase<Derivedother_V> & other_V,
+  const Eigen::MatrixXi & other_Ele, 
+  Eigen::PlainObjectBase<DerivedsqrD> & sqrD,
+  Eigen::PlainObjectBase<DerivedI> & I,
+  Eigen::PlainObjectBase<DerivedC> & C) const
+{
+  assert(other_Ele.cols() == 1 && 
+    "Only implemented for other as list of points");
+  assert(other_V.cols() == V.cols() && "other must match this dimension");
+  sqrD.setConstant(other_Ele.rows(),1,std::numeric_limits<double>::infinity());
+  I.resize(other_Ele.rows(),1);
+  C.resize(other_Ele.rows(),other_V.cols());
+  // All points in other_V currently think they need to check against root of
+  // this. The point of using another AABB is to quickly prune chunks of
+  // other_V so that most points just check some subtree of this.
+
+  // This holds a conservative estimate of max(sqr_D) where sqr_D is the
+  // current best minimum squared distance for all points in this subtree
+  double min_sqr_d = std::numeric_limits<double>::infinity();
+  squared_distance_helper(
+    V,Ele,&other,other_V,other_Ele,min_sqr_d,sqrD,I,C);
+}
+
+template <typename DerivedV, int DIM>
+template < 
+  typename Derivedother_V,
+  typename DerivedsqrD, 
+  typename DerivedI, 
+  typename DerivedC>
+inline typename igl::AABB<DerivedV,DIM>::Scalar igl::AABB<DerivedV,DIM>::squared_distance_helper(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::MatrixXi & Ele, 
+  const AABB<Derivedother_V,DIM> * other,
+  const Eigen::PlainObjectBase<Derivedother_V> & other_V,
+  const Eigen::MatrixXi & other_Ele, 
+  const Scalar /*min_sqr_d*/,
+  Eigen::PlainObjectBase<DerivedsqrD> & sqrD,
+  Eigen::PlainObjectBase<DerivedI> & I,
+  Eigen::PlainObjectBase<DerivedC> & C) const
+{
+  using namespace std;
+  using namespace Eigen;
+
+  // This implementation is a bit disappointing. There's no major speed up. Any
+  // performance gains seem to come from accidental cache coherency and
+  // diminish for larger "other" (the opposite of what was intended).
+
+  // Base case
+  if(other->is_leaf() && this->is_leaf())
+  {
+    Scalar sqr_d = sqrD(other->m_primitive);
+    int i = I(other->m_primitive);
+    RowVectorDIMS c = C.row(      other->m_primitive);
+    RowVectorDIMS p = other_V.row(other->m_primitive);
+    leaf_squared_distance(V,Ele,p,sqr_d,i,c);
+    sqrD( other->m_primitive) = sqr_d;
+    I(    other->m_primitive) = i;
+    C.row(other->m_primitive) = c;
+    //cout<<"leaf: "<<sqr_d<<endl;
+    //other->m_max_sqr_d = sqr_d;
+    return sqr_d;
+  }
+
+  if(other->is_leaf())
+  {
+    Scalar sqr_d = sqrD(other->m_primitive);
+    int i = I(other->m_primitive);
+    RowVectorDIMS c = C.row(      other->m_primitive);
+    RowVectorDIMS p = other_V.row(other->m_primitive);
+    sqr_d = squared_distance(V,Ele,p,sqr_d,i,c);
+    sqrD( other->m_primitive) = sqr_d;
+    I(    other->m_primitive) = i;
+    C.row(other->m_primitive) = c;
+    //other->m_max_sqr_d = sqr_d;
+    return sqr_d;
+  }
+
+  //// Exact minimum squared distance between arbitary primitives inside this and
+  //// othre's bounding boxes
+  //const auto & min_squared_distance = [&](
+  //  const AABB<DerivedV,DIM> * A,
+  //  const AABB<Derivedother_V,DIM> * B)->Scalar
+  //{
+  //  return A->m_box.squaredExteriorDistance(B->m_box);
+  //};
+
+  if(this->is_leaf())
+  {
+    //if(min_squared_distance(this,other) < other->m_max_sqr_d)
+    if(true)
+    {
+      this->squared_distance_helper(
+        V,Ele,other->m_left,other_V,other_Ele,0,sqrD,I,C);
+      this->squared_distance_helper(
+        V,Ele,other->m_right,other_V,other_Ele,0,sqrD,I,C);
+    }else
+    {
+      // This is never reached...
+    }
+    //// we know other is not a leaf
+    //other->m_max_sqr_d = std::max(other->m_left->m_max_sqr_d,other->m_right->m_max_sqr_d);
+    return 0;
+  }
+
+  // FORCE DOWN TO OTHER LEAF EVAL
+  //if(min_squared_distance(this,other) < other->m_max_sqr_d)
+  if(true)
+  {
+    if(true)
+    {
+      this->squared_distance_helper(
+        V,Ele,other->m_left,other_V,other_Ele,0,sqrD,I,C);
+      this->squared_distance_helper(
+        V,Ele,other->m_right,other_V,other_Ele,0,sqrD,I,C);
+    }else // this direction never seems to be faster
+    {
+      this->m_left->squared_distance_helper(
+        V,Ele,other,other_V,other_Ele,0,sqrD,I,C);
+      this->m_right->squared_distance_helper(
+        V,Ele,other,other_V,other_Ele,0,sqrD,I,C);
+    }
+  }else
+  {
+    // this is never reached ... :-(
+  }
+  //// we know other is not a leaf
+  //other->m_max_sqr_d = std::max(other->m_left->m_max_sqr_d,other->m_right->m_max_sqr_d);
+
+  return 0;
+#if 0 // False
+
+  // _Very_ conservative approximation of maximum squared distance between
+  // primitives inside this and other's bounding boxes
+  const auto & max_squared_distance = [](
+    const AABB<DerivedV,DIM> * A,
+    const AABB<Derivedother_V,DIM> * B)->Scalar
+  {
+    AlignedBox<Scalar,DIM> combo = A->m_box;
+    combo.extend(B->m_box);
+    return combo.diagonal().squaredNorm();
+  };
+
+  //// other base-case
+  //if(other->is_leaf())
+  //{
+  //  double sqr_d = sqrD(other->m_primitive);
+  //  int i = I(other->m_primitive);
+  //  RowVectorDIMS c = C.row(m_primitive);
+  //  RowVectorDIMS p = other_V.row(m_primitive);
+  //  leaf_squared_distance(V,Ele,p,sqr_d,i,c);
+  //  sqrD(other->m_primitive) = sqr_d;
+  //  I(other->m_primitive) = i;
+  //  C.row(m_primitive) = c;
+  //  return;
+  //}
+  std::vector<const AABB<DerivedV,DIM> * > this_list;
+  if(this->is_leaf())
+  {
+    this_list.push_back(this);
+  }else
+  {
+    assert(this->m_left);
+    this_list.push_back(this->m_left);
+    assert(this->m_right);
+    this_list.push_back(this->m_right);
+  }
+  std::vector<AABB<Derivedother_V,DIM> *> other_list;
+  if(other->is_leaf())
+  {
+    other_list.push_back(other);
+  }else
+  {
+    assert(other->m_left);
+    other_list.push_back(other->m_left);
+    assert(other->m_right);
+    other_list.push_back(other->m_right);
+  }
+
+  //const std::function<Scalar(
+  //  const AABB<Derivedother_V,DIM> * other)
+  //    > max_sqr_d = [&sqrD,&max_sqr_d](const AABB<Derivedother_V,DIM> * other)->Scalar
+  //  {
+  //    if(other->is_leaf())
+  //    {
+  //      return sqrD(other->m_primitive);
+  //    }else
+  //    {
+  //      return std::max(max_sqr_d(other->m_left),max_sqr_d(other->m_right));
+  //    }
+  //  };
+
+  //// Potentially recurse on all pairs, if minimum distance is less than running
+  //// bound
+  //Eigen::Matrix<Scalar,Eigen::Dynamic,1> other_max_sqr_d =
+  //  Eigen::Matrix<Scalar,Eigen::Dynamic,1>::Constant(other_list.size(),1,min_sqr_d);
+  for(size_t child = 0;child<other_list.size();child++)
+  {
+    auto other_tree = other_list[child];
+
+    Eigen::Matrix<Scalar,Eigen::Dynamic,1> this_max_sqr_d(this_list.size(),1);
+    for(size_t t = 0;t<this_list.size();t++)
+    {
+      const auto this_tree = this_list[t];
+      this_max_sqr_d(t) = max_squared_distance(this_tree,other_tree);
+    }
+    if(this_list.size() ==2 &&
+      ( this_max_sqr_d(0) > this_max_sqr_d(1))
+      )
+    {
+      std::swap(this_list[0],this_list[1]);
+      //std::swap(this_max_sqr_d(0),this_max_sqr_d(1));
+    }
+    const Scalar sqr_d = this_max_sqr_d.minCoeff();
+
+
+    for(size_t t = 0;t<this_list.size();t++)
+    {
+      const auto this_tree = this_list[t];
+
+      //const auto mm = max_sqr_d(other_tree);
+      //const Scalar mc = other_max_sqr_d(child);
+      //assert(mc == mm);
+      // Only look left/right in this_list if can possible decrease somebody's
+      // distance in this_tree.
+      const Scalar min_this_other = min_squared_distance(this_tree,other_tree); 
+      if(
+          min_this_other < sqr_d && 
+          min_this_other < other_tree->m_max_sqr_d)
+      {
+        //cout<<"before: "<<other_max_sqr_d(child)<<endl;
+        //other_max_sqr_d(child) = std::min(
+        //  other_max_sqr_d(child),
+        //  this_tree->squared_distance_helper(
+        //    V,Ele,other_tree,other_V,other_Ele,other_max_sqr_d(child),sqrD,I,C));
+        //cout<<"after: "<<other_max_sqr_d(child)<<endl;
+          this_tree->squared_distance_helper(
+            V,Ele,other_tree,other_V,other_Ele,0,sqrD,I,C);
+      }
+    }
+  }
+  //const Scalar ret = other_max_sqr_d.maxCoeff();
+  //const auto mm = max_sqr_d(other);
+  //assert(mm == ret);
+  //cout<<"non-leaf: "<<ret<<endl;
+  //return ret;
+  if(!other->is_leaf())
+  {
+    other->m_max_sqr_d = std::max(other->m_left->m_max_sqr_d,other->m_right->m_max_sqr_d);
+  }
+  return 0;
+#endif
+}
+
+template <typename DerivedV, int DIM>
+inline void igl::AABB<DerivedV,DIM>::leaf_squared_distance(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::MatrixXi & Ele, 
+  const RowVectorDIMS & p,
+  Scalar & sqr_d,
+  int & i,
+  RowVectorDIMS & c) const
+{
+  using namespace Eigen;
+  using namespace std;
+  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);
+}
+
+
+template <typename DerivedV, int DIM>
+inline void igl::AABB<DerivedV,DIM>::set_min(
+  const RowVectorDIMS & 
+#ifndef NDEBUG
+  p
+#endif
+  ,
+  const Scalar sqr_d_candidate,
+  const int i_candidate,
+  const RowVectorDIMS & c_candidate,
+  Scalar & sqr_d,
+  int & i,
+  RowVectorDIMS & c) const
+{
+#ifndef NDEBUG
+  //std::cout<<matlab_format(c_candidate,"c_candidate")<<std::endl;
+  const Scalar pc_norm = (p-c_candidate).squaredNorm();
+  const Scalar diff = fabs(sqr_d_candidate - pc_norm);
+  assert(diff<=1e-10 && "distance should match norm of difference");
+#endif
+  if(sqr_d_candidate < sqr_d)
+  {
+    i = i_candidate;
+    c = c_candidate;
+    sqr_d = sqr_d_candidate;
+  }
+}
+
+
+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),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),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),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;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+template void igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 3>::init(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&);
+template void igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 2>::init(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&);
+template void igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 3>::squared_distance<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -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::Matrix<int, -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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&) const;
+template void igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 2>::squared_distance<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -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::Matrix<int, -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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&) const;
+template double igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 3>::squared_distance(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, 1, 3, 1, 1, 3> const&, int&, Eigen::Matrix<double, 1, 3, 1, 1, 3>&) const;
+template double igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 2>::squared_distance(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, 1, 2, 1, 1, 2> const&, int&, Eigen::Matrix<double, 1, 2, 1, 1, 2>&) const;
+template void igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 2>::squared_distance<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -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::Matrix<int, -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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&) const;
+template void igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 3>::squared_distance<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -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::Matrix<int, -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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&) const;
+
+
+#endif
+
+

include/igl/AABB.h

@@ -289,935 +289,9 @@ 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"
-#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>
-inline void igl::AABB<DerivedV,DIM>::init(
-    const Eigen::PlainObjectBase<DerivedV> & V,
-    const Eigen::MatrixXi & Ele, 
-    const Eigen::PlainObjectBase<Derivedbb_mins> & bb_mins,
-    const Eigen::PlainObjectBase<Derivedbb_maxs> & bb_maxs,
-    const Eigen::VectorXi & elements,
-    const int i)
-{
-  using namespace std;
-  using namespace Eigen;
-  deinit();
-  if(bb_mins.size() > 0)
-  {
-    assert(bb_mins.rows() == bb_maxs.rows() && "Serial tree arrays must match");
-    assert(bb_mins.cols() == V.cols() && "Serial tree array dim must match V");
-    assert(bb_mins.cols() == bb_maxs.cols() && "Serial tree arrays must match");
-    assert(bb_mins.rows() == elements.rows() &&
-        "Serial tree arrays must match");
-    // construct from serialization
-    m_box.extend(bb_mins.row(i).transpose());
-    m_box.extend(bb_maxs.row(i).transpose());
-    m_primitive = elements(i);
-    // Not leaf then recurse
-    if(m_primitive == -1)
-    {
-      m_left = new AABB();
-      m_left->init( V,Ele,bb_mins,bb_maxs,elements,2*i+1);
-      m_right = new AABB();
-      m_right->init( V,Ele,bb_mins,bb_maxs,elements,2*i+2);
-      //m_depth = std::max( m_left->m_depth, m_right->m_depth)+1;
-    }
-  }else
-  {
-    VectorXi allI = colon<int>(0,Ele.rows()-1);
-    MatrixXDIMS BC;
-    if(Ele.cols() == 1)
-    {
-      // points
-      BC = V;
-    }else
-    {
-      // Simplices
-      barycenter(V,Ele,BC);
-    }
-    MatrixXi SI(BC.rows(),BC.cols());
-    {
-      MatrixXDIMS _;
-      MatrixXi IS;
-      igl::sort(BC,1,true,_,IS);
-      // Need SI(i) to tell which place i would be sorted into
-      const int dim = IS.cols();
-      for(int i = 0;i<IS.rows();i++)
-      {
-        for(int d = 0;d<dim;d++)
-        {
-          SI(IS(i,d),d) = i;
-        }
-      }
-    }
-    init(V,Ele,SI,allI);
-  }
-}
-
-  template <typename DerivedV, int DIM>
-inline void igl::AABB<DerivedV,DIM>::init(
-    const Eigen::PlainObjectBase<DerivedV> & V,
-    const Eigen::MatrixXi & Ele)
-{
-  using namespace Eigen;
-  // deinit will be immediately called...
-  return init(V,Ele,MatrixXDIMS(),MatrixXDIMS(),VectorXi(),0);
-}
-
-  template <typename DerivedV, int DIM>
-inline void igl::AABB<DerivedV,DIM>::init(
-    const Eigen::PlainObjectBase<DerivedV> & V,
-    const Eigen::MatrixXi & Ele, 
-    const Eigen::MatrixXi & SI,
-    const Eigen::VectorXi & I)
-{
-  using namespace Eigen;
-  using namespace std;
-  deinit();
-  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>();
-  // Compute bounding box
-  for(int i = 0;i<I.rows();i++)
-  {
-    for(int c = 0;c<Ele.cols();c++)
-    {
-      m_box.extend(V.row(Ele(I(i),c)).transpose());
-      m_box.extend(V.row(Ele(I(i),c)).transpose());
-    }
-  }
-  switch(I.size())
-  {
-    case 0:
-      {
-        assert(false);
-      }
-    case 1:
-      {
-        m_primitive = I(0);
-        break;
-      }
-    default:
-      {
-        // Compute longest direction
-        int max_d = -1;
-        m_box.diagonal().maxCoeff(&max_d);
-        // Can't use median on BC directly because many may have same value,
-        // but can use median on sorted BC indices
-        VectorXi SIdI(I.rows());
-        for(int i = 0;i<I.rows();i++)
-        {
-          SIdI(i) = SI(I(i),max_d);
-        }
-        // Since later I use <= I think I don't need to worry about odd/even
-        // Pass by copy to avoid changing input
-        const auto median = [](VectorXi A)->Scalar
-        {
-          size_t n = A.size()/2;
-          nth_element(A.data(),A.data()+n,A.data()+A.size());
-          if(A.rows() % 2 == 1)
-          {
-            return A(n);
-          }else
-          {
-            nth_element(A.data(),A.data()+n-1,A.data()+A.size());
-            return 0.5*(A(n)+A(n-1));
-          }
-        };
-        const Scalar med = median(SIdI);
-        VectorXi LI((I.rows()+1)/2),RI(I.rows()/2);
-        assert(LI.rows()+RI.rows() == I.rows());
-        // Distribute left and right
-        {
-          int li = 0;
-          int ri = 0;
-          for(int i = 0;i<I.rows();i++)
-          {
-            if(SIdI(i)<=med)
-            {
-              LI(li++) = I(i);
-            }else
-            {
-              RI(ri++) = I(i);
-            }
-          }
-        }
-        //m_depth = 0;
-        if(LI.rows()>0)
-        {
-          m_left = new AABB();
-          m_left->init(V,Ele,SI,LI);
-          //m_depth = std::max(m_depth, m_left->m_depth+1);
-        }
-        if(RI.rows()>0)
-        {
-          m_right = new AABB();
-          m_right->init(V,Ele,SI,RI);
-          //m_depth = std::max(m_depth, m_right->m_depth+1);
-        }
-      }
-  }
-}
-
-template <typename DerivedV, int DIM>
-inline bool igl::AABB<DerivedV,DIM>::is_leaf() const
-{
-  return m_primitive != -1;
-}
-
-template <typename DerivedV, int DIM>
-template <typename Derivedq>
-inline std::vector<int> igl::AABB<DerivedV,DIM>::find(
-    const Eigen::PlainObjectBase<DerivedV> & V,
-    const Eigen::MatrixXi & Ele, 
-    const Eigen::PlainObjectBase<Derivedq> & q,
-    const bool first) const
-{
-  using namespace std;
-  using namespace Eigen;
-  assert(q.size() == DIM && 
-      "Query dimension should match aabb dimension");
-  assert(Ele.cols() == V.cols()+1 && 
-      "AABB::find only makes sense for (d+1)-simplices");
-  const Scalar epsilon = igl::EPS<Scalar>();
-  // Check if outside bounding box
-  bool inside = m_box.contains(q.transpose());
-  if(!inside)
-  {
-    return std::vector<int>();
-  }
-  assert(m_primitive==-1 || (m_left == NULL && m_right == NULL));
-  if(is_leaf())
-  {
-    // Initialize to some value > -epsilon
-    Scalar a1=0,a2=0,a3=0,a4=0;
-    switch(DIM)
-    {
-      case 3:
-        {
-          // Barycentric coordinates
-          typedef Eigen::Matrix<Scalar,1,3> RowVector3S;
-          const RowVector3S V1 = V.row(Ele(m_primitive,0));
-          const RowVector3S V2 = V.row(Ele(m_primitive,1));
-          const RowVector3S V3 = V.row(Ele(m_primitive,2));
-          const RowVector3S V4 = V.row(Ele(m_primitive,3));
-          a1 = volume_single(V2,V4,V3,(RowVector3S)q);
-          a2 = volume_single(V1,V3,V4,(RowVector3S)q);
-          a3 = volume_single(V1,V4,V2,(RowVector3S)q);
-          a4 = volume_single(V1,V2,V3,(RowVector3S)q);
-          break;
-        }
-      case 2:
-        {
-          // Barycentric coordinates
-          typedef Eigen::Matrix<Scalar,2,1> Vector2S;
-          const Vector2S V1 = V.row(Ele(m_primitive,0));
-          const Vector2S V2 = V.row(Ele(m_primitive,1));
-          const Vector2S V3 = V.row(Ele(m_primitive,2));
-          // Hack for now to keep templates simple. If becomes bottleneck
-          // consider using std::enable_if_t 
-          const Vector2S q2 = q.head(2);
-          a1 = doublearea_single(V1,V2,q2);
-          a2 = doublearea_single(V2,V3,q2);
-          a3 = doublearea_single(V3,V1,q2);
-          break;
-        }
-      default:assert(false);
-    }
-    // Normalization is important for correcting sign
-    Scalar sum = a1+a2+a3+a4;
-    a1 /= sum;
-    a2 /= sum;
-    a3 /= sum;
-    a4 /= sum;
-    if(
-        a1>=-epsilon && 
-        a2>=-epsilon && 
-        a3>=-epsilon && 
-        a4>=-epsilon)
-    {
-      return std::vector<int>(1,m_primitive);
-    }else
-    {
-      return std::vector<int>();
-    }
-  }
-  std::vector<int> left = m_left->find(V,Ele,q,first);
-  if(first && !left.empty())
-  {
-    return left;
-  }
-  std::vector<int> right = m_right->find(V,Ele,q,first);
-  if(first)
-  {
-    return right;
-  }
-  left.insert(left.end(),right.begin(),right.end());
-  return left;
-}
-
-template <typename DerivedV, int DIM>
-inline int igl::AABB<DerivedV,DIM>::subtree_size() const
-{
-  // 1 for self
-  int n = 1;
-  int n_left = 0,n_right = 0;
-  if(m_left != NULL)
-  {
-    n_left = m_left->subtree_size();
-  }
-  if(m_right != NULL)
-  {
-    n_right = m_right->subtree_size();
-  }
-  n += 2*std::max(n_left,n_right);
-  return n;
-}
-
-
-template <typename DerivedV, int DIM>
-template <typename Derivedbb_mins, typename Derivedbb_maxs>
-inline void igl::AABB<DerivedV,DIM>::serialize(
-    Eigen::PlainObjectBase<Derivedbb_mins> & bb_mins,
-    Eigen::PlainObjectBase<Derivedbb_maxs> & bb_maxs,
-    Eigen::VectorXi & elements,
-    const int i) const
-{
-  using namespace std;
-  using namespace Eigen;
-  // Calling for root then resize output
-  if(i==0)
-  {
-    const int m = subtree_size();
-    //cout<<"m: "<<m<<endl;
-    bb_mins.resize(m,DIM);
-    bb_maxs.resize(m,DIM);
-    elements.resize(m,1);
-  }
-  //cout<<i<<" ";
-  bb_mins.row(i) = m_box.min();
-  bb_maxs.row(i) = m_box.max();
-  elements(i) = m_primitive;
-  if(m_left != NULL)
-  {
-    m_left->serialize(bb_mins,bb_maxs,elements,2*i+1);
-  }
-  if(m_right != NULL)
-  {
-    m_right->serialize(bb_mins,bb_maxs,elements,2*i+2);
-  }
-}
-
-template <typename DerivedV, int DIM>
-inline typename igl::AABB<DerivedV,DIM>::Scalar 
-igl::AABB<DerivedV,DIM>::squared_distance(
-  const Eigen::PlainObjectBase<DerivedV> & V,
-  const Eigen::MatrixXi & Ele, 
-  const RowVectorDIMS & p,
-  int & i,
-  RowVectorDIMS & c) const
-{
-  return squared_distance(V,Ele,p,std::numeric_limits<Scalar>::infinity(),i,c);
-}
-
-
-template <typename DerivedV, int DIM>
-inline typename igl::AABB<DerivedV,DIM>::Scalar 
-igl::AABB<DerivedV,DIM>::squared_distance(
-  const Eigen::PlainObjectBase<DerivedV> & V,
-  const Eigen::MatrixXi & Ele, 
-  const RowVectorDIMS & p,
-  Scalar min_sqr_d,
-  int & i,
-  RowVectorDIMS & c) const
-{
-  using namespace Eigen;
-  using namespace std;
-  Scalar sqr_d = min_sqr_d;
-  //assert(DIM == 3 && "Code has only been tested for DIM == 3");
-  assert((Ele.cols() == 3 || Ele.cols() == 2 || Ele.cols() == 1)
-    && "Code has only been tested for simplex sizes 3,2,1");
-
-  assert(m_primitive==-1 || (m_left == NULL && m_right == NULL));
-  if(is_leaf())
-  {
-    leaf_squared_distance(V,Ele,p,sqr_d,i,c);
-  }else
-  {
-    bool looked_left = false;
-    bool looked_right = false;
-    const auto & look_left = [&]()
-    {
-      int i_left;
-      RowVectorDIMS c_left = c;
-      Scalar sqr_d_left = m_left->squared_distance(V,Ele,p,sqr_d,i_left,c_left);
-      set_min(p,sqr_d_left,i_left,c_left,sqr_d,i,c);
-      looked_left = true;
-    };
-    const auto & look_right = [&]()
-    {
-      int i_right;
-      RowVectorDIMS c_right = c;
-      Scalar sqr_d_right = m_right->squared_distance(V,Ele,p,sqr_d,i_right,c_right);
-      set_min(p,sqr_d_right,i_right,c_right,sqr_d,i,c);
-      looked_right = true;
-    };
-
-    // must look left or right if in box
-    if(m_left->m_box.contains(p.transpose()))
-    {
-      look_left();
-    }
-    if(m_right->m_box.contains(p.transpose()))
-    {
-      look_right();
-    }
-    // if haven't looked left and could be less than current min, then look
-    Scalar  left_min_sqr_d = m_left->m_box.squaredExteriorDistance(p.transpose());
-    Scalar right_min_sqr_d = m_right->m_box.squaredExteriorDistance(p.transpose());
-    if(left_min_sqr_d < right_min_sqr_d)
-    {
-      if(!looked_left && left_min_sqr_d<sqr_d)
-      {
-        look_left();
-      }
-      if( !looked_right && right_min_sqr_d<sqr_d)
-      {
-        look_right();
-      }
-    }else
-    {
-      if( !looked_right && right_min_sqr_d<sqr_d)
-      {
-        look_right();
-      }
-      if(!looked_left && left_min_sqr_d<sqr_d)
-      {
-        look_left();
-      }
-    }
-  }
-  return sqr_d;
-}
-
-template <typename DerivedV, int DIM>
-template <
-  typename DerivedP, 
-  typename DerivedsqrD, 
-  typename DerivedI, 
-  typename DerivedC>
-inline void igl::AABB<DerivedV,DIM>::squared_distance(
-  const Eigen::PlainObjectBase<DerivedV> & V,
-  const Eigen::MatrixXi & Ele, 
-  const Eigen::PlainObjectBase<DerivedP> & P,
-  Eigen::PlainObjectBase<DerivedsqrD> & sqrD,
-  Eigen::PlainObjectBase<DerivedI> & I,
-  Eigen::PlainObjectBase<DerivedC> & C) const
-{
-  assert(P.cols() == V.cols() && "cols in P should match dim of cols in V");
-  sqrD.resize(P.rows(),1);
-  I.resize(P.rows(),1);
-  C.resize(P.rows(),P.cols());
-  for(int p = 0;p<P.rows();p++)
-  {
-    RowVectorDIMS Pp = P.row(p), c;
-    int Ip;
-    sqrD(p) = squared_distance(V,Ele,Pp,Ip,c);
-    I(p) = Ip;
-    C.row(p).head(DIM) = c;
-  }
-}
-
-template <typename DerivedV, int DIM>
-template < 
-  typename Derivedother_V,
-  typename DerivedsqrD, 
-  typename DerivedI, 
-  typename DerivedC>
-inline void igl::AABB<DerivedV,DIM>::squared_distance(
-  const Eigen::PlainObjectBase<DerivedV> & V,
-  const Eigen::MatrixXi & Ele, 
-  const AABB<Derivedother_V,DIM> & other,
-  const Eigen::PlainObjectBase<Derivedother_V> & other_V,
-  const Eigen::MatrixXi & other_Ele, 
-  Eigen::PlainObjectBase<DerivedsqrD> & sqrD,
-  Eigen::PlainObjectBase<DerivedI> & I,
-  Eigen::PlainObjectBase<DerivedC> & C) const
-{
-  assert(other_Ele.cols() == 1 && 
-    "Only implemented for other as list of points");
-  assert(other_V.cols() == V.cols() && "other must match this dimension");
-  sqrD.setConstant(other_Ele.rows(),1,std::numeric_limits<double>::infinity());
-  I.resize(other_Ele.rows(),1);
-  C.resize(other_Ele.rows(),other_V.cols());
-  // All points in other_V currently think they need to check against root of
-  // this. The point of using another AABB is to quickly prune chunks of
-  // other_V so that most points just check some subtree of this.
-
-  // This holds a conservative estimate of max(sqr_D) where sqr_D is the
-  // current best minimum squared distance for all points in this subtree
-  double min_sqr_d = std::numeric_limits<double>::infinity();
-  squared_distance_helper(
-    V,Ele,&other,other_V,other_Ele,min_sqr_d,sqrD,I,C);
-}
-
-template <typename DerivedV, int DIM>
-template < 
-  typename Derivedother_V,
-  typename DerivedsqrD, 
-  typename DerivedI, 
-  typename DerivedC>
-inline typename igl::AABB<DerivedV,DIM>::Scalar igl::AABB<DerivedV,DIM>::squared_distance_helper(
-  const Eigen::PlainObjectBase<DerivedV> & V,
-  const Eigen::MatrixXi & Ele, 
-  const AABB<Derivedother_V,DIM> * other,
-  const Eigen::PlainObjectBase<Derivedother_V> & other_V,
-  const Eigen::MatrixXi & other_Ele, 
-  const Scalar /*min_sqr_d*/,
-  Eigen::PlainObjectBase<DerivedsqrD> & sqrD,
-  Eigen::PlainObjectBase<DerivedI> & I,
-  Eigen::PlainObjectBase<DerivedC> & C) const
-{
-  using namespace std;
-  using namespace Eigen;
-
-  // This implementation is a bit disappointing. There's no major speed up. Any
-  // performance gains seem to come from accidental cache coherency and
-  // diminish for larger "other" (the opposite of what was intended).
-
-  // Base case
-  if(other->is_leaf() && this->is_leaf())
-  {
-    Scalar sqr_d = sqrD(other->m_primitive);
-    int i = I(other->m_primitive);
-    RowVectorDIMS c = C.row(      other->m_primitive);
-    RowVectorDIMS p = other_V.row(other->m_primitive);
-    leaf_squared_distance(V,Ele,p,sqr_d,i,c);
-    sqrD( other->m_primitive) = sqr_d;
-    I(    other->m_primitive) = i;
-    C.row(other->m_primitive) = c;
-    //cout<<"leaf: "<<sqr_d<<endl;
-    //other->m_max_sqr_d = sqr_d;
-    return sqr_d;
-  }
-
-  if(other->is_leaf())
-  {
-    Scalar sqr_d = sqrD(other->m_primitive);
-    int i = I(other->m_primitive);
-    RowVectorDIMS c = C.row(      other->m_primitive);
-    RowVectorDIMS p = other_V.row(other->m_primitive);
-    sqr_d = squared_distance(V,Ele,p,sqr_d,i,c);
-    sqrD( other->m_primitive) = sqr_d;
-    I(    other->m_primitive) = i;
-    C.row(other->m_primitive) = c;
-    //other->m_max_sqr_d = sqr_d;
-    return sqr_d;
-  }
-
-  //// Exact minimum squared distance between arbitary primitives inside this and
-  //// othre's bounding boxes
-  //const auto & min_squared_distance = [&](
-  //  const AABB<DerivedV,DIM> * A,
-  //  const AABB<Derivedother_V,DIM> * B)->Scalar
-  //{
-  //  return A->m_box.squaredExteriorDistance(B->m_box);
-  //};
-
-  if(this->is_leaf())
-  {
-    //if(min_squared_distance(this,other) < other->m_max_sqr_d)
-    if(true)
-    {
-      this->squared_distance_helper(
-        V,Ele,other->m_left,other_V,other_Ele,0,sqrD,I,C);
-      this->squared_distance_helper(
-        V,Ele,other->m_right,other_V,other_Ele,0,sqrD,I,C);
-    }else
-    {
-      // This is never reached...
-    }
-    //// we know other is not a leaf
-    //other->m_max_sqr_d = std::max(other->m_left->m_max_sqr_d,other->m_right->m_max_sqr_d);
-    return 0;
-  }
-
-  // FORCE DOWN TO OTHER LEAF EVAL
-  //if(min_squared_distance(this,other) < other->m_max_sqr_d)
-  if(true)
-  {
-    if(true)
-    {
-      this->squared_distance_helper(
-        V,Ele,other->m_left,other_V,other_Ele,0,sqrD,I,C);
-      this->squared_distance_helper(
-        V,Ele,other->m_right,other_V,other_Ele,0,sqrD,I,C);
-    }else // this direction never seems to be faster
-    {
-      this->m_left->squared_distance_helper(
-        V,Ele,other,other_V,other_Ele,0,sqrD,I,C);
-      this->m_right->squared_distance_helper(
-        V,Ele,other,other_V,other_Ele,0,sqrD,I,C);
-    }
-  }else
-  {
-    // this is never reached ... :-(
-  }
-  //// we know other is not a leaf
-  //other->m_max_sqr_d = std::max(other->m_left->m_max_sqr_d,other->m_right->m_max_sqr_d);
-
-  return 0;
-#if 0 // False
-
-  // _Very_ conservative approximation of maximum squared distance between
-  // primitives inside this and other's bounding boxes
-  const auto & max_squared_distance = [](
-    const AABB<DerivedV,DIM> * A,
-    const AABB<Derivedother_V,DIM> * B)->Scalar
-  {
-    AlignedBox<Scalar,DIM> combo = A->m_box;
-    combo.extend(B->m_box);
-    return combo.diagonal().squaredNorm();
-  };
-
-  //// other base-case
-  //if(other->is_leaf())
-  //{
-  //  double sqr_d = sqrD(other->m_primitive);
-  //  int i = I(other->m_primitive);
-  //  RowVectorDIMS c = C.row(m_primitive);
-  //  RowVectorDIMS p = other_V.row(m_primitive);
-  //  leaf_squared_distance(V,Ele,p,sqr_d,i,c);
-  //  sqrD(other->m_primitive) = sqr_d;
-  //  I(other->m_primitive) = i;
-  //  C.row(m_primitive) = c;
-  //  return;
-  //}
-  std::vector<const AABB<DerivedV,DIM> * > this_list;
-  if(this->is_leaf())
-  {
-    this_list.push_back(this);
-  }else
-  {
-    assert(this->m_left);
-    this_list.push_back(this->m_left);
-    assert(this->m_right);
-    this_list.push_back(this->m_right);
-  }
-  std::vector<AABB<Derivedother_V,DIM> *> other_list;
-  if(other->is_leaf())
-  {
-    other_list.push_back(other);
-  }else
-  {
-    assert(other->m_left);
-    other_list.push_back(other->m_left);
-    assert(other->m_right);
-    other_list.push_back(other->m_right);
-  }
-
-  //const std::function<Scalar(
-  //  const AABB<Derivedother_V,DIM> * other)
-  //    > max_sqr_d = [&sqrD,&max_sqr_d](const AABB<Derivedother_V,DIM> * other)->Scalar
-  //  {
-  //    if(other->is_leaf())
-  //    {
-  //      return sqrD(other->m_primitive);
-  //    }else
-  //    {
-  //      return std::max(max_sqr_d(other->m_left),max_sqr_d(other->m_right));
-  //    }
-  //  };
-
-  //// Potentially recurse on all pairs, if minimum distance is less than running
-  //// bound
-  //Eigen::Matrix<Scalar,Eigen::Dynamic,1> other_max_sqr_d =
-  //  Eigen::Matrix<Scalar,Eigen::Dynamic,1>::Constant(other_list.size(),1,min_sqr_d);
-  for(size_t child = 0;child<other_list.size();child++)
-  {
-    auto other_tree = other_list[child];
-
-    Eigen::Matrix<Scalar,Eigen::Dynamic,1> this_max_sqr_d(this_list.size(),1);
-    for(size_t t = 0;t<this_list.size();t++)
-    {
-      const auto this_tree = this_list[t];
-      this_max_sqr_d(t) = max_squared_distance(this_tree,other_tree);
-    }
-    if(this_list.size() ==2 &&
-      ( this_max_sqr_d(0) > this_max_sqr_d(1))
-      )
-    {
-      std::swap(this_list[0],this_list[1]);
-      //std::swap(this_max_sqr_d(0),this_max_sqr_d(1));
-    }
-    const Scalar sqr_d = this_max_sqr_d.minCoeff();
-
-
-    for(size_t t = 0;t<this_list.size();t++)
-    {
-      const auto this_tree = this_list[t];
-
-      //const auto mm = max_sqr_d(other_tree);
-      //const Scalar mc = other_max_sqr_d(child);
-      //assert(mc == mm);
-      // Only look left/right in this_list if can possible decrease somebody's
-      // distance in this_tree.
-      const Scalar min_this_other = min_squared_distance(this_tree,other_tree); 
-      if(
-          min_this_other < sqr_d && 
-          min_this_other < other_tree->m_max_sqr_d)
-      {
-        //cout<<"before: "<<other_max_sqr_d(child)<<endl;
-        //other_max_sqr_d(child) = std::min(
-        //  other_max_sqr_d(child),
-        //  this_tree->squared_distance_helper(
-        //    V,Ele,other_tree,other_V,other_Ele,other_max_sqr_d(child),sqrD,I,C));
-        //cout<<"after: "<<other_max_sqr_d(child)<<endl;
-          this_tree->squared_distance_helper(
-            V,Ele,other_tree,other_V,other_Ele,0,sqrD,I,C);
-      }
-    }
-  }
-  //const Scalar ret = other_max_sqr_d.maxCoeff();
-  //const auto mm = max_sqr_d(other);
-  //assert(mm == ret);
-  //cout<<"non-leaf: "<<ret<<endl;
-  //return ret;
-  if(!other->is_leaf())
-  {
-    other->m_max_sqr_d = std::max(other->m_left->m_max_sqr_d,other->m_right->m_max_sqr_d);
-  }
-  return 0;
-#endif
-}
-
-template <typename DerivedV, int DIM>
-inline void igl::AABB<DerivedV,DIM>::leaf_squared_distance(
-  const Eigen::PlainObjectBase<DerivedV> & V,
-  const Eigen::MatrixXi & Ele, 
-  const RowVectorDIMS & p,
-  Scalar & sqr_d,
-  int & i,
-  RowVectorDIMS & c) const
-{
-  using namespace Eigen;
-  using namespace std;
-  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);
-}
-
-
-template <typename DerivedV, int DIM>
-inline void igl::AABB<DerivedV,DIM>::set_min(
-  const RowVectorDIMS & 
-#ifndef NDEBUG
-  p
-#endif
-  ,
-  const Scalar sqr_d_candidate,
-  const int i_candidate,
-  const RowVectorDIMS & c_candidate,
-  Scalar & sqr_d,
-  int & i,
-  RowVectorDIMS & c) const
-{
-#ifndef NDEBUG
-  //std::cout<<matlab_format(c_candidate,"c_candidate")<<std::endl;
-  const Scalar pc_norm = (p-c_candidate).squaredNorm();
-  const Scalar diff = fabs(sqr_d_candidate - pc_norm);
-  assert(diff<=1e-10 && "distance should match norm of difference");
+#ifndef IGL_STATIC_LIBRARY
+#  include "AABB.cpp"
 #endif
-  if(sqr_d_candidate < sqr_d)
-  {
-    i = i_candidate;
-    c = c_candidate;
-    sqr_d = sqr_d_candidate;
-  }
-}
-
-
-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),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),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),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

@@ -147,4 +147,5 @@ template void igl::ambient_occlusion<Eigen::Matrix<double, -1, -1, 0, -1, -1>, E
 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::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::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> >&);
+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::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> >&);
 #endif

include/igl/embree/ambient_occlusion.cpp

@@ -58,4 +58,5 @@ template void igl::embree::ambient_occlusion<Eigen::Matrix<double, -1, -1, 0, -1
 template void igl::embree::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> >(igl::embree::EmbreeIntersector 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> >&);
 template void igl::embree::ambient_occlusion<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -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<int, -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&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template void igl::embree::ambient_occlusion<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > 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> >&);
+template void igl::embree::ambient_occlusion<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -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<int, -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&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/embree/ambient_occlusion.h

@@ -7,7 +7,7 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #ifndef IGL_EMBREE_AMBIENT_OCCLUSION_H
 #define IGL_EMBREE_AMBIENT_OCCLUSION_H
-#include <igl/igl_inline.h>
+#include "../igl_inline.h"
 #include <Eigen/Core>
 namespace igl
 {

include/igl/embree/reorient_facets_raycast.cpp

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2013 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 
+//
+// 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 "reorient_facets_raycast.h"
 #include "../per_face_normals.h"
@@ -17,8 +17,8 @@
 #include <limits>
 
 template <
-  typename DerivedV, 
-  typename DerivedF, 
+  typename DerivedV,
+  typename DerivedF,
   typename DerivedI,
   typename DerivedC>
 IGL_INLINE void igl::embree::reorient_facets_raycast(
@@ -36,37 +36,37 @@ IGL_INLINE void igl::embree::reorient_facets_raycast(
   using namespace std;
   assert(F.cols() == 3);
   assert(V.cols() == 3);
-  
+
   // number of faces
   const int m = F.rows();
-  
+
   MatrixXi FF = F;
   if (facet_wise) {
     C.resize(m);
     for (int i = 0; i < m; ++i) C(i) = i;
-  
+
   } else {
     if (is_verbose) cout << "extracting patches... ";
     bfs_orient(F,FF,C);
   }
   if (is_verbose) cout << (C.maxCoeff() + 1)  << " components. ";
-  
+
   // number of patches
   const int num_cc = C.maxCoeff()+1;
-  
+
   // Init Embree
   EmbreeIntersector ei;
   ei.init(V.template cast<float>(),FF);
-  
+
   // face normal
   MatrixXd N;
   per_face_normals(V,FF,N);
-  
+
   // face area
   Matrix<typename DerivedV::Scalar,Dynamic,1> A;
   doublearea(V,FF,A);
   double area_total = A.sum();
-  
+
   // determine number of rays per component according to its area
   VectorXd area_per_component;
   area_per_component.setZero(num_cc);
@@ -80,7 +80,7 @@ IGL_INLINE void igl::embree::reorient_facets_raycast(
     num_rays_per_component(c) = max<int>(static_cast<int>(rays_total * area_per_component(c) / area_total), rays_minimum);
   }
   rays_total = num_rays_per_component.sum();
-  
+
   // generate all the rays
   if (is_verbose) cout << "generating rays... ";
   uniform_real_distribution<float> rdist;
@@ -147,12 +147,12 @@ IGL_INLINE void igl::embree::reorient_facets_raycast(
     }
   }
   if (is_verbose) cout << ray_face.size()  << " rays. ";
-  
+
   // per component voting: first=front, second=back
   vector<pair<float, float>> C_vote_distance(num_cc, make_pair(0, 0));      // sum of distance between ray origin and intersection
   vector<pair<int  , int  >> C_vote_infinity(num_cc, make_pair(0, 0));      // number of rays reaching infinity
   vector<pair<int  , int  >> C_vote_parity(num_cc, make_pair(0, 0));        // sum of parity count for each ray
-  
+
   if (is_verbose) cout << "shooting rays... ";
 #pragma omp parallel for
   for (int i = 0; i < (int)ray_face.size(); ++i)
@@ -161,7 +161,7 @@ IGL_INLINE void igl::embree::reorient_facets_raycast(
     Vector3f o = ray_ori [i];
     Vector3f d = ray_dir [i];
     int c = C(f);
-    
+
     // shoot ray toward front & back
     vector<Hit> hits_front;
     vector<Hit> hits_back;
@@ -171,13 +171,13 @@ IGL_INLINE void igl::embree::reorient_facets_raycast(
     ei.intersectRay(o, -d, hits_back , num_rays_back );
     if (!hits_front.empty() && hits_front[0].id == f) hits_front.erase(hits_front.begin());
     if (!hits_back .empty() && hits_back [0].id == f) hits_back .erase(hits_back .begin());
-    
+
     if (use_parity) {
 #pragma omp atomic
       C_vote_parity[c].first  += hits_front.size() % 2;
 #pragma omp atomic
       C_vote_parity[c].second += hits_back .size() % 2;
-    
+
     } else {
       if (hits_front.empty())
       {
@@ -187,7 +187,7 @@ IGL_INLINE void igl::embree::reorient_facets_raycast(
 #pragma omp atomic
         C_vote_distance[c].first += hits_front[0].t;
       }
-    
+
       if (hits_back.empty())
       {
 #pragma omp atomic
@@ -198,14 +198,14 @@ IGL_INLINE void igl::embree::reorient_facets_raycast(
       }
     }
   }
-  
+
   I.resize(m);
   for(int f = 0; f < m; ++f)
   {
     int c = C(f);
     if (use_parity) {
       I(f) = C_vote_parity[c].first > C_vote_parity[c].second ? 1 : 0;      // Ideally, parity for the front/back side should be 1/0 (i.e., parity sum for all rays should be smaller on the front side)
-    
+
     } else {
       I(f) = (C_vote_infinity[c].first == C_vote_infinity[c].second && C_vote_distance[c].first <  C_vote_distance[c].second) ||
               C_vote_infinity[c].first <  C_vote_infinity[c].second
@@ -219,8 +219,8 @@ IGL_INLINE void igl::embree::reorient_facets_raycast(
 }
 
 template <
-  typename DerivedV, 
-  typename DerivedF, 
+  typename DerivedV,
+  typename DerivedF,
   typename DerivedFF,
   typename DerivedI>
 IGL_INLINE void igl::embree::reorient_facets_raycast(
@@ -255,4 +255,5 @@ IGL_INLINE void igl::embree::reorient_facets_raycast(
 // Explicit template specialization
 template void igl::embree::reorient_facets_raycast<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::embree::reorient_facets_raycast<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<bool, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, int, bool, bool, bool, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::embree::reorient_facets_raycast<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, int, bool, bool, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/list_to_matrix.cpp

@@ -146,5 +146,6 @@ template bool igl::list_to_matrix<double, Eigen::Matrix<double, -1, -1, 0, -1, -
 template bool igl::list_to_matrix<double, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::vector<double, std::allocator<double> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template bool igl::list_to_matrix<double, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::vector<double, std::allocator<double> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template bool igl::list_to_matrix<unsigned long, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::vector<unsigned long, std::allocator<unsigned long> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template bool igl::list_to_matrix<unsigned __int64, class Eigen::Matrix<int, -1, -1, 0, -1, -1> >(class std::vector<unsigned __int64, class std::allocator<unsigned __int64> > const &, class Eigen::PlainObjectBase<class Eigen::Matrix<int, -1, -1, 0, -1, -1> > &);
 #endif
 #endif

include/igl/png/readPNG.cpp

@@ -0,0 +1,47 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Daniele Panozzo <daniele.panozzo@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 "readPNG.h"
+#include <stb_image.h>
+
+IGL_INLINE bool igl::png::readPNG(
+  const std::string png_file,
+  Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic>& R,
+  Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic>& G,
+  Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic>& B,
+  Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic>& A
+)
+{
+  int cols,rows,n;
+  unsigned char *data = stbi_load(png_file.c_str(), &cols, &rows, &n, 4);
+  if(data == NULL) {
+    return false;
+  }
+
+  R.resize(cols,rows);
+  G.resize(cols,rows);
+  B.resize(cols,rows);
+  A.resize(cols,rows);
+
+  for (unsigned i=0; i<rows; ++i) {
+    for (unsigned j=0; j<cols; ++j) {
+      R(j,rows-1-i) = data[4*(j + cols * i) + 0];
+      G(j,rows-1-i) = data[4*(j + cols * i) + 1];
+      B(j,rows-1-i) = data[4*(j + cols * i) + 2];
+      A(j,rows-1-i) = data[4*(j + cols * i) + 3];
+    }
+  }
+
+  stbi_image_free(data);
+
+  return true;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+// generated by autoexplicit.sh
+#endif

include/igl/png/readPNG.h

@@ -0,0 +1,39 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Daniele Panozzo <daniele.panozzo@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_PNG_READ_PNG_H
+#define IGL_PNG_READ_PNG_H
+#include "../igl_inline.h"
+#include <Eigen/Core>
+#include <string>
+
+namespace igl
+{
+  namespace png
+  {
+    // Read an image from a .png file into 4 memory buffers
+    //
+    // Input:
+    //  png_file  path to .png file
+    // Output:
+    //  R,G,B,A texture channels
+    // Returns true on success, false on failure
+    //
+    IGL_INLINE bool readPNG(const std::string png_file,
+    Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic>& R,
+    Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic>& G,
+    Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic>& B,
+    Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic>& A
+    );
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "readPNG.cpp"
+#endif
+
+#endif

include/igl/png/writePNG.cpp

@@ -0,0 +1,47 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Daniele Panozzo <daniele.panozzo@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 "writePNG.h"
+#include <stb_image_write.h>
+#include <vector>
+
+IGL_INLINE bool igl::png::writePNG(
+  const Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic>& R,
+  const Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic>& G,
+  const Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic>& B,
+  const Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic>& A,
+  const std::string png_file
+)
+{
+  assert((R.rows() == G.rows()) && (G.rows() == B.rows()) && (B.rows() == A.rows()));
+  assert((R.cols() == G.cols()) && (G.cols() == B.cols()) && (B.cols() == A.cols()));
+
+  const int comp = 4;                                  // 4 Channels Red, Green, Blue, Alpha
+  const int stride_in_bytes = R.rows()*comp;           // Lenght of one row in bytes
+  std::vector<unsigned char> data(R.size()*comp,0);     // The image itself;
+
+  for (unsigned i = 0; i<R.rows();++i)
+  {
+    for (unsigned j = 0; j < R.cols(); ++j)
+    {
+        data[(j * R.rows() * comp) + (i * comp) + 0] = R(i,R.cols()-1-j);
+        data[(j * R.rows() * comp) + (i * comp) + 1] = G(i,R.cols()-1-j);
+        data[(j * R.rows() * comp) + (i * comp) + 2] = B(i,R.cols()-1-j);
+        data[(j * R.rows() * comp) + (i * comp) + 3] = A(i,R.cols()-1-j);
+    }
+  }
+
+  stbi_write_png(png_file.c_str(), R.rows(), R.cols(), comp, data.data(), stride_in_bytes);
+
+  return true;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+// generated by autoexplicit.sh
+
+#endif

include/igl/png/writePNG.h

@@ -0,0 +1,41 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Daniele Panozzo <daniele.panozzo@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_PNG_WRITE_PNG_H
+#define IGL_PNG_WRITE_PNG_H
+#include "../igl_inline.h"
+#include <Eigen/Core>
+#include <string>
+
+namespace igl
+{
+  namespace png
+  {
+    // Writes an image to a png file
+    //
+    // Input:
+    //  R,G,B,A texture channels
+    // Output:
+    //  png_file  path to .png file
+    // Returns true on success, false on failure
+    //
+    IGL_INLINE bool writePNG
+    (
+      const Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic>& R,
+      const Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic>& G,
+      const Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic>& B,
+      const Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic>& A,
+      const std::string png_file
+    );
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "writePNG.cpp"
+#endif
+
+#endif

include/igl/random_dir.cpp

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2013 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 
+//
+// 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 "random_dir.h"
 #include <igl/PI.h>
@@ -18,14 +18,14 @@ IGL_INLINE Eigen::Vector3d igl::random_dir()
   double r = sqrt(1.0-z*z);
   double x = r * cos(t);
   double y = r * sin(t);
-  return Vector3d(x,y,z); 
+  return Vector3d(x,y,z);
 }
 
 IGL_INLINE Eigen::MatrixXd igl::random_dir_stratified(const int n)
 {
   using namespace Eigen;
   using namespace std;
-  const double m = floor(sqrt(double(n)));
+  const double m = std::floor(sqrt(double(n)));
   MatrixXd N(n,3);
   int row = 0;
   for(int i = 0;i<m;i++)

include/igl/round.cpp

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2013 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 
+//
+// 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 "round.h"
 #include <cmath>
@@ -13,7 +13,7 @@
 template <typename DerivedX >
 IGL_INLINE DerivedX igl::round(const DerivedX r)
 {
-  return (r > 0.0) ? floor(r + 0.5) : ceil(r - 0.5);
+  return (r > 0.0) ? std::floor(r + 0.5) : std::ceil(r - 0.5);
 }
 
 template < typename DerivedX, typename DerivedY>

include/igl/slice.cpp

@@ -279,4 +279,5 @@ template Eigen::Matrix<double, -1, -1, 0, -1, -1> igl::slice<Eigen::Matrix<doubl
 template void igl::slice<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -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&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::slice<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::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::slice<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::slice<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
 #endif

include/igl/slice_mask.cpp

@@ -117,4 +117,8 @@ IGL_INLINE DerivedX igl::slice_mask(
 
 #ifdef IGL_STATIC_LIBRARY
 template void igl::slice_mask<Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Array<bool, -1, 1, 0, -1, 1> const&, int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::slice_mask<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Array<bool, -1, 1, 0, -1, 1> const&, Eigen::Array<bool, -1, 1, 0, -1, 1> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::slice_mask<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Array<bool, -1, 1, 0, -1, 1> const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::slice_mask<Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Array<bool, -1, 1, 0, -1, 1> const&, Eigen::Array<bool, -1, 1, 0, -1, 1> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+
 #endif

include/igl/triangle/triangulate.h

@@ -7,7 +7,7 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #ifndef IGL_TRIANGLE_TRIANGULATE_H
 #define IGL_TRIANGLE_TRIANGULATE_H
-#include <igl/igl_inline.h>
+#include "../igl_inline.h"
 #include <string>
 #include <Eigen/Core>
 

include/igl/unproject_onto_mesh.cpp

@@ -73,5 +73,6 @@ IGL_INLINE bool igl::unproject_onto_mesh(
 #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> >&);
+template bool igl::unproject_onto_mesh<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 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<double, -1, 1, 0, -1, 1> >&);
 #endif
 

include/igl/viewer/Viewer.cpp

@@ -346,13 +346,13 @@ namespace viewer
   T,t     Toggle filled faces
   Z       Snap to canonical view
   [,]     Toggle between rotation control types (e.g. trackball, two-axis
-          valuator with fixed up)
+          valuator with fixed up))"
 #ifdef IGL_VIEWER_WITH_NANOGUI
+		R"(
   ;       Toggle vertex labels
-  :       Toggle face labels
+  :       Toggle face labels)"
 #endif
-
-)");
+);
     std::cout<<usage;
 #endif
   }

index.html

@@ -14,9 +14,9 @@
 
 <h1 id="libigl-asimplecgeometryprocessinglibrary">libigl - A simple C++ geometry processing library</h1>
 
-<figure>
-<img src="libigl-teaser.png" alt="" />
-</figure>
+<p><a href="https://travis-ci.org/libigl/libigl"><img src="https://travis-ci.org/libigl/libigl.svg?branch=master" alt="Build Status" /></a>
+<a href="https://ci.appveyor.com/project/danielepanozzo/libigl-6hjk1"><img src="https://ci.appveyor.com/api/projects/status/mf3t9rnhco0vhly8?svg=true" alt="Build status" /></a>
+<img src="libigl-teaser.png" alt="" /></p>
 
 <p><a href="https://github.com/libigl/libigl/">https://github.com/libigl/libigl/</a></p>
 
@@ -51,6 +51,9 @@ the header-only default mode: (i.e. just include the headers you want to use).</
 group prototypes a lot in MATLAB, and we have a useful <a href="matlab-to-eigen.html">MATLAB to libigl+Eigen
 conversion table</a>.</p>
 
+<p>We regularly test compiling our library on Mac OS X with clang, Linux with gcc
+and Windows with Visual Studio 2015 Community Edition.</p>
+
 <h2 id="tutorial">Tutorial</h2>
 
 <p>As of version 1.0, libigl includes an introductory
@@ -168,7 +171,7 @@ and updating (i.e. committing) that change to the hash.</p>
 subrepos:</p>
 
 <pre><code class="bash">git pull
-git submodule update -- recursive
+git submodule update --recursive
 </code></pre>
 
 <h2 id="unittesting">Unit testing</h2>
@@ -199,7 +202,7 @@ BibTeX entry:</p>
   title = {{libigl}: A simple {C++} geometry processing library},
   author = {Alec Jacobson and Daniele Panozzo and others},
   note = {http://libigl.github.io/libigl/},
-  year = {2015},
+  year = {2016},
 }
 </code></pre>
 
@@ -210,23 +213,35 @@ Eurographics/ACM Symposium on Geometry Processing software award. Here are a
 few labs/companies/institutions using libigl:</p>
 
 <ul>
+<li><a href="http://www.adobe.com/technology/">Adobe Research</a></li>
+<li><a href="http://www.ea.com">Electronic Arts, Inc</a></li>
+<li><a href="http://meshconsultants.ca/">Mesh</a>, consultants, Canada</li>
+<li><a href="http://graphics.pixar.com/research/">Pixar Research</a></li>
 <li><a href="http://esotericsoftware.com/">Spine by Esoteric Software</a> is an animation tool dedicated to 2D characters.</li>
 <li>Columbia University, <a href="http://www.cs.columbia.edu/cg/">Columbia Computer Graphics Group</a>, USA</li>
 <li><a href="http://www.graphics.cornell.edu/">Cornell University</a>, USA</li>
+<li><a href="http://dcgi.felk.cvut.cz/">Czech Technical University in Prague</a>, Czech</li>
 <li>EPF Lausanne, <a href="http://lgg.epfl.ch/people.php">Computer Graphics and Geometry Laboratory</a>, Switzerland</li>
 <li>ETH Zurich, <a href="http://igl.ethz.ch/">Interactive Geometry Lab</a> and <a href="http://ait.inf.ethz.ch/">Advanced Technologies Lab</a>, Swizterland</li>
 <li>George Mason University, <a href="http://cs.gmu.edu/~ygingold/">CraGL</a>, USA</li>
-<li><a href="http://www.ust.hk/">Hong Kong University of Science and Technology</a>, USA</li>
+<li><a href="http://www.ust.hk/">Hong Kong University of Science and Technology</a>, Hong Kong</li>
 <li><a href="http://www.nii.ac.jp/en/">National Institute of Informatics</a>, Japan</li>
 <li>New York University, <a href="http://mrl.nyu.edu/">Media Research Lab</a>, USA</li>
 <li>NYUPoly, <a href="http://game.engineering.nyu.edu/">Game Innovation Lab</a>, USA</li>
-<li><a href="http://www.telecom-paristech.fr/en/formation-et-innovation-dans-le-numerique.html">Telecom ParisTech</a>, Paris, France</li>
+<li><a href="https://www.cg.tu-berlin.de">TU Berlin</a>, Germany</li>
 <li><a href="http://www.tudelft.nl/en/">TU Delft</a>, Netherlands</li>
+<li><a href="https://www.tuwien.ac.at/en/tuwien_home/">TU Wien</a>, Austria</li>
+<li><a href="http://www.telecom-paristech.fr/en/formation-et-innovation-dans-le-numerique.html">Telecom ParisTech</a>, Paris, France</li>
 <li><a href="http://mtm.ufsc.br/~leo/">Universidade Federal de Santa Catarina</a>, Brazil</li>
-<li><a href="http://www.usi.ch/en">Università della Svizzera Italiana</a>, Switzerland</li>
 <li><a href="http://vecg.cs.ucl.ac.uk/">University College London</a>, England</li>
+<li><a href="http://vis.berkeley.edu/">University of California Berkeley</a>, USA</li>
 <li><a href="http://www.cam.ac.uk/">University of Cambridge</a>, England</li>
 <li><a href="http://cg.cis.upenn.edu/">University of Pennsylvania</a>, USA</li>
+<li><a href="http://www.cs.utexas.edu/users/evouga/">University of Texas at Austin</a>, USA</li>
+<li><a href="https://www.csc.uvic.ca/Research/graphics/">University of Victoria</a>, Canada</li>
+<li><a href="http://www.usi.ch/en">Università della Svizzera Italiana</a>, Switzerland</li>
+<li><a href="http://www.univ-tlse3.fr/">Université Toulouse III Paul Sabatier</a>, France</li>
+<li><a href="http://www.math.zju.edu.cn/cagd/">Zhejiang University</a>, China</li>
 </ul>
 
 <h2 id="contact">Contact</h2>
@@ -248,8 +263,8 @@ page</a>.</p>
 
 <h2 id="copyright">Copyright</h2>
 
-<p>2015 Alec Jacobson, Daniele Panozzo, Christian Schüller, Olga Diamanti, Qingnan
-Zhou, Nico Pietroni, Stefan Brugger, Kenshi Takayama, Wenzel Jakob, Nikolas De
+<p>2016 Alec Jacobson, Daniele Panozzo, Christian Schüller, Olga Diamanti, Qingnan
+Zhou, Sebastian Koch, Nico Pietroni, Stefan Brugger, Kenshi Takayama, Wenzel Jakob, Nikolas De
 Giorgis, Luigi Rocca, Leonardo Sacht, Kevin Walliman, Olga Sorkine-Hornung, and others.</p>
 
 <p>Please see individual files for appropriate copyright notices.</p>

optional/index.html

@@ -50,14 +50,6 @@ make
 <p>These are (admittedly unpopular) functions that have never been used by us
 statically so we haven&#8217;t explicit instantiations (yet).</p>
 
-<h4 id="examples">Examples</h4>
-
-<p>You can make a slew of examples by issuing:</p>
-
-<pre><code>cd ../examples
-make
-</code></pre>
-
 <h4 id="external">External</h4>
 
 <p>Finally there are a number of external libraries that we include in
@@ -316,7 +308,7 @@ int main(int argc, char * argv[])
 {
 Eigen::MatrixXd V;
 Eigen::MatrixXi F;
-return (argc&gt;=2 &amp;amp;&amp;amp; igl::read_triangle_mesh(argv[1],V,F)?0:1);
+return (argc&gt;=2 &amp;&amp; igl::read_triangle_mesh(argv[1],V,F)?0:1);
 }
 </code></pre>
 

python/CMakeLists.txt

@@ -57,7 +57,7 @@ option(LIBIGL_WITH_NANOGUI          "Use Nanogui menu"   OFF)
 option(LIBIGL_WITH_CGAL             "Use CGAL"           OFF)
 option(LIBIGL_WITH_BOOLEAN          "Use Cork boolean"   OFF)
 option(LIBIGL_WITH_COMISO           "Use CoMiso"         ON)
-option(LIBIGL_WITH_EMBREE           "Use Embree"         OFF)
+option(LIBIGL_WITH_EMBREE           "Use Embree"         ON)
 option(LIBIGL_WITH_LIM              "Use LIM"            ON)
 option(LIBIGL_WITH_MATLAB           "Use Matlab"         OFF)
 option(LIBIGL_WITH_MOSEK            "Use MOSEK"          OFF)
@@ -142,10 +142,10 @@ elseif (UNIX)
   #Enable flag if undefined symbols appear on pyigl module import to get notified about the missing symbols at link time
   option(CHECK_UNDEFINED        "Check for undefined symbols"    OFF)
 
-  # Strip unnecessary sections of the binary on Linux/Mac OS 
+  # Strip unnecessary sections of the binary on Linux/Mac OS
   if(APPLE)
     set_target_properties(pyigl PROPERTIES MACOSX_RPATH ".")
-    
+
     if (NOT CHECK_UNDEFINED)
       set_target_properties(pyigl PROPERTIES LINK_FLAGS "-undefined dynamic_lookup -dead_strip")
     endif()

python/README.md

@@ -5,7 +5,7 @@
 <span style="color:#F62217">
 Everything in this folder is currently being developed and it is likely to be
 changed radically in the next couple of months, breaking compatibility between
-different version. We plan to stabilize the python API by the end of 2015.
+different version. We plan to stabilize the python API by the end of 2016.
 </span>
 
 ## Introduction
@@ -153,7 +153,7 @@ page](https://github.com/libigl/libigl/issues).
 
 ## Copyright
 2015 Alec Jacobson, Daniele Panozzo, Christian Schüller, Olga Diamanti, Qingnan
-Zhou, Nico Pietroni, Stefan Brugger, Kenshi Takayama, Wenzel Jakob, Nikolas De
+Zhou, Sebastian Koch, Nico Pietroni, Stefan Brugger, Kenshi Takayama, Wenzel Jakob, Nikolas De
 Giorgis, Luigi Rocca, Leonardo Sacht, Olga Sorkine-Hornung, and others.
 
 Please see individual files for appropriate copyright notices.

python/iglhelpers.py

@@ -8,7 +8,9 @@ def p2e(m):
             return igl.eigen.MatrixXi(m)
         elif m.dtype.type == np.float64:
             return igl.eigen.MatrixXd(m)
-        raise TypeError("p2e only support dtype float64 or int32")
+        elif m.dtype.type == np.bool:
+            return igl.eigen.MatrixXb(m)
+        raise TypeError("p2e only support dtype float64, int32 and bool")
     if sparse.issparse(m):
         # convert in a dense matrix with triples
         coo = m.tocoo()
@@ -34,6 +36,8 @@ def e2p(m):
         return np.array(m, dtype='float64')
     elif isinstance(m, igl.eigen.MatrixXi):
         return np.array(m, dtype='int32')
+    elif isinstance(m, igl.eigen.MatrixXb):
+        return np.array(m, dtype='bool')
     elif isinstance(m, igl.eigen.SparseMatrixd):
         coo = np.array(m.toCOO())
         I = coo[:, 0]

python/py_doc.cpp

@@ -32,6 +32,63 @@ const char *__doc_igl_local_basis = R"igl_Qu8mg5v7(// Compute a local orthogonal
   //   B3 eigen matrix #F by 3, normal of the triangle
   //
   // See also: adjacency_matrix)igl_Qu8mg5v7";
+const char *__doc_igl_signed_distance = R"igl_Qu8mg5v7(// Computes signed distance to a mesh
+  //
+  // Inputs:
+  //   P  #P by 3 list of query point positions
+  //   V  #V by 3 list of vertex positions
+  //   F  #F by ss list of triangle indices, ss should be 3 unless sign_type ==
+  //     SIGNED_DISTANCE_TYPE_UNSIGNED
+  //   sign_type  method for computing distance _sign_ S
+  // Outputs:
+  //   S  #P list of smallest signed distances
+  //   I  #P list of facet indices corresponding to smallest distances
+  //   C  #P by 3 list of closest points
+  //   N  #P by 3 list of closest normals (only set if
+  //     sign_type=SIGNED_DISTANCE_TYPE_PSEUDONORMAL)
+  //
+  // Known bugs: This only computes distances to triangles. So unreferenced
+  // vertices and degenerate triangles are ignored.)igl_Qu8mg5v7";
+const char *__doc_igl_signed_distance_pseudonormal = R"igl_Qu8mg5v7(// Computes signed distance to mesh
+  //
+  // Inputs:
+  //   tree  AABB acceleration tree (see AABB.h)
+  //   F  #F by 3 list of triangle indices
+  //   FN  #F by 3 list of triangle normals 
+  //   VN  #V by 3 list of vertex normals (ANGLE WEIGHTING)
+  //   EN  #E by 3 list of edge normals (UNIFORM WEIGHTING)
+  //   EMAP  #F*3 mapping edges in F to E
+  //   q  Query point
+  // Returns signed distance to mesh
+  //)igl_Qu8mg5v7";
+const char *__doc_igl_signed_distance_winding_number = R"igl_Qu8mg5v7(// Inputs:
+  //   tree  AABB acceleration tree (see cgal/point_mesh_squared_distance.h)
+  //   hier  Winding number evaluation hierarchy
+  //   q  Query point
+  // Returns signed distance to mesh)igl_Qu8mg5v7";
+const char *__doc_igl_triangle_triangulate = R"igl_Qu8mg5v7(// Triangulate the interior of a polygon using the triangle library.
+    //
+    // Inputs:
+    //   V #V by 2 list of 2D vertex positions
+    //   E #E by 2 list of vertex ids forming unoriented edges of the boundary of the polygon
+    //   H #H by 2 coordinates of points contained inside holes of the polygon
+    //   flags  string of options pass to triangle (see triangle documentation)
+    // Outputs:
+    //   V2  #V2 by 2  coordinates of the vertives of the generated triangulation
+    //   F2  #F2 by 3  list of indices forming the faces of the generated triangulation
+    //
+    // TODO: expose the option to prevent Steiner points on the boundary
+    //)igl_Qu8mg5v7";
+const char *__doc_igl_embree_ambient_occlusion = R"igl_Qu8mg5v7(// Compute ambient occlusion per given point
+    //
+    // Inputs:
+    //    ei  EmbreeIntersector containing (V,F)
+    //    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)
+    //)igl_Qu8mg5v7";
 const char *__doc_igl_cotmatrix = R"igl_Qu8mg5v7(// Constructs the cotangent stiffness matrix (discrete laplacian) for a given
   // mesh (V,F).
   //
@@ -86,6 +143,13 @@ const char *__doc_igl_per_face_normals = R"igl_Qu8mg5v7(// Compute face normals
   //   per_face_normals(V,F,Vector3d(1,1,1).normalized(),N);)igl_Qu8mg5v7";
 const char *__doc_igl_per_face_normals_stable = R"igl_Qu8mg5v7(// Special version where order of face indices is guaranteed not to effect
   // output.)igl_Qu8mg5v7";
+const char *__doc_igl_quad_planarity = R"igl_Qu8mg5v7(// Compute planarity of the faces of a quad mesh
+  // Inputs:
+  //   V  #V by 3 eigen Matrix of mesh vertex 3D positions
+  //   F  #F by 4 eigen Matrix of face (quad) indices
+  // Output:
+  //   P  #F by 1 eigen Matrix of mesh face (quad) planarities
+  //)igl_Qu8mg5v7";
 const char *__doc_igl_readOFF = R"igl_Qu8mg5v7(// Read a mesh from an ascii obj file, filling in vertex positions, normals
   // and texture coordinates. Mesh may have faces of any number of degree
   //
@@ -194,6 +258,21 @@ const char *__doc_igl_massmatrix = R"igl_Qu8mg5v7(// Constructs the mass (area)
   //
   // See also: adjacency_matrix
   //)igl_Qu8mg5v7";
+const char *__doc_igl_unproject_onto_mesh = R"igl_Qu8mg5v7(// 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)igl_Qu8mg5v7";
 const char *__doc_igl_colon = R"igl_Qu8mg5v7(// Colon operator like matlab's colon operator. Enumerats values between low
   // and hi with step step.
   // Templates:
@@ -271,6 +350,14 @@ const char *__doc_igl_gaussian_curvature = R"igl_Qu8mg5v7(// Compute discrete lo
   // Output:
   //   K  #V by 1 eigen Matrix of discrete gaussian curvature values
   //)igl_Qu8mg5v7";
+const char *__doc_igl_planarize_quad_mesh = R"igl_Qu8mg5v7(// Inputs:
+  //   Vin        #V by 3 eigen Matrix of mesh vertex 3D positions
+  //   F          #F by 4 eigen Matrix of face (quad) indices
+  //   maxIter    maximum numbers of iterations
+  //   threshold  minimum allowed threshold for non-planarity
+  // Output:
+  //   Vout       #V by 3 eigen Matrix of planar mesh vertex 3D positions
+  //)igl_Qu8mg5v7";
 const char *__doc_igl_avg_edge_length = R"igl_Qu8mg5v7(// Compute the average edge length for the given triangle mesh
   // Templates:
   //   DerivedV derived from vertex positions matrix type: i.e. MatrixXd
@@ -340,6 +427,17 @@ const char *__doc_igl_upsample = R"igl_Qu8mg5v7(// Subdivide a mesh without movi
   //
   // Known issues:
   //   - assumes (V,F) is edge-manifold.)igl_Qu8mg5v7";
+const char *__doc_igl_slice_mask = R"igl_Qu8mg5v7(// Act like the matlab X(row_mask,col_mask) operator, where
+  // row_mask, col_mask are non-negative integer indices.
+  // 
+  // Inputs:
+  //   X  m by n matrix
+  //   R  m list of row bools
+  //   C  n list of column bools
+  // Output:
+  //   Y  #trues-in-R by #trues-in-C matrix
+  //
+  // See also: slice_mask)igl_Qu8mg5v7";
 const char *__doc_igl_point_mesh_squared_distance = R"igl_Qu8mg5v7(// Compute distances from a set of points P to a triangle mesh (V,F)
   //
   // Inputs:
@@ -382,6 +480,25 @@ const char *__doc_igl_setdiff = R"igl_Qu8mg5v7(// Set difference of elements of
   //   C  (k<=m)-long vector of unique elements appearing in A but not in B
   //   IA  (k<=m)-long list of indices into A so that C = A(IA)
   //)igl_Qu8mg5v7";
+const char *__doc_igl_copyleft_tetgen_tetrahedralize = R"igl_Qu8mg5v7(// Mesh the interior of a surface mesh (V,F) using tetgen
+      //
+      // Inputs:
+      //   V  #V by 3 vertex position list
+      //   F  #F list of polygon face indices into V (0-indexed)
+      //   switches  string of tetgen options (See tetgen documentation) e.g.
+      //     "pq1.414a0.01" tries to mesh the interior of a given surface with
+      //       quality and area constraints
+      //     "" will mesh the convex hull constrained to pass through V (ignores F)
+      // Outputs:
+      //   TV  #V by 3 vertex position list
+      //   TT  #T by 4 list of tet face indices
+      //   TF  #F by 3 list of triangle face indices
+      // Returns status:
+      //   0 success
+      //   1 tetgen threw exception
+      //   2 tetgen did not crash but could not create any tets (probably there are
+      //     holes, duplicate faces etc.)
+      //   -1 other error)igl_Qu8mg5v7";
 const char *__doc_igl_comb_frame_field = R"igl_Qu8mg5v7(// Inputs:
   //   V            #V by 3 eigen Matrix of mesh vertex 3D positions
   //   F            #F by 4 eigen Matrix of face (quad) indices
@@ -672,12 +789,12 @@ const char *__doc_igl_cross_field_missmatch = R"igl_Qu8mg5v7(// Inputs:
   //   Handle_MMatch    #F by 3 eigen Matrix containing the integer missmatch of the cross field
   //                    across all face edges
   //)igl_Qu8mg5v7";
-const char *__doc_igl_grad = R"igl_Qu8mg5v7(// Gradient of a scalar function defined on piecewise linear elements (mesh)
-  // is constant on each triangle i,j,k:
-  // grad(Xijk) = (Xj-Xi) * (Vi - Vk)^R90 / 2A + (Xk-Xi) * (Vj - Vi)^R90 / 2A
-  // where Xi is the scalar value at vertex i, Vi is the 3D position of vertex
-  // i, and A is the area of triangle (i,j,k). ^R90 represent a rotation of
-  // 90 degrees
+const char *__doc_igl_grad = R"igl_Qu8mg5v7(// Gradient of a scalar function defined on piecewise linear elements (mesh)
+  // is constant on each triangle i,j,k:
+  // grad(Xijk) = (Xj-Xi) * (Vi - Vk)^R90 / 2A + (Xk-Xi) * (Vj - Vi)^R90 / 2A
+  // where Xi is the scalar value at vertex i, Vi is the 3D position of vertex
+  // i, and A is the area of triangle (i,j,k). ^R90 represent a rotation of
+  // 90 degrees
   //)igl_Qu8mg5v7";
 const char *__doc_igl_slice_into = R"igl_Qu8mg5v7(// Act like the matlab Y(row_indices,col_indices) = X
   // 

python/py_doc.h

@@ -1,10 +1,16 @@
 extern const char *__doc_igl_principal_curvature;
 extern const char *__doc_igl_local_basis;
+extern const char *__doc_igl_signed_distance;
+extern const char *__doc_igl_signed_distance_pseudonormal;
+extern const char *__doc_igl_signed_distance_winding_number;
+extern const char *__doc_igl_triangle_triangulate;
+extern const char *__doc_igl_embree_ambient_occlusion;
 extern const char *__doc_igl_cotmatrix;
 extern const char *__doc_igl_floor;
 extern const char *__doc_igl_slice;
 extern const char *__doc_igl_per_face_normals;
 extern const char *__doc_igl_per_face_normals_stable;
+extern const char *__doc_igl_quad_planarity;
 extern const char *__doc_igl_readOFF;
 extern const char *__doc_igl_per_vertex_normals;
 extern const char *__doc_igl_sortrows;
@@ -14,6 +20,7 @@ extern const char *__doc_igl_cat;
 extern const char *__doc_igl_eigs;
 extern const char *__doc_igl_per_corner_normals;
 extern const char *__doc_igl_massmatrix;
+extern const char *__doc_igl_unproject_onto_mesh;
 extern const char *__doc_igl_colon;
 extern const char *__doc_igl_fit_rotations;
 extern const char *__doc_igl_fit_rotations_planar;
@@ -21,14 +28,17 @@ extern const char *__doc_igl_fit_rotations_SSE;
 extern const char *__doc_igl_rotate_vectors;
 extern const char *__doc_igl_read_triangle_mesh;
 extern const char *__doc_igl_gaussian_curvature;
+extern const char *__doc_igl_planarize_quad_mesh;
 extern const char *__doc_igl_avg_edge_length;
 extern const char *__doc_igl_barycentric_coordinates;
 extern const char *__doc_igl_lscm;
 extern const char *__doc_igl_find_cross_field_singularities;
 extern const char *__doc_igl_upsample;
+extern const char *__doc_igl_slice_mask;
 extern const char *__doc_igl_point_mesh_squared_distance;
 extern const char *__doc_igl_parula;
 extern const char *__doc_igl_setdiff;
+extern const char *__doc_igl_copyleft_tetgen_tetrahedralize;
 extern const char *__doc_igl_comb_frame_field;
 extern const char *__doc_igl_map_vertices_to_circle;
 extern const char *__doc_igl_writeOBJ;

python/py_igl.cpp

@@ -62,6 +62,15 @@
 #include <igl/fit_rotations.h>
 #include <igl/polar_svd.h>
 #include <igl/covariance_scatter_matrix.h>
+#include <igl/slice_mask.h>
+#include <igl/signed_distance.h>
+#include <igl/quad_planarity.h>
+#include <igl/planarize_quad_mesh.h>
+#include <igl/triangle/triangulate.h>
+#include <igl/copyleft/tetgen/tetrahedralize.h>
+#include <igl/embree/ambient_occlusion.h>
+#include <igl/unproject_onto_mesh.h>
+//#include <igl/.h>
 
 
 void python_export_igl(py::module &m)
@@ -126,5 +135,14 @@ void python_export_igl(py::module &m)
 #include "py_igl/py_fit_rotations.cpp"
 #include "py_igl/py_polar_svd.cpp"
 #include "py_igl/py_covariance_scatter_matrix.cpp"
+#include "py_igl/py_slice_mask.cpp"
+#include "py_igl/py_signed_distance.cpp"
+#include "py_igl/py_quad_planarity.cpp"
+#include "py_igl/py_planarize_quad_mesh.cpp"
+#include "py_igl/py_triangle_triangulate.cpp"
+#include "py_igl/py_copyleft_tetgen_tetrahedralize.cpp"
+#include "py_igl/py_embree_ambient_occlusion.cpp"
+#include "py_igl/py_unproject_onto_mesh.cpp"
+//#include "py_igl/py_.cpp"
 
 }

python/py_igl/py_AABB.cpp

@@ -5,7 +5,7 @@ AABB
 .def(py::init<const igl::AABB<Eigen::MatrixXd,3>& >())
 .def("init",[](igl::AABB<Eigen::MatrixXd,3>& tree, const Eigen::MatrixXd& V, const Eigen::MatrixXi& Ele)
 {
-    return tree.init(V, Ele, Eigen::Matrix<double, Eigen::Dynamic, 3>(), Eigen::Matrix<double, Eigen::Dynamic, 3>(), Eigen::MatrixXi(), 0); 
+    return tree.init(V, Ele, Eigen::Matrix<double, Eigen::Dynamic, 3>(), Eigen::Matrix<double, Eigen::Dynamic, 3>(), Eigen::VectorXi(), 0); 
 })
 .def("squared_distance", [](const igl::AABB<Eigen::MatrixXd,3>& tree, const Eigen::MatrixXd& V, const Eigen::MatrixXi& Ele, const Eigen::MatrixXd& P, Eigen::MatrixXd& sqrD, Eigen::MatrixXi& I, Eigen::MatrixXd& C)
 {

python/py_igl/py_copyleft_tetgen_tetrahedralize.cpp

@@ -0,0 +1,16 @@
+
+m.def("copyleft_tetgen_tetrahedralize", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  const std::string switches,
+  Eigen::MatrixXd& TV,
+  Eigen::MatrixXi& TT,
+  Eigen::MatrixXi& TF
+)
+{
+  return igl::copyleft::tetgen::tetrahedralize(V, F, switches, TV, TT, TF);
+}, __doc_igl_copyleft_tetgen_tetrahedralize,
+py::arg("V"), py::arg("F"), py::arg("switches"), py::arg("TV"), py::arg("TT"), py::arg("TF"));
+
+

python/py_igl/py_embree_ambient_occlusion.cpp

@@ -0,0 +1,16 @@
+
+
+m.def("embree_ambient_occlusion", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  const Eigen::MatrixXd& P,
+  const Eigen::MatrixXd& N,
+  const int num_samples,
+  Eigen::MatrixXd& S
+)
+{
+  return igl::embree::ambient_occlusion(V, F, P, N, num_samples, S);
+}, __doc_igl_embree_ambient_occlusion,
+py::arg("V"), py::arg("F"), py::arg("P"), py::arg("N"), py::arg("num_samples"), py::arg("S"));
+

python/py_igl/py_planarize_quad_mesh.cpp

@@ -0,0 +1,14 @@
+
+m.def("planarize_quad_mesh", []
+(
+  const Eigen::MatrixXd& Vin,
+  const Eigen::MatrixXi& F,
+  const int maxIter,
+  const double & threshold,
+  Eigen::MatrixXd& Vout
+)
+{
+  return igl::planarize_quad_mesh(Vin, F, maxIter, threshold, Vout);
+}, __doc_igl_planarize_quad_mesh,
+py::arg("Vin"), py::arg("F"), py::arg("maxIter"), py::arg("threshold"), py::arg("Vout"));
+

python/py_igl/py_point_mesh_squared_distance.cpp

@@ -8,7 +8,7 @@ m.def("point_mesh_squared_distance", []
   Eigen::MatrixXd& C
 )
 {
-  assert_is_VectorX("I",I);
+//  assert_is_VectorX("I",I);
   return igl::point_mesh_squared_distance(P, V, Ele, sqrD, I, C);
 }, __doc_igl_point_mesh_squared_distance,
 py::arg("P"), py::arg("V"), py::arg("Ele"), py::arg("sqrD"), py::arg("I"), py::arg("C"));

python/py_igl/py_quad_planarity.cpp

@@ -0,0 +1,15 @@
+
+
+m.def("quad_planarity", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  Eigen::MatrixXd& P
+)
+{
+  Eigen::VectorXd Pv;
+  igl::quad_planarity(V, F, Pv);
+  P = Pv;
+}, __doc_igl_quad_planarity,
+py::arg("V"), py::arg("F"), py::arg("P"));
+

python/py_igl/py_signed_distance.cpp

@@ -0,0 +1,132 @@
+
+py::enum_<igl::SignedDistanceType>(m, "SignedDistanceType")
+    .value("SIGNED_DISTANCE_TYPE_PSEUDONORMAL", igl::SIGNED_DISTANCE_TYPE_PSEUDONORMAL)
+    .value("SIGNED_DISTANCE_TYPE_WINDING_NUMBER", igl::SIGNED_DISTANCE_TYPE_WINDING_NUMBER)
+    .value("SIGNED_DISTANCE_TYPE_DEFAULT", igl::SIGNED_DISTANCE_TYPE_DEFAULT)
+    .value("SIGNED_DISTANCE_TYPE_UNSIGNED", igl::SIGNED_DISTANCE_TYPE_UNSIGNED)
+    .value("NUM_SIGNED_DISTANCE_TYPE", igl::NUM_SIGNED_DISTANCE_TYPE)
+    .export_values();
+
+
+m.def("signed_distance", []
+(
+  const Eigen::MatrixXd& P,
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  const igl::SignedDistanceType sign_type,
+  Eigen::MatrixXd& S,
+  Eigen::MatrixXi& I,
+  Eigen::MatrixXd& C,
+  Eigen::MatrixXd& N
+)
+{
+  Eigen::VectorXd Sv;
+  Eigen::VectorXi Iv;
+  igl::signed_distance(P, V, F, sign_type, Sv, Iv, C, N);
+  S = Sv;
+  I = Iv;
+}, __doc_igl_signed_distance,
+py::arg("P"), py::arg("V"), py::arg("F"), py::arg("sign_type"), py::arg("S"), py::arg("I"), py::arg("C"), py::arg("N"));
+
+//m.def("signed_distance_pseudonormal", []
+//(
+//  const AABB<Eigen::MatrixXd, 3> & tree,
+//  const Eigen::MatrixXd& V,
+//  const Eigen::MatrixXi& F,
+//  const Eigen::MatrixXd& FN,
+//  const Eigen::MatrixXd& VN,
+//  const Eigen::MatrixXd& EN,
+//  const Eigen::MatrixXi& EMAP,
+//  const Eigen::MatrixXd& q
+//)
+//{
+//  assert_is_VectorX("q", q);
+//  assert_is_VectorX("EMAP",EMAP);
+//  return igl::signed_distance_pseudonormal(tree, V, F, FN, VN, EN, EMAP, q);
+//}, __doc_igl_signed_distance_pseudonormal,
+//py::arg("tree"), py::arg("V"), py::arg("F"), py::arg("FN"), py::arg("VN"), py::arg("EN"), py::arg("EMAP"), py::arg("q"));
+
+m.def("signed_distance_pseudonormal", []
+(
+  const Eigen::MatrixXd& P,
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  const igl::AABB<Eigen::MatrixXd, 3> & tree,
+  const Eigen::MatrixXd& FN,
+  const Eigen::MatrixXd& VN,
+  const Eigen::MatrixXd& EN,
+  const Eigen::MatrixXi& EMAP,
+  Eigen::MatrixXd& S,
+  Eigen::MatrixXi& I,
+  Eigen::MatrixXd& C,
+  Eigen::MatrixXd& N
+)
+{
+  assert_is_VectorX("EMAP", EMAP);
+  Eigen::VectorXi EMAPv;
+  if (EMAP.size() != 0)
+    EMAPv = EMAP;
+  Eigen::VectorXd Sv;
+  Eigen::VectorXi Iv;
+  igl::signed_distance_pseudonormal(P, V, F, tree, FN, VN, EN, EMAPv, Sv, Iv, C, N);
+  S = Sv;
+  I = Iv;
+}, __doc_igl_signed_distance_pseudonormal,
+py::arg("P"), py::arg("V"), py::arg("F"), py::arg("tree"), py::arg("FN"), py::arg("VN"), py::arg("EN"), py::arg("EMAP"), py::arg("S"), py::arg("I"), py::arg("C"), py::arg("N"));
+
+//m.def("signed_distance_pseudonormal", []
+//(
+//  const AABB<Eigen::MatrixXd, 3> & tree,
+//  const Eigen::MatrixXd& V,
+//  const Eigen::MatrixXi& F,
+//  const Eigen::MatrixXd& FN,
+//  const Eigen::MatrixXd& VN,
+//  const Eigen::MatrixXd& EN,
+//  const Eigen::MatrixXi & EMAP,
+//  const Eigen::MatrixXd & q,
+//  double & s,
+//  double & sqrd,
+//  int & i,
+//  Eigen::MatrixXd & c,
+//  Eigen::MatrixXd & n
+//)
+//{
+//  assert_is_VectorX("EMAP",EMAP);
+//  assert_is_VectorX("q",q);
+//  return igl::signed_distance_pseudonormal(tree, V, F, FN, VN, EN, EMAP, q, s, sqrd, i, c, n);
+//}, __doc_igl_signed_distance_pseudonormal,
+//py::arg("tree"), py::arg("V"), py::arg("F"), py::arg("FN"), py::arg("VN"), py::arg("EN"), py::arg("EMAP"), py::arg("q"), py::arg("s"), py::arg("sqrd"), py::arg("i"), py::arg("c"), py::arg("n"));
+
+//m.def("signed_distance_pseudonormal", []
+//(
+//  const AABB<Eigen::MatrixXd, 2> & tree,
+//  const Eigen::MatrixXd& V,
+//  const Eigen::MatrixXi& F,
+//  const Eigen::MatrixXd& FN,
+//  const Eigen::MatrixXd& VN,
+//  const Eigen::MatrixXd & q,
+//  double & s,
+//  double & sqrd,
+//  int & i,
+//  Eigen::MatrixXd & c,
+//  Eigen::MatrixXd & n
+//)
+//{
+//  assert_is_VectorX("q",q);
+//  return igl::signed_distance_pseudonormal(tree, V, F, FN, VN, q, s, sqrd, i, c, n);
+//}, __doc_igl_signed_distance_pseudonormal,
+//py::arg("tree"), py::arg("V"), py::arg("F"), py::arg("FN"), py::arg("VN"), py::arg("q"), py::arg("s"), py::arg("sqrd"), py::arg("i"), py::arg("c"), py::arg("n"));
+
+//m.def("signed_distance_winding_number", []
+//(
+//  AABB<Eigen::MatrixXd, 3> & tree,
+//  const Eigen::MatrixXd& V,
+//  const Eigen::MatrixXi& F,
+//  igl::WindingNumberAABB<Eigen::Vector3d> & hier,
+//  Eigen::RowVector3d & q
+//)
+//{
+//  return igl::signed_distance_winding_number(tree, V, F, hier, q);
+//}, __doc_igl_signed_distance_winding_number,
+//py::arg("tree"), py::arg("V"), py::arg("F"), py::arg("hier"), py::arg("q"));
+

python/py_igl/py_slice.cpp

@@ -41,6 +41,19 @@ m.def("slice", []
 }, __doc_igl_slice,
 py::arg("X"), py::arg("R"), py::arg("C"), py::arg("Y"));
 
+m.def("slice", []
+(
+  const Eigen::MatrixXd& X,
+  const Eigen::MatrixXi& R,
+  const int dim,
+  Eigen::MatrixXd& Y
+)
+{
+  assert_is_VectorX("R",R);
+  return igl::slice(X,R,dim,Y);
+}, __doc_igl_slice,
+py::arg("X"), py::arg("R"), py::arg("dim"), py::arg("Y"));
+
 m.def("slice", []
 (
   const Eigen::MatrixXd& X,

python/py_igl/py_slice_mask.cpp

@@ -0,0 +1,76 @@
+m.def("slice_mask", []
+(
+  const Eigen::MatrixXd& X,
+  const Eigen::Matrix<bool, Eigen::Dynamic, Eigen::Dynamic> & R,
+  const Eigen::Matrix<bool, Eigen::Dynamic, Eigen::Dynamic> & C,
+  Eigen::MatrixXd& Y
+)
+{
+  assert_is_VectorX("R",R);
+  assert_is_VectorX("C",C);
+  return igl::slice_mask(X, R, C, Y);
+}, __doc_igl_slice_mask,
+py::arg("X"), py::arg("R"), py::arg("C"), py::arg("Y"));
+
+m.def("slice_mask", []
+(
+  const Eigen::MatrixXd& X,
+  const Eigen::Matrix<bool, Eigen::Dynamic, Eigen::Dynamic> & R,
+  const int dim,
+  Eigen::MatrixXd& Y
+)
+{
+  assert_is_VectorX("R",R);
+  return igl::slice_mask(X, R, dim, Y);
+}, __doc_igl_slice_mask,
+py::arg("X"), py::arg("R"), py::arg("dim"), py::arg("Y"));
+
+m.def("slice_mask", []
+(
+  const Eigen::MatrixXi& X,
+  const Eigen::Matrix<bool, Eigen::Dynamic, Eigen::Dynamic> & R,
+  const Eigen::Matrix<bool, Eigen::Dynamic, Eigen::Dynamic> & C,
+  Eigen::MatrixXi& Y
+)
+{
+  assert_is_VectorX("R",R);
+  assert_is_VectorX("C",C);
+  return igl::slice_mask(X, R, C, Y);
+}, __doc_igl_slice_mask,
+py::arg("X"), py::arg("R"), py::arg("C"), py::arg("Y"));
+
+m.def("slice_mask", []
+(
+  const Eigen::MatrixXi& X,
+  const Eigen::Matrix<bool, Eigen::Dynamic, Eigen::Dynamic> & R,
+  const int dim,
+  Eigen::MatrixXi& Y
+)
+{
+  assert_is_VectorX("R",R);
+  return igl::slice_mask(X, R, dim, Y);
+}, __doc_igl_slice_mask,
+py::arg("X"), py::arg("R"), py::arg("dim"), py::arg("Y"));
+
+//m.def("slice_mask", []
+//(
+//  const Eigen::MatrixXd& X,
+//  Eigen::Array<bool, Eigen::Dynamic, 1> & R,
+//  Eigen::Array<bool, Eigen::Dynamic, 1> & C
+//)
+//{
+//  return igl::slice_mask(X, R, C);
+//}, __doc_igl_slice_mask,
+//py::arg("X"), py::arg("R"), py::arg("C"));
+
+//m.def("slice_mask", []
+//(
+//  const Eigen::MatrixXd& X,
+//  Eigen::Array<bool, Eigen::Dynamic, 1> & R,
+//  int dim
+//)
+//{
+//  return igl::slice_mask(X, R, dim);
+//}, __doc_igl_slice_mask,
+//py::arg("X"), py::arg("R"), py::arg("dim"));
+

python/py_igl/py_slice_tets.cpp

@@ -13,11 +13,9 @@ m.def("slice_tets", []
   Eigen::VectorXd pl;
   if (plane.size() != 0)
     pl = plane;
-  assert_is_VectorX("J", J);
   Eigen::VectorXi Jv;
-  if (J.size() != 0)
-    Jv = J;
-  return igl::slice_tets(V, T, pl, U, G, Jv, BC);
+  igl::slice_tets(V, T, pl, U, G, Jv, BC);
+  J = Jv;
 }, __doc_igl_slice_tets,
 py::arg("V"), py::arg("T"), py::arg("plane"), py::arg("U"), py::arg("G"), py::arg("J"), py::arg("BC"));
 

python/py_igl/py_triangle_triangulate.cpp

@@ -0,0 +1,16 @@
+
+
+m.def("triangle_triangulate", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& E,
+  const Eigen::MatrixXd& H,
+  const std::string flags,
+  Eigen::MatrixXd& V2,
+  Eigen::MatrixXi& F2
+)
+{
+  return igl::triangle::triangulate(V, E, H, flags, V2, F2);
+}, __doc_igl_triangle_triangulate,
+py::arg("V"), py::arg("E"), py::arg("H"), py::arg("flags"), py::arg("V2"), py::arg("F2"));
+

python/py_igl/py_unproject_onto_mesh.cpp

@@ -0,0 +1,57 @@
+// COMPLETE BINDINGS ========================
+
+m.def("unproject_onto_mesh", []
+(
+  const Eigen::MatrixXd & pos,
+  const Eigen::MatrixXd & model,
+  const Eigen::MatrixXd & proj,
+  const Eigen::MatrixXd & viewport,
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  Eigen::MatrixXi& fid, // TODO: Can we replace this with integer object reference?
+  Eigen::MatrixXd& bc
+)
+{
+  assert_is_Vector2("pos", pos);
+  Eigen::Vector2f posv;
+  if (pos.size() != 0)
+    posv = Eigen::Vector2f(pos.cast<float>());
+  assert_is_Matrix4("model", model);
+  Eigen::Matrix4f modelm;
+  if (model.size() != 0)
+    modelm = model.cast<float>();
+  assert_is_Matrix4("proj", proj);
+  Eigen::Matrix4f projm;
+  if (proj.size() != 0)
+    projm = proj.cast<float>();
+  assert_is_Vector4("viewport", viewport);
+  Eigen::Vector4f viewportv;
+  if (viewport.size() != 0)
+    viewportv = Eigen::Vector4f(viewport.cast<float>());
+
+  Eigen::VectorXd bcv;
+  int fidi;
+  bool ret = igl::unproject_onto_mesh(posv, modelm, projm, viewportv, V, F, fidi, bcv);
+  fid(0, 0) = fidi;
+  bc = bcv;
+  return ret;
+}, __doc_igl_unproject_onto_mesh,
+py::arg("pos"), py::arg("model"), py::arg("proj"), py::arg("viewport"), py::arg("V"), py::arg("F"), py::arg("fid"), py::arg("bc"));
+
+// INCOMPLETE BINDINGS ========================
+
+//m.def("unproject_onto_mesh", []
+//(
+//  Eigen::Vector2f & pos,
+//  Eigen::Matrix4f & model,
+//  Eigen::Matrix4f & proj,
+//  Eigen::Vector4f & viewport,
+//  std::function<bool (const Eigen::Vector3f &, const Eigen::Vector3f &, igl::Hit &)> & shoot_ray,
+//  int & fid,
+//  Eigen::MatrixXd& bc
+//)
+//{
+//  return igl::unproject_onto_mesh(pos, model, proj, viewport, shoot_ray, fid, bc);
+//}, __doc_igl_unproject_onto_mesh,
+//py::arg("pos"), py::arg("model"), py::arg("proj"), py::arg("viewport"), py::arg("shoot_ray"), py::arg("fid"), py::arg("bc"));
+

python/py_igl_viewer.cpp

@@ -356,6 +356,10 @@ py::class_<igl::viewer::ViewerCore> viewercore_class(me, "ViewerCore");
     .def("open_dialog_load_mesh", &igl::viewer::Viewer::open_dialog_load_mesh)
     .def("open_dialog_save_mesh", &igl::viewer::Viewer::open_dialog_save_mesh)
 
+    // Input handling
+    .def_readwrite("current_mouse_x", &igl::viewer::Viewer::current_mouse_x)
+    .def_readwrite("current_mouse_y", &igl::viewer::Viewer::current_mouse_y)
+
     // Callbacks
     .def_readwrite("callback_init", &igl::viewer::Viewer::callback_init)
     .def_readwrite("callback_pre_draw", &igl::viewer::Viewer::callback_pre_draw)

python/py_vector.cpp

@@ -126,9 +126,15 @@ py::class_<Type> bind_eigen_2(py::module &m, const char *name,
         .def("rightCols", [](Type &m, const int& k) { return Type(m.rightCols(k)); })
         .def("leftCols", [](Type &m, const int& k) { return Type(m.leftCols(k)); })
 
+        .def("setLeftCols", [](Type &m, const int& k, const Type& v) { return Type(m.leftCols(k) = v); })
+        .def("setRightCols", [](Type &m, const int& k, const Type& v) { return Type(m.rightCols(k) = v); })
+
         .def("topRows", [](Type &m, const int& k) { return Type(m.topRows(k)); })
         .def("bottomRows", [](Type &m, const int& k) { return Type(m.bottomRows(k)); })
 
+        .def("setTopRows", [](Type &m, const int& k, const Type& v) { return Type(m.topRows(k) = v); })
+        .def("setBottomRows", [](Type &m, const int& k, const Type& v) { return Type(m.bottomRows(k) = v); })
+
         .def("topLeftCorner", [](Type &m, const int& p, const int&q) { return Type(m.topLeftCorner(p,q)); })
         .def("bottomLeftCorner", [](Type &m, const int& p, const int&q) { return Type(m.bottomLeftCorner(p,q)); })
         .def("topRightCorner", [](Type &m, const int& p, const int&q) { return Type(m.topRightCorner(p,q)); })
@@ -167,6 +173,7 @@ py::class_<Type> bind_eigen_2(py::module &m, const char *name,
         .def("cwiseQuotient", [](const Type &m1, const Type &m2) -> Type { return m1.cwiseQuotient(m2); })
 
         /* Row and column-wise operations */
+        .def("rowwiseSet", [](Type &m, const Type &m2) {return Type(m.rowwise() = Eigen::Matrix<Scalar, 1, Eigen::Dynamic>(m2));} )
         .def("rowwiseSum", [](const Type &m) {return Type(m.rowwise().sum());} )
         .def("rowwiseProd", [](const Type &m) {return Type(m.rowwise().prod());} )
         .def("rowwiseMean", [](const Type &m) {return Type(m.rowwise().mean());} )
@@ -175,6 +182,7 @@ py::class_<Type> bind_eigen_2(py::module &m, const char *name,
         .def("rowwiseMinCoeff", [](const Type &m) {return Type(m.rowwise().minCoeff());} )
         .def("rowwiseMaxCoeff", [](const Type &m) {return Type(m.rowwise().maxCoeff());} )
 
+        .def("colwiseSet", [](Type &m, const Type &m2) {return Type(m.colwise() = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>(m2));} )
         .def("colwiseSum", [](const Type &m) {return Type(m.colwise().sum());} )
         .def("colwiseProd", [](const Type &m) {return Type(m.colwise().prod());} )
         .def("colwiseMean", [](const Type &m) {return Type(m.colwise().mean());} )
@@ -249,6 +257,26 @@ py::class_<Type> bind_eigen_2(py::module &m, const char *name,
         /* Comparison operators */
         .def(py::self == py::self)
         .def(py::self != py::self)
+        .def("__lt__", []
+        (const Type &a, const Scalar& b) -> Eigen::Matrix<bool,Eigen::Dynamic,Eigen::Dynamic>
+        {
+          return Eigen::Matrix<bool, Eigen::Dynamic, Eigen::Dynamic>(a.array() < b);
+        })
+        .def("__gt__", []
+        (const Type &a, const Scalar& b) -> Eigen::Matrix<bool,Eigen::Dynamic,Eigen::Dynamic>
+        {
+          return Eigen::Matrix<bool, Eigen::Dynamic, Eigen::Dynamic>(a.array() > b);
+        })
+        .def("__le__", []
+        (const Type &a, const Scalar& b) -> Eigen::Matrix<bool,Eigen::Dynamic,Eigen::Dynamic>
+        {
+          return Eigen::Matrix<bool, Eigen::Dynamic, Eigen::Dynamic>(a.array() <= b);
+        })
+        .def("__ge__", []
+        (const Type &a, const Scalar& b) -> Eigen::Matrix<bool,Eigen::Dynamic,Eigen::Dynamic>
+        {
+          return Eigen::Matrix<bool, Eigen::Dynamic, Eigen::Dynamic>(a.array() >= b);
+        })
 
         .def("transposeInPlace", [](Type &m) { m.transposeInPlace(); })
         /* Other transformations */
@@ -603,6 +631,9 @@ void python_export_vector(py::module &m) {
 //    py::implicitly_convertible<py::buffer, Eigen::MatrixXi>();
     //py::implicitly_convertible<double, Eigen::MatrixXi>();
 
+    /* Bindings for MatrixXb */
+    bind_eigen_2<Eigen::Matrix<bool,Eigen::Dynamic,Eigen::Dynamic> > (me, "MatrixXb");
+
     /* Bindings for MatrixXuc */
     bind_eigen_2<Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> > (me, "MatrixXuc");
     // py::implicitly_convertible<py::buffer, Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> >();

python/python_shared.cpp

@@ -24,6 +24,7 @@ PYBIND11_PLUGIN(pyigl) {
         .. autosummary::
            :toctree: _generate
 
+
     )pyigldoc");
 
     python_export_vector(m);

python/python_shared.h

@@ -31,6 +31,26 @@ void assert_is_RowVectorX(const std::string name, const Eigen::PlainObjectBase<S
     throw std::runtime_error(name + " must be a row vector.");
 }
 
+template<typename Scalar>
+void assert_is_Vector2(const std::string name, const Eigen::PlainObjectBase<Scalar>& v)
+{
+  if (v.size() == 0)
+    return;
+
+  if ((v.cols() != 1) || (v.rows() != 2))
+    throw std::runtime_error(name + " must be a column vector with 2 entries.");
+}
+
+template<typename Scalar>
+void assert_is_RowVector2(const std::string name, const Eigen::PlainObjectBase<Scalar>& v)
+{
+  if (v.size() == 0)
+    return;
+
+  if ((v.cols() != 2) || (v.rows() != 1))
+    throw std::runtime_error(name + " must be a row vector with 2 entries.");
+}
+
 template<typename Scalar>
 void assert_is_Vector3(const std::string name, const Eigen::PlainObjectBase<Scalar>& v)
 {

python/scripts/generate_bindings.py

@@ -100,7 +100,7 @@ def map_parameter_types(name, cpp_type, parsed_types, errors, enum_types):
                 else:
                     result.append("MatrixXd&")
                 break
-            if t == "MatrixXi":
+            if t == "MatrixXi" or t == "VectorXi":
                 result.append("MatrixXi&")
                 break
             if t == "MatrixXd" or t == "VectorXd":
@@ -197,7 +197,10 @@ if __name__ == '__main__':
     path = os.path.dirname(__file__)
     if path != "":
         os.chdir(path)
-    shutil.rmtree("generated")
+    try:
+        shutil.rmtree("generated")
+    except:
+        pass # Ignore missing generated directory
     os.makedirs("generated/complete")
     os.mkdir("generated/partial")
 

python/tutorial/105_Overlays.py

@@ -65,7 +65,7 @@ l1 = 'x: ' + str(m[0,0]) + ' y: ' + str(m[0,1]) + ' z: ' + str(m[0,2])
 viewer.data.add_label(m.transpose(),l1)
 
 l2 = 'x: ' + str(M[0,0]) + ' y: ' + str(M[0,1]) + ' z: ' + str(M[0,2])
-viewer.data.add_label(M.transpose(),l2);
+viewer.data.add_label(M.transpose(),l2)
 
 # Launch the viewer
-viewer.launch();
+viewer.launch()

python/tutorial/509_Planarization.py

@@ -0,0 +1,117 @@
+# Add the igl library to the modules search path
+import sys, os
+
+sys.path.insert(0, os.getcwd() + "/../")
+
+import pyigl as igl
+import random
+from math import cos, sin, pi
+
+TUTORIAL_SHARED_PATH = "../../tutorial/shared/"
+
+viewer = igl.viewer.Viewer()
+
+# Quad mesh generated from conjugate field
+VQC = igl.eigen.MatrixXd()
+FQC = igl.eigen.MatrixXi()
+FQCtri = igl.eigen.MatrixXi()
+PQC0 = igl.eigen.MatrixXd()
+PQC1 = igl.eigen.MatrixXd()
+PQC2 = igl.eigen.MatrixXd()
+PQC3 = igl.eigen.MatrixXd()
+
+# Planarized quad mesh
+VQCplan = igl.eigen.MatrixXd()
+FQCtriplan = igl.eigen.MatrixXi()
+PQC0plan = igl.eigen.MatrixXd()
+PQC1plan = igl.eigen.MatrixXd()
+PQC2plan = igl.eigen.MatrixXd()
+PQC3plan = igl.eigen.MatrixXd()
+
+
+def key_down(viewer, key, modifier):
+    if key == ord('1'):
+        # Draw the triangulated quad mesh
+        viewer.data.set_mesh(VQC, FQCtri)
+
+        # Assign a color to each quad that corresponds to its planarity
+        planarity = igl.eigen.MatrixXd()
+        igl.quad_planarity(VQC, FQC, planarity)
+        Ct = igl.eigen.MatrixXd()
+        igl.jet(planarity, 0, 0.01, Ct)
+        C = igl.eigen.MatrixXd(FQCtri.rows(), 3)
+        C.setTopRows(Ct.rows(), Ct)
+        C.setBottomRows(Ct.rows(), Ct)
+        viewer.data.set_colors(C)
+
+        # Plot a line for each edge of the quad mesh
+        viewer.data.add_edges(PQC0, PQC1, igl.eigen.MatrixXd([[0, 0, 0]]))
+        viewer.data.add_edges(PQC1, PQC2, igl.eigen.MatrixXd([[0, 0, 0]]))
+        viewer.data.add_edges(PQC2, PQC3, igl.eigen.MatrixXd([[0, 0, 0]]))
+        viewer.data.add_edges(PQC3, PQC0, igl.eigen.MatrixXd([[0, 0, 0]]))
+
+    elif key == ord('2'):
+        # Draw the planar quad mesh
+        viewer.data.set_mesh(VQCplan, FQCtri)
+
+        # Assign a color to each quad that corresponds to its planarity
+        planarity = igl.eigen.MatrixXd()
+        igl.quad_planarity(VQCplan, FQC, planarity)
+        Ct = igl.eigen.MatrixXd()
+        igl.jet(planarity, 0, 0.01, Ct)
+        C = igl.eigen.MatrixXd(FQCtri.rows(), 3)
+        C.setTopRows(Ct.rows(), Ct)
+        C.setBottomRows(Ct.rows(), Ct)
+        viewer.data.set_colors(C)
+
+        # Plot a line for each edge of the quad mesh
+        viewer.data.add_edges(PQC0plan, PQC1plan, igl.eigen.MatrixXd([[0, 0, 0]]))
+        viewer.data.add_edges(PQC1plan, PQC2plan, igl.eigen.MatrixXd([[0, 0, 0]]))
+        viewer.data.add_edges(PQC2plan, PQC3plan, igl.eigen.MatrixXd([[0, 0, 0]]))
+        viewer.data.add_edges(PQC3plan, PQC0plan, igl.eigen.MatrixXd([[0, 0, 0]]))
+
+    else:
+        return False
+
+    return True
+
+
+# Load a quad mesh generated by a conjugate field
+igl.readOFF(TUTORIAL_SHARED_PATH + "inspired_mesh_quads_Conjugate.off", VQC, FQC)
+
+# Convert it to a triangle mesh
+FQCtri.resize(2 * FQC.rows(), 3)
+
+FQCtriUpper = igl.eigen.MatrixXi(FQC.rows(), 3)
+FQCtriLower = igl.eigen.MatrixXi(FQC.rows(), 3)
+
+FQCtriUpper.setCol(0, FQC.col(0))
+FQCtriUpper.setCol(1, FQC.col(1))
+FQCtriUpper.setCol(2, FQC.col(2))
+FQCtriLower.setCol(0, FQC.col(2))
+FQCtriLower.setCol(1, FQC.col(3))
+FQCtriLower.setCol(2, FQC.col(0))
+
+FQCtri.setTopRows(FQCtriUpper.rows(), FQCtriUpper)
+FQCtri.setBottomRows(FQCtriLower.rows(), FQCtriLower)
+
+igl.slice(VQC, FQC.col(0), 1, PQC0)
+igl.slice(VQC, FQC.col(1), 1, PQC1)
+igl.slice(VQC, FQC.col(2), 1, PQC2)
+igl.slice(VQC, FQC.col(3), 1, PQC3)
+
+# Planarize it
+igl.planarize_quad_mesh(VQC, FQC, 100, 0.005, VQCplan)
+
+# Convert the planarized mesh to triangles
+igl.slice(VQCplan, FQC.col(0), 1, PQC0plan)
+igl.slice(VQCplan, FQC.col(1), 1, PQC1plan)
+igl.slice(VQCplan, FQC.col(2), 1, PQC2plan)
+igl.slice(VQCplan, FQC.col(3), 1, PQC3plan)
+
+# Launch the viewer
+key_down(viewer, ord('2'), 0)
+viewer.core.invert_normals = True
+viewer.core.show_lines = False
+viewer.callback_key_down = key_down
+viewer.launch()

python/tutorial/604_Triangle.py

@@ -0,0 +1,21 @@
+# Add the igl library to the modules search path
+import sys, os
+sys.path.insert(0, os.getcwd() + "/../")
+
+import pyigl as igl
+
+# Input polygon
+V = igl.eigen.MatrixXd([[-1, -1], [1, -1], [1, 1], [-1, 1], [-2, -2], [2, -2], [2, 2], [-2, 2]])
+E = igl.eigen.MatrixXi([[0, 1], [1, 2], [2, 3], [3, 0], [4, 5], [5, 6], [6,7], [7,4]])
+H = igl.eigen.MatrixXd([[0, 0]])
+
+# Triangulated Interior
+V2 = igl.eigen.MatrixXd()
+F2 = igl.eigen.MatrixXi()
+
+igl.triangle_triangulate(V, E, H, "a0.005q", V2, F2)
+
+# Plot the mesh
+viewer = igl.viewer.Viewer()
+viewer.data.set_mesh(V2, F2)
+viewer.launch()

python/tutorial/605_Tetgen.py

@@ -0,0 +1,71 @@
+# Add the igl library to the modules search path
+import sys, os
+sys.path.insert(0, os.getcwd() + "/../")
+
+import pyigl as igl
+
+TUTORIAL_SHARED_PATH = "../../tutorial/shared/"
+
+
+# Input polygon
+V = igl.eigen.MatrixXd()
+F = igl.eigen.MatrixXi()
+B = igl.eigen.MatrixXd()
+
+# Tetrahedralized interior
+TV = igl.eigen.MatrixXd()
+TT = igl.eigen.MatrixXi()
+TF = igl.eigen.MatrixXi()
+
+viewer = igl.viewer.Viewer()
+
+
+def key_down(viewer, key, modifier):
+    if key >= ord('1') and key <= ord('9'):
+        t = float((key - ord('1')) + 1) / 9.0
+        v = igl.eigen.MatrixXd()
+        v = B.col(2) - B.col(2).minCoeff()
+        v /= v.col(0).maxCoeff()
+
+        s = []
+        for i in range(v.size()):
+            if v[i, 0] < t:
+                s.append(i)
+
+        V_temp = igl.eigen.MatrixXd(len(s) * 4, 3)
+        F_temp = igl.eigen.MatrixXi(len(s) * 4, 3)
+
+        for i in range(len(s)):
+            V_temp.setRow(i * 4 + 0, TV.row(TT[s[i], 0]))
+            V_temp.setRow(i * 4 + 1, TV.row(TT[s[i], 1]))
+            V_temp.setRow(i * 4 + 2, TV.row(TT[s[i], 2]))
+            V_temp.setRow(i * 4 + 3, TV.row(TT[s[i], 3]))
+
+            F_temp.setRow(i * 4 + 0, igl.eigen.MatrixXi([[(i*4)+0, (i*4)+1, (i*4)+3]]))
+            F_temp.setRow(i * 4 + 1, igl.eigen.MatrixXi([[(i*4)+0, (i*4)+2, (i*4)+1]]))
+            F_temp.setRow(i * 4 + 2, igl.eigen.MatrixXi([[(i*4)+3, (i*4)+2, (i*4)+0]]))
+            F_temp.setRow(i * 4 + 3, igl.eigen.MatrixXi([[(i*4)+1, (i*4)+2, (i*4)+3]]))
+
+        viewer.data.clear()
+        viewer.data.set_mesh(V_temp, F_temp)
+        viewer.data.set_face_based(True)
+
+    else:
+        return False
+
+    return True
+
+
+# Load a surface mesh
+igl.readOFF(TUTORIAL_SHARED_PATH + "fertility.off", V, F)
+
+# Tetrahedralize the interior
+igl.copyleft_tetgen_tetrahedralize(V, F, "pq1.414Y", TV, TT, TF)
+
+# Compute barycenters
+igl.barycenter(TV, TT, B)
+
+# Plot the generated mesh
+key_down(viewer, ord('5'), 0)
+viewer.callback_key_down = key_down
+viewer.launch()

python/tutorial/606_AmbientOcclusion.py

@@ -0,0 +1,65 @@
+# Add the igl library to the modules search path
+import sys, os
+
+import math
+
+sys.path.insert(0, os.getcwd() + "/../")
+
+import pyigl as igl
+
+TUTORIAL_SHARED_PATH = "../../tutorial/shared/"
+
+
+# Mesh + AO values + Normals
+V = igl.eigen.MatrixXd()
+F = igl.eigen.MatrixXi()
+AO = igl.eigen.MatrixXd()
+N = igl.eigen.MatrixXd()
+
+
+viewer = igl.viewer.Viewer()
+
+
+def key_down(viewer, key, modifier):
+
+    color = igl.eigen.MatrixXd([[0.9, 0.85, 0.9]])
+
+    if key == ord('1'):
+        # Show the mesh without the ambient occlusion factor
+        viewer.data.set_colors(color)
+    elif key == ord('2'):
+        # Show the mesh with the ambient occlusion factor
+        C = color.replicate(V.rows(), 1)
+        for i in range(C.rows()):
+            C.setRow(i, C.row(i) * AO[i, 0])
+        viewer.data.set_colors(C)
+    elif key == ord('.'):
+        viewer.core.lighting_factor += 0.1
+    elif key == ord(','):
+        viewer.core.lighting_factor -= 0.1
+    else:
+        return False
+
+    viewer.core.lighting_factor = min(max(viewer.core.lighting_factor, 0.0), 1.0)
+    return True
+
+
+print("Press 1 to turn off Ambient Occlusion\nPress 2 to turn on Ambient Occlusion\nPress . to turn up lighting\nPress , to turn down lighting")
+
+# Load a surface mesh
+igl.readOFF(TUTORIAL_SHARED_PATH + "fertility.off", V, F)
+
+# Calculate vertex normals
+igl.per_vertex_normals(V, F, N)
+
+# Compute ambient occlusion factor using embree
+igl.embree_ambient_occlusion(V, F, V, N, 500, AO)
+AO = 1.0 - AO
+
+# Plot the generated mesh
+viewer.data.set_mesh(V, F)
+key_down(viewer, ord('2'), 0)
+viewer.callback_key_down = key_down
+viewer.core.show_lines = False
+viewer.core.lighting_factor = 0.0
+viewer.launch()

python/tutorial/607_Picking.py

@@ -0,0 +1,45 @@
+# Add the igl library to the modules search path
+import sys, os
+sys.path.insert(0, os.getcwd() + "/../")
+
+import pyigl as igl
+
+TUTORIAL_SHARED_PATH = "../../tutorial/shared/"
+
+
+# Mesh with per-face color
+V = igl.eigen.MatrixXd()
+F = igl.eigen.MatrixXi()
+C = igl.eigen.MatrixXd()
+
+viewer = igl.viewer.Viewer()
+
+def mouse_down(viewer, a, b):
+    bc = igl.eigen.MatrixXd()
+
+    # Cast a ray in the view direction starting from the mouse position
+    fid = igl.eigen.MatrixXi([-1])
+    coord = igl.eigen.MatrixXd([viewer.current_mouse_x, viewer.core.viewport[3] - viewer.current_mouse_y])
+    hit = igl.unproject_onto_mesh(coord, viewer.core.view * viewer.core.model,
+      viewer.core.proj, viewer.core.viewport, V, F, fid, bc)
+    if hit:
+        C.setRow(fid[0, 0], igl.eigen.MatrixXd([[1, 0, 0]]))
+        viewer.data.set_colors(C)
+        return True
+
+    return False
+
+
+print("Usage: [LeftMouseClick] to select a face")
+
+# Load a mesh in OFF format
+igl.readOFF(TUTORIAL_SHARED_PATH + "fertility.off", V, F)
+
+# Initialize white
+C.setConstant(F.rows(), 3, 1.0)
+
+viewer.data.set_mesh(V, F)
+viewer.data.set_colors(C)
+viewer.core.show_lines = False
+viewer.callback_mouse_down = mouse_down
+viewer.launch()

python/tutorial/704_SignedDistance.py

@@ -2,11 +2,13 @@ from __future__ import print_function
 
 # Add the igl library to the modules search path
 import sys, os
+
 sys.path.insert(0, os.getcwd() + "/../")
 
 import pyigl as igl
 from iglhelpers import e2p
 import math
+
 TUTORIAL_SHARED_PATH = "../../tutorial/shared/"
 
 global V, F, T, tree, FN, VN, EN, E, EMAP, max_distance, slice_z, overlay
@@ -27,9 +29,21 @@ overlay = False
 
 viewer = igl.viewer.Viewer()
 
+
+def append_mesh(C_vis, F_vis, V_vis, V, F, color):
+    F_vis.conservativeResize(F_vis.rows() + F.rows(), 3)
+    F_vis.setBottomRows(F.rows(), F + V_vis.rows())
+    V_vis.conservativeResize(V_vis.rows() + V.rows(), 3)
+    V_vis.setBottomRows(V.rows(), V)
+    C_vis.conservativeResize(C_vis.rows() + V.rows(), 3)
+    colorM = igl.eigen.MatrixXd(V.rows(), C_vis.cols())
+    colorM.rowwiseSet(color)
+    C_vis.setBottomRows(V.rows(), colorM)
+
+
 def update_visualization(viewer):
     global V, F, T, tree, FN, VN, EN, E, EMAP, max_distance, slice_z, overlay
-    plane = igl.eigen.MatrixXd([0.0, 0.0, 1.0, -((1-slice_z) * V.col(2).minCoeff() + slice_z * V.col(2).maxCoeff())])
+    plane = igl.eigen.MatrixXd([0.0, 0.0, 1.0, -((1 - slice_z) * V.col(2).minCoeff() + slice_z * V.col(2).maxCoeff())])
     V_vis = igl.eigen.MatrixXd()
     F_vis = igl.eigen.MatrixXi()
 
@@ -38,60 +52,47 @@ def update_visualization(viewer):
     bary = igl.eigen.SparseMatrixd()
     igl.slice_tets(V, T, plane, V_vis, F_vis, J, bary)
     max_l = 0.03
-#    while True:
-#        l = igl.eigen.MatrixXd()
-#        igl.edge_lengths(V_vis, F_vis, l)
-#        l /= (V_vis.colwise().maxCoeff() - V_vis.colwise().minCoeff()).norm()
-#        
-#        if l.maxCoeff() < max_l:
-#            break
-#        
-#        bad = e2p(l.rowwiseMaxCoeff())
-#        bad = bad > max_l
-#        F_vis_bad = igl.eigen.MatrixXi()
-#        F_vis_good = igl.eigen.MatrixXi()
-#        igl::slice_mask(F_vis, bad, 1, F_vis_bad);
-#        igl::slice_mask(F_vis, (bad!=true).eval(), 1, F_vis_good);
-#        igl.upsample(V_vis, F_vis_bad)
-#        F_vis = igl.cat(1, F_vis_bad, F_vis_good)
-
-
-#    #Compute signed distance
-#    S_vis = igl.eigen.MatrixXd()
-#    I = igl.eigen.MatrixXi()
-#    N = igl.eigen.MatrixXd()
-#    C = igl.eigen.MatrixXd()
-
-#    # Bunny is a watertight mesh so use pseudonormal for signing
-#    igl.signed_distance_pseudonormal(V_vis, V, F, tree, FN, VN, EN, EMAP, S_vis, I, C, N)
-
-#    # push to [0,1] range
-#    S_vis.array() = 0.5*(S_vis.array()/max_distance)+0.5;
-#    C_vis = igl.eigen.MatrixXi()
-#    # color without normalizing
-#    igl.parula(S_vis, False, C_vis)
-
-
-#    const auto & append_mesh = [&C_vis,&F_vis,&V_vis](const Eigen::MatrixXd & V, const Eigen::MatrixXi & F, const RowVector3d & color)
-
-#    F_vis.conservativeResize(F_vis.rows() + F.rows(), 3)
-#    F_vis.bottomRows(F.rows()) = F.array() + V_vis.rows()
-#    V_vis.conservativeResize(V_vis.rows() + V.rows(), 3)
-#    V_vis.bottomRows(V.rows()) = V
-#    C_vis.conservativeResize(C_vis.rows() + V.rows(), 3)
-#    C_vis.bottomRows(V.rows()).rowwise() = color
-
-#    if overlay:
-#        append_mesh(V, F, RowVector3d(0.8,0.8,0.8))
+    while True:
+        l = igl.eigen.MatrixXd()
+        igl.edge_lengths(V_vis, F_vis, l)
+        l /= (V_vis.colwiseMaxCoeff() - V_vis.colwiseMinCoeff()).norm()
+
+        if l.maxCoeff() < max_l:
+            break
+
+        bad = l.rowwiseMaxCoeff() > max_l
+        notbad = l.rowwiseMaxCoeff() <= max_l  # TODO replace by ~ operator
+        F_vis_bad = igl.eigen.MatrixXi()
+        F_vis_good = igl.eigen.MatrixXi()
+        igl.slice_mask(F_vis, bad, 1, F_vis_bad)
+        igl.slice_mask(F_vis, notbad, 1, F_vis_good)
+        igl.upsample(V_vis, F_vis_bad)
+        F_vis = igl.cat(1, F_vis_bad, F_vis_good)
+
+    # Compute signed distance
+    S_vis = igl.eigen.MatrixXd()
+    I = igl.eigen.MatrixXi()
+    N = igl.eigen.MatrixXd()
+    C = igl.eigen.MatrixXd()
+
+    # Bunny is a watertight mesh so use pseudonormal for signing
+    igl.signed_distance_pseudonormal(V_vis, V, F, tree, FN, VN, EN, EMAP, S_vis, I, C, N)
+
+    # push to [0,1] range
+    S_vis = 0.5 * (S_vis / max_distance) + 0.5
+    C_vis = igl.eigen.MatrixXd()
+    # color without normalizing
+    igl.parula(S_vis, False, C_vis)
+
+    if overlay:
+        append_mesh(C_vis, F_vis, V_vis, V, F, igl.eigen.MatrixXd([[0.8, 0.8, 0.8]]))
 
     viewer.data.clear()
     viewer.data.set_mesh(V_vis, F_vis)
-#    viewer.data.set_colors(C_vis)
+    viewer.data.set_colors(C_vis)
     viewer.core.lighting_factor = overlay
 
 
-
-
 def key_down(viewer, key, modifier):
     global slice_z, overlay
 
@@ -111,30 +112,26 @@ def key_down(viewer, key, modifier):
 print("Press [space] to toggle showing surface.")
 print("Press '.'/',' to push back/pull forward slicing plane.")
 
-#Load mesh: (V,T) tet-mesh of convex hull, F contains original surface triangles
-igl.readMESH(TUTORIAL_SHARED_PATH + "bunny.mesh", V, T, F);
+# Load mesh: (V,T) tet-mesh of convex hull, F contains original surface triangles
+igl.readMESH(TUTORIAL_SHARED_PATH + "bunny.mesh", V, T, F)
 
-#Call to point_mesh_squared_distance to determine bounds
+# Call to point_mesh_squared_distance to determine bounds
 sqrD = igl.eigen.MatrixXd()
 I = igl.eigen.MatrixXi()
 C = igl.eigen.MatrixXd()
 igl.point_mesh_squared_distance(V, V, F, sqrD, I, C)
 max_distance = math.sqrt(sqrD.maxCoeff())
 
-#Precompute signed distance AABB tree
+# Precompute signed distance AABB tree
 tree.init(V, F)
 
-#Precompute vertex, edge and face normals
+# Precompute vertex, edge and face normals
 igl.per_face_normals(V, F, FN)
 igl.per_vertex_normals(V, F, igl.PER_VERTEX_NORMALS_WEIGHTING_TYPE_ANGLE, FN, VN)
 igl.per_edge_normals(V, F, igl.PER_EDGE_NORMALS_WEIGHTING_TYPE_UNIFORM, FN, EN, E, EMAP)
 
-#Plot the generated mesh
+# Plot the generated mesh
 update_visualization(viewer);
 viewer.callback_key_down = key_down
 viewer.core.show_lines = False
 viewer.launch()
-
-
-
-

style-guidelines.html

@@ -1,3 +1,17 @@
+<!DOCTYPE html>
+<html>
+<head>
+	<meta charset="utf-8"/>
+	<title>libigl</title>
+	<meta name="author" content="Alec Jacobson and Daniele Panozzo and others"/>
+	<link type="text/css" rel="stylesheet" href="tutorial/style.css"/>
+<script type='text/javascript' src='http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML'></script>
+<link rel='stylesheet' href='http://yandex.st/highlightjs/7.3/styles/default.min.css'>
+<script src='http://yandex.st/highlightjs/7.3/highlight.min.js'></script>
+<script>hljs.initHighlightingOnLoad();</script>
+</head>
+<body>
+
 <h1 id="libiglstyleguidelines">Libigl Style Guidelines</h1>
 
 <p>Libigl is used and developed by many people. This document highlights some
@@ -39,8 +53,6 @@ namespace igl
   // This is an example of a function, it takes a templated parameter and
   // shovels it into cout
   //
-  // Templates:
-  //   T  type that supports
   // Input:
   //   input  some input of a Printable type
   // Returns true for the sake of returning something
@@ -94,8 +106,10 @@ new pair of .h/.cpp files with this sub-function.</p>
 then avoid crowding the namespace by creating lambda functions within the
 function implementation.</p>
 
-<p>These lambda functions must still be documented with clear <a href="#header-documentation">input and output
-arguments</a>.</p>
+<p>These lambda functions must still be documented with clear <a href="#headerdocumentation">input and output
+arguments</a>. Avoid using full capturing of all automatic
+variables: do not use <code>[&amp;]</code> or <code>[=]</code>. Rather specify each captured variable
+individually.</p>
 
 <h3 id="avoidhelperclasses">Avoid &#8220;helper&#8221; classes</h3>
 
@@ -157,7 +171,7 @@ than pointers (e.g. <code>Matrix * mat</code>) or value (e.g. <code>Matrix mat</
 <p>All functions should be implemented with at least one overload that has a
 <code>void</code> or simple return type (e.g. <code>bool</code> on success/failure). With this
 implementation its then possible to write an overload that returns a single
-output.</p>
+output. Please see <a href="#templatingwitheigen">Templating with Eigen</a>.</p>
 
 <p>For example:</p>
 
@@ -168,6 +182,77 @@ template &lt;typename Atype&gt;
 Eigen::SparseMatrix&lt;Atype&gt; adjacency_matrix(const ... &amp; F);
 </code></pre>
 
+<h2 id="templatingwitheigen">Templating with Eigen</h2>
+
+<p>Functions taking Eigen dense matrices/arrays as inputs and outputs (but <strong>not</strong>
+return arguments), should template on top of <code>Eigen::PlainObjectBase</code>. **Each
+parameter** should be derived using its own template.</p>
+
+<p>For example,</p>
+
+<pre><code class="cpp">template &lt;typename DerivedV, typename DerivedF, typename DerivedBC&gt;
+void barycenter(
+  const Eigen::PlainObjectBase&lt;DerivedV&gt; &amp; V,
+  const Eigen::PlainObjectBase&lt;DerivedF&gt; &amp; F,
+  const Eigen::PlainObjectBase&lt;DerivedBC&gt; &amp; BC);
+</code></pre>
+
+<p>The <code>Derived*</code> template encodes the scalar type (e.g. <code>double</code>, <code>int</code>), the
+number of rows and cols at compile time, and the data storage (Row-major vs.
+column-major). </p>
+
+<p>Returning Eigen types is discouraged. In cases where the size and scalar type
+are a fixed <strong>and matching</strong> function of an input <code>Derived*</code> template, then
+return that <code>Derived*</code> type. <strong>Do not</strong> return
+<code>Eigen::PlainObjectBase&lt;...&gt;</code> types. For example, this function scales fits a
+given set of points to the unit cube. The return is a new set of vertex
+positions so its type should <em>match</em> that of the input points:</p>
+
+<pre><code class="cpp">template &lt;typename DerivedV&gt;
+void DerivedV fit_to_unit_cube(const Eigen::PlainObjectBase&lt;DerivedV&gt; &amp; V);
+</code></pre>
+
+<p>To implement this function, it is <strong>required</strong> to implement a more generic
+output-argument version and call that. So a full implementation looks like:</p>
+
+<p>In <code>igl/fit_in_unit_cube.h</code>:</p>
+
+<pre><code class="cpp">template &lt;typename DerivedV, typename DerivedW&gt;
+void fit_to_unit_cube(
+  const Eigen::PlainObjectBase&lt;DerivedV&gt; &amp; V,
+  Eigen::PlainObjectBase&lt;DerivedW&gt; &amp; W);
+template &lt;typename DerivedV&gt;
+void DerivedV fit_to_unit_cube(const Eigen::PlainObjectBase&lt;DerivedV&gt; &amp; V);
+</code></pre>
+
+<p>In <code>igl/fit_in_unit_cube.cpp</code>:</p>
+
+<pre><code>template &lt;typename DerivedV, typename DerivedW&gt;
+void fit_to_unit_cube(
+  const Eigen::PlainObjectBase&lt;DerivedV&gt; &amp; V,
+  Eigen::PlainObjectBase&lt;DerivedW&gt; &amp; W)
+{
+  W = (V.rowwise()-V.colwise().minCoeff()).array() /
+    (V.maxCoeff()-V.minCoeff());
+}
+
+template &lt;typename DerivedV&gt;
+void DerivedV fit_to_unit_cube(const Eigen::PlainObjectBase&lt;DerivedV&gt; &amp; V)
+{
+  DerivedV W;
+  fit_to_unit_cube(V,W);
+  return W;
+}
+</code></pre>
+
+<p>Notice that <code>W</code> is declared as a <code>DerivedV</code> type and <strong>not</strong>
+<code>Eigen::PlainObjectBase&lt;DerivedV&gt;</code> type.</p>
+
+<p><strong>Note:</strong> Not all functions are suitable for returning Eigen types. For example
+<code>igl::barycenter</code> above outputs a #F by dim list of barycenters. Returning a
+<code>DerivedV</code> type would be inappropriate since the number of rows in <code>DerivedV</code>
+will be #V and may not match the number of rows in <code>DerivedF</code> (#F).</p>
+
 <h2 id="functionnamingconventions">Function naming conventions</h2>
 
 <p>Functions (and <a href="#filefunction">thus also files</a>) should have simple,
@@ -179,8 +264,8 @@ unnecessary prefaces. For example, instead of <code>compute_adjacency_matrix</co
 
 <h2 id="variablenamingconventions">Variable naming conventions</h2>
 
-<p>Libigl prefers short (even single character) variable names <em>with heavy
-documentation</em> in the comments in the header file or above the declaration of
+<p>Libigl prefers short (even single character) variable names _with heavy
+documentation_ in the comments in the header file or above the declaration of
 the function. When possible use <code>V</code> to mean a list of vertex positions and <code>F</code>
 to mean a list of faces/triangles.</p>
 
@@ -293,3 +378,17 @@ edited by you first. This means for
 
 <p>Whenever possible <code>#include</code> directives should be placed in the <code>.cpp</code>
 implementation file rather than the <code>.h</code> header file.</p>
+
+<h2 id="warnings">Warnings</h2>
+
+<p>Code should compile without firing any warnings.</p>
+
+<h3 id="anexception">An Exception</h3>
+
+<p>The only exception is for the use of the deprecated
+<code>Eigen::DynamicSparseMatrix</code> in core sub-routines (e.g. <code>igl::cat</code>). This class
+is still supported and faster than the standard, non-deprecated Eigen
+implementation so we&#8217;re keeping it as long as possible and profitable.</p>
+
+</body>
+</html>

style-guidelines.md

@@ -379,3 +379,14 @@ edited by you first. This means for
 
 Whenever possible `#include` directives should be placed in the `.cpp`
 implementation file rather than the `.h` header file.
+
+## Warnings
+
+Code should compile without firing any warnings.
+
+### An Exception
+
+The only exception is for the use of the deprecated
+`Eigen::DynamicSparseMatrix` in core sub-routines (e.g. `igl::cat`). This class
+is still supported and faster than the standard, non-deprecated Eigen
+implementation so we're keeping it as long as possible and profitable.

tutorial/107_ScreenCapture/CMakeLists.txt

@@ -0,0 +1,8 @@
+cmake_minimum_required(VERSION 2.8.12)
+project(107_ScreenCapture)
+
+add_executable(${PROJECT_NAME}_bin
+  main.cpp)
+target_include_directories(${PROJECT_NAME}_bin PRIVATE ${LIBIGL_INCLUDE_DIRS})
+target_compile_definitions(${PROJECT_NAME}_bin PRIVATE ${LIBIGL_DEFINITIONS})
+target_link_libraries(${PROJECT_NAME}_bin ${LIBIGL_LIBRARIES} ${LIBIGL_EXTRA_LIBRARIES})

tutorial/107_ScreenCapture/main.cpp

@@ -0,0 +1,82 @@
+#include <igl/readOFF.h>
+#include <igl/viewer/Viewer.h>
+#include <iostream>
+#include "tutorial_shared_path.h"
+#include <igl/png/writePNG.h>
+#include <igl/png/readPNG.h>
+
+// This function is called every time a keyboard button is pressed
+bool key_down(igl::viewer::Viewer& viewer, unsigned char key, int modifier)
+{
+  if (key == '1')
+  {
+    // Allocate temporary buffers
+    Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> R(1280,800);
+    Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> G(1280,800);
+    Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> B(1280,800);
+    Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> A(1280,800);
+
+    // Draw the scene in the buffers
+    viewer.core.draw_buffer(viewer.data,viewer.opengl,false,R,G,B,A);
+
+    // Save it to a PNG
+    igl::png::writePNG(R,G,B,A,"out.png");
+  }
+
+  if (key == '2')
+  {
+    // Allocate temporary buffers
+    Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> R,G,B,A;
+
+    // Read the PNG
+    igl::png::readPNG("out.png",R,G,B,A);
+
+    // Replace the mesh with a triangulated square
+    Eigen::MatrixXd V(4,3);
+    V <<
+      -0.5,-0.5,0,
+       0.5,-0.5,0,
+       0.5, 0.5,0,
+      -0.5, 0.5,0;
+    Eigen::MatrixXi F(2,3);
+    F <<
+      0,1,2,
+      2,3,0;
+    Eigen::MatrixXd UV(4,2);
+    UV <<
+      0,0,
+      1,0,
+      1,1,
+      0,1;
+
+    viewer.data.clear();
+    viewer.data.set_mesh(V,F);
+    viewer.data.set_uv(UV);
+    viewer.core.align_camera_center(V);
+    viewer.core.show_texture = true;
+
+    // Use the image as a texture
+    viewer.data.set_texture(R,G,B);
+
+  }
+
+
+  return false;
+}
+
+int main(int argc, char *argv[])
+{
+  // Load a mesh in OFF format
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  igl::readOFF(TUTORIAL_SHARED_PATH "/bunny.off", V, F);
+
+  std::cerr << "Press 1 to render the scene and save it in a png." << std::endl;
+  std::cerr << "Press 2 to load the saved png and use it as a texture." << std::endl;
+
+  // Plot the mesh and register the callback
+  igl::viewer::Viewer viewer;
+  viewer.callback_key_down = &key_down;
+  viewer.data.set_mesh(V, F);
+  viewer.launch();
+}

tutorial/509_Planarization/main.cpp

@@ -25,7 +25,7 @@ Eigen::MatrixXd PQC0plan, PQC1plan, PQC2plan, PQC3plan;
 
 
 // Scale for visualizing the fields
-double global_scale;
+double global_scale; //TODO: not used
 
 
 bool key_down(igl::viewer::Viewer& viewer, unsigned char key, int modifier)

tutorial/607_Picking/main.cpp

@@ -9,8 +9,10 @@ int main(int argc, char *argv[])
   // 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);
+
   // Initialize white
   C = Eigen::MatrixXd::Constant(F.rows(),3,1);
   igl::viewer::Viewer viewer;
@@ -22,15 +24,8 @@ int main(int argc, char *argv[])
     // 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))
+    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;

tutorial/CMakeLists.txt

@@ -89,6 +89,9 @@ if(TUTORIALS_CHAPTER1)
   add_subdirectory("104_Colors")
   add_subdirectory("105_Overlays")
   add_subdirectory("106_ViewerMenu")
+  if(LIBIGL_WITH_PNG)
+    add_subdirectory("107_ScreenCapture")
+  endif()
 endif()
 
 # Chapter 2

tutorial/tutorial.html.REMOVED.git-id

@@ -1 +1 @@
-f46cce260ca12390184cb5fef8205ea121b3cb97
+0ba0c22d44c641478413576281dc8e745316e8c7

tutorial/tutorial.md.REMOVED.git-id

@@ -1 +1 @@
-f87bfabc63b2facf8ce8a778e86028722a7af22f
+6c7992b923fba06ca0986754535c23cb220875da