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

Former-commit-id: 54e518632b681ff5e54eddf46c2cdebd6c40b317

README.md

@@ -195,6 +195,7 @@ Eurographics/ACM Symposium on Geometry Processing software award. Here are a
 few labs/companies/institutions using libigl:
 
  - [Adobe Research](http://www.adobe.com/technology/)  
+ - [Electronic Arts, Inc](http://www.ea.com)
  - [Mesh](http://meshconsultants.ca/), consultants, Canada
  - [Pixar Research](http://graphics.pixar.com/research/)
  - [Spine by Esoteric Software](http://esotericsoftware.com/) is an animation tool dedicated to 2D characters.

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/barycentric_coordinates.cpp

@@ -106,4 +106,5 @@ template void igl::barycentric_coordinates<Eigen::Matrix<double, -1, -1, 0, -1,
 template void igl::barycentric_coordinates<Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&);
 template void igl::barycentric_coordinates<Eigen::Matrix<double, 1, 2, 1, 1, 2>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&);
 template void igl::barycentric_coordinates<Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&);
+template void igl::barycentric_coordinates<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/circulation.cpp

@@ -44,8 +44,8 @@ IGL_INLINE std::vector<int> igl::circulation(
   int ei = e;
   while(true)
   {
-    N.push_back(fi);
     step(ei,fi,ei,fi);
+    N.push_back(fi);
     // back to start?
     if(fi == f0)
     {

include/igl/collapse_edge.cpp

@@ -111,8 +111,8 @@ IGL_INLINE bool igl::collapse_edge(
   // finally, reindex faces and edges incident on d. Do this last so asserts
   // make sense.
   //
-  // Could actually skip first two, since those are always the two collpased
-  // faces.
+  // Could actually skip first and last, since those are always the two
+  // collpased faces.
   for(auto f : nV2Fd)
   {
     for(int v = 0;v<3;v++)

include/igl/copyleft/cgal/outer_element.cpp

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2015 Qingnan Zhou <qnzhou@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 "outer_element.h"
 #include <iostream>
@@ -23,7 +23,7 @@ IGL_INLINE void igl::copyleft::cgal::outer_vertex(
         IndexType & v_index,
         Eigen::PlainObjectBase<DerivedA> & A)
 {
-    // Algorithm: 
+    // Algorithm:
     //    Find an outer vertex (i.e. vertex reachable from infinity)
     //    Return the vertex with the largest X value.
     //    If there is a tie, pick the one with largest Y value.
@@ -112,7 +112,7 @@ IGL_INLINE void igl::copyleft::cgal::outer_edge(
     const ScalarArray3& outer_v = V.row(outer_vid);
     assert(candidate_faces.size() > 0);
 
-    auto get_vertex_index = [&](const IndexArray3& f, Index vid) -> Index 
+    auto get_vertex_index = [&](const IndexArray3& f, Index vid) -> Index
     {
         if (f[0] == vid) return 0;
         if (f[1] == vid) return 1;
@@ -204,80 +204,14 @@ IGL_INLINE void igl::copyleft::cgal::outer_edge(
     std::copy(incident_faces.begin(), incident_faces.end(), A.data());
 }
 
-template<
-    typename DerivedV,
-    typename DerivedF,
-    typename DerivedN,
-    typename DerivedI,
-    typename IndexType
-    >
-IGL_INLINE void igl::copyleft::cgal::outer_facet(
-        const Eigen::PlainObjectBase<DerivedV> & V,
-        const Eigen::PlainObjectBase<DerivedF> & F,
-        const Eigen::PlainObjectBase<DerivedN> & N,
-        const Eigen::PlainObjectBase<DerivedI> & I,
-        IndexType & f,
-        bool & flipped) {
-    // Algorithm:
-    //    Find an outer edge.
-    //    Find the incident facet with the largest absolute X normal component.
-    //    If there is a tie, keep the one with positive X component.
-    //    If there is still a tie, pick the face with the larger signed index
-    //    (flipped face has negative index).
-    typedef typename DerivedV::Scalar Scalar;
-    typedef typename DerivedV::Index Index;
-    const size_t INVALID = std::numeric_limits<size_t>::max();
-
-    Index v1,v2;
-    Eigen::Matrix<Index,Eigen::Dynamic,1> incident_faces;
-    outer_edge(V, F, I, v1, v2, incident_faces);
-    assert(incident_faces.size() > 0);
-
-    auto generic_fabs = [&](const Scalar& val) -> const Scalar {
-        if (val >= 0) return val;
-        else return -val;
-    };
-
-    Scalar max_nx = 0;
-    size_t outer_fid = INVALID;
-    const size_t num_incident_faces = incident_faces.size();
-    for (size_t i=0; i<num_incident_faces; i++) 
-    {
-        const auto& fid = incident_faces(i);
-        const Scalar nx = N(fid, 0);
-        if (outer_fid == INVALID) {
-            max_nx = nx;
-            outer_fid = fid;
-        } else {
-            if (generic_fabs(nx) > generic_fabs(max_nx)) {
-                max_nx = nx;
-                outer_fid = fid;
-            } else if (nx == -max_nx && nx > 0) {
-                max_nx = nx;
-                outer_fid = fid;
-            } else if (nx == max_nx) {
-                if ((max_nx >= 0 && outer_fid < fid) ||
-                    (max_nx <  0 && outer_fid > fid)) {
-                    max_nx = nx;
-                    outer_fid = fid;
-                }
-            }
-        }
-    }
-
-    assert(outer_fid != INVALID);
-    f = outer_fid;
-    flipped = max_nx < 0;
-}
 
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
-template void igl::copyleft::cgal::outer_facet<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<long, -1, 1, 0, -1, 1>, int>(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<long, -1, 1, 0, -1, 1> > const&, int&, bool&);
-template void igl::copyleft::cgal::outer_facet<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 1, -1, -1>, unsigned long>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 1, -1, -1> > const&, unsigned long&, bool&);
-template void igl::copyleft::cgal::outer_facet<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, int>(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, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int&, bool&);
-template void igl::copyleft::cgal::outer_facet<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<int, -1, 1, 0, -1, 1>, int>(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<int, -1, 1, 0, -1, 1> > const&, int&, bool&);
 template void igl::copyleft::cgal::outer_edge<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, long, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> > const&, long&, long&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
 template void igl::copyleft::cgal::outer_edge<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, long, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, long&, long&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
+template void igl::copyleft::cgal::outer_edge<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, long, Eigen::Matrix<long, -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<long, -1, 1, 0, -1, 1> > const&, long&, long&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
+template void igl::copyleft::cgal::outer_edge<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, long, Eigen::Matrix<long, -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> > const&, long&, long&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
+template void igl::copyleft::cgal::outer_edge<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 1, -1, -1>, long, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 1, -1, -1> > const&, long&, long&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/copyleft/cgal/outer_element.h

@@ -73,38 +73,6 @@ namespace igl
               IndexType & v2,
               Eigen::PlainObjectBase<DerivedA> & A);
 
-
-      // Find a facet that is reachable from infinity without crossing any faces.
-      // Such facet is called "outer facet."
-      //
-      // Precondition: The input mesh must have all self-intersection resolved.
-      // I.e there is no duplicated vertices, no overlapping edge and no
-      // intersecting faces (the only exception is there could be topologically
-      // duplicated faces).  See cgal::remesh_self_intersections.h for how to
-      // obtain such input.
-      //
-      // Inputs:
-      //   V  #V by 3 list of vertex positions
-      //   F  #F by 3 list of triangle indices into V
-      //   N  #N by 3 list of face normals
-      //   I  #I list of facets to consider
-      // Outputs:
-      //   f  Index of the outer facet.
-      //   flipped  true iff the normal of f points inwards.
-      template<
-          typename DerivedV,
-          typename DerivedF,
-          typename DerivedN,
-          typename DerivedI,
-          typename IndexType
-          >
-      IGL_INLINE void outer_facet(
-              const Eigen::PlainObjectBase<DerivedV> & V,
-              const Eigen::PlainObjectBase<DerivedF> & F,
-              const Eigen::PlainObjectBase<DerivedN> & N,
-              const Eigen::PlainObjectBase<DerivedI> & I,
-              IndexType & f,
-              bool & flipped);
     }
   }
 }

include/igl/copyleft/cgal/outer_facet.cpp

@@ -79,10 +79,82 @@ IGL_INLINE void igl::copyleft::cgal::outer_facet(
     flipped = adj_faces[order[0]] > 0;
 }
 
+
+
+template<
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedN,
+    typename DerivedI,
+    typename IndexType
+    >
+IGL_INLINE void igl::copyleft::cgal::outer_facet(
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::PlainObjectBase<DerivedF> & F,
+        const Eigen::PlainObjectBase<DerivedN> & N,
+        const Eigen::PlainObjectBase<DerivedI> & I,
+        IndexType & f,
+        bool & flipped) {
+    // Algorithm:
+    //    Find an outer edge.
+    //    Find the incident facet with the largest absolute X normal component.
+    //    If there is a tie, keep the one with positive X component.
+    //    If there is still a tie, pick the face with the larger signed index
+    //    (flipped face has negative index).
+    typedef typename DerivedV::Scalar Scalar;
+    typedef typename DerivedV::Index Index;
+    const size_t INVALID = std::numeric_limits<size_t>::max();
+
+    Index v1,v2;
+    Eigen::Matrix<Index,Eigen::Dynamic,1> incident_faces;
+    outer_edge(V, F, I, v1, v2, incident_faces);
+    assert(incident_faces.size() > 0);
+
+    auto generic_fabs = [&](const Scalar& val) -> const Scalar {
+        if (val >= 0) return val;
+        else return -val;
+    };
+
+    Scalar max_nx = 0;
+    size_t outer_fid = INVALID;
+    const size_t num_incident_faces = incident_faces.size();
+    for (size_t i=0; i<num_incident_faces; i++) 
+    {
+        const auto& fid = incident_faces(i);
+        const Scalar nx = N(fid, 0);
+        if (outer_fid == INVALID) {
+            max_nx = nx;
+            outer_fid = fid;
+        } else {
+            if (generic_fabs(nx) > generic_fabs(max_nx)) {
+                max_nx = nx;
+                outer_fid = fid;
+            } else if (nx == -max_nx && nx > 0) {
+                max_nx = nx;
+                outer_fid = fid;
+            } else if (nx == max_nx) {
+                if ((max_nx >= 0 && outer_fid < fid) ||
+                    (max_nx <  0 && outer_fid > fid)) {
+                    max_nx = nx;
+                    outer_fid = fid;
+                }
+            }
+        }
+    }
+
+    assert(outer_fid != INVALID);
+    f = outer_fid;
+    flipped = max_nx < 0;
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+template void igl::copyleft::cgal::outer_facet<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<long, -1, 1, 0, -1, 1>, int>(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<long, -1, 1, 0, -1, 1> > const&, int&, bool&);
+template void igl::copyleft::cgal::outer_facet<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 1, -1, -1>, unsigned long>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 1, -1, -1> > const&, unsigned long&, bool&);
+template void igl::copyleft::cgal::outer_facet<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, int>(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, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int&, bool&);
+template void igl::copyleft::cgal::outer_facet<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<int, -1, 1, 0, -1, 1>, int>(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<int, -1, 1, 0, -1, 1> > const&, int&, bool&);
 template void igl::copyleft::cgal::outer_facet<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, int>(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> > const&, int&, bool&);
 template void igl::copyleft::cgal::outer_facet<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, int>(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<long, -1, 1, 0, -1, 1> > const&, int&, bool&);
 template void igl::copyleft::cgal::outer_facet<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, int>(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> > const&, int&, bool&);

include/igl/copyleft/cgal/outer_facet.h

@@ -47,6 +47,39 @@ namespace igl
                 const Eigen::PlainObjectBase<DerivedI> & I,
                 IndexType & f,
                 bool & flipped);
+
+      // Find a facet that is reachable from infinity without crossing any faces.
+      // Such facet is called "outer facet."
+      //
+      // Precondition: The input mesh must have all self-intersection resolved.
+      // I.e there is no duplicated vertices, no overlapping edge and no
+      // intersecting faces (the only exception is there could be topologically
+      // duplicated faces).  See cgal::remesh_self_intersections.h for how to
+      // obtain such input.
+      //
+      // Inputs:
+      //   V  #V by 3 list of vertex positions
+      //   F  #F by 3 list of triangle indices into V
+      //   N  #N by 3 list of face normals
+      //   I  #I list of facets to consider
+      // Outputs:
+      //   f  Index of the outer facet.
+      //   flipped  true iff the normal of f points inwards.
+      template<
+          typename DerivedV,
+          typename DerivedF,
+          typename DerivedN,
+          typename DerivedI,
+          typename IndexType
+          >
+      IGL_INLINE void outer_facet(
+              const Eigen::PlainObjectBase<DerivedV> & V,
+              const Eigen::PlainObjectBase<DerivedF> & F,
+              const Eigen::PlainObjectBase<DerivedN> & N,
+              const Eigen::PlainObjectBase<DerivedI> & I,
+              IndexType & f,
+              bool & flipped);
+
     }
 
   }

include/igl/copyleft/progressive_hulls.cpp

@@ -14,7 +14,8 @@ IGL_INLINE bool igl::copyleft::progressive_hulls(
   const Eigen::MatrixXi & F,
   const size_t max_m,
   Eigen::MatrixXd & U,
-  Eigen::MatrixXi & G)
+  Eigen::MatrixXi & G,
+  Eigen::VectorXi & J)
 {
   int m = F.rows();
   return decimate(
@@ -23,5 +24,6 @@ IGL_INLINE bool igl::copyleft::progressive_hulls(
     progressive_hulls_cost_and_placement,
     max_faces_stopping_condition(m,max_m),
     U,
-    G);
+    G,
+    J);
 }

include/igl/copyleft/progressive_hulls.h

@@ -26,13 +26,15 @@ namespace igl
     // Outputs:
     //   U  #U by dim list of output vertex posistions (can be same ref as V)
     //   G  #G by 3 list of output face indices into U (can be same ref as G)
+    //   J  #G list of indices into F of birth faces
     // Returns true if m was reached (otherwise #G > m)
     IGL_INLINE bool progressive_hulls(
       const Eigen::MatrixXd & V,
       const Eigen::MatrixXi & F,
       const size_t max_m,
       Eigen::MatrixXd & U,
-      Eigen::MatrixXi & G);
+      Eigen::MatrixXi & G,
+      Eigen::VectorXi & J);
   }
 }
 #ifndef IGL_STATIC_LIBRARY

include/igl/decimate.cpp

@@ -17,7 +17,8 @@ IGL_INLINE bool igl::decimate(
   const Eigen::MatrixXi & F,
   const size_t max_m,
   Eigen::MatrixXd & U,
-  Eigen::MatrixXi & G)
+  Eigen::MatrixXi & G,
+  Eigen::VectorXi & J)
 {
   int m = F.rows();
   const auto & shortest_edge_and_midpoint = [](
@@ -40,7 +41,8 @@ IGL_INLINE bool igl::decimate(
     shortest_edge_and_midpoint,
     max_faces_stopping_condition(m,max_m),
     U,
-    G);
+    G,
+    J);
 }
 
 IGL_INLINE bool igl::decimate(
@@ -72,7 +74,8 @@ IGL_INLINE bool igl::decimate(
       const int,
       const int)> & stopping_condition,
   Eigen::MatrixXd & U,
-  Eigen::MatrixXi & G)
+  Eigen::MatrixXi & G,
+  Eigen::VectorXi & J)
 {
   using namespace Eigen;
   using namespace std;
@@ -132,6 +135,7 @@ IGL_INLINE bool igl::decimate(
   }
   // remove all IGL_COLLAPSE_EDGE_NULL faces
   MatrixXi F2(F.rows(),3);
+  J.resize(F.rows());
   int m = 0;
   for(int f = 0;f<F.rows();f++)
   {
@@ -140,10 +144,13 @@ IGL_INLINE bool igl::decimate(
       F(f,1) != IGL_COLLAPSE_EDGE_NULL || 
       F(f,2) != IGL_COLLAPSE_EDGE_NULL)
     {
-      F2.row(m++) = F.row(f);
+      F2.row(m) = F.row(f);
+      J(m) = f;
+      m++;
     }
   }
   F2.conservativeResize(m,F2.cols());
+  J.conservativeResize(m);
   VectorXi _1;
   remove_unreferenced(V,F2,U,G,_1);
   return clean_finish;

include/igl/decimate.h

@@ -25,13 +25,15 @@ namespace igl
   // Outputs:
   //   U  #U by dim list of output vertex posistions (can be same ref as V)
   //   G  #G by 3 list of output face indices into U (can be same ref as G)
+  //   J  #G list of indices into F of birth face
   // Returns true if m was reached (otherwise #G > m)
   IGL_INLINE bool decimate(
     const Eigen::MatrixXd & V,
     const Eigen::MatrixXi & F,
     const size_t max_m,
     Eigen::MatrixXd & U,
-    Eigen::MatrixXi & G);
+    Eigen::MatrixXi & G,
+    Eigen::VectorXi & J);
   // Inputs:
   //   cost_and_placement  function computing cost of collapsing an edge and 3d
   //     position where it should be placed:
@@ -70,7 +72,8 @@ namespace igl
       const int,
       const int)> & stopping_condition,
     Eigen::MatrixXd & U,
-    Eigen::MatrixXi & G);
+    Eigen::MatrixXi & G,
+    Eigen::VectorXi & J);
 
 }
 

include/igl/in_element.h

@@ -11,6 +11,7 @@
 #include "igl_inline.h"
 #include "AABB.h"
 #include <Eigen/Core>
+#include <Eigen/Sparse>
 
 namespace igl
 {

include/igl/infinite_cost_stopping_condition.cpp

@@ -0,0 +1,108 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "infinite_cost_stopping_condition.h"
+
+IGL_INLINE void igl::infinite_cost_stopping_condition(
+  const std::function<void(
+    const int,
+    const Eigen::MatrixXd &,
+    const Eigen::MatrixXi &,
+    const Eigen::MatrixXi &,
+    const Eigen::VectorXi &,
+    const Eigen::MatrixXi &,
+    const Eigen::MatrixXi &,
+    double &,
+    Eigen::RowVectorXd &)> & cost_and_placement,
+  std::function<bool(
+    const Eigen::MatrixXd &,
+    const Eigen::MatrixXi &,
+    const Eigen::MatrixXi &,
+    const Eigen::VectorXi &,
+    const Eigen::MatrixXi &,
+    const Eigen::MatrixXi &,
+    const std::set<std::pair<double,int> > &,
+    const std::vector<std::set<std::pair<double,int> >::iterator > &,
+    const Eigen::MatrixXd &,
+    const int,
+    const int,
+    const int,
+    const int,
+    const int)> & stopping_condition)
+{
+  stopping_condition = 
+    [&cost_and_placement]
+    (
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    const Eigen::MatrixXi & E,
+    const Eigen::VectorXi & EMAP,
+    const Eigen::MatrixXi & EF,
+    const Eigen::MatrixXi & EI,
+    const std::set<std::pair<double,int> > & Q,
+    const std::vector<std::set<std::pair<double,int> >::iterator > & Qit,
+    const Eigen::MatrixXd & C,
+    const int e,
+    const int /*e1*/,
+    const int /*e2*/,
+    const int /*f1*/,
+    const int /*f2*/)->bool
+    {
+      Eigen::RowVectorXd p;
+      double cost;
+      cost_and_placement(e,V,F,E,EMAP,EF,EI,cost,p);
+      return isinf(cost);
+    };
+}
+
+IGL_INLINE 
+  std::function<bool(
+    const Eigen::MatrixXd &,
+    const Eigen::MatrixXi &,
+    const Eigen::MatrixXi &,
+    const Eigen::VectorXi &,
+    const Eigen::MatrixXi &,
+    const Eigen::MatrixXi &,
+    const std::set<std::pair<double,int> > &,
+    const std::vector<std::set<std::pair<double,int> >::iterator > &,
+    const Eigen::MatrixXd &,
+    const int,
+    const int,
+    const int,
+    const int,
+    const int)> 
+  igl::infinite_cost_stopping_condition(
+    const std::function<void(
+      const int,
+      const Eigen::MatrixXd &,
+      const Eigen::MatrixXi &,
+      const Eigen::MatrixXi &,
+      const Eigen::VectorXi &,
+      const Eigen::MatrixXi &,
+      const Eigen::MatrixXi &,
+      double &,
+      Eigen::RowVectorXd &)> & cost_and_placement)
+{
+  std::function<bool(
+    const Eigen::MatrixXd &,
+    const Eigen::MatrixXi &,
+    const Eigen::MatrixXi &,
+    const Eigen::VectorXi &,
+    const Eigen::MatrixXi &,
+    const Eigen::MatrixXi &,
+    const std::set<std::pair<double,int> > &,
+    const std::vector<std::set<std::pair<double,int> >::iterator > &,
+    const Eigen::MatrixXd &,
+    const int,
+    const int,
+    const int,
+    const int,
+    const int)> stopping_condition;
+  infinite_cost_stopping_condition(cost_and_placement,stopping_condition);
+  return stopping_condition;
+}
+

include/igl/infinite_cost_stopping_condition.h

@@ -0,0 +1,85 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_INFINITE_COST_STOPPING_CONDITION_H
+#define IGL_INFINITE_COST_STOPPING_CONDITION_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+#include <set>
+#include <functional>
+namespace igl
+{
+  // Stopping condition function compatible with igl::decimate. The output
+  // function handle will return true if cost of next edge is infinite.
+  //
+  // Inputs:
+  //   cost_and_placement  handle being used by igl::collapse_edge
+  // Outputs:
+  //   stopping_condition
+  //
+  IGL_INLINE void infinite_cost_stopping_condition(
+    const std::function<void(
+      const int,
+      const Eigen::MatrixXd &,
+      const Eigen::MatrixXi &,
+      const Eigen::MatrixXi &,
+      const Eigen::VectorXi &,
+      const Eigen::MatrixXi &,
+      const Eigen::MatrixXi &,
+      double &,
+      Eigen::RowVectorXd &)> & cost_and_placement,
+    std::function<bool(
+      const Eigen::MatrixXd &,
+      const Eigen::MatrixXi &,
+      const Eigen::MatrixXi &,
+      const Eigen::VectorXi &,
+      const Eigen::MatrixXi &,
+      const Eigen::MatrixXi &,
+      const std::set<std::pair<double,int> > &,
+      const std::vector<std::set<std::pair<double,int> >::iterator > &,
+      const Eigen::MatrixXd &,
+      const int,
+      const int,
+      const int,
+      const int,
+      const int)> & stopping_condition);
+  IGL_INLINE 
+    std::function<bool(
+      const Eigen::MatrixXd &,
+      const Eigen::MatrixXi &,
+      const Eigen::MatrixXi &,
+      const Eigen::VectorXi &,
+      const Eigen::MatrixXi &,
+      const Eigen::MatrixXi &,
+      const std::set<std::pair<double,int> > &,
+      const std::vector<std::set<std::pair<double,int> >::iterator > &,
+      const Eigen::MatrixXd &,
+      const int,
+      const int,
+      const int,
+      const int,
+      const int)> 
+    infinite_cost_stopping_condition(
+      const std::function<void(
+        const int,
+        const Eigen::MatrixXd &,
+        const Eigen::MatrixXi &,
+        const Eigen::MatrixXi &,
+        const Eigen::VectorXi &,
+        const Eigen::MatrixXi &,
+        const Eigen::MatrixXi &,
+        double &,
+        Eigen::RowVectorXd &)> & cost_and_placement);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "infinite_cost_stopping_condition.cpp"
+#endif
+#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/parula.cpp

@@ -66,6 +66,7 @@ IGL_INLINE void igl::parula(
 // generated by autoexplicit.sh
 template void igl::parula<Eigen::Array<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Array<int, -1, 1, 0, -1, 1> > const&, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::parula<double>(double, double*);
+template void igl::parula<double>(double, double&, double&, double&);
 template void igl::parula<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&, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::parula<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&, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/parula.h

@@ -44,8 +44,8 @@ namespace igl
     const bool normalize,
     Eigen::PlainObjectBase<DerivedC> & C);
   // Inputs:
-  //   min_z  value at blue
-  //   max_z  value at red
+  //   min_Z  value at blue
+  //   max_Z  value at red
   template <typename DerivedZ, typename DerivedC>
   IGL_INLINE void parula(
     const Eigen::PlainObjectBase<DerivedZ> & Z,

include/igl/per_edge_normals.cpp

@@ -116,4 +116,7 @@ IGL_INLINE void igl::per_edge_normals(
 // Explicit template specialization
 template void igl::per_edge_normals<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<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<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::per_edge_normals<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<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&, igl::PerEdgeNormalsWeightingType, 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<int, -1, 1, 0, -1, 1> >&);
+template void igl::per_edge_normals<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<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&, igl::PerEdgeNormalsWeightingType, 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<int, -1, -1, 0, -1, -1> >&);
+template void igl::per_edge_normals<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<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&, igl::PerEdgeNormalsWeightingType, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::per_edge_normals<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<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<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/per_edge_normals.h

@@ -25,7 +25,8 @@ namespace igl
   // Inputs:
   //   V  #V by 3 eigen Matrix of mesh vertex 3D positions
   //   F  #F by 3 eigen Matrix of face (triangle) indices
-  //   weighting  weighting type
+  //   weight  weighting type
+  //   FN  #F by 3 matrix of 3D face normals per face
   // Output:
   //   N  #2 by 3 matrix of mesh edge 3D normals per row
   //   E  #E by 2 matrix of edge indices per row

include/igl/point_mesh_squared_distance.cpp

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2014 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 "point_mesh_squared_distance.h"
 #include "AABB.h"
@@ -17,7 +17,7 @@ template <
 IGL_INLINE void igl::point_mesh_squared_distance(
   const Eigen::PlainObjectBase<DerivedP> & P,
   const Eigen::PlainObjectBase<DerivedV> & V,
-  const Eigen::MatrixXi & Ele, 
+  const Eigen::MatrixXi & Ele,
   Eigen::PlainObjectBase<DerivedsqrD> & sqrD,
   Eigen::PlainObjectBase<DerivedI> & I,
   Eigen::PlainObjectBase<DerivedC> & C)
@@ -46,5 +46,6 @@ IGL_INLINE void igl::point_mesh_squared_distance(
 }
 
 #ifdef IGL_STATIC_LIBRARY
+template void igl::point_mesh_squared_distance<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<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::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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::point_mesh_squared_distance<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<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::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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/point_mesh_squared_distance.h

@@ -10,6 +10,7 @@
 #include "igl_inline.h"
 #include <Eigen/Core>
 #include <vector>
+
 namespace igl
 {
   // Compute distances from a set of points P to a triangle mesh (V,F)
@@ -17,7 +18,7 @@ namespace igl
   // Inputs:
   //   P  #P by 3 list of query point positions
   //   V  #V by 3 list of vertex positions
-  //   F  #F by (3|2|1) list of (triangle|edge|point) indices
+  //   Ele  #Ele by (3|2|1) list of (triangle|edge|point) indices
   // Outputs:
   //   sqrD  #P list of smallest squared distances
   //   I  #P list of primitive indices corresponding to smallest distances

include/igl/polar_svd.cpp

@@ -75,4 +75,6 @@ template void igl::polar_svd<Eigen::Matrix<float,3,3,0,3,3>,Eigen::Matrix<float,
 template void igl::polar_svd<Eigen::Matrix<double, 3, 3, 0, 3, 3>, Eigen::Matrix<double, 3, 3, 0, 3, 3>, Eigen::Matrix<double, 3, 3, 0, 3, 3>, Eigen::Matrix<double, 3, 3, 0, 3, 3>, Eigen::Matrix<double, 3, 1, 0, 3, 1>, Eigen::Matrix<double, 3, 3, 0, 3, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 3, 0, 3, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 3, 0, 3, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 3, 0, 3, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 3, 0, 3, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 1, 0, 3, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 3, 0, 3, 3> >&);
 template void igl::polar_svd<Eigen::Matrix<float, 2, 2, 0, 2, 2>, Eigen::Matrix<float, 2, 2, 0, 2, 2>, Eigen::Matrix<float, 2, 2, 0, 2, 2>, Eigen::Matrix<float, 2, 2, 0, 2, 2>, Eigen::Matrix<float, 2, 1, 0, 2, 1>, Eigen::Matrix<float, 2, 2, 0, 2, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<float, 2, 2, 0, 2, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, 2, 2, 0, 2, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<float, 2, 2, 0, 2, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<float, 2, 2, 0, 2, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<float, 2, 1, 0, 2, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<float, 2, 2, 0, 2, 2> >&);
 template void igl::polar_svd<Eigen::Matrix<double, 2, 2, 0, 2, 2>, Eigen::Matrix<double, 2, 2, 0, 2, 2>, Eigen::Matrix<double, 2, 2, 0, 2, 2>, Eigen::Matrix<double, 2, 2, 0, 2, 2>, Eigen::Matrix<double, 2, 1, 0, 2, 1>, Eigen::Matrix<double, 2, 2, 0, 2, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 2, 2, 0, 2, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 2, 2, 0, 2, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, 2, 2, 0, 2, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, 2, 2, 0, 2, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, 2, 1, 0, 2, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, 2, 2, 0, 2, 2> >&);
+template void igl::polar_svd<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::polar_svd<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/readMESH.cpp

@@ -470,4 +470,6 @@ template bool igl::readMESH<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix
 template bool igl::readMESH<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
 // generated by autoexplicit.sh
 template bool igl::readMESH<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template bool igl::readMESH<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template bool igl::readMESH<double, int>(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&);
 #endif

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/simplify_polyhedron.cpp

@@ -0,0 +1,103 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "simplify_polyhedron.h"
+#include "decimate.h"
+#include "circulation.h"
+#include "per_face_normals.h"
+#include "infinite_cost_stopping_condition.h"
+#include <functional>
+
+IGL_INLINE void igl::simplify_polyhedron(
+  const Eigen::MatrixXd & OV,
+  const Eigen::MatrixXi & OF,
+  Eigen::MatrixXd & V,
+  Eigen::MatrixXi & F,
+  Eigen::VectorXi & J)
+{
+  // TODO: to generalize to open meshes, 0-cost should keep all incident
+  // boundary edges on their original lines. (for non-manifold meshes,
+  // igl::decimate needs to be generalized)
+
+  Eigen::MatrixXd N;
+  // Function for computing cost of collapsing edge (0 if at least one
+  // direction doesn't change pointset, inf otherwise) and placement (in lowest
+  // cost direction).
+  const auto & perfect= [&N](
+    const int e,
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    const Eigen::MatrixXi & E,
+    const Eigen::VectorXi & EMAP,
+    const Eigen::MatrixXi & EF,
+    const Eigen::MatrixXi & EI,
+    double & cost,
+    Eigen::RowVectorXd & p)
+  {
+    // Function for ocmputing cost (0 or inf) of collapsing edge by placing
+    // vertex at `positive` end of edge.
+    const auto & perfect_directed = [&N](
+      const int e,
+      const bool positive,
+      const Eigen::MatrixXd & V,
+      const Eigen::MatrixXi & F,
+      const Eigen::MatrixXi & E,
+      const Eigen::VectorXi & EMAP,
+      const Eigen::MatrixXi & EF,
+      const Eigen::MatrixXi & EI,
+      double & cost,
+      Eigen::RowVectorXd & p)
+    {
+      const auto vi = E(e,positive);
+      const auto vj = E(e,!positive);
+      p = V.row(vj);
+      std::vector<int> faces = igl::circulation(e,positive,F,E,EMAP,EF,EI);
+      cost = 0;
+      for(auto f : faces)
+      {
+        // Skip the faces being collapsed
+        if(f == EF(e,0) || f == EF(e,1))
+        {
+          continue;
+        }
+        const Eigen::RowVectorXd nbefore = N.row(f);
+        // Face with vi replaced with vj
+        const Eigen::RowVector3i fafter(
+            F(f,0) == vi ? vj : F(f,0),
+            F(f,1) == vi ? vj : F(f,1),
+            F(f,2) == vi ? vj : F(f,2));
+        Eigen::RowVectorXd nafter;
+        igl::per_face_normals(V,fafter,nafter);
+        const double epsilon = 1e-10;
+        // if normal changed then not feasible, break
+        if((nbefore-nafter).norm() > epsilon)
+        {
+          cost = std::numeric_limits<double>::infinity();
+          break;
+        }
+      }
+    }; 
+    p.resize(3);
+    double cost0, cost1;
+    Eigen::RowVectorXd p0, p1;
+    perfect_directed(e,false,V,F,E,EMAP,EF,EI,cost0,p0);
+    perfect_directed(e,true,V,F,E,EMAP,EF,EI,cost1,p1);
+    if(cost0 < cost1)
+    {
+      cost = cost0;
+      p = p0;
+    }else
+    {
+      cost = cost1;
+      p = p1;
+    }
+  };
+  igl::per_face_normals(OV,OF,N);
+  igl::decimate(
+    OV,OF,perfect,igl::infinite_cost_stopping_condition(perfect),V,F,J);
+}
+

include/igl/simplify_polyhedron.h

@@ -0,0 +1,37 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_SIMPLIFY_POLYHEDRON_H
+#define IGL_SIMPLIFY_POLYHEDRON_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Simplify a polyhedron represented as a triangle mesh (OV,OF) by collapsing
+  // any edge that doesn't contribute to defining surface's pointset. This
+  // _would_ also make sense for open and non-manifold meshes, but the current
+  // implementation only works with closed manifold surfaces with well defined
+  // triangle normals.
+  //
+  // Inputs:
+  //   OV  #OV by 3 list of input mesh vertex positions
+  //   OF  #OF by 3 list of input mesh triangle indices into OV
+  // Outputs:
+  //   V  #V by 3 list of output mesh vertex positions
+  //   F  #F by 3 list of input mesh triangle indices into V
+  //   J  #F list of indices into OF of birth parents
+  IGL_INLINE void simplify_polyhedron(
+    const Eigen::MatrixXd & OV,
+    const Eigen::MatrixXi & OF,
+    Eigen::MatrixXd & V,
+    Eigen::MatrixXi & F,
+    Eigen::VectorXi & J);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "simplify_polyhedron.cpp"
+#endif
+#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/slice_tets.cpp

@@ -226,4 +226,5 @@ IGL_INLINE void igl::slice_tets(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 template void igl::slice_tets<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, 4, 1, 0, 4, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, double>(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, 4, 1, 0, 4, 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<int, -1, 1, 0, -1, 1> >&, Eigen::SparseMatrix<double, 0, int>&);
+template void igl::slice_tets<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<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, double>(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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::SparseMatrix<double, 0, int>&);
 #endif

include/igl/unique.cpp

@@ -267,4 +267,5 @@ template void igl::unique<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<i
 template void igl::unique<double>(std::vector<double, std::allocator<double> > const&, std::vector<double, std::allocator<double> >&);
 template void igl::unique<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> >&);
 template void igl::unique<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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::unique_rows<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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/unique_simplices.cpp

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

include/igl/viewer/Viewer.cpp

@@ -6,6 +6,8 @@
 // v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 
+// Must defined this before including Viewer.h
+#define IGL_VIEWER_VIEWER_CPP
 #include "Viewer.h"
 
 #ifdef _WIN32
@@ -273,6 +275,33 @@ namespace viewer
 
   IGL_INLINE Viewer::Viewer()
   {
+    // This mess is to help debug problems arising when compiling
+    // libiglviewer.a with(without) IGL_STATIC_LIBRARY defined and including
+    // Viewer.h without(with) IGL_STATIC_LIBRARY defined.
+#ifdef IGL_STATIC_LIBRARY
+    std::cout<<"igl_with_nanogui_defined_consistently: "<<igl_with_nanogui_defined_consistently()<<std::endl;
+    std::cout<<"igl_with_nanogui_defined_at_compile: "<<  igl_with_nanogui_defined_at_compile()  <<std::endl;
+    std::cout<<"igl_with_nanogui_defined_at_include: "<<  igl_with_nanogui_defined_at_include()  <<std::endl;
+    // First try to first assert
+    assert(igl_with_nanogui_defined_consistently() && 
+      "Must compile and include with IGL_VIEWER_WITH_NANOGUI defined consistently");
+#ifdef NDEBUG
+    // Second print warning since it's hopeless that this will run if wrong.
+    if(!igl_with_nanogui_defined_consistently())
+    {
+      std::cerr<<
+        "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<std::endl<<std::endl<<
+        "WARNING: Seems that IGL_WITH_NANOGUI " <<
+        (igl_with_nanogui_defined_at_compile() ? "was" : "was not") <<
+        " defined"<<std::endl<<"during compilation of Viewer.cpp and "<<
+        (igl_with_nanogui_defined_at_include() ? "was" : "was not") <<
+        " defined"<<std::endl<<"during inclusion of Viewer.h"<<std::endl <<
+        "You're about to get some nasty memory errors."<<std::endl<<std::endl<<
+        "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<std::endl;
+    }
+#endif
+#endif
+
 #ifdef IGL_VIEWER_WITH_NANOGUI
     ngui = nullptr;
     screen = nullptr;
@@ -317,9 +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)"
+#endif
+);
     std::cout<<usage;
 #endif
   }
@@ -530,6 +563,14 @@ namespace viewer
         }
         return true;
       }
+#ifdef IGL_VIEWER_WITH_NANOGUI
+      case ';':
+        core.show_vertid = !core.show_vertid;
+        return true;
+      case ':':
+        core.show_faceid = !core.show_faceid;
+        return true;
+#endif
       default: break;//do nothing
     }
     return false;
@@ -1003,5 +1044,20 @@ namespace viewer
     launch_shut();
     return EXIT_SUCCESS;
   }
+  IGL_INLINE bool Viewer::igl_with_nanogui_defined_at_compile()
+  {
+#ifdef IGL_VIEWER_WITH_NANOGUI
+    return true;
+#else
+    return false;
+#endif
+  }
+  IGL_INLINE bool Viewer::igl_with_nanogui_defined_consistently()
+  {
+    return 
+      igl_with_nanogui_defined_at_compile() == 
+      igl_with_nanogui_defined_at_include();
+  }
 } // end namespace
 }
+

include/igl/viewer/Viewer.h

@@ -153,10 +153,28 @@ namespace viewer
     void* callback_key_down_data;
     void* callback_key_up_data;
 
+  public:
+
+    static IGL_INLINE bool igl_with_nanogui_defined_at_compile();
+    static IGL_INLINE bool igl_with_nanogui_defined_consistently();
+
   public:
       EIGEN_MAKE_ALIGNED_OPERATOR_NEW
   };
 
+  bool igl_with_nanogui_defined_at_include();
+#ifndef IGL_VIEWER_VIEWER_CPP
+  bool igl_with_nanogui_defined_at_include()
+  {
+    // this must be inlined here.
+#ifdef IGL_VIEWER_WITH_NANOGUI
+    return true;
+#else
+    return false;
+#endif
+  }
+#endif
+
 } // end namespace
 } // end namespace
 

include/igl/writeMESH.cpp

@@ -8,7 +8,7 @@
 #include "writeMESH.h"
 
 #include "verbose.h"
-#include "matrix_to_list.h"
+#include "list_to_matrix.h"
 #include <Eigen/Core>
 
 #include <iostream>
@@ -140,4 +140,7 @@ IGL_INLINE bool igl::writeMESH(
 template bool igl::writeMESH<Eigen::Matrix<double, 8, 3, 0, 8, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, 12, 3, 0, 12, 3> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<double, 8, 3, 0, 8, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, 12, 3, 0, 12, 3> > const&);
 //template bool igl::writeMESH<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&);
 template bool igl::writeMESH<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, 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> > const&);
+
+template bool igl::writeMESH<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&);
+template bool igl::writeMESH<double, int>(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&);
 #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,34 @@ 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://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,7 +262,7 @@ page</a>.</p>
 
 <h2 id="copyright">Copyright</h2>
 
-<p>2015 Alec Jacobson, Daniele Panozzo, Christian Schüller, Olga Diamanti, Qingnan
+<p>2016 Alec Jacobson, Daniele Panozzo, Christian Schüller, Olga Diamanti, Qingnan
 Zhou, Nico Pietroni, Stefan Brugger, Kenshi Takayama, Wenzel Jakob, Nikolas De
 Giorgis, Luigi Rocca, Leonardo Sacht, Kevin Walliman, Olga Sorkine-Hornung, and others.</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/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,40 @@ 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_cotmatrix = R"igl_Qu8mg5v7(// Constructs the cotangent stiffness matrix (discrete laplacian) for a given
   // mesh (V,F).
   //
@@ -150,6 +184,21 @@ const char *__doc_igl_jet = R"igl_Qu8mg5v7(// JET like MATLAB's jet
   //   r  red value
   //   g  green value
   //   b  blue value)igl_Qu8mg5v7";
+const char *__doc_igl_cat = R"igl_Qu8mg5v7(// Perform concatenation of a two matrices along a single dimension
+  // If dim == 1, then C = [A;B]. If dim == 2 then C = [A B]
+  // 
+  // Template:
+  //   Scalar  scalar data type for sparse matrices like double or int
+  //   Mat  matrix type for all matrices (e.g. MatrixXd, SparseMatrix)
+  //   MatC  matrix type for ouput matrix (e.g. MatrixXd) needs to support
+  //     resize
+  // Inputs:
+  //   A  first input matrix
+  //   B  second input matrix
+  //   dim  dimension along which to concatenate, 0 or 1
+  // Outputs:
+  //   C  output matrix
+  //   )igl_Qu8mg5v7";
 const char *__doc_igl_eigs = R"igl_Qu8mg5v7(See eigs for the documentation.)igl_Qu8mg5v7";
 const char *__doc_igl_per_corner_normals = R"igl_Qu8mg5v7(// Compute vertex normals via vertex position list, face list
   // Inputs:
@@ -198,6 +247,29 @@ const char *__doc_igl_colon = R"igl_Qu8mg5v7(// Colon operator like matlab's col
   //     than hi, vice versa if hi<low
   // Output:
   //   I  list of values from low to hi with step size step)igl_Qu8mg5v7";
+const char *__doc_igl_fit_rotations = R"igl_Qu8mg5v7(// Known issues: This seems to be implemented in Eigen/Geometry:
+  // Eigen::umeyama
+  //
+  // FIT_ROTATIONS Given an input mesh and new positions find rotations for
+  // every covariance matrix in a stack of covariance matrices
+  // 
+  // Inputs:
+  //   S  nr*dim by dim stack of covariance matrices
+  //   single_precision  whether to use single precision (faster)
+  // Outputs:
+  //   R  dim by dim * nr list of rotations
+  //)igl_Qu8mg5v7";
+const char *__doc_igl_fit_rotations_planar = R"igl_Qu8mg5v7(// FIT_ROTATIONS Given an input mesh and new positions find 2D rotations for
+  // every vertex that best maps its one ring to the new one ring
+  // 
+  // Inputs:
+  //   S  nr*dim by dim stack of covariance matrices, third column and every
+  //   third row will be ignored
+  // Outputs:
+  //   R  dim by dim * nr list of rotations, third row and third column of each
+  //   rotation will just be identity
+  //)igl_Qu8mg5v7";
+const char *__doc_igl_fit_rotations_SSE = R"igl_Qu8mg5v7(See fit_rotations_SSE for the documentation.)igl_Qu8mg5v7";
 const char *__doc_igl_rotate_vectors = R"igl_Qu8mg5v7(// Rotate the vectors V by A radiants on the tangent plane spanned by B1 and
   // B2
   //
@@ -245,6 +317,17 @@ const char *__doc_igl_avg_edge_length = R"igl_Qu8mg5v7(// Compute the average ed
   //   l  average edge length
   //
   // See also: adjacency_matrix)igl_Qu8mg5v7";
+const char *__doc_igl_barycentric_coordinates = R"igl_Qu8mg5v7(// Compute barycentric coordinates in a tet
+  //
+  // Inputs:
+  //   P  #P by 3 Query points in 3d
+  //   A  #P by 3 Tet corners in 3d
+  //   B  #P by 3 Tet corners in 3d
+  //   C  #P by 3 Tet corners in 3d
+  //   D  #P by 3 Tet corners in 3d
+  // Outputs:
+  //   L  #P by 4 list of barycentric coordinates
+  //   )igl_Qu8mg5v7";
 const char *__doc_igl_lscm = R"igl_Qu8mg5v7(// Compute a Least-squares conformal map parametrization (equivalently
   // derived in "Intrinsic Parameterizations of Surface Meshes" [Desbrun et al.
   // 2002] and "Least Squares Conformal Maps for Automatic Texture Atlas
@@ -273,6 +356,52 @@ const char *__doc_igl_find_cross_field_singularities = R"igl_Qu8mg5v7(// Inputs:
   //   isSingularity    #V by 1 boolean eigen Vector indicating the presence of a singularity on a vertex
   //   singularityIndex #V by 1 integer eigen Vector containing the singularity indices
   //)igl_Qu8mg5v7";
+const char *__doc_igl_upsample = R"igl_Qu8mg5v7(// Subdivide a mesh without moving vertices: loop subdivision but odd
+  // vertices stay put and even vertices are just edge midpoints
+  // 
+  // Templates:
+  //   MatV  matrix for vertex positions, e.g. MatrixXd
+  //   MatF  matrix for vertex positions, e.g. MatrixXi
+  // Inputs:
+  //   V  #V by dim  mesh vertices
+  //   F  #F by 3  mesh triangles
+  // Outputs:
+  //   NV new vertex positions, V is guaranteed to be at top
+  //   NF new list of face indices
+  //
+  // NOTE: V should not be the same as NV,
+  // NOTE: F should not be the same as NF, use other proto
+  //
+  // 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:
+  //   P  #P by 3 list of query point positions
+  //   V  #V by 3 list of vertex positions
+  //   Ele  #Ele by (3|2|1) list of (triangle|edge|point) indices
+  // Outputs:
+  //   sqrD  #P list of smallest squared distances
+  //   I  #P list of primitive indices corresponding to smallest distances
+  //   C  #P by 3 list of closest points
+  //
+  // Known bugs: This only computes distances to given primitivess. So
+  // unreferenced vertices are ignored. However, degenerate primitives are
+  // handled correctly: triangle [1 2 2] is treated as a segment [1 2], and
+  // triangle [1 1 1] is treated as a point. So one _could_ add extra
+  // combinatorially degenerate rows to Ele for all unreferenced vertices to
+  // also get distances to points.)igl_Qu8mg5v7";
 const char *__doc_igl_parula = R"igl_Qu8mg5v7(// PARULA like MATLAB's parula
   //
   // Inputs:
@@ -361,6 +490,27 @@ const char *__doc_igl_active_set = R"igl_Qu8mg5v7(// Known Bugs: rows of [Aeq;Ai
   // Benchmark: For a harmonic solve on a mesh with 325K facets, matlab 2.2
   // secs, igl/min_quad_with_fixed.h 7.1 secs
   //)igl_Qu8mg5v7";
+const char *__doc_igl_per_edge_normals = R"igl_Qu8mg5v7(// Compute face normals via vertex position list, face list
+  // Inputs:
+  //   V  #V by 3 eigen Matrix of mesh vertex 3D positions
+  //   F  #F by 3 eigen Matrix of face (triangle) indices
+  //   weight  weighting type
+  //   FN  #F by 3 matrix of 3D face normals per face
+  // Output:
+  //   N  #2 by 3 matrix of mesh edge 3D normals per row
+  //   E  #E by 2 matrix of edge indices per row
+  //   EMAP  #E by 1 matrix of indices from all edges to E
+  //)igl_Qu8mg5v7";
+const char *__doc_igl_covariance_scatter_matrix = R"igl_Qu8mg5v7(// Construct the covariance scatter matrix for a given arap energy
+  // Inputs:
+  //   V  #V by Vdim list of initial domain positions
+  //   F  #F by 3 list of triangle indices into V
+  //   energy  ARAPEnergyType enum value defining which energy is being used.
+  //     See ARAPEnergyType.h for valid options and explanations.
+  // Outputs:
+  //   CSM dim*#V/#F by dim*#V sparse matrix containing special laplacians along
+  //     the diagonal so that when multiplied by V gives covariance matrix
+  //     elements, can be used to speed up covariance matrix computation)igl_Qu8mg5v7";
 const char *__doc_igl_boundary_facets = R"igl_Qu8mg5v7(// BOUNDARY_FACETS Determine boundary faces (edges) of tetrahedra (triangles)
   // stored in T (analogous to qptoolbox's `outline` and `boundary_faces`).
   //
@@ -385,6 +535,27 @@ const char *__doc_igl_compute_frame_field_bisectors = R"igl_Qu8mg5v7(// Compute
   //   BIS1  #F by 3 eigen Matrix of the first per face frame field bisector
   //   BIS2  #F by 3 eigen Matrix of the second per face frame field bisector
   //)igl_Qu8mg5v7";
+const char *__doc_igl_edge_lengths = R"igl_Qu8mg5v7(// Constructs a list of lengths of edges opposite each index in a face
+  // (triangle/tet) list
+  //
+  // Templates:
+  //   DerivedV derived from vertex positions matrix type: i.e. MatrixXd
+  //   DerivedF derived from face indices matrix type: i.e. MatrixXi
+  //   DerivedL derived from edge lengths matrix type: i.e. MatrixXd
+  // Inputs:
+  //   V  eigen matrix #V by 3
+  //   F  #F by 2 list of mesh edges
+  //    or
+  //   F  #F by 3 list of mesh faces (must be triangles)
+  //    or
+  //   T  #T by 4 list of mesh elements (must be tets)
+  // Outputs:
+  //   L  #F by {1|3|6} list of edge lengths 
+  //     for edges, column of lengths
+  //     for triangles, columns correspond to edges [1,2],[2,0],[0,1]
+  //     for tets, columns correspond to edges
+  //     [3 0],[3 1],[3 2],[1 2],[2 0],[0 1]
+  //)igl_Qu8mg5v7";
 const char *__doc_igl_readOBJ = 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
   //
@@ -487,6 +658,18 @@ const char *__doc_igl_min_quad_with_fixed_solve = R"igl_Qu8mg5v7(// Solves a sys
   //   sol  #unknowns+#lagrange by cols solution to linear system
   // Returns true on success, false on error)igl_Qu8mg5v7";
 const char *__doc_igl_min_quad_with_fixed = R"igl_Qu8mg5v7(See min_quad_with_fixed for the documentation.)igl_Qu8mg5v7";
+const char *__doc_igl_writeMESH = R"igl_Qu8mg5v7(// save a tetrahedral volume mesh to a .mesh file
+  //
+  // Templates:
+  //   Scalar  type for positions and vectors (will be cast as double)
+  //   Index  type for indices (will be cast to int)
+  // Input:
+  //   mesh_file_name  path of .mesh file
+  //   V  double matrix of vertex positions  #V by 3
+  //   T  #T list of tet indices into vertex positions
+  //   F  #F list of face indices into vertex positions
+  //
+  // Known bugs: Holes and regions are not supported)igl_Qu8mg5v7";
 const char *__doc_igl_unique = R"igl_Qu8mg5v7(// Act like matlab's [C,IA,IC] = unique(X)
   //
   // Templates:
@@ -534,12 +717,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
   // 
@@ -550,6 +733,23 @@ const char *__doc_igl_slice_into = R"igl_Qu8mg5v7(// Act like the matlab Y(row_i
   //   Y  ym by yn lhs matrix
   // Output:
   //   Y  ym by yn lhs matrix, same as input but Y(R,C) = X)igl_Qu8mg5v7";
+const char *__doc_igl_slice_tets = R"igl_Qu8mg5v7(// SLICE_TETS Slice through a tet mesh (V,T) along a given plane (via its
+  // implicit equation).
+  //
+  // Inputs:
+  //   V  #V by 3 list of tet mesh vertices
+  //   T  #T by 4 list of tet indices into V 
+  //   plane  list of 4 coefficients in the plane equation: [x y z 1]'*plane = 0
+  //   Optional:
+  //     'Manifold' followed by whether to stitch together triangles into a
+  //       manifold mesh {true}: results in more compact U but slightly slower.
+  // Outputs:
+  //   U  #U by 3 list of triangle mesh vertices along slice
+  //   G  #G by 3 list of triangles indices into U
+  //   J  #G list of indices into T revealing from which tet each faces comes
+  //   BC  #U by #V list of barycentric coordinates (or more generally: linear
+  //     interpolation coordinates) so that U = BC*V
+  // )igl_Qu8mg5v7";
 const char *__doc_igl_n_polyvector = R"igl_Qu8mg5v7(// Inputs:
   //   v0, v1         the two #3 by 1 vectors
   //   normalized     boolean, if false, then the vectors are normalized prior to the calculation
@@ -577,6 +777,19 @@ const char *__doc_igl_boundary_loop = R"igl_Qu8mg5v7(// Compute list of ordered
   // Outputs:
   //   L  list of loops where L[i] = ordered list of boundary vertices in loop i
   //)igl_Qu8mg5v7";
+const char *__doc_igl_polar_svd = R"igl_Qu8mg5v7(// Computes the polar decomposition (R,T) of a matrix A using SVD singular
+  // value decomposition
+  //
+  // Inputs:
+  //   A  3 by 3 matrix to be decomposed
+  // Outputs:
+  //   R  3 by 3 rotation matrix part of decomposition (**always rotataion**)
+  //   T  3 by 3 stretch matrix part of decomposition
+  //   U  3 by 3 left-singular vectors
+  //   S  3 by 1 singular values
+  //   V  3 by 3 right-singular vectors
+  //
+  //)igl_Qu8mg5v7";
 const char *__doc_igl_comb_cross_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
@@ -592,6 +805,20 @@ const char *__doc_igl_invert_diag = R"igl_Qu8mg5v7(// Templates:
   //   X  an m by n sparse matrix
   // Outputs:
   //   Y  an m by n sparse matrix)igl_Qu8mg5v7";
+const char *__doc_igl_readMESH = R"igl_Qu8mg5v7(// load a tetrahedral volume mesh from a .mesh file
+  //
+  // Templates:
+  //   Scalar  type for positions and vectors (will be read as double and cast
+  //     to Scalar)
+  //   Index  type for indices (will be read as int and cast to Index)
+  // Input:
+  //   mesh_file_name  path of .mesh file
+  // Outputs:
+  //   V  double matrix of vertex positions  #V by 3
+  //   T  #T list of tet indices into vertex positions
+  //   F  #F list of face indices into vertex positions
+  //
+  // Known bugs: Holes and regions are not supported)igl_Qu8mg5v7";
 const char *__doc_igl_copyleft_comiso_miq = R"igl_Qu8mg5v7(// Inputs:
     //   V              #V by 3 list of mesh vertex 3D positions
     //   F              #F by 3 list of faces indices in V

python/py_doc.h

@@ -1,5 +1,8 @@
 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_cotmatrix;
 extern const char *__doc_igl_floor;
 extern const char *__doc_igl_slice;
@@ -10,24 +13,35 @@ extern const char *__doc_igl_per_vertex_normals;
 extern const char *__doc_igl_sortrows;
 extern const char *__doc_igl_barycenter;
 extern const char *__doc_igl_jet;
+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_colon;
+extern const char *__doc_igl_fit_rotations;
+extern const char *__doc_igl_fit_rotations_planar;
+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_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_comb_frame_field;
 extern const char *__doc_igl_map_vertices_to_circle;
 extern const char *__doc_igl_writeOBJ;
 extern const char *__doc_igl_active_set;
+extern const char *__doc_igl_per_edge_normals;
+extern const char *__doc_igl_covariance_scatter_matrix;
 extern const char *__doc_igl_boundary_facets;
 extern const char *__doc_igl_compute_frame_field_bisectors;
+extern const char *__doc_igl_edge_lengths;
 extern const char *__doc_igl_readOBJ;
 extern const char *__doc_igl_cut_mesh_from_singularities;
 extern const char *__doc_igl_readDMAT;
@@ -37,6 +51,7 @@ extern const char *__doc_igl_doublearea_quad;
 extern const char *__doc_igl_min_quad_with_fixed_precompute;
 extern const char *__doc_igl_min_quad_with_fixed_solve;
 extern const char *__doc_igl_min_quad_with_fixed;
+extern const char *__doc_igl_writeMESH;
 extern const char *__doc_igl_unique;
 extern const char *__doc_igl_unique_rows;
 extern const char *__doc_igl_arap_precomputation;
@@ -44,10 +59,13 @@ extern const char *__doc_igl_arap_solve;
 extern const char *__doc_igl_cross_field_missmatch;
 extern const char *__doc_igl_grad;
 extern const char *__doc_igl_slice_into;
+extern const char *__doc_igl_slice_tets;
 extern const char *__doc_igl_n_polyvector;
 extern const char *__doc_igl_harmonic;
 extern const char *__doc_igl_boundary_loop;
+extern const char *__doc_igl_polar_svd;
 extern const char *__doc_igl_comb_cross_field;
 extern const char *__doc_igl_invert_diag;
+extern const char *__doc_igl_readMESH;
 extern const char *__doc_igl_copyleft_comiso_miq;
 extern const char *__doc_igl_copyleft_comiso_nrosy;

python/py_igl.cpp

@@ -49,6 +49,23 @@
 #include <igl/comb_frame_field.h>
 #include <igl/n_polyvector.h>
 
+#include <igl/point_mesh_squared_distance.h>
+#include <igl/AABB.h>
+#include <igl/readMESH.h>
+#include <igl/writeMESH.h>
+#include <igl/slice_tets.h>
+#include <igl/edge_lengths.h>
+#include <igl/upsample.h>
+#include <igl/cat.h>
+#include <igl/per_edge_normals.h>
+#include <igl/barycentric_coordinates.h>
+#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/.h>
+
 
 void python_export_igl(py::module &m)
 {
@@ -99,4 +116,22 @@ void python_export_igl(py::module &m)
 #include "py_igl/py_comb_frame_field.cpp"
 #include "py_igl/py_n_polyvector.cpp"
 
+#include "py_igl/py_point_mesh_squared_distance.cpp"
+#include "py_igl/py_AABB.cpp"
+#include "py_igl/py_readMESH.cpp"
+#include "py_igl/py_writeMESH.cpp"
+#include "py_igl/py_slice_tets.cpp"
+#include "py_igl/py_edge_lengths.cpp"
+#include "py_igl/py_upsample.cpp"
+#include "py_igl/py_cat.cpp"
+#include "py_igl/py_per_edge_normals.cpp"
+#include "py_igl/py_barycentric_coordinates.cpp"
+#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_.cpp"
+
 }

python/py_igl/py_AABB.cpp

@@ -0,0 +1,14 @@
+py::class_<igl::AABB<Eigen::MatrixXd,3> > AABB(m, "AABB");
+
+AABB
+.def(py::init<>())
+.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::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)
+{
+    return tree.squared_distance(V, Ele, P, sqrD, I, C);
+})
+;

python/py_igl/py_barycentric_coordinates.cpp

@@ -0,0 +1,29 @@
+
+
+m.def("barycentric_coordinates", []
+(
+  const Eigen::MatrixXd& P,
+  const Eigen::MatrixXd& A,
+  const Eigen::MatrixXd& B,
+  const Eigen::MatrixXd& C,
+  const Eigen::MatrixXd& D,
+  Eigen::MatrixXd& L
+)
+{
+  return igl::barycentric_coordinates(P, A, B, C, D, L);
+}, __doc_igl_barycentric_coordinates,
+py::arg("P"), py::arg("A"), py::arg("B"), py::arg("C"), py::arg("D"), py::arg("L"));
+
+m.def("barycentric_coordinates", []
+(
+  const Eigen::MatrixXd& P,
+  const Eigen::MatrixXd& A,
+  const Eigen::MatrixXd& B,
+  const Eigen::MatrixXd& C,
+  Eigen::MatrixXd& L
+)
+{
+  return igl::barycentric_coordinates(P, A, B, C, L);
+}, __doc_igl_barycentric_coordinates,
+py::arg("P"), py::arg("A"), py::arg("B"), py::arg("C"), py::arg("L"));
+

python/py_igl/py_cat.cpp

@@ -0,0 +1,56 @@
+m.def("cat", []
+(
+  const int dim,
+  const Eigen::MatrixXd& A,
+  const Eigen::MatrixXd& B,
+  Eigen::MatrixXd& C
+)
+{
+  return igl::cat(dim, A, B, C);
+}, __doc_igl_cat,
+py::arg("dim"), py::arg("A"), py::arg("B"), py::arg("C"));
+
+m.def("cat", []
+(
+  const int dim,
+  Eigen::MatrixXd& A,
+  Eigen::MatrixXd& B
+)
+{
+  return igl::cat(dim, A, B);
+}, __doc_igl_cat,
+py::arg("dim"), py::arg("A"), py::arg("B"));
+
+m.def("cat", []
+(
+  const int dim,
+  Eigen::MatrixXi& A,
+  Eigen::MatrixXi& B
+)
+{
+  return igl::cat(dim, A, B);
+}, __doc_igl_cat,
+py::arg("dim"), py::arg("A"), py::arg("B"));
+
+//m.def("cat", []
+//(
+//  const std::vector<std::vector<Eigen::MatrixXd > > & A, 
+//  Eigen::MatrixXd & C
+//)
+//{
+//  return igl::cat(A, C);
+//}, __doc_igl_cat,
+//py::arg("A"), py::arg("C"));
+
+m.def("cat", []
+(
+  const int dim,
+  const Eigen::SparseMatrix<double>& A,
+  const Eigen::SparseMatrix<double>& B,
+  Eigen::SparseMatrix<double>& C
+)
+{
+  return igl::cat(dim, A, B, C);
+}, __doc_igl_cat,
+py::arg("dim"), py::arg("A"), py::arg("B"), py::arg("C"));
+

python/py_igl/py_covariance_scatter_matrix.cpp

@@ -0,0 +1,11 @@
+m.def("covariance_scatter_matrix", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  const igl::ARAPEnergyType energy,
+  Eigen::SparseMatrix<double>& CSM
+)
+{
+  return igl::covariance_scatter_matrix(V,F,energy,CSM);
+}, __doc_igl_covariance_scatter_matrix,
+py::arg("V"), py::arg("F"), py::arg("energy"), py::arg("CSM"));

python/py_igl/py_edge_lengths.cpp

@@ -0,0 +1,11 @@
+m.def("edge_lengths", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  Eigen::MatrixXd& L
+)
+{
+  return igl::edge_lengths(V, F, L);
+}, __doc_igl_edge_lengths,
+py::arg("V"), py::arg("F"), py::arg("L"));
+

python/py_igl/py_fit_rotations.cpp

@@ -0,0 +1,22 @@
+m.def("fit_rotations", []
+(
+  const Eigen::MatrixXd& S,
+  const bool single_precision,
+  Eigen::MatrixXd& R
+)
+{
+  return igl::fit_rotations(S, single_precision, R);
+}, __doc_igl_fit_rotations,
+py::arg("S"), py::arg("single_precision"), py::arg("R"));
+
+
+m.def("fit_rotations_planar", []
+(
+  const Eigen::MatrixXd& S,
+  Eigen::MatrixXd& R
+)
+{
+  return igl::fit_rotations_planar(S, R);
+}, __doc_igl_fit_rotations_planar,
+py::arg("S"), py::arg("R"));
+

python/py_igl/py_parula.cpp

@@ -1,3 +1,25 @@
+//m.def("parula", []
+//(
+//  const double f,
+//  T * rgb
+//)
+//{
+//  return igl::parula(f, rgb);
+//}, __doc_igl_parula,
+//py::arg("f"), py::arg("rgb"));
+
+m.def("parula", []
+(
+  const double f,
+  double & r,
+  double & g,
+  double & b
+)
+{
+  return igl::parula(f, r, g, b);
+}, __doc_igl_parula,
+py::arg("f"), py::arg("r"), py::arg("g"), py::arg("b"));
+
 m.def("parula", []
 (
   const Eigen::MatrixXd& Z,
@@ -6,7 +28,7 @@ m.def("parula", []
 )
 {
   assert_is_VectorX("Z",Z);
-  return igl::parula(Z,normalize,C);
+  return igl::parula(Z, normalize, C);
 }, __doc_igl_parula,
 py::arg("Z"), py::arg("normalize"), py::arg("C"));
 
@@ -19,6 +41,7 @@ m.def("parula", []
 )
 {
   assert_is_VectorX("Z",Z);
-  return igl::parula(Z,min_Z,max_Z,C);
+  return igl::parula(Z, min_Z, max_Z, C);
 }, __doc_igl_parula,
 py::arg("Z"), py::arg("min_Z"), py::arg("max_Z"), py::arg("C"));
+

python/py_igl/py_per_edge_normals.cpp

@@ -0,0 +1,50 @@
+py::enum_<igl::PerEdgeNormalsWeightingType>(m, "PerEdgeNormalsWeightingType")
+    .value("PER_EDGE_NORMALS_WEIGHTING_TYPE_UNIFORM", igl::PER_EDGE_NORMALS_WEIGHTING_TYPE_UNIFORM)
+    .value("PER_EDGE_NORMALS_WEIGHTING_TYPE_AREA", igl::PER_EDGE_NORMALS_WEIGHTING_TYPE_AREA)
+    .value("PER_EDGE_NORMALS_WEIGHTING_TYPE_DEFAULT", igl::PER_EDGE_NORMALS_WEIGHTING_TYPE_DEFAULT)
+    .value("NUM_PER_EDGE_NORMALS_WEIGHTING_TYPE", igl::NUM_PER_EDGE_NORMALS_WEIGHTING_TYPE)
+    .export_values();
+
+
+m.def("per_edge_normals", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  const igl::PerEdgeNormalsWeightingType weight,
+  const Eigen::MatrixXd& FN,
+  Eigen::MatrixXd& N,
+  Eigen::MatrixXi& E,
+  Eigen::MatrixXi& EMAP
+)
+{
+  return igl::per_edge_normals(V, F, weight, FN, N, E, EMAP);
+}, __doc_igl_per_edge_normals,
+py::arg("V"), py::arg("F"), py::arg("weight"), py::arg("FN"), py::arg("N"), py::arg("E"), py::arg("EMAP"));
+
+m.def("per_edge_normals", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  const igl::PerEdgeNormalsWeightingType weight,
+  Eigen::MatrixXd& N,
+  Eigen::MatrixXi& E,
+  Eigen::MatrixXi& EMAP
+)
+{
+  return igl::per_edge_normals(V, F, weight, N, E, EMAP);
+}, __doc_igl_per_edge_normals,
+py::arg("V"), py::arg("F"), py::arg("weight"), py::arg("N"), py::arg("E"), py::arg("EMAP"));
+
+m.def("per_edge_normals", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  Eigen::MatrixXd& N,
+  Eigen::MatrixXi& E,
+  Eigen::MatrixXi& EMAP
+)
+{
+  return igl::per_edge_normals(V, F, N, E, EMAP);
+}, __doc_igl_per_edge_normals,
+py::arg("V"), py::arg("F"), py::arg("N"), py::arg("E"), py::arg("EMAP"));
+

python/py_igl/py_point_mesh_squared_distance.cpp

@@ -0,0 +1,15 @@
+m.def("point_mesh_squared_distance", []
+(
+  const Eigen::MatrixXd& P,
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& Ele,
+  Eigen::MatrixXd& sqrD,
+  Eigen::MatrixXi& I,
+  Eigen::MatrixXd& C
+)
+{
+//  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_polar_svd.cpp

@@ -0,0 +1,26 @@
+m.def("polar_svd", []
+(
+  const Eigen::MatrixXd& A,
+  Eigen::MatrixXd& R,
+  Eigen::MatrixXd& T,
+  Eigen::MatrixXd& U,
+  Eigen::MatrixXd& S,
+  Eigen::MatrixXd& V
+)
+{
+  return igl::polar_svd(A, R, T, U, S, V);
+}, __doc_igl_polar_svd,
+py::arg("A"), py::arg("R"), py::arg("T"), py::arg("U"), py::arg("S"), py::arg("V"));
+
+
+m.def("polar_svd", []
+(
+  const Eigen::MatrixXd& A,
+  Eigen::MatrixXd& R,
+  Eigen::MatrixXd& T
+)
+{
+  return igl::polar_svd(A, R, T);
+}, __doc_igl_polar_svd,
+py::arg("A"), py::arg("R"), py::arg("T"));
+

python/py_igl/py_readMESH.cpp

@@ -0,0 +1,24 @@
+m.def("readMESH", []
+(
+  const std::string mesh_file_name,
+  std::vector<std::vector<double> > & V,
+  std::vector<std::vector<int> > & T,
+  std::vector<std::vector<int> > & F
+)
+{
+  return igl::readMESH(mesh_file_name, V, T, F);
+}, __doc_igl_readMESH,
+py::arg("mesh_file_name"), py::arg("V"), py::arg("T"), py::arg("F"));
+
+m.def("readMESH", []
+(
+  const std::string mesh_file_name,
+  Eigen::MatrixXd& V,
+  Eigen::MatrixXi& T,
+  Eigen::MatrixXi& F
+)
+{
+  return igl::readMESH(mesh_file_name, V, T, F);
+}, __doc_igl_readMESH,
+py::arg("mesh_file_name"), py::arg("V"), py::arg("T"), py::arg("F"));
+

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_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

@@ -0,0 +1,21 @@
+m.def("slice_tets", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& T,
+  const Eigen::MatrixXd& plane,
+  Eigen::MatrixXd& U,
+  Eigen::MatrixXi& G,
+  Eigen::MatrixXi& J,
+  Eigen::SparseMatrix<double> & BC
+)
+{
+  assert_is_VectorX("plane", plane);
+  Eigen::VectorXd pl;
+  if (plane.size() != 0)
+    pl = plane;
+  Eigen::VectorXi Jv;
+  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_upsample.cpp

@@ -0,0 +1,22 @@
+m.def("upsample", []
+(
+  Eigen::MatrixXd& V,
+  Eigen::MatrixXi& F
+)
+{
+  return igl::upsample(V, F);
+}, __doc_igl_upsample,
+py::arg("V"), py::arg("F"));
+
+m.def("upsample", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  Eigen::MatrixXd& NV,
+  Eigen::MatrixXi& NF
+)
+{
+  return igl::upsample(V, F, NV, NF);
+}, __doc_igl_upsample,
+py::arg("V"), py::arg("F"), py::arg("NV"), py::arg("NF"));
+

python/py_igl/py_writeMESH.cpp

@@ -0,0 +1,24 @@
+m.def("writeMESH", []
+(
+  const std::string mesh_file_name,
+  const std::vector<std::vector<double> > & V,
+  const std::vector<std::vector<int> > & T,
+  const std::vector<std::vector<int> > & F
+)
+{
+  return igl::writeMESH(mesh_file_name, V, T, F);
+}, __doc_igl_writeMESH,
+py::arg("mesh_file_name"), py::arg("V"), py::arg("T"), py::arg("F"));
+
+m.def("writeMESH", []
+(
+  const std::string str,
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXd& T,
+  const Eigen::MatrixXi& F
+)
+{
+  return igl::writeMESH(str, V, T, F);
+}, __doc_igl_writeMESH,
+py::arg("str"), py::arg("V"), py::arg("T"), py::arg("F"));
+

python/py_vector.cpp

@@ -91,6 +91,7 @@ py::class_<Type> bind_eigen_2(py::module &m, const char *name,
         .def("size", [](const Type &m) { return m.size(); })
         .def("cols", [](const Type &m) { return m.cols(); })
         .def("rows", [](const Type &m) { return m.rows(); })
+        .def("shape", [](const Type &m) { return std::tuple<int,int>(m.rows(), m.cols()); })
 
         /* Extract rows and colums */
         .def("col", [](const Type &m, int i) {
@@ -128,6 +129,9 @@ py::class_<Type> bind_eigen_2(py::module &m, const char *name,
         .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)); })
@@ -248,6 +252,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 */
@@ -356,6 +380,8 @@ py::class_<Type> bind_eigen_sparse_2(py::module &m, const char *name,
         .def("size", [](const Type &m) { return m.size(); })
         .def("cols", [](const Type &m) { return m.cols(); })
         .def("rows", [](const Type &m) { return m.rows(); })
+        .def("shape", [](const Type &m) { return std::tuple<int,int>(m.rows(), m.cols()); })
+
 
         /* Initialization */
         .def("setZero", [](Type &m) { m.setZero(); })
@@ -515,6 +541,7 @@ py::class_<Type> bind_eigen_diagonal_2(py::module &m, const char *name,
         .def("size", [](const Type &m) { return m.size(); })
         .def("cols", [](const Type &m) { return m.cols(); })
         .def("rows", [](const Type &m) { return m.rows(); })
+        .def("shape", [](const Type &m) { return std::tuple<int,int>(m.rows(), m.cols()); })
 
         /* Initialization */
         .def("setZero", [](Type &m) { m.setZero(); })
@@ -599,6 +626,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> >();
@@ -648,6 +678,8 @@ void python_export_vector(py::module &m) {
        else
           return "Numerical Issue";
     })
+    .def("analyzePattern",[](Eigen::SimplicialLLT<Eigen::SparseMatrix<double > >& s, const Eigen::SparseMatrix<double>& a) { return s.analyzePattern(a); })
+    .def("factorize",[](Eigen::SimplicialLLT<Eigen::SparseMatrix<double > >& s, const Eigen::SparseMatrix<double>& a) { return s.factorize(a); })
     .def("solve",[](const Eigen::SimplicialLLT<Eigen::SparseMatrix<double > >& s, const Eigen::MatrixXd& rhs) { return Eigen::MatrixXd(s.solve(rhs)); })
     ;
 

python/scripts/basic_function.mako

@@ -0,0 +1,46 @@
+% for enum in enums:
+py::enum_<\
+% for n in enum['namespaces']:
+${n}::\
+% endfor
+${enum['name']}>(m, "${enum['name']}")
+% for c in enum['constants']:
+    .value("${c}", \
+% for n in enum['namespaces']:
+${n}::\
+% endfor
+${c})
+% endfor
+    .export_values();
+% endfor
+
+
+% for func in functions:
+m.def("${func['name']}", []
+(
+  % for p in func['parameters'][:-1]:
+  ${p['type']} ${p['name']},
+  % endfor
+  ${func['parameters'][-1]['type']} ${func['parameters'][-1]['name']}
+)
+{
+  return \
+% for n in func['namespaces']:
+${n}::\
+% endfor
+${func['name']}(\
+% for p in func['parameters'][:-1]:
+${p['name']}, \
+% endfor
+${func['parameters'][-1]['name']});
+}, __doc_\
+% for n in func['namespaces']:
+${n}_\
+% endfor
+${func['name']},
+% for p in func['parameters'][:-1]:
+py::arg("${p['name']}"), \
+% endfor
+py::arg("${func['parameters'][-1]['name']}"));
+
+% endfor

python/scripts/generate_bindings.py

@@ -0,0 +1,301 @@
+#!/usr/bin/env python3
+#
+#  Syntax: generate_docstrings.py <path_to_c++_header_files> <path_to_python_files>
+#
+#  Extract documentation from C++ header files to use it in libiglPython bindings
+#
+
+import os, sys, glob
+import pickle
+
+import shutil
+from joblib import Parallel, delayed
+from multiprocessing import cpu_count
+from mako.template import Template
+from parser import parse
+
+
+# http://stackoverflow.com/questions/3207219/how-to-list-all-files-of-a-directory-in-python
+def get_filepaths(directory):
+    """
+    This function will generate the file names in a directory
+    tree by walking the tree either top-down or bottom-up. For each
+    directory in the tree rooted at directory top (including top itself),
+    it yields a 3-tuple (dirpath, dirnames, filenames).
+    """
+    file_paths = []  # List which will store all of the full filepaths.
+
+    # Walk the tree.
+    for root, directories, files in os.walk(directory):
+        for filename in files:
+            # Join the two strings in order to form the full filepath.
+            filepath = os.path.join(root, filename)
+            file_paths.append(filepath)  # Add it to the list.
+
+    return file_paths  # Self-explanatory.
+
+
+def get_name_from_path(path, basepath, prefix, postfix):
+    f_clean = path[len(basepath):]
+    f_clean = f_clean.replace(basepath, "")
+    f_clean = f_clean.replace(postfix, "")
+    f_clean = f_clean.replace(prefix, "")
+    f_clean = f_clean.replace("/", "_")
+    f_clean = f_clean.replace("\\", "_")
+    f_clean = f_clean.replace(" ", "_")
+    f_clean = f_clean.replace(".", "_")
+    return f_clean
+
+
+def map_parameter_types(name, cpp_type, parsed_types, errors, enum_types):
+    # TODO Replace with proper regex matching and derive types from templates, comment parsing, names in cpp files
+    # CAUTION: This is work in progress mapping code to get a grip of the problem
+    # Types to map
+    #    const int dim -> const int& dim ?
+    result = []
+
+    if cpp_type.startswith("const"):
+        result.append("const ")
+        cpp_type = cpp_type[6:]  # Strip const part
+
+    # Handle special types
+    skip_parsing = False
+    if cpp_type.startswith("MatY"):
+        result.append("Eigen::SparseMatrix<double>&")
+        skip_parsing = True
+    if cpp_type == "std::vector<std::vector<Scalar> > &":
+        result.append("std::vector<std::vector<double> > &")
+        skip_parsing = True
+    if cpp_type == "std::vector<std::vector<Index> > &":
+        result.append("std::vector<std::vector<int> > &")
+        skip_parsing = True
+    for constant in enum_types:
+        if cpp_type.endswith(constant):
+            result.append(cpp_type)
+            skip_parsing = True
+
+    if len(parsed_types) == 0:
+        errors.append("Empty typechain: %s" % cpp_type)
+        if cpp_type == "int" or cpp_type == "bool":
+            return cpp_type, True
+        else:
+            return cpp_type, False
+
+    # print(parsed_types, cpp_type)
+    if not skip_parsing:
+        for i, t in enumerate(parsed_types):
+
+            if t == "Eigen":
+                result.append("Eigen::")
+                continue
+            if t == "std":
+                result.append("std::")
+                continue
+
+            if t == "PlainObjectBase" or t == "MatrixBase":
+                if name == "F":
+                    result.append("MatrixXi&")
+                elif name == "V":
+                    result.append("MatrixXd&")
+                else:
+                    result.append("MatrixXd&")
+                break
+            if t == "MatrixXi" or t == "VectorXi":
+                result.append("MatrixXi&")
+                break
+            if t == "MatrixXd" or t == "VectorXd":
+                result.append("MatrixXd&")
+                break
+            if t == "SparseMatrix" and len(parsed_types) >= i + 2 and (
+                    parsed_types[i + 1] == "Scalar" or parsed_types[i + 1] == "T"):
+                result.append("SparseMatrix<double>&")
+                break
+            if t == "SparseVector" and len(parsed_types) >= i + 2 and (parsed_types[i + 1] == "Scalar" or parsed_types[
+                    i + 1] == "T"):
+                result.append("SparseMatrix<double>&")
+                break
+
+            if t == "bool" or t == "int" or t == "double" or t == "unsigned" or t == "string":
+                if cpp_type.endswith("&"):
+                    result.append(t + " &")
+                else:
+                    result.append(t)
+                break
+
+            else:
+                errors.append("Unknown typechain: %s" % cpp_type)
+                return cpp_type, False
+
+
+    return "".join(result), True
+
+
+if __name__ == '__main__':
+
+    if len(sys.argv) != 2:
+        print('Syntax: %s <path_to_c++_files>' % sys.argv[0])
+        exit(-1)
+
+    errors = {"missing": [], "empty": [], "others": [], "incorrect": [], "render": [], "various": []}
+    files = {"complete": [], "partial": [], "errors": [], "others": [], "empty": []}
+
+    # List all files in the given folder and subfolders
+    cpp_base_path = sys.argv[1]
+    cpp_file_paths = get_filepaths(cpp_base_path)
+
+    # Add all the .h filepaths to a dict
+    print("Collecting cpp files for parsing...")
+    mapping = {}
+    cppmapping = {}
+    for f in cpp_file_paths:
+        if f.endswith(".h"):
+            name = get_name_from_path(f, cpp_base_path, "", ".h")
+            mapping[name] = f
+
+        if f.endswith(".cpp"):
+            name = get_name_from_path(f, cpp_base_path, "", ".cpp")
+            cppmapping[name] = f
+
+    # Add all python binding files to a list
+    implemented_names = list(mapping.keys())  # ["point_mesh_squared_distance"]
+    implemented_names.sort()
+    single_postfix = ""
+    single_prefix = ""
+
+    # Create a list of all cpp header files
+    files_to_parse = []
+    cppfiles_to_parse = []
+    for n in implemented_names:
+        files_to_parse.append(mapping[n])
+
+        if n not in cppmapping:
+            errors["missing"].append("No cpp source file for function %s found." % n)
+        else:
+            cppfiles_to_parse.append(cppmapping[n])
+
+    # Parse c++ header files
+    print("Parsing header files...")
+    load_headers = False
+    if load_headers:
+        with open("headers.dat", 'rb') as fs:
+            dicts = pickle.load(fs)
+    else:
+        job_count = cpu_count()
+        dicts = Parallel(n_jobs=job_count)(delayed(parse)(path) for path in files_to_parse)
+
+    if not load_headers:
+        print("Saving parsed header files...")
+        with open("headers.dat", 'wb') as fs:
+            pickle.dump(dicts, fs)
+
+    # Not yet needed, as explicit template parsing does not seem to be supported in clang
+    # Parse c++ source files
+    # cppdicts = Parallel(n_jobs=job_count)(delayed(parse)(path) for path in cppfiles_to_parse)
+
+    # Change directory to become independent of execution directory
+    print("Generating directory tree for binding files...")
+    path = os.path.dirname(__file__)
+    if path != "":
+        os.chdir(path)
+    try:
+        shutil.rmtree("generated")
+    except:
+        pass # Ignore missing generated directory
+    os.makedirs("generated/complete")
+    os.mkdir("generated/partial")
+
+    print("Generating and writing binding files...")
+    for idx, n in enumerate(implemented_names):
+        d = dicts[idx]
+        contained_elements = sum(map(lambda x: len(x), d.values()))
+
+        # Skip files that don't contain functions/enums/classes
+        if contained_elements == 0:
+            errors["empty"].append("Function %s contains no parseable content in cpp header. Something might be wrong." % n)
+            files["empty"].append(n)
+            continue
+
+        # Add functions with classes to others
+        if len(d["classes"]) != 0 or len(d["structs"]) != 0:
+            errors["others"].append("Function %s contains classes/structs in cpp header. Skipping" % n)
+            files["others"].append(n)
+            continue
+
+        # Work on files that contain only functions/enums and namespaces
+        if len(d["functions"]) + len(d["namespaces"]) + len(d["enums"]) == contained_elements:
+            correct_functions = []
+            incorrect_functions = []
+
+            # Collect enums to generate binding files
+            enums = []
+            enum_types = []
+            for e in d["enums"]:
+                enums.append({"name": e.name, "namespaces": d["namespaces"], "constants": e.constants})
+                enum_types.append(e.name)
+
+            # Collect functions to generate binding files
+            for f in d["functions"]:
+                parameters = []
+                correct_function = True
+                f_errors = []
+                for p in f.parameters:
+                    typ, correct = map_parameter_types(p[0], p[1], p[2], f_errors, enum_types)
+                    correct_function &= correct
+                    parameters.append({"name": p[0], "type": typ})
+
+                if correct_function and len(parameters) > 0: #TODO add constants like EPS
+                    correct_functions.append({"parameters": parameters, "namespaces": d["namespaces"], "name": f.name})
+                elif len(parameters) > 0:
+                    incorrect_functions.append({"parameters": parameters, "namespaces": d["namespaces"], "name": f.name})
+                    errors["incorrect"].append("Incorrect function in %s: %s, %s\n" % (n, f.name, ",".join(f_errors)))
+                else:
+                    errors["various"].append("Function without pars in %s: %s, %s\n" % (n, f.name, ","
+                                                                                                     "".join(f_errors)))
+
+            # Write binding files
+            try:
+                tpl = Template(filename='basic_function.mako')
+                rendered = tpl.render(functions=correct_functions, enums=enums)
+                tpl1 = Template(filename='basic_function.mako')
+                rendered1 = tpl.render(functions=incorrect_functions, enums=enums)
+                path = "generated/"
+                if len(incorrect_functions) == 0 and (len(correct_functions) != 0 or len(enums) != 0):
+                    path += "complete/"
+                    with open(path + single_prefix + "py_" + n + ".cpp", 'w') as fs:
+                        fs.write(rendered)
+                    files["complete"].append(n)
+                else:
+                    path += "partial/"
+                    with open(path + single_prefix + "py_" + n + ".cpp", 'w') as fs:
+                        fs.write("// COMPLETE BINDINGS ========================\n")
+                        fs.write(rendered)
+                        fs.write("\n\n\n\n// INCOMPLETE BINDINGS ========================\n")
+                        fs.write(rendered1)
+
+                    if len(correct_functions) != 0:
+                        files["partial"].append(n)
+                    else:
+                        files["errors"].append(n)
+
+            except Exception as e:
+                files["errors"].append(n)
+                errors["render"].append("Template rendering failed:" + n + " " + str(correct_functions) + ", incorrect "
+                                                                                                "functions are " + str(
+                    incorrect_functions) + str(e) + "\n")
+
+    print("Writing error and overview files...")
+    with open("errors.txt" + single_postfix, 'w') as fs:
+        l = list(errors.keys())
+        l.sort()
+        for k in l:
+            fs.write("%s: %i \n" %(k, len(errors[k])))
+            fs.writelines("\n".join(errors[k]))
+            fs.write("\n\n\n")
+
+    with open("files.txt" + single_postfix, 'w') as fs:
+        l = list(files.keys())
+        l.sort()
+        for k in l:
+            fs.write("%s: %i \n" %(k, len(files[k])))
+            fs.writelines("\n".join(files[k]))
+            fs.write("\n\n\n")

python/scripts/generate_docstrings.py

@@ -87,16 +87,16 @@ if __name__ == '__main__':
     cpplines = []
 
     for idx, n in enumerate(implemented_names):
-        dict = dicts[idx]
-        contained_elements = sum(map(lambda x: len(x), dict.values()))
+        d = dicts[idx]
+        contained_elements = sum(map(lambda x: len(x), d.values()))
         # Check for files that don't contain functions/enums/classes
         if contained_elements == 0:
             print("Function %s contains no parseable content in cpp header. Something might be wrong." % n)
             continue
         else:
             names = []
-            namespaces = "_".join(dict["namespaces"])  # Assumption that all entities lie in deepest namespace
-            for f in dict["functions"]:
+            namespaces = "_".join(d["namespaces"])  # Assumption that all entities lie in deepest namespace
+            for f in d["functions"]:
                 h_string = "extern const char *__doc_" + namespaces + "_" + f.name + ";\n"
                 docu_string = "See " + f.name + " for the documentation."
                 if f.documentation != "":

python/scripts/parser.py

@@ -3,10 +3,16 @@ import os
 from threading import Thread
 
 import clang.cindex
+import ccsyspath
 import itertools
 from mako.template import Template
 
 
+
+
+
+
+
 def get_annotations(node):
     return [c.displayname for c in node.get_children()
             if c.kind == clang.cindex.CursorKind.ANNOTATE_ATTR]
@@ -17,33 +23,33 @@ class Function(object):
         self.name = cursor.spelling
         self.annotations = get_annotations(cursor)
         self.access = cursor.access_specifier
-        #        template_pars = [c.extent for c in cursor.get_children() if c.kind == clang.cindex.CursorKind.TEMPLATE_TYPE_PARAMETER]
-        #        parameter_dec = [c for c in cursor.get_children() if c.kind == clang.cindex.CursorKind.PARM_DECL]
-        #        print(parameter_dec, template_pars)
-        #        print(cursor.get_num_template_arguments(), cursor.get_template_argument_type(0), cursor.get_template_argument_value(0), template_pars, parameter_dec)
-        self.parameters = []
-        self.parnames = []
+#        template_pars = [c.extent for c in cursor.get_children() if c.kind == clang.cindex.CursorKind.TEMPLATE_TYPE_PARAMETER]
+        parameter_dec = [c for c in cursor.get_children() if c.kind == clang.cindex.CursorKind.PARM_DECL]
+        
+        parameters = []
+        for p in parameter_dec:
+            children = []
+            for c in p.get_children():
+#                print(c.spelling)
+                children.append(c.spelling)
+            parameters.append((p.spelling, p.type.spelling, children))
+
+        self.parameters = parameters
         self.documentation = cursor.raw_comment
 
 
 class Enum(object):
     def __init__(self, cursor):
         self.name = cursor.spelling
-        self.annotations = get_annotations(cursor)
-        self.access = cursor.access_specifier
-        #        template_pars = [c.extent for c in cursor.get_children() if c.kind == clang.cindex.CursorKind.TEMPLATE_TYPE_PARAMETER]
-        #        parameter_dec = [c for c in cursor.get_children() if c.kind == clang.cindex.CursorKind.PARM_DECL]
-        #        print(parameter_dec, template_pars)
-        #        print(cursor.get_num_template_arguments(), cursor.get_template_argument_type(0), cursor.get_template_argument_value(0), template_pars, parameter_dec)
-        self.parameters = []
-        self.parnames = []
+        self.constants = [c.spelling for c in cursor.get_children() if c.kind ==
+                       clang.cindex.CursorKind.ENUM_CONSTANT_DECL]
         self.documentation = cursor.raw_comment
 
-# class Class(object):
-#    def __init__(self, cursor):
-#        self.name = cursor.spelling
+class Class(object):
+    def __init__(self, cursor):
+        self.name = cursor.spelling
 #        self.functions = []
-#        self.annotations = get_annotations(cursor)
+        self.annotations = get_annotations(cursor)
 
 #        for c in cursor.get_children():
 #            if (c.kind == clang.cindex.CursorKind.CXX_METHOD and
@@ -51,19 +57,16 @@ class Enum(object):
 #                f = Function(c)
 #                self.functions.append(f)
 
-def find_namespace_node(c):
-    if (c.kind == clang.cindex.CursorKind.NAMESPACE and c.spelling == "igl"):
-        return c
-    else:
-        for child_node in c.get_children():
-            return find_namespace_node(child_node)
 
 
 def traverse(c, path, objects):
     if c.location.file and not c.location.file.name.endswith(path):
         return
 
-    # print(c.kind, c.spelling)
+    if c.spelling == "PARULA_COLOR_MAP": # Fix to prevent python stack overflow from infinite recursion
+        return
+
+#    print(c.kind, c.spelling)
 
 
     if c.kind == clang.cindex.CursorKind.TRANSLATION_UNIT or c.kind == clang.cindex.CursorKind.UNEXPOSED_DECL:
@@ -76,7 +79,7 @@ def traverse(c, path, objects):
         pass
 
     elif c.kind == clang.cindex.CursorKind.FUNCTION_TEMPLATE:
-        # print("Function Template", c.spelling, c.raw_comment)
+#        print("Function Template", c.spelling, c.raw_comment)
         objects["functions"].append(Function(c))
         return
 
@@ -90,6 +93,18 @@ def traverse(c, path, objects):
         objects["enums"].append(Enum(c))
         return
 
+    elif c.kind == clang.cindex.CursorKind.CLASS_DECL:
+        objects["classes"].append(Class(c))
+        return
+
+    elif c.kind == clang.cindex.CursorKind.CLASS_TEMPLATE:
+        objects["classes"].append(Class(c))
+        return
+
+    elif c.kind == clang.cindex.CursorKind.STRUCT_DECL:
+        objects["structs"].append(Class(c))
+        return
+
     else:
         # print("Unknown", c.kind, c.spelling)
         pass
@@ -100,20 +115,16 @@ def traverse(c, path, objects):
 
 def parse(path):
     index = clang.cindex.Index.create()
-    tu = index.parse(path, ['-x', 'c++', '-std=c++11', '-fparse-all-comments', '-DIGL_STATIC_LIBRARY'])
-    # Clang can't parse files with missing definitions, add static library definition
-    objects = {"functions": [], "enums": [], "namespaces": [], "classes": []}
+    # Clang can't parse files with missing definitions, add static library definition or not?
+    args = ['-x', 'c++', '-std=c++11', '-fparse-all-comments', '-DIGL_STATIC_LIBRARY']
+    args.append('-I/usr/include/eigen3/') # TODO Properly add all needed includes
+    syspath = ccsyspath.system_include_paths('clang++-3.7') # Add the system libraries
+    incargs = [(b'-I' + inc).decode("utf-8") for inc in syspath]
+    args.extend(incargs)
+
+    tu = index.parse(path, args)
+    objects = {"functions": [], "enums": [], "namespaces": [], "classes": [], "structs": []}
     traverse(tu.cursor, path, objects)
-
-    #    tpl = Template(filename='bind.mako')
-    #    rendered = tpl.render(functions=functions)
-
-    #    OUTPUT_DIR = 'generated'
-
-    #    if not os.path.isdir(OUTPUT_DIR): os.mkdir(OUTPUT_DIR)
-
-    #    with open("generated/{}.bind.cc".format(sys.argv[1]), "w") as f:
-    #        f.write(rendered)
     return objects
 
 if __name__ == '__main__':

python/tutorial/704_SignedDistance.py

@@ -0,0 +1,134 @@
+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
+
+V = igl.eigen.MatrixXd()
+F = igl.eigen.MatrixXi()
+T = igl.eigen.MatrixXi()
+tree = igl.AABB()
+FN = igl.eigen.MatrixXd()
+VN = igl.eigen.MatrixXd()
+EN = igl.eigen.MatrixXd()
+E = igl.eigen.MatrixXi()
+EMAP = igl.eigen.MatrixXi()
+
+max_distance = 1
+slice_z = 0.5
+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)
+    # C_vis.setBottomRows(V.rows(), color)
+
+
+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())])
+    V_vis = igl.eigen.MatrixXd()
+    F_vis = igl.eigen.MatrixXi()
+
+    # Extract triangle mesh slice through volume mesh and subdivide nasty triangles
+    J = igl.eigen.MatrixXi()
+    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.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.core.lighting_factor = overlay
+
+
+def key_down(viewer, key, modifier):
+    global slice_z, overlay
+
+    if key == ord(' '):
+        overlay = not overlay
+    elif key == ord('.'):
+        slice_z = min(slice_z + 0.01, 0.99)
+    elif key == ord(','):
+        slice_z = max(slice_z - 0.01, 0.01)
+    else:
+        return False
+
+    update_visualization(viewer)
+    return True
+
+
+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);
+
+# 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
+tree.init(V, F)
+
+# 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
+update_visualization(viewer);
+viewer.callback_key_down = key_down
+viewer.core.show_lines = False
+viewer.launch()

shared/cmake/CMakeLists.txt

@@ -115,6 +115,8 @@ if(LIBIGL_WITH_BBW)
   if(LIBIGL_USE_STATIC_LIBRARY)
     CompileIGL_Module("bbw")
     if(LIBIGL_WITH_MOSEK)
+      find_package(MOSEK REQUIRED)
+      list(APPEND LIBIGL_INCLUDE_DIRS ${MOSEK_INCLUDE_DIR})
       target_include_directories(iglbbw PRIVATE ${MOSEK_INCLUDE_DIR})
     else()
       target_compile_definitions(iglbbw PRIVATE -DIGL_NO_MOSEK)

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>. <strong>Each
+parameter</strong> 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,
@@ -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.