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

Former-commit-id: 53f243d14bd8207ca17dc8e3a546fc246097a222

examples/skeleton-posing/example.cpp

@@ -49,6 +49,9 @@
 
 #ifdef __APPLE__
 #include <GLUT/glut.h>
+#ifndef GLUT_ACTIVE_COMMAND
+#  define GLUT_ACTIVE_COMMAND 9
+#endif
 #include <Carbon/Carbon.h>
 #else
 #include <GL/glut.h>

include/igl/all_edges.cpp

@@ -7,9 +7,10 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "all_edges.h"
 
+template <typename DerivedF, typename DerivedE>
 IGL_INLINE void igl::all_edges(
-  const Eigen::MatrixXi & F,
-  Eigen::MatrixXi & E)
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedE> & E)
 {
   E.resize(F.rows()*F.cols(),F.cols()-1);
   switch(F.cols())
@@ -41,3 +42,8 @@ IGL_INLINE void igl::all_edges(
       return;
   }
 }
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+template void igl::all_edges<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> >&);
+#endif

include/igl/all_edges.h

@@ -22,9 +22,10 @@ namespace igl
   // directed edge including repeats (meaning interior edges on a surface will
   // show up once for each direction and non-manifold edges may appear more than
   // once for each direction).
+  template <typename DerivedF, typename DerivedE>
   IGL_INLINE void all_edges(
-    const Eigen::MatrixXi & F,
-    Eigen::MatrixXi & E);
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    Eigen::PlainObjectBase<DerivedE> & E);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/boost/components.h

@@ -23,7 +23,8 @@ namespace igl
   IGL_INLINE void components(
     const Eigen::SparseMatrix<AScalar> & A,
     Eigen::PlainObjectBase<DerivedC> & C);
-  // Ditto but for mesh faces as input
+  // Ditto but for mesh faces as input. This computes connected components of
+  // **vertices** where **edges** establish connectivity.
   //
   // Inputs:
   //   F  n by 3 list of triangle indices

include/igl/cgal/CGAL_includes.hpp

@@ -31,9 +31,10 @@
 #include <CGAL/AABB_traits.h>
 #include <CGAL/AABB_triangle_primitive.h>
 
-// Boolean operations
-#include <CGAL/Polyhedron_3.h>
-#include <CGAL/Nef_polyhedron_3.h>
+// ALEC: IS ANY IGL FUNCTION USING THESE?
+//// Boolean operations
+//#include <CGAL/Polyhedron_3.h>
+//#include <CGAL/Nef_polyhedron_3.h>
 
 // Delaunay Triangulation in 3D
 #include <CGAL/Delaunay_triangulation_3.h>

include/igl/cgal/SelfIntersectMesh.h

@@ -29,9 +29,27 @@ namespace igl
   // or 
   //     CGAL::Exact_predicates_exact_constructions_kernel
 
-  template <typename Kernel>
+  template <
+    typename Kernel,
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedVV,
+    typename DerivedFF,
+    typename DerivedIF,
+    typename DerivedJ,
+    typename DerivedIM>
   class SelfIntersectMesh
   {
+    typedef 
+      SelfIntersectMesh<
+      Kernel,
+      DerivedV,
+      DerivedF,
+      DerivedVV,
+      DerivedFF,
+      DerivedIF,
+      DerivedJ,
+      DerivedIM> Self;
     public:
       // 3D Primitives
       typedef CGAL::Point_3<Kernel>    Point_3;
@@ -39,8 +57,8 @@ namespace igl
       typedef CGAL::Triangle_3<Kernel> Triangle_3; 
       typedef CGAL::Plane_3<Kernel>    Plane_3;
       typedef CGAL::Tetrahedron_3<Kernel> Tetrahedron_3; 
-      typedef CGAL::Polyhedron_3<Kernel> Polyhedron_3; 
-      typedef CGAL::Nef_polyhedron_3<Kernel> Nef_polyhedron_3; 
+      //typedef CGAL::Polyhedron_3<Kernel> Polyhedron_3; 
+      //typedef CGAL::Nef_polyhedron_3<Kernel> Nef_polyhedron_3; 
       // 2D Primitives
       typedef CGAL::Point_2<Kernel>    Point_2;
       typedef CGAL::Segment_2<Kernel>  Segment_2; 
@@ -62,21 +80,22 @@ namespace igl
         Box;
 
       // Input mesh
-      const Eigen::MatrixXd & V;
-      const Eigen::MatrixXi & F;
+      const Eigen::PlainObjectBase<DerivedV> & V;
+      const Eigen::PlainObjectBase<DerivedF> & F;
       // Number of self-intersecting triangle pairs
-      int count;
+      typedef typename DerivedF::Index Index;
+      Index count;
       std::vector<std::list<CGAL::Object> > F_objects;
       Triangles T;
-      std::list<int> lIF;
+      std::list<Index> lIF;
       std::vector<bool> offensive;
-      std::vector<int> offending_index;
-      std::vector<int> offending;
+      std::vector<Index> offending_index;
+      std::vector<Index> offending;
       // Make a short name for the edge map's key
-      typedef std::pair<int,int> EMK;
+      typedef std::pair<Index,Index> EMK;
       // Make a short name for the type stored at each edge, the edge map's
       // value
-      typedef std::list<int> EMV;
+      typedef std::list<Index> EMV;
       // Make a short name for the edge map
       typedef std::map<EMK,EMV> EdgeMap;
       EdgeMap edge2faces;
@@ -88,27 +107,26 @@ namespace igl
       //
       // See also: remesh_self_intersections.h
       inline SelfIntersectMesh(
-        const Eigen::MatrixXd & V,
-        const Eigen::MatrixXi & F,
-        const RemeshSelfIntersectionsParam & params,
-        Eigen::MatrixXd & VV,
-        Eigen::MatrixXi & FF,
-        Eigen::MatrixXi & IF,
-        Eigen::VectorXi & J,
-        Eigen::VectorXi & IM
-        );
+          const Eigen::PlainObjectBase<DerivedV> & V,
+          const Eigen::PlainObjectBase<DerivedF> & F,
+          const RemeshSelfIntersectionsParam & params,
+          Eigen::PlainObjectBase<DerivedVV> & VV,
+          Eigen::PlainObjectBase<DerivedFF> & FF,
+          Eigen::PlainObjectBase<DerivedIF> & IF,
+          Eigen::PlainObjectBase<DerivedJ> & J,
+          Eigen::PlainObjectBase<DerivedIM> & IM);
     private:
       // Helper function to mark a face as offensive
       //
       // Inputs:
       //   f  index of face in F
-      inline void mark_offensive(const int f);
+      inline void mark_offensive(const Index f);
       // Helper function to count intersections between faces
       //
       // Input:
       //   fa  index of face A in F
       //   fb  index of face B in F
-      inline void count_intersection(const int fa,const int fb);
+      inline void count_intersection( const Index fa, const Index fb);
       // Helper function for box_intersect. Intersect two triangles A and B,
       // append the intersection object (point,segment,triangle) to a running
       // list for A and B
@@ -123,8 +141,8 @@ namespace igl
       inline bool intersect(
           const Triangle_3 & A, 
           const Triangle_3 & B, 
-          const int fa,
-          const int fb);
+          const Index fa,
+          const Index fb);
       // Helper function for box_intersect. In the case where A and B have
       // already been identified to share a vertex, then we only want to add
       // possible segment intersections. Assumes truly duplicate triangles are
@@ -143,17 +161,17 @@ namespace igl
       inline bool single_shared_vertex(
           const Triangle_3 & A,
           const Triangle_3 & B,
-          const int fa,
-          const int fb,
-          const int va,
-          const int vb);
+          const Index fa,
+          const Index fb,
+          const Index va,
+          const Index vb);
       // Helper handling one direction
       inline bool single_shared_vertex(
           const Triangle_3 & A,
           const Triangle_3 & B,
-          const int fa,
-          const int fb,
-          const int va);
+          const Index fa,
+          const Index fb,
+          const Index va);
       // Helper function for box_intersect. In the case where A and B have
       // already been identified to share two vertices, then we only want to add
       // a possible coplanar (Triangle) intersection. Assumes truly degenerate
@@ -161,8 +179,8 @@ namespace igl
       inline bool double_shared_vertex(
           const Triangle_3 & A,
           const Triangle_3 & B,
-          const int fa,
-          const int fb);
+          const Index fa,
+          const Index fb);
 
     public:
       // Callback function called during box self intersections test. Means
@@ -239,25 +257,57 @@ namespace igl
 //  boost::function<void(const Box &a,const Box &b)> cb
 //    = boost::bind(&::box_intersect, this, _1,_2);
 //   
-template <typename Kernel>
-inline void igl::SelfIntersectMesh<Kernel>::box_intersect(
-  igl::SelfIntersectMesh<Kernel> * SIM, 
-  const typename igl::SelfIntersectMesh<Kernel>::Box &a, 
-  const typename igl::SelfIntersectMesh<Kernel>::Box &b)
+template <
+  typename Kernel,
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedVV,
+  typename DerivedFF,
+  typename DerivedIF,
+  typename DerivedJ,
+  typename DerivedIM>
+inline void igl::SelfIntersectMesh<
+  Kernel,
+  DerivedV,
+  DerivedF,
+  DerivedVV,
+  DerivedFF,
+  DerivedIF,
+  DerivedJ,
+  DerivedIM>::box_intersect(
+  Self * SIM, 
+  const typename Self::Box &a, 
+  const typename Self::Box &b)
 {
   SIM->box_intersect(a,b);
 }
 
-template <typename Kernel>
-inline igl::SelfIntersectMesh<Kernel>::SelfIntersectMesh(
-  const Eigen::MatrixXd & V,
-  const Eigen::MatrixXi & F,
+template <
+  typename Kernel,
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedVV,
+  typename DerivedFF,
+  typename DerivedIF,
+  typename DerivedJ,
+  typename DerivedIM>
+inline igl::SelfIntersectMesh<
+  Kernel,
+  DerivedV,
+  DerivedF,
+  DerivedVV,
+  DerivedFF,
+  DerivedIF,
+  DerivedJ,
+  DerivedIM>::SelfIntersectMesh(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
   const RemeshSelfIntersectionsParam & params,
-  Eigen::MatrixXd & VV,
-  Eigen::MatrixXi & FF,
-  Eigen::MatrixXi & IF,
-  Eigen::VectorXi & J,
-  Eigen::VectorXi & IM):
+  Eigen::PlainObjectBase<DerivedVV> & VV,
+  Eigen::PlainObjectBase<DerivedFF> & FF,
+  Eigen::PlainObjectBase<DerivedIF> & IF,
+  Eigen::PlainObjectBase<DerivedJ> & J,
+  Eigen::PlainObjectBase<DerivedIM> & IM):
   V(V),
   F(F),
   count(0),
@@ -305,9 +355,9 @@ inline igl::SelfIntersectMesh<Kernel>::SelfIntersectMesh(
   assert(lIF.size()%2 == 0);
   IF.resize(lIF.size()/2,2);
   {
-    int i=0;
+    Index i=0;
     for(
-      typename list<int>::const_iterator ifit = lIF.begin();
+      typename list<Index>::const_iterator ifit = lIF.begin();
       ifit!=lIF.end();
       )
     {
@@ -326,19 +376,21 @@ inline igl::SelfIntersectMesh<Kernel>::SelfIntersectMesh(
 
   int NF_count = 0;
   // list of new faces, we'll fix F later
-  vector<MatrixXi> NF(offending.size());
+  vector<
+    typename Eigen::Matrix<typename DerivedFF::Scalar,Dynamic,Dynamic>
+    > NF(offending.size());
   // list of new vertices
   list<Point_3> NV;
-  int NV_count = 0;
+  Index NV_count = 0;
   vector<CDT_plus_2> cdt(offending.size());
   vector<Plane_3> P(offending.size());
   // Use map for *all* faces
-  map<typename CDT_plus_2::Vertex_handle,int> v2i;
+  map<typename CDT_plus_2::Vertex_handle,Index> v2i;
   // Loop over offending triangles
-  for(int o = 0;o<(int)offending.size();o++)
+  for(Index o = 0;o<(Index)offending.size();o++)
   {
     // index in F
-    const int f = offending[o];
+    const Index f = offending[o];
     projected_delaunay(T[f],F_objects[f],cdt[o]);
     // Q: Is this also delaunay in 3D?
     // A: No, because the projection is affine and delaunay is not affine
@@ -349,7 +401,7 @@ inline igl::SelfIntersectMesh<Kernel>::SelfIntersectMesh(
     P[o] = Plane_3(T[f].vertex(0),T[f].vertex(1),T[f].vertex(2));
     // Build index map
     {
-      int i=0;
+      Index i=0;
       for(
         typename CDT_plus_2::Finite_vertices_iterator vit = cdt[o].finite_vertices_begin();
         vit != cdt[o].finite_vertices_end();
@@ -377,8 +429,8 @@ inline igl::SelfIntersectMesh<Kernel>::SelfIntersectMesh(
           for(int e = 0; e<3 && !found;e++)
           {
             // Index of F's eth edge in V
-            int i = F(f,(e+1)%3);
-            int j = F(f,(e+2)%3);
+            Index i = F(f,(e+1)%3);
+            Index j = F(f,(e+2)%3);
             // Be sure that i<j
             if(i>j)
             {
@@ -387,7 +439,7 @@ inline igl::SelfIntersectMesh<Kernel>::SelfIntersectMesh(
             assert(edge2faces.count(EMK(i,j))==1);
             // loop over neighbors
             for(
-              list<int>::const_iterator nit = edge2faces[EMK(i,j)].begin();
+              typename list<Index>::const_iterator nit = edge2faces[EMK(i,j)].begin();
               nit != edge2faces[EMK(i,j)].end() && !found;
               nit++)
             {
@@ -397,7 +449,7 @@ inline igl::SelfIntersectMesh<Kernel>::SelfIntersectMesh(
                 continue;
               }
               // index of neighbor in offending (to find its cdt)
-              int no = offending_index[*nit];
+              Index no = offending_index[*nit];
               // Loop over vertices of that neighbor's cdt (might not have been
               // processed yet, but then it's OK because it'll just be empty)
               for(
@@ -425,7 +477,7 @@ inline igl::SelfIntersectMesh<Kernel>::SelfIntersectMesh(
       }
     }
     {
-      int i = 0;
+      Index i = 0;
       // Resize to fit new number of triangles
       NF[o].resize(cdt[o].number_of_faces(),3);
       NF_count+=NF[o].rows();
@@ -442,16 +494,16 @@ inline igl::SelfIntersectMesh<Kernel>::SelfIntersectMesh(
       }
     }
   }
-  assert(NV_count == (int)NV.size());
+  assert(NV_count == (Index)NV.size());
   // Build output
 #ifndef NDEBUG
   {
-    int off_count = 0;
-    for(int f = 0;f<F.rows();f++)
+    Index off_count = 0;
+    for(Index f = 0;f<F.rows();f++)
     {
       off_count+= (offensive[f]?1:0);
     }
-    assert(off_count==(int)offending.size());
+    assert(off_count==(Index)offending.size());
   }
 #endif
   // Append faces
@@ -459,8 +511,8 @@ inline igl::SelfIntersectMesh<Kernel>::SelfIntersectMesh(
   J.resize(FF.rows());
   // First append non-offending original faces
   // There's an Eigen way to do this in one line but I forget
-  int off = 0;
-  for(int f = 0;f<F.rows();f++)
+  Index off = 0;
+  for(Index f = 0;f<F.rows();f++)
   {
     if(!offensive[f])
     {
@@ -469,9 +521,9 @@ inline igl::SelfIntersectMesh<Kernel>::SelfIntersectMesh(
       off++;
     }
   }
-  assert(off == (int)(F.rows()-offending.size()));
+  assert(off == (Index)(F.rows()-offending.size()));
   // Now append replacement faces for offending faces
-  for(int o = 0;o<(int)offending.size();o++)
+  for(Index o = 0;o<(Index)offending.size();o++)
   {
     FF.block(off,0,NF[o].rows(),3) = NF[o];
     J.block(off,0,NF[o].rows(),1).setConstant(offending[o]);
@@ -481,13 +533,13 @@ inline igl::SelfIntersectMesh<Kernel>::SelfIntersectMesh(
   VV.resize(V.rows()+NV_count,3);
   VV.block(0,0,V.rows(),3) = V;
   {
-    int i = 0;
+    Index i = 0;
     for(
       typename list<Point_3>::const_iterator nvit = NV.begin();
       nvit != NV.end();
       nvit++)
     {
-      for(int d = 0;d<3;d++)
+      for(Index d = 0;d<3;d++)
       {
         const Point_3 & p = *nvit;
         VV(V.rows()+i,d) = CGAL::to_double(p[d]);
@@ -502,10 +554,10 @@ inline igl::SelfIntersectMesh<Kernel>::SelfIntersectMesh(
     }
   }
   IM.resize(VV.rows(),1);
-  map<Point_3,int> vv2i;
+  map<Point_3,Index> vv2i;
   // Safe to check for duplicates using double for original vertices: if
   // incoming reps are different then the points are unique.
-  for(int v = 0;v<V.rows();v++)
+  for(Index v = 0;v<V.rows();v++)
   {
     const Point_3 p(V(v,0),V(v,1),V(v,2));
     if(vv2i.count(p)==0)
@@ -517,7 +569,7 @@ inline igl::SelfIntersectMesh<Kernel>::SelfIntersectMesh(
   }
   // Must check for duplicates of new vertices using exact.
   {
-    int v = V.rows();
+    Index v = V.rows();
     for(
       typename list<Point_3>::const_iterator nvit = NV.begin();
       nvit != NV.end();
@@ -547,8 +599,24 @@ inline igl::SelfIntersectMesh<Kernel>::SelfIntersectMesh(
 }
 
 
-template <typename Kernel>
-inline void igl::SelfIntersectMesh<Kernel>::mark_offensive(const int f)
+template <
+  typename Kernel,
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedVV,
+  typename DerivedFF,
+  typename DerivedIF,
+  typename DerivedJ,
+  typename DerivedIM>
+inline void igl::SelfIntersectMesh<
+  Kernel,
+  DerivedV,
+  DerivedF,
+  DerivedVV,
+  DerivedFF,
+  DerivedIF,
+  DerivedJ,
+  DerivedIM>::mark_offensive(const Index f)
 {
   using namespace std;
   lIF.push_back(f);
@@ -558,11 +626,11 @@ inline void igl::SelfIntersectMesh<Kernel>::mark_offensive(const int f)
     offending_index[f]=offending.size();
     offending.push_back(f);
     // Add to edge map
-    for(int e = 0; e<3;e++)
+    for(Index e = 0; e<3;e++)
     {
       // Index of F's eth edge in V
-      int i = F(f,(e+1)%3);
-      int j = F(f,(e+2)%3);
+      Index i = F(f,(e+1)%3);
+      Index j = F(f,(e+2)%3);
       // Be sure that i<j
       if(i>j)
       {
@@ -571,7 +639,7 @@ inline void igl::SelfIntersectMesh<Kernel>::mark_offensive(const int f)
       // Create new list if there is no entry
       if(edge2faces.count(EMK(i,j))==0)
       {
-        edge2faces[EMK(i,j)] = list<int>();
+        edge2faces[EMK(i,j)] = list<Index>();
       }
       // append to list
       edge2faces[EMK(i,j)].push_back(f);
@@ -579,10 +647,26 @@ inline void igl::SelfIntersectMesh<Kernel>::mark_offensive(const int f)
   }
 }
 
-template <typename Kernel>
-inline void igl::SelfIntersectMesh<Kernel>::count_intersection(
-  const int fa,
-  const int fb)
+template <
+  typename Kernel,
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedVV,
+  typename DerivedFF,
+  typename DerivedIF,
+  typename DerivedJ,
+  typename DerivedIM>
+inline void igl::SelfIntersectMesh<
+  Kernel,
+  DerivedV,
+  DerivedF,
+  DerivedVV,
+  DerivedFF,
+  DerivedIF,
+  DerivedJ,
+  DerivedIM>::count_intersection(
+  const Index fa,
+  const Index fb)
 {
   mark_offensive(fa);
   mark_offensive(fb);
@@ -594,12 +678,28 @@ inline void igl::SelfIntersectMesh<Kernel>::count_intersection(
   }
 }
 
-template <typename Kernel>
-inline bool igl::SelfIntersectMesh<Kernel>::intersect(
+template <
+  typename Kernel,
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedVV,
+  typename DerivedFF,
+  typename DerivedIF,
+  typename DerivedJ,
+  typename DerivedIM>
+inline bool igl::SelfIntersectMesh<
+  Kernel,
+  DerivedV,
+  DerivedF,
+  DerivedVV,
+  DerivedFF,
+  DerivedIF,
+  DerivedJ,
+  DerivedIM>::intersect(
   const Triangle_3 & A, 
   const Triangle_3 & B, 
-  const int fa,
-  const int fb)
+  const Index fa,
+  const Index fb)
 {
   // Determine whether there is an intersection
   if(!CGAL::do_intersect(A,B))
@@ -617,14 +717,30 @@ inline bool igl::SelfIntersectMesh<Kernel>::intersect(
   return true;
 }
 
-template <typename Kernel>
-inline bool igl::SelfIntersectMesh<Kernel>::single_shared_vertex(
+template <
+  typename Kernel,
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedVV,
+  typename DerivedFF,
+  typename DerivedIF,
+  typename DerivedJ,
+  typename DerivedIM>
+inline bool igl::SelfIntersectMesh<
+  Kernel,
+  DerivedV,
+  DerivedF,
+  DerivedVV,
+  DerivedFF,
+  DerivedIF,
+  DerivedJ,
+  DerivedIM>::single_shared_vertex(
   const Triangle_3 & A,
   const Triangle_3 & B,
-  const int fa,
-  const int fb,
-  const int va,
-  const int vb)
+  const Index fa,
+  const Index fb,
+  const Index va,
+  const Index vb)
 {
   ////using namespace std;
   //CGAL::Object result = CGAL::intersection(A,B);
@@ -647,13 +763,29 @@ inline bool igl::SelfIntersectMesh<Kernel>::single_shared_vertex(
   return single_shared_vertex(B,A,fb,fa,vb);
 }
 
-template <typename Kernel>
-inline bool igl::SelfIntersectMesh<Kernel>::single_shared_vertex(
+template <
+  typename Kernel,
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedVV,
+  typename DerivedFF,
+  typename DerivedIF,
+  typename DerivedJ,
+  typename DerivedIM>
+inline bool igl::SelfIntersectMesh<
+  Kernel,
+  DerivedV,
+  DerivedF,
+  DerivedVV,
+  DerivedFF,
+  DerivedIF,
+  DerivedJ,
+  DerivedIM>::single_shared_vertex(
   const Triangle_3 & A,
   const Triangle_3 & B,
-  const int fa,
-  const int fb,
-  const int va)
+  const Index fa,
+  const Index fb,
+  const Index va)
 {
   // This was not a good idea. It will not handle coplanar triangles well.
   using namespace std;
@@ -712,12 +844,28 @@ inline bool igl::SelfIntersectMesh<Kernel>::single_shared_vertex(
 }
 
 
-template <typename Kernel>
-inline bool igl::SelfIntersectMesh<Kernel>::double_shared_vertex(
+template <
+  typename Kernel,
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedVV,
+  typename DerivedFF,
+  typename DerivedIF,
+  typename DerivedJ,
+  typename DerivedIM>
+inline bool igl::SelfIntersectMesh<
+  Kernel,
+  DerivedV,
+  DerivedF,
+  DerivedVV,
+  DerivedFF,
+  DerivedIF,
+  DerivedJ,
+  DerivedIM>::double_shared_vertex(
   const Triangle_3 & A,
   const Triangle_3 & B,
-  const int fa,
-  const int fb)
+  const Index fa,
+  const Index fb)
 {
   using namespace std;
   // Cheaper way to do this than calling do_intersect?
@@ -775,15 +923,38 @@ inline bool igl::SelfIntersectMesh<Kernel>::double_shared_vertex(
   return false;
 }
 
-template <typename Kernel>
-inline void igl::SelfIntersectMesh<Kernel>::box_intersect(
+template <
+  typename Kernel,
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedVV,
+  typename DerivedFF,
+  typename DerivedIF,
+  typename DerivedJ,
+  typename DerivedIM>
+inline void igl::SelfIntersectMesh<
+  Kernel,
+  DerivedV,
+  DerivedF,
+  DerivedVV,
+  DerivedFF,
+  DerivedIF,
+  DerivedJ,
+  DerivedIM>::box_intersect(
   const Box& a, 
   const Box& b)
 {
   using namespace std;
+  // Could we write this as a static function of:
+  //
+  // F.row(fa)
+  // F.row(fb)
+  // A
+  // B
+
   // index in F and T
-  int fa = a.handle()-T.begin();
-  int fb = b.handle()-T.begin();
+  Index fa = a.handle()-T.begin();
+  Index fb = b.handle()-T.begin();
   const Triangle_3 & A = *a.handle();
   const Triangle_3 & B = *b.handle();
   // I'm not going to deal with degenerate triangles, though at some point we
@@ -791,14 +962,14 @@ inline void igl::SelfIntersectMesh<Kernel>::box_intersect(
   assert(!a.handle()->is_degenerate());
   assert(!b.handle()->is_degenerate());
   // Number of combinatorially shared vertices
-  int comb_shared_vertices = 0;
+  Index comb_shared_vertices = 0;
   // Number of geometrically shared vertices (*not* including combinatorially
   // shared)
-  int geo_shared_vertices = 0;
+  Index geo_shared_vertices = 0;
   // Keep track of shared vertex indices (we only handles single shared
   // vertices as a special case, so just need last/first/only ones)
-  int va=-1,vb=-1;
-  int ea,eb;
+  Index va=-1,vb=-1;
+  Index ea,eb;
   for(ea=0;ea<3;ea++)
   {
     for(eb=0;eb<3;eb++)
@@ -816,7 +987,7 @@ inline void igl::SelfIntersectMesh<Kernel>::box_intersect(
       }
     }
   }
-  const int total_shared_vertices = comb_shared_vertices + geo_shared_vertices;
+  const Index total_shared_vertices = comb_shared_vertices + geo_shared_vertices;
   if(comb_shared_vertices== 3)
   {
     // Combinatorially duplicate face, these should be removed by preprocessing
@@ -906,8 +1077,24 @@ done:
 //   A_objects_3  updated list of intersection objects for A
 // Outputs:
 //   cdt  Contrained delaunay triangulation in projected 2D plane
-template <typename Kernel>
-inline void igl::SelfIntersectMesh<Kernel>::projected_delaunay(
+template <
+  typename Kernel,
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedVV,
+  typename DerivedFF,
+  typename DerivedIF,
+  typename DerivedJ,
+  typename DerivedIM>
+inline void igl::SelfIntersectMesh<
+  Kernel,
+  DerivedV,
+  DerivedF,
+  DerivedVV,
+  DerivedFF,
+  DerivedIF,
+  DerivedJ,
+  DerivedIM>::projected_delaunay(
   const Triangle_3 & A,
   const std::list<CGAL::Object> & A_objects_3,
   CDT_plus_2 & cdt)
@@ -918,12 +1105,12 @@ inline void igl::SelfIntersectMesh<Kernel>::projected_delaunay(
   Plane_3 P(A.vertex(0),A.vertex(1),A.vertex(2));
   // Insert triangle into vertices
   typename CDT_plus_2::Vertex_handle corners[3];
-  for(int i = 0;i<3;i++)
+  for(Index i = 0;i<3;i++)
   {
     corners[i] = cdt.insert(P.to_2d(A.vertex(i)));
   }
   // Insert triangle edges as constraints
-  for(int i = 0;i<3;i++)
+  for(Index i = 0;i<3;i++)
   {
     cdt.insert_constraint( corners[(i+1)%3], corners[(i+2)%3]);
   }
@@ -953,11 +1140,11 @@ inline void igl::SelfIntersectMesh<Kernel>::projected_delaunay(
     {
       //cerr<<REDRUM("Poly...")<<endl;
       const std::vector<Point_3 > & poly = *polyp;
-      const int m = poly.size();
+      const Index m = poly.size();
       assert(m>=2);
-      for(int p = 0;p<m;p++)
+      for(Index p = 0;p<m;p++)
       {
-        const int np = (p+1)%m;
+        const Index np = (p+1)%m;
         cdt.insert_constraint(P.to_2d(poly[p]),P.to_2d(poly[np]));
       }
     }else

include/igl/cgal/remesh_self_intersections.cpp

@@ -10,15 +10,23 @@
 #include <igl/C_STR.h>
 #include <list>
 
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedVV,
+  typename DerivedFF,
+  typename DerivedIF,
+  typename DerivedJ,
+  typename DerivedIM>
 IGL_INLINE void igl::remesh_self_intersections(
-  const Eigen::MatrixXd & V,
-  const Eigen::MatrixXi & F,
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
   const RemeshSelfIntersectionsParam & params,
-  Eigen::MatrixXd & VV,
-  Eigen::MatrixXi & FF,
-  Eigen::MatrixXi & IF,
-  Eigen::VectorXi & J,
-  Eigen::VectorXi & IM)
+  Eigen::PlainObjectBase<DerivedVV> & VV,
+  Eigen::PlainObjectBase<DerivedFF> & FF,
+  Eigen::PlainObjectBase<DerivedIF> & IF,
+  Eigen::PlainObjectBase<DerivedJ> & J,
+  Eigen::PlainObjectBase<DerivedIM> & IM)
 {
   using namespace std;
   if(params.detect_only)
@@ -36,7 +44,18 @@ IGL_INLINE void igl::remesh_self_intersections(
 //#endif
 
     typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
-    SelfIntersectMesh<Kernel> SIM = SelfIntersectMesh<Kernel>(V,F,params,VV,FF,IF,J,IM);
+    typedef 
+      SelfIntersectMesh<
+        Kernel,
+        DerivedV,
+        DerivedF,
+        DerivedVV,
+        DerivedFF,
+        DerivedIF,
+        DerivedJ,
+        DerivedIM> 
+      SelfIntersectMeshK;
+    SelfIntersectMeshK SIM = SelfIntersectMeshK(V,F,params,VV,FF,IF,J,IM);
 
 //#ifdef __APPLE__
 //    signal(SIGFPE,SIG_DFL);
@@ -45,6 +64,22 @@ IGL_INLINE void igl::remesh_self_intersections(
   }else
   {
     typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
-    SelfIntersectMesh<Kernel> SIM = SelfIntersectMesh<Kernel>(V,F,params,VV,FF,IF,J,IM);
+    typedef 
+      SelfIntersectMesh<
+        Kernel,
+        DerivedV,
+        DerivedF,
+        DerivedVV,
+        DerivedFF,
+        DerivedIF,
+        DerivedJ,
+        DerivedIM> 
+      SelfIntersectMeshK;
+    SelfIntersectMeshK SIM = SelfIntersectMeshK(V,F,params,VV,FF,IF,J,IM);
   }
 }
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+template void igl::remesh_self_intersections<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::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::RemeshSelfIntersectionsParam 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::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/cgal/remesh_self_intersections.h

@@ -53,15 +53,23 @@ namespace igl
   // any resulting additional vertices along that edge may not get properly
   // connected so that the output mesh has the same global topology. This is
   // because 
+  template <
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedVV,
+    typename DerivedFF,
+    typename DerivedIF,
+    typename DerivedJ,
+    typename DerivedIM>
   IGL_INLINE void remesh_self_intersections(
-    const Eigen::MatrixXd & V,
-    const Eigen::MatrixXi & F,
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedF> & F,
     const RemeshSelfIntersectionsParam & params,
-    Eigen::MatrixXd & VV,
-    Eigen::MatrixXi & FF,
-    Eigen::MatrixXi & IF,
-    Eigen::VectorXi & J,
-    Eigen::VectorXi & IM);
+    Eigen::PlainObjectBase<DerivedVV> & VV,
+    Eigen::PlainObjectBase<DerivedFF> & FF,
+    Eigen::PlainObjectBase<DerivedIF> & IF,
+    Eigen::PlainObjectBase<DerivedJ> & J,
+    Eigen::PlainObjectBase<DerivedIM> & IM);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/is_edge_manifold.cpp

@@ -6,12 +6,56 @@
 // 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 "is_edge_manifold.h"
+#include "all_edges.h"
+#include "unique_simplices.h"
 
 #include <algorithm>
+#include <vector>
+
+template <
+  typename DerivedF, 
+  typename DerivedBF,
+  typename DerivedE,
+  typename DerivedEMAP,
+  typename DerivedBE>
+IGL_INLINE bool igl::is_edge_manifold(
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedBF>& BF,
+  Eigen::PlainObjectBase<DerivedE>& E,
+  Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+  Eigen::PlainObjectBase<DerivedBE>& BE)
+{
+  using namespace Eigen;
+  typedef typename DerivedF::Index Index;
+  typedef Matrix<typename DerivedF::Scalar,Dynamic,1> VectorXF;
+  typedef Matrix<typename DerivedF::Scalar,Dynamic,2> MatrixXF2;
+  MatrixXF2 allE;
+  all_edges(F,allE);
+  // Find unique undirected edges and mapping
+  VectorXF _;
+  unique_simplices(allE,E,_,EMAP);
+  std::vector<typename DerivedF::Index> count(E.rows(),0);
+  for(Index e = 0;e<EMAP.rows();e++)
+  {
+    count[EMAP[e]]++;
+  }
+  const Index m = F.rows();
+  BF.resize(m,3);
+  BE.resize(E.rows(),1);
+  bool all = true;
+  for(Index e = 0;e<EMAP.rows();e++)
+  {
+    const bool manifold = count[EMAP[e]] <= 2;
+    all &= BF(e%m,e/m) = manifold;
+    BE(EMAP[e]) = manifold;
+  }
+  return all;
+}
 
 template <typename DerivedV, typename DerivedF>
-IGL_INLINE bool igl::is_edge_manifold(const Eigen::PlainObjectBase<DerivedV>& /*V*/,
-                                 const Eigen::PlainObjectBase<DerivedF>& F)
+IGL_INLINE bool igl::is_edge_manifold(
+  const Eigen::PlainObjectBase<DerivedV>& /*V*/,
+  const Eigen::PlainObjectBase<DerivedF>& F)
 {
   // List of edges (i,j,f,c) where edge i<j is associated with corner i of face
   // f

include/igl/is_edge_manifold.h

@@ -10,7 +10,6 @@
 #include "igl_inline.h"
 
 #include <Eigen/Core>
-#include <vector>
 
 namespace igl 
 {
@@ -25,6 +24,18 @@ namespace igl
   //  Does not check for non-manifold vertices
   //
   // See also: is_vertex_manifold
+  template <
+    typename DerivedF, 
+    typename DerivedBF,
+    typename DerivedE,
+    typename DerivedEMAP,
+    typename DerivedBE>
+  IGL_INLINE bool is_edge_manifold(
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    Eigen::PlainObjectBase<DerivedBF>& BF,
+    Eigen::PlainObjectBase<DerivedE>& E,
+    Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+    Eigen::PlainObjectBase<DerivedBE>& BE);
   template <typename DerivedV, typename DerivedF>
   IGL_INLINE bool is_edge_manifold(
     const Eigen::PlainObjectBase<DerivedV>& V,

include/igl/is_vertex_manifold.cpp

@@ -0,0 +1,89 @@
+#include "is_vertex_manifold.h"
+#include "triangle_triangle_adjacency.h"
+#include "vertex_triangle_adjacency.h"
+#include "unique.h"
+#include <vector>
+#include <cassert>
+#include <map>
+#include <queue>
+#include <iostream>
+
+template <typename DerivedF,typename DerivedB>
+IGL_INLINE bool igl::is_vertex_manifold(
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedB>& B)
+{
+  using namespace std;
+  using namespace Eigen;
+  assert(F.cols() == 3 && "F must contain triangles");
+  typedef typename DerivedF::Scalar Index;
+  typedef typename DerivedF::Index FIndex;
+  const FIndex m = F.rows();
+  const Index n = F.maxCoeff()+1;
+  vector<vector<vector<FIndex > > > TT;
+  vector<vector<vector<FIndex > > > TTi;
+  triangle_triangle_adjacency(F,TT,TTi);
+
+  vector<vector<FIndex > > V2F,_1;
+  vertex_triangle_adjacency(n,F,V2F,_1);
+
+  const auto & check_vertex = [&](const Index v)->bool
+  {
+    vector<FIndex> uV2Fv;
+    {
+      vector<size_t> _1,_2;
+      unique(V2F[v],uV2Fv,_1,_2);
+    }
+    const FIndex one_ring_size = uV2Fv.size();
+    if(one_ring_size == 0)
+    {
+      return false;
+    }
+    const FIndex g = uV2Fv[0];
+    queue<Index> Q;
+    Q.push(g);
+    map<FIndex,bool> seen;
+    while(!Q.empty())
+    {
+      const FIndex f = Q.front();
+      Q.pop();
+      if(seen.count(f)==1)
+      {
+        continue;
+      }
+      seen[f] = true;
+      // Face f's neighbor lists opposite opposite each corner
+      for(const auto & c : TT[f])
+      {
+        // Each neighbor
+        for(const auto & n : c)
+        {
+          bool contains_v = false;
+          for(Index nc = 0;nc<F.cols();nc++)
+          {
+            if(F(n,nc) == v)
+            {
+              contains_v = true;
+              break;
+            }
+          }
+          if(seen.count(n)==0 && contains_v)
+          {
+            Q.push(n);
+          }
+        }
+      }
+    }
+    return one_ring_size == seen.size();
+  };
+
+  // Unreferenced vertices are considered non-manifold
+  B.setConstant(n,1,false);
+  // Loop over all vertices touched by F
+  bool all = true;
+  for(Index v = 0;v<n;v++)
+  {
+    all &= B(v) = check_vertex(v);
+  }
+  return all;
+}

include/igl/is_vertex_manifold.h

@@ -0,0 +1,39 @@
+// 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 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_IS_VERTEX_MANIFOLD_H
+#define IGL_IS_VERTEX_MANIFOLD_H
+#include "igl_inline.h"
+
+#include <Eigen/Core>
+
+namespace igl 
+{
+  // Check if a mesh is vertex-manifold. This only checks whether the faces
+  // incident on each vertex form exactly one connected component. Vertices
+  // incident on non-manifold edges are not consider non-manifold by this
+  // function (see is_edge_manifold.h). Unreferenced verties are considered
+  // non-manifold (zero components).
+  //
+  // Inputs:
+  //   F  #F by 3 list of triangle indices
+  // Outputs:
+  //   B  #V list whether 
+  // Returns whether mesh is vertex manifold.
+  //
+  // See also: is_edge_manifold
+  template <typename DerivedF,typename DerivedB>
+  IGL_INLINE bool is_vertex_manifold(
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    Eigen::PlainObjectBase<DerivedB>& B);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "is_vertex_manifold.cpp"
+#endif
+
+#endif

include/igl/next_filename.cpp

@@ -24,7 +24,7 @@ bool igl::next_filename(
   while(true)
   {
     next = STR(prefix << setfill('0') << setw(zeros)<<i<<suffix);
-    if(!file_exists(next.c_str()))
+    if(!file_exists(next))
     {
       return true;
     }

include/igl/per_vertex_normals.cpp

@@ -32,13 +32,13 @@ IGL_INLINE void igl::per_vertex_normals(
   return per_vertex_normals(V,F,PER_VERTEX_NORMALS_WEIGHTING_TYPE_DEFAULT,N);
 }
 
-template <typename DerivedV, typename DerivedF>
+template <typename DerivedV, typename DerivedF, typename DerivedFN, typename DerivedN>
 IGL_INLINE void igl::per_vertex_normals(
   const Eigen::PlainObjectBase<DerivedV>& V,
   const Eigen::PlainObjectBase<DerivedF>& F,
   const igl::PerVertexNormalsWeightingType weighting,
-  const Eigen::PlainObjectBase<DerivedV>& FN,
-  Eigen::PlainObjectBase<DerivedV> & N)
+  const Eigen::PlainObjectBase<DerivedFN>& FN,
+  Eigen::PlainObjectBase<DerivedN> & N)
 {
   // Resize for output
   N.setZero(V.rows(),3);
@@ -99,9 +99,4 @@ IGL_INLINE void igl::per_vertex_normals(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 // generated by autoexplicit.sh
-template void igl::per_vertex_normals<Eigen::Matrix<double, -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::PerVertexNormalsWeightingType, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
-template void igl::per_vertex_normals<Eigen::Matrix<double, -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> >&);
-template void igl::per_vertex_normals<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&);
-template void igl::per_vertex_normals<Eigen::Matrix<double, -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> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/per_vertex_normals.h

@@ -46,13 +46,13 @@ namespace igl
     Eigen::PlainObjectBase<DerivedV> & N);
   // Inputs:
   //   FN  #F by 3 matrix of face (triangle) normals
-  template <typename DerivedV, typename DerivedF>
+  template <typename DerivedV, typename DerivedF, typename DerivedFN, typename DerivedN>
   IGL_INLINE void per_vertex_normals(
     const Eigen::PlainObjectBase<DerivedV>& V,
     const Eigen::PlainObjectBase<DerivedF>& F,
     const PerVertexNormalsWeightingType weighting,
-    const Eigen::PlainObjectBase<DerivedV>& FN,
-    Eigen::PlainObjectBase<DerivedV> & N);
+    const Eigen::PlainObjectBase<DerivedFN>& FN,
+    Eigen::PlainObjectBase<DerivedN> & N);
   // Without weighting
   template <typename DerivedV, typename DerivedF>
   IGL_INLINE void per_vertex_normals(

include/igl/sort.cpp

@@ -16,12 +16,12 @@
 #include <algorithm>
 #include <iostream>
 
-template <typename DerivedX, typename DerivedIX>
+template <typename DerivedX, typename DerivedY, typename DerivedIX>
 IGL_INLINE void igl::sort(
   const Eigen::PlainObjectBase<DerivedX>& X,
   const int dim,
   const bool ascending,
-  Eigen::PlainObjectBase<DerivedX>& Y,
+  Eigen::PlainObjectBase<DerivedY>& Y,
   Eigen::PlainObjectBase<DerivedIX>& IX)
 {
   // get number of rows (or columns)
@@ -74,12 +74,12 @@ IGL_INLINE void igl::sort(
   }
 }
 
-template <typename DerivedX, typename DerivedIX>
+template <typename DerivedX, typename DerivedY, typename DerivedIX>
 IGL_INLINE void igl::sort_new(
   const Eigen::PlainObjectBase<DerivedX>& X,
   const int dim,
   const bool ascending,
-  Eigen::PlainObjectBase<DerivedX>& Y,
+  Eigen::PlainObjectBase<DerivedY>& Y,
   Eigen::PlainObjectBase<DerivedIX>& IX)
 {
   // get number of rows (or columns)
@@ -137,12 +137,12 @@ IGL_INLINE void igl::sort_new(
   }
 }
 
-template <typename DerivedX, typename DerivedIX>
+template <typename DerivedX, typename DerivedY, typename DerivedIX>
 IGL_INLINE void igl::sort2(
   const Eigen::PlainObjectBase<DerivedX>& X,
   const int dim,
   const bool ascending,
-  Eigen::PlainObjectBase<DerivedX>& Y,
+  Eigen::PlainObjectBase<DerivedY>& Y,
   Eigen::PlainObjectBase<DerivedIX>& IX)
 {
   using namespace Eigen;
@@ -213,24 +213,10 @@ IGL_INLINE void igl::sort(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 // generated by autoexplicit.sh
-template void igl::sort<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::sort<int>(std::vector<int, std::allocator<int> > const&, bool, std::vector<int, std::allocator<int> >&, std::vector<unsigned long, std::allocator<unsigned long> >&);
 // generated by autoexplicit.sh
-template void igl::sort<Eigen::Matrix<float, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::sort<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::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 // generated by autoexplicit.sh
-template void igl::sort<Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::sort<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 // generated by autoexplicit.sh
-template void igl::sort<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-template void igl::sort<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&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-template void igl::sort<igl::SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> > >(std::vector<igl::SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> >, std::allocator<igl::SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> > > > const&, bool, std::vector<igl::SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> >, std::allocator<igl::SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> > > >&, std::vector<unsigned long, std::allocator<unsigned long> >&);
-template void igl::sort<int>(std::vector<int, std::allocator<int> > const&, bool, std::vector<int, std::allocator<int> >&, std::vector<unsigned long, std::allocator<unsigned long> >&);
-template void igl::sort<igl::SortableRow<Eigen::Matrix<double, -1, 1, 0, -1, 1> > >(std::vector<igl::SortableRow<Eigen::Matrix<double, -1, 1, 0, -1, 1> >, std::allocator<igl::SortableRow<Eigen::Matrix<double, -1, 1, 0, -1, 1> > > > const&, bool, std::vector<igl::SortableRow<Eigen::Matrix<double, -1, 1, 0, -1, 1> >, std::allocator<igl::SortableRow<Eigen::Matrix<double, -1, 1, 0, -1, 1> > > >&, std::vector<unsigned long, std::allocator<unsigned long> >&);
-template void igl::sort2<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&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-template void igl::sort2<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-template void igl::sort_new<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-template void igl::sort<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template void igl::sort<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template void igl::sort<Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template void igl::sort<Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 2, 0, -1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&);
-template void igl::sort<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-template void igl::sort<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&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/sort.h

@@ -30,28 +30,28 @@ namespace igl
   //   Y  m by n matrix whose entries are sorted
   //   IX  m by n matrix of indices so that if dim = 1, then in matlab notation
   //     for j = 1:n, Y(:,j) = X(I(:,j),j); end
-  template <typename DerivedX, typename DerivedIX>
+  template <typename DerivedX, typename DerivedY, typename DerivedIX>
   IGL_INLINE void sort(
     const Eigen::PlainObjectBase<DerivedX>& X,
     const int dim,
     const bool ascending,
-    Eigen::PlainObjectBase<DerivedX>& Y,
+    Eigen::PlainObjectBase<DerivedY>& Y,
     Eigen::PlainObjectBase<DerivedIX>& IX);
-  template <typename DerivedX, typename DerivedIX>
+  template <typename DerivedX, typename DerivedY, typename DerivedIX>
   // Only better if size(X,dim) is small
   IGL_INLINE void sort_new(
     const Eigen::PlainObjectBase<DerivedX>& X,
     const int dim,
     const bool ascending,
-    Eigen::PlainObjectBase<DerivedX>& Y,
+    Eigen::PlainObjectBase<DerivedY>& Y,
     Eigen::PlainObjectBase<DerivedIX>& IX);
   // Special case if size(X,dim) == 2
-  template <typename DerivedX, typename DerivedIX>
+  template <typename DerivedX, typename DerivedY, typename DerivedIX>
   IGL_INLINE void sort2(
     const Eigen::PlainObjectBase<DerivedX>& X,
     const int dim,
     const bool ascending,
-    Eigen::PlainObjectBase<DerivedX>& Y,
+    Eigen::PlainObjectBase<DerivedY>& Y,
     Eigen::PlainObjectBase<DerivedIX>& IX);
 
   // Act like matlab's [Y,I] = SORT(X) for std library vectors

include/igl/triangle_triangle_adjacency.cpp

@@ -6,8 +6,9 @@
 // v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "triangle_triangle_adjacency.h"
-
-#include <igl/is_edge_manifold.h>
+#include "is_edge_manifold.h"
+#include "all_edges.h"
+#include <map>
 #include <algorithm>
 
 template <typename Scalar, typename Index>
@@ -98,6 +99,63 @@ IGL_INLINE void igl::triangle_triangle_adjacency(const Eigen::PlainObjectBase<Sc
   triangle_triangle_adjacency_extractTTi(F,TTT,TTi);
 }
 
+template <
+  typename DerivedF, 
+  typename TTIndex, 
+  typename TTiIndex>
+  IGL_INLINE void igl::triangle_triangle_adjacency(
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    std::vector<std::vector<std::vector<TTIndex> > > & TT,
+    std::vector<std::vector<std::vector<TTiIndex> > > & TTi)
+{
+  using namespace Eigen;
+  using namespace std;
+  using namespace igl;
+  assert(F.cols() == 3 && "Faces must be triangles");
+  typedef typename DerivedF::Index Index;
+  // number of faces
+  const int m = F.rows();
+  // All occurances of directed edges
+  MatrixXi E;
+  all_edges(F,E);
+  assert(E.rows() == 3*m);
+  // uE2E[i] --> {j,k,...} means unique edge i corresponds to face edges j and
+  // k (where j-edge comes is the j/m edge of face j%m)
+  map<pair<Index,Index>,vector<Index> > uE2E;
+  for(int e = 0;e<E.rows();e++)
+  {
+    Index i = E(e,0);
+    Index j = E(e,1);
+    if(i<j)
+    {
+      uE2E[pair<Index,Index>(i,j)].push_back(e);
+    }else
+    {
+      uE2E[pair<Index,Index>(j,i)].push_back(e);
+    }
+  }
+  // E2E[i] --> {j,k,...} means face edge i corresponds to other faces edges j
+  // and k
+  TT.resize (m,vector<vector<TTIndex> >(F.cols()));
+  TTi.resize(m,vector<vector<TTiIndex> >(F.cols()));
+  for(int e = 0;e<E.rows();e++)
+  {
+    const Index i = E(e,0);
+    const Index j = E(e,1);
+    const Index f = e%m;
+    const Index c = e/m;
+    const vector<Index> & N = 
+      i<j ? uE2E[pair<Index,Index>(i,j)] : uE2E[pair<Index,Index>(j,i)];
+    for(const auto & ne : N)
+    {
+      const Index nf = ne%m;
+      const Index nc = ne/m;
+      TT[f][c].push_back(nf);
+      TTi[f][c].push_back(nc);
+    }
+  }
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 template void igl::triangle_triangle_adjacency<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);

include/igl/triangle_triangle_adjacency.h

@@ -29,6 +29,7 @@ namespace igl
   //   TTi  #F by #3 adjacent matrix, the element i,j is the id of edge of the triangle TT(i,j) that is adjacent with triangle i
   // NOTE: the first edge of a triangle is [0,1] the second [1,2] and the third [2,3].
   //       this convention is DIFFERENT from cotmatrix_entries.h
+  // Known bug: this should not need to take V as input.
 
   template <typename Scalar, typename Index>
   IGL_INLINE void triangle_triangle_adjacency(const Eigen::PlainObjectBase<Scalar>& V,
@@ -57,6 +58,24 @@ namespace igl
   IGL_INLINE void triangle_triangle_adjacency_extractTTi(const Eigen::PlainObjectBase<Index>& F,
                                 std::vector<std::vector<int> >& TTT,
                                 Eigen::PlainObjectBase<Index>& TTi);
+  // Adjacency list version, which works with non-manifold meshes
+  //
+  // Inputs:
+  //   F  #F by 3 list of triangle indices
+  // Outputs:
+  //   TT  #F by 3 list of lists so that TT[i][c] --> {j,k,...} means that faces j and
+  //     k etc. are edge-neighbors of face i on face i's edge opposite corner c
+  //   TTj  #F list of lists so that TTj[i][c] --> {j,k,...} means that face
+  //     TT[i][c][0] is an edge-neighbor of face i incident on the edge of face
+  //     TT[i][c][0] opposite corner j, and TT[i][c][1] " corner k, etc.
+  template <
+    typename DerivedF, 
+    typename TTIndex, 
+    typename TTiIndex>
+    IGL_INLINE void triangle_triangle_adjacency(
+      const Eigen::PlainObjectBase<DerivedF> & F,
+      std::vector<std::vector<std::vector<TTIndex> > > & TT,
+      std::vector<std::vector<std::vector<TTiIndex> > > & TTi);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/vertex_triangle_adjacency.cpp

@@ -7,24 +7,23 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "vertex_triangle_adjacency.h"
 
-#include "verbose.h"
-
-template <typename DerivedV, typename DerivedF, typename IndexType>
+template <typename DerivedF, typename VFType, typename VFiType>
 IGL_INLINE void igl::vertex_triangle_adjacency(
-                   const Eigen::PlainObjectBase<DerivedV>& V,
-                   const Eigen::PlainObjectBase<DerivedF>& F,
-                   std::vector<std::vector<IndexType> >& VF,
-                   std::vector<std::vector<IndexType> >& VFi)
+  const typename DerivedF::Scalar n,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  std::vector<std::vector<VFType> >& VF,
+  std::vector<std::vector<VFiType> >& VFi)
 {
   VF.clear();
   VFi.clear();
 
-  VF.resize(V.rows());
-  VFi.resize(V.rows());
+  VF.resize(n);
+  VFi.resize(n);
 
-  for(int fi=0; fi<F.rows(); ++fi)
+  typedef typename DerivedF::Index Index;
+  for(Index fi=0; fi<F.rows(); ++fi)
   {
-    for(int i = 0; i < F.cols(); ++i)
+    for(Index i = 0; i < F.cols(); ++i)
     {
       VF[F(fi,i)].push_back(fi);
       VFi[F(fi,i)].push_back(i);
@@ -32,6 +31,17 @@ IGL_INLINE void igl::vertex_triangle_adjacency(
   }
 }
 
+
+template <typename DerivedV, typename DerivedF, typename IndexType>
+IGL_INLINE void igl::vertex_triangle_adjacency(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  std::vector<std::vector<IndexType> >& VF,
+  std::vector<std::vector<IndexType> >& VFi)
+{
+  return vertex_triangle_adjacency(V.rows(),F,VF,VFi);
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 // generated by autoexplicit.sh

include/igl/vertex_triangle_adjacency.h

@@ -17,7 +17,8 @@ namespace igl
   // vertex_face_adjacency constructs the vertex-face topology of a given mesh (V,F)
   //
   // Inputs:
-  //   V  #V by 3 list of vertex coordinates
+  //   //V  #V by 3 list of vertex coordinates
+  //   n  number of vertices #V (e.g. `F.maxCoeff()+1` or `V.rows()`)
   //   F  #F by dim list of mesh faces (must be triangles)
   // Outputs:
   //   VF  #V list of lists of incident faces (adjacency list)
@@ -27,13 +28,21 @@ namespace igl
   // See also: edges, cotmatrix, diag, vv
   //
   // Known bugs: this should not take V as an input parameter.
-
+  // Known bugs/features: if a facet is combinatorially degenerate then faces
+  // will appear multiple times in VF and correpondingly in VFI (j appears
+  // twice in F.row(i) then i will appear twice in VF[j])
+  template <typename DerivedF, typename VFType, typename VFiType>
+  IGL_INLINE void vertex_triangle_adjacency(
+    const typename DerivedF::Scalar n,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    std::vector<std::vector<VFType> >& VF,
+    std::vector<std::vector<VFiType> >& VFi);
   template <typename DerivedV, typename DerivedF, typename IndexType>
   IGL_INLINE void vertex_triangle_adjacency(
-                     const Eigen::PlainObjectBase<DerivedV>& V,
-                     const Eigen::PlainObjectBase<DerivedF>& F,
-                     std::vector<std::vector<IndexType> >& VF,
-                     std::vector<std::vector<IndexType> >& VFi);
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    std::vector<std::vector<IndexType> >& VF,
+    std::vector<std::vector<IndexType> >& VFi);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/writeDMAT.cpp

@@ -94,6 +94,34 @@ IGL_INLINE bool igl::writeDMAT(
   return true;
 }
 
+template <typename Scalar>
+IGL_INLINE bool igl::writeDMAT(
+  const std::string file_name, 
+  const std::vector<Scalar > W)
+{
+  FILE * fp = fopen(file_name.c_str(),"w");
+  if(fp == NULL)
+  {
+    fprintf(stderr,"IOError: writeDMAT() could not open %s...",file_name.c_str());
+    return false; 
+  }
+  int num_rows = (int)W.size();
+  int num_cols = 0;
+  if(num_rows > 0)
+  {
+    num_cols = 1;
+  }
+  // first line contains number of columns and number of rows
+  fprintf(fp,"%d %d\n",num_cols,num_rows);
+  // loop over rows (down columns) quickly
+  for(int i = 0;i < num_rows;i++)
+  {
+    fprintf(fp,"%0.15lf\n",(double)W[i]);
+  }
+  fclose(fp);
+  return true;
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 // generated by autoexplicit.sh

include/igl/writeDMAT.h

@@ -32,6 +32,10 @@ namespace igl
   IGL_INLINE bool writeDMAT(
     const std::string file_name, 
     const std::vector<std::vector<Scalar> > W);
+  template <typename Scalar>
+  IGL_INLINE bool writeDMAT(
+    const std::string file_name, 
+    const std::vector<Scalar > W);
 }
 
 #ifndef IGL_STATIC_LIBRARY