Browse Source

Fixed minor bugs in bindings, moved AABB implementation to own file, added boolean matrix and comparison functions to python bindings

Former-commit-id: aa6e79dc36e4b86c402d1da2eced1c1a27b73372
Sebastian Koch 9 years ago
parent
commit
6e4a3738aa

+ 955 - 0
include/igl/AABB.cpp

@@ -0,0 +1,955 @@
+// 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 void 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 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&;
+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
+
+

+ 0 - 932
include/igl/AABB.h

@@ -288,936 +288,4 @@ public:
       EIGEN_MAKE_ALIGNED_OPERATOR_NEW
     };
 }
-
-// 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");
-#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

+ 1 - 0
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

+ 9 - 75
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

+ 0 - 32
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);
     }
   }
 }

+ 72 - 0
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&);

+ 33 - 0
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);
+
     }
 
   }

+ 1 - 0
include/igl/in_element.h

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

+ 1 - 0
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

+ 2 - 2
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,

+ 3 - 0
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

+ 2 - 1
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

+ 6 - 5
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

+ 2 - 1
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

+ 2 - 0
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

+ 2 - 0
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

+ 4 - 0
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

+ 1 - 0
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

+ 1 - 0
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

+ 1 - 0
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

+ 4 - 1
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

+ 5 - 1
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]

+ 51 - 6
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).
   //
@@ -340,6 +374,17 @@ const char *__doc_igl_upsample = R"igl_Qu8mg5v7(// Subdivide a mesh without movi
   //
   // Known issues:
   //   - assumes (V,F) is edge-manifold.)igl_Qu8mg5v7";
+const char *__doc_igl_slice_mask = R"igl_Qu8mg5v7(// Act like the matlab X(row_mask,col_mask) operator, where
+  // row_mask, col_mask are non-negative integer indices.
+  // 
+  // Inputs:
+  //   X  m by n matrix
+  //   R  m list of row bools
+  //   C  n list of column bools
+  // Output:
+  //   Y  #trues-in-R by #trues-in-C matrix
+  //
+  // See also: slice_mask)igl_Qu8mg5v7";
 const char *__doc_igl_point_mesh_squared_distance = R"igl_Qu8mg5v7(// Compute distances from a set of points P to a triangle mesh (V,F)
   //
   // Inputs:
@@ -672,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
   // 

+ 4 - 0
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;
@@ -26,6 +29,7 @@ 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;

+ 7 - 0
python/py_igl.cpp

@@ -62,6 +62,9 @@
 #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)
@@ -126,5 +129,9 @@ void python_export_igl(py::module &m)
 #include "py_igl/py_fit_rotations.cpp"
 #include "py_igl/py_polar_svd.cpp"
 #include "py_igl/py_covariance_scatter_matrix.cpp"
+#include "py_igl/py_slice_mask.cpp"
+#include "py_igl/py_signed_distance.cpp"
+
+//#include "py_igl/py_.cpp"
 
 }

+ 134 - 0
python/py_igl/py_signed_distance.cpp

@@ -0,0 +1,134 @@
+
+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;
+  if (S.size() != 0)
+    Sv = S;
+  Eigen::VectorXi Iv;
+  if (I.size() != 0)
+    Iv = I;
+  return igl::signed_distance(P, V, F, sign_type, Sv, Iv, C, N);
+}, __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"));
+

+ 76 - 0
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"));
+

+ 2 - 4
python/py_igl/py_slice_tets.cpp

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

+ 26 - 0
python/py_vector.cpp

@@ -129,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)); })
@@ -249,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 */
@@ -603,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> >();

+ 46 - 0
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

+ 301 - 0
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")

+ 4 - 4
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 != "":

+ 50 - 39
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__':

+ 52 - 58
python/tutorial/704_SignedDistance.py

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