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

Conflicts:
	include/igl/cut_mesh_from_singularities.h

Former-commit-id: 00538ec2c16dd56468ecd5e799084fd86fa6b392

README.md

@@ -406,6 +406,6 @@ BibTeX entry:
   title = {{libigl}: A simple {C++} geometry processing library},
   author = {Alec Jacobson and Daniele Panozzo and others},
   note = {http://igl.ethz.ch/projects/libigl/},
-  year = {2013},
+  year = {2014},
 }
 ```

VERSION.txt

@@ -3,4 +3,4 @@
 # Anyone may increment Minor to indicate a small change.
 # Major indicates a large change or large number of changes (upload to website)
 # World indicates a substantial change or release
-1.0.3
+1.0.4

build/Makefile.conf

@@ -113,6 +113,7 @@ endif
 ifeq ($(IGL_USERNAME),olkido)
     IGL_WITH_MATLAB=0
 		IGL_WITH_COMISO=1
+		IGL_WITH_XML=1
 		COMISO=/Users/olkido/Documents/igl/MIQ/src/CoMISo
 	  IGL_WITH_XML=1
     AFLAGS= -m64

build/Makefile_xml

@@ -1,7 +1,7 @@
 include Makefile.conf
 
 .PHONY: all
-all: 
+all:
 debug:
 #all: libiglxml
 #debug: libiglxml
@@ -20,7 +20,7 @@ OBJ_FILES=$(addprefix obj/,$(notdir $(CPP_FILES:.cpp=.o)))
 # include igl headers
 INC+=-I../include/
 
-# EXPECTS THAT CFLAGS IS ALREADY SET APPROPRIATELY 
+# EXPECTS THAT CFLAGS IS ALREADY SET APPROPRIATELY
 
 # Eigen dependency
 EIGEN3_INC=-I$(DEFAULT_PREFIX)/include/eigen3 -I$(DEFAULT_PREFIX)/include/eigen3/unsupported
@@ -42,7 +42,7 @@ INC+=$(ANTTWEAKBAR_INC)
 
 #CFLAGS+=-std=c++11
 
-obj: 
+obj:
 	mkdir -p obj
 
 ../lib/libiglxml.a: $(OBJ_FILES)

examples/multi-viewport/Makefile

@@ -13,7 +13,7 @@ LIBIGL_INC=-I$(LIBIGL)/include
 LIBIGL_LIB=-L$(LIBIGL)/lib -ligl -liglembree
 
 EIGEN3_INC=-I/opt/local/include/eigen3 -I/opt/local/include/eigen3/unsupported
-CFLAGS+=-std=c++11 -g
+CFLAGS+=-std=c++11 -g -Wno-deprecated-declarations
 
 ANTTWEAKBAR_INC=-I$(LIBIGL)/external/AntTweakBar/include
 ANTTWEAKBAR_LIB=-L$(LIBIGL)/external/AntTweakBar/lib -lAntTweakBar -framework AppKit

examples/multi-viewport/example.cpp

@@ -49,7 +49,6 @@ int width,height;
 igl::Camera down_camera;
 bool trackball_on = false;
 int down_mouse_x,down_mouse_y,move_x,move_y;
-int down_vp;
 // Position of light
 float light_pos[4] = {0.1,0.1,-0.9,0};
 // Vertex positions, normals, colors and centroid
@@ -222,7 +221,7 @@ void display()
   // All smooth points
   glEnable( GL_POINT_SMOOTH );
 
-  // Flash lights
+  // Flashlights
   lights();
   for(int vp = 0;vp<NUM_VIEWPORTS;vp++)
   {
@@ -321,7 +320,6 @@ void display()
   report_gl_error();
 
   glutSwapBuffers();
-  glutPostRedisplay();
 }
 
 // Initialize colors to a boring green
@@ -332,24 +330,6 @@ void init_C()
   C.col(2).setConstant(0.3);
 }
 
-int in_viewport(const int x, const int y)
-{
-  int down_vp = -1;
-  for(int vp = 0;vp<NUM_VIEWPORTS;vp++)
-  {
-    if(
-      x >= viewports[vp].x &&
-      y >= viewports[vp].y &&
-      x <  viewports[vp].x+viewports[vp].width &&
-      y <  viewports[vp].y+viewports[vp].height)
-    {
-      down_vp = vp;
-      break;
-    }
-  }
-  return down_vp;
-}
-
 const double SNAP_DIST = 10;
 void mouse_move(int mouse_x, int mouse_y)
 {
@@ -359,8 +339,10 @@ void mouse_move(int mouse_x, int mouse_y)
   using namespace std;
   move_x = mouse_x;
   move_y = mouse_y;
-  const int in_vp = in_viewport(mouse_x,height-mouse_y);
-  if(in_vp >= 0)
+  const int in_vp = find_if(viewports,viewports+NUM_VIEWPORTS,
+    [&mouse_x,&mouse_y](const AugViewport & vp) -> bool
+    {return vp.inside(mouse_x,height-mouse_y);})-viewports;
+  if(in_vp < NUM_VIEWPORTS)
   {
     if(
       viewports[in_vp].width > 0 &&
@@ -391,6 +373,7 @@ void mouse_move(int mouse_x, int mouse_y)
   {
     glutSetCursor(GLUT_CURSOR_LEFT_ARROW);
   }
+  glutPostRedisplay();
 }
 
 void mouse(int glutButton, int glutState, int mouse_x, int mouse_y)
@@ -422,9 +405,11 @@ void mouse(int glutButton, int glutState, int mouse_x, int mouse_y)
         vert_on = true;
       } else
       {
-        down_vp = in_viewport(mouse_x,height-mouse_y);
+        down_vp = find_if(viewports,viewports+NUM_VIEWPORTS,
+          [&mouse_x,&mouse_y](const AugViewport & vp) -> bool
+          {return vp.inside(mouse_x,height-mouse_y);})-viewports;
         // down
-        if(down_vp >= 0)
+        if(down_vp < NUM_VIEWPORTS)
         {
           glutSetCursor(GLUT_CURSOR_CYCLE);
           // collect information for trackball
@@ -436,6 +421,7 @@ void mouse(int glutButton, int glutState, int mouse_x, int mouse_y)
       }
     break;
   }
+  glutPostRedisplay();
 }
 
 void mouse_drag(int mouse_x, int mouse_y)
@@ -476,6 +462,7 @@ void mouse_drag(int mouse_x, int mouse_y)
       viewports[down_vp].mouse_y(mouse_y,height),
       viewports[down_vp].camera.m_rotation_conj.coeffs().data());
   }
+  glutPostRedisplay();
 }
 
 
@@ -490,8 +477,10 @@ void key(unsigned char key, int mouse_x, int mouse_y)
       exit(0);
     case 'Z':
     {
-      const int in_vp = in_viewport(mouse_x,height-mouse_y);
-      if(in_vp >= 0)
+      const int in_vp = find_if(viewports,viewports+NUM_VIEWPORTS,
+          [&mouse_x,&mouse_y](const AugViewport & vp) -> bool
+          {return vp.inside(mouse_x,height-mouse_y);})-viewports;
+      if(in_vp < NUM_VIEWPORTS)
       {
         igl::snap_to_canonical_view_quat(
           viewports[in_vp].camera.m_rotation_conj.coeffs().data(),
@@ -504,6 +493,7 @@ void key(unsigned char key, int mouse_x, int mouse_y)
       cout<<"Unknown key command: "<<key<<" "<<int(key)<<endl;
   }
 
+  glutPostRedisplay();
 }
 
 int main(int argc, char * argv[])

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/HalfEdgeIterator.h

@@ -12,6 +12,7 @@
 #include <Eigen/Core>
 
 #include <vector>
+#include <igl/igl_inline.h>
 
 namespace igl
 {

include/igl/RotateWidget.h

@@ -70,6 +70,8 @@ namespace igl
       inline bool is_down() const;
       inline void draw() const;
       inline void draw_guide() const;
+    public:
+        EIGEN_MAKE_ALIGNED_OPERATOR_NEW
   };
 }
 

include/igl/Viewport.h

@@ -55,6 +55,13 @@ namespace igl
     {
       return mx - x;
     }
+    // Returns whether point (mx,my) is in extend of Viewport
+    bool inside(const int mx, const int my) const
+    {
+      return 
+        mx >= x && my >= y && 
+        mx < x+width && my < y+height;
+    }
   };
 }
 

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

@@ -0,0 +1,132 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2014 Nico Pietroni <nico.pietroni@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 "comb_line_field.h"
+
+#include <vector>
+#include <deque>
+#include "per_face_normals.h"
+#include "is_border_vertex.h"
+#include "rotation_matrix_from_directions.h"
+
+#include "triangle_triangle_adjacency.h"
+
+namespace igl {
+template <typename DerivedV, typename DerivedF>
+class CombLine
+{
+public:
+
+    const Eigen::PlainObjectBase<DerivedV> &V;
+    const Eigen::PlainObjectBase<DerivedF> &F;
+    const Eigen::PlainObjectBase<DerivedV> &PD1;
+    Eigen::PlainObjectBase<DerivedV> N;
+
+private:
+    // internal
+    Eigen::PlainObjectBase<DerivedF> TT;
+    Eigen::PlainObjectBase<DerivedF> TTi;
+
+
+private:
+
+
+    static inline double Sign(double a){return (double)((a>0)?+1:-1);}
+
+
+private:
+
+    // returns the 180 deg rotation of a (around n) most similar to target b
+    // a and b should be in the same plane orthogonal to N
+    static inline Eigen::Matrix<typename DerivedV::Scalar, 3, 1> K_PI_line(const Eigen::Matrix<typename DerivedV::Scalar, 3, 1>& a,
+                                                                           const Eigen::Matrix<typename DerivedV::Scalar, 3, 1>& b)
+    {
+        typename DerivedV::Scalar scorea = a.dot(b);
+        if (scorea<0)
+            return -a;
+        else
+            return a;
+    }
+
+
+
+public:
+
+    inline CombLine(const Eigen::PlainObjectBase<DerivedV> &_V,
+                    const Eigen::PlainObjectBase<DerivedF> &_F,
+                    const Eigen::PlainObjectBase<DerivedV> &_PD1):
+        V(_V),
+        F(_F),
+        PD1(_PD1)
+    {
+        igl::per_face_normals(V,F,N);
+        igl::triangle_triangle_adjacency(V,F,TT,TTi);
+    }
+
+    inline void comb(Eigen::PlainObjectBase<DerivedV> &PD1out)
+    {
+        PD1out.setZero(F.rows(),3);PD1out<<PD1;
+
+        Eigen::VectorXi mark = Eigen::VectorXi::Constant(F.rows(),false);
+
+        std::deque<int> d;
+
+        d.push_back(0);
+        mark(0) = true;
+
+        while (!d.empty())
+        {
+            int f0 = d.at(0);
+            d.pop_front();
+            for (int k=0; k<3; k++)
+            {
+                int f1 = TT(f0,k);
+                if (f1==-1) continue;
+                if (mark(f1)) continue;
+
+                Eigen::Matrix<typename DerivedV::Scalar, 3, 1> dir0  = PD1out.row(f0);
+                Eigen::Matrix<typename DerivedV::Scalar, 3, 1> dir1  = PD1out.row(f1);
+                Eigen::Matrix<typename DerivedV::Scalar, 3, 1> n0    = N.row(f0);
+                Eigen::Matrix<typename DerivedV::Scalar, 3, 1> n1    = N.row(f1);
+
+                Eigen::Matrix<typename DerivedV::Scalar, 3, 1> dir0Rot = igl::rotation_matrix_from_directions(n0, n1)*dir0;
+                dir0Rot.normalize();
+                Eigen::Matrix<typename DerivedV::Scalar, 3, 1> targD   = K_PI_line(dir1,dir0Rot);
+
+                PD1out.row(f1)  = targD;
+                //PD2out.row(f1)  = n1.cross(targD).normalized();
+
+                mark(f1) = true;
+                d.push_back(f1);
+
+            }
+        }
+
+        // everything should be marked
+        for (int i=0; i<F.rows(); i++)
+        {
+            assert(mark(i));
+        }
+    }
+
+};
+}
+
+template <typename DerivedV, typename DerivedF>
+IGL_INLINE void igl::comb_line_field(const Eigen::PlainObjectBase<DerivedV> &V,
+                                     const Eigen::PlainObjectBase<DerivedF> &F,
+                                     const Eigen::PlainObjectBase<DerivedV> &PD1,
+                                     Eigen::PlainObjectBase<DerivedV> &PD1out)
+{
+    igl::CombLine<DerivedV, DerivedF> cmb(V, F, PD1);
+    cmb.comb(PD1out);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+#endif

include/igl/comb_line_field.h

@@ -0,0 +1,39 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2014 Nico Pietroni <nico.pietroni@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_COMB_LINE_FIELD_H
+#define IGL_COMB_LINE_FIELD_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Computes principal matchings of the vectors of a cross field across face edges,
+  // and generates a combed cross field defined on the mesh faces
+
+  // Inputs:
+  //   V          #V by 3 eigen Matrix of mesh vertex 3D positions
+  //   F          #F by 4 eigen Matrix of face (quad) indices
+  //   PD1in      #F by 3 eigen Matrix of the first per face cross field vector
+  //   PD2in      #F by 3 eigen Matrix of the second per face cross field vector
+  // Output:
+  //   PD1out      #F by 3 eigen Matrix of the first combed cross field vector
+  //   PD2out      #F by 3 eigen Matrix of the second combed cross field vector
+  //
+
+
+  template <typename DerivedV, typename DerivedF>
+  IGL_INLINE void comb_line_field(const Eigen::PlainObjectBase<DerivedV> &V,
+                                  const Eigen::PlainObjectBase<DerivedF> &F,
+                                  const Eigen::PlainObjectBase<DerivedV> &PD1in,
+                                  Eigen::PlainObjectBase<DerivedV> &PD1out);
+}
+#ifndef IGL_STATIC_LIBRARY
+#include "comb_line_field.cpp"
+#endif
+
+#endif

include/igl/comiso/miq.cpp

@@ -975,7 +975,7 @@ IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::SolvePoisson(Eigen::Vect
 
   clearUserConstraint();
   // copy the user constraints number
-  for (int i = 0; i < hardFeatures.size(); ++i)
+  for (size_t i = 0; i < hardFeatures.size(); ++i)
   {
     addSharpEdgeConstraint(hardFeatures[i][0],hardFeatures[i][1]);
   }
@@ -1446,7 +1446,7 @@ IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::AddSingularityRound()
 template <typename DerivedV, typename DerivedF>
 IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::AddToRoundVertices(std::vector<int> ids)
 {
-  for (int i = 0; i < ids.size(); ++i)
+  for (size_t i = 0; i < ids.size(); ++i)
   {
     if (ids[i] < 0 || ids[i] >= V.rows())
       std::cerr << "WARNING: Ignored round vertex constraint, vertex " << ids[i] << " does not exist in the mesh." << std::endl;
@@ -1790,7 +1790,7 @@ IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::addSharpEdgeConstraint(i
   // prepare constraint
   std::vector<int> c(Handle_SystemInfo.num_vert_variables*2 + 1);
 
-  for (int i = 0; i < c.size(); ++i)
+  for (size_t i = 0; i < c.size(); ++i)
   {
     c[i] = 0;
   }
@@ -2272,7 +2272,9 @@ IGL_INLINE void igl::miq(const Eigen::PlainObjectBase<DerivedV> &V,
            DirectRound,
            iter,
            localIter,
-           DoRound);
+           DoRound,
+           roundVertices,
+           hardFeatures);
 
 }
 

include/igl/cut_mesh_from_singularities.cpp

@@ -9,6 +9,7 @@
 #include "cut_mesh_from_singularities.h"
 
 #include <igl/triangle_triangle_adjacency.h>
+#include <igl/edge_topology.h>
 
 #include <vector>
 #include <deque>
@@ -32,6 +33,7 @@ namespace igl {
     Eigen::PlainObjectBase<DerivedF> TT;
     Eigen::PlainObjectBase<DerivedF> TTi;
 
+    Eigen::MatrixXi E, F2E, E2F;
   protected:
 
     inline bool IsRotSeam(const int f0,const int edge)
@@ -75,17 +77,32 @@ namespace igl {
     inline void Retract(Eigen::PlainObjectBase<DerivedO> &Handle_Seams)
     {
       std::vector<int> e(V.rows(),0); // number of edges per vert
-
-      for (unsigned f=0; f<F.rows(); f++)
+      // for (unsigned f=0; f<F.rows(); f++)
+      // {
+      //   for (int s = 0; s<3; s++)
+      //   {
+      //     if (Handle_Seams(f,s))
+      //       if (TT(f,s)<=f)
+      //       {
+      //         e[ F(f,s) ] ++;
+      //         e[ F(f,(s+1)%3) ] ++;
+      //       }
+      //   }
+      // }
+      for (int ei=0; ei<E.rows(); ++ei)
       {
-        for (int s = 0; s<3; s++)
+        //only need one face
+        int f0 = E2F(ei,0);
+        if (f0==-1)
+          f0 = E2F(ei,1);
+        int k=0;
+        for (k=0; k<3; ++k)
+          if (F2E(f0,k)==ei)
+            break;
+        if (Handle_Seams(f0,k))
         {
-          if (Handle_Seams(f,s))
-            if (TT(f,s)<=f)
-            {
-              e[ F(f,s) ] ++;
-              e[ F(f,(s+1)%3) ] ++;
-            }
+          e[ F(f0,k) ] ++;
+          e[ F(f0,(k+1)%3) ] ++;
         }
       }
 
@@ -115,7 +132,7 @@ namespace igl {
           }
         }
 
-        if (guard++>10000)
+        if (guard++>0)
           over = true;
 
       } while (!over);
@@ -131,6 +148,7 @@ namespace igl {
     Handle_MMatch(Handle_MMatch_)
     {
       triangle_triangle_adjacency(V,F,TT,TTi);
+      edge_topology(V,F,E,F2E,E2F);
     };
 
     inline void cut(Eigen::PlainObjectBase<DerivedO> &Handle_Seams)

include/igl/cut_mesh_from_singularities.h

@@ -22,8 +22,6 @@ namespace igl
   // MMatch:            #F by 3 list of per corner integer mismatch
   // seams:             #F by 3 list of per corner booleans that denotes if an edge is a seam or not
   //
-  // TODO: make the name of the variables consistent in the cpp
-
   template <typename DerivedV, typename DerivedF, typename DerivedM, typename DerivedO>
   IGL_INLINE void cut_mesh_from_singularities(const Eigen::PlainObjectBase<DerivedV> &V,
                                                    const Eigen::PlainObjectBase<DerivedF> &F,

include/igl/draw_point.cpp

@@ -30,7 +30,9 @@ IGL_INLINE void igl::draw_point(
 
   float f;
   glGetFloatv(GL_POINT_SIZE_MAX,&f);
-  assert(requested_r<=0.5*f);
+  // THIS IS OVERZEALOUS on Mac OS X: OpenGL reports a smaller point size than
+  // possible.
+  //assert(requested_r<=0.5*f);
   double r = (requested_r<0.5*f?requested_r:0.5*f);
 
   //glDisable(GL_DEPTH_TEST);
@@ -87,11 +89,18 @@ IGL_INLINE void igl::draw_point(
   const double requested_r,
   const bool selected)
 {
-  return draw_point(P(0),P(1),P(2),requested_r,selected);
+  switch(P.size())
+  {
+    case 2:
+      return draw_point(P(0),P(1),0,requested_r,selected);
+    default:
+      return draw_point(P(0),P(1),P(2),requested_r,selected);
+  }
 }
 
 #ifdef IGL_STATIC_LIBRARY
 template void igl::draw_point<Eigen::Matrix<double, 3, 1, 0, 3, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 1, 0, 3, 1> > const&, double, bool);
+template void igl::draw_point<Eigen::Matrix<double, 2, 1, 0, 2, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 2, 1, 0, 2, 1> > const&, double, bool); 
 #endif
 
 #endif

include/igl/edge_topology.cpp

@@ -18,6 +18,13 @@ IGL_INLINE void igl::edge_topology(
                                    Eigen::MatrixXi& EF)
 {
   // Only needs to be edge-manifold
+  if (V.rows() ==0 || F.rows()==0)
+  {
+    EV = Eigen::MatrixXi::Constant(0,2,-1);
+    FE = Eigen::MatrixXi::Constant(0,3,-1);
+    EF = Eigen::MatrixXi::Constant(0,2,-1);
+    return;
+  }
   assert(igl::is_edge_manifold(V,F));
   std::vector<std::vector<int> > ETT;
   for(int f=0;f<F.rows();++f)
@@ -36,7 +43,7 @@ IGL_INLINE void igl::edge_topology(
 
   // count the number of edges (assume manifoldness)
   int En = 1; // the last is always counted
-  for(unsigned i=0;i<ETT.size()-1;++i)
+  for(int i=0;i<int(ETT.size())-1;++i)
     if (!((ETT[i][0] == ETT[i+1][0]) && (ETT[i][1] == ETT[i+1][1])))
       ++En;
 
@@ -75,7 +82,6 @@ IGL_INLINE void igl::edge_topology(
 
   // Sort the relation EF, accordingly to EV
   // the first one is the face on the left of the edge
-
   for(unsigned i=0; i<EF.rows(); ++i)
   {
     int fid = EF(i,0);

include/igl/file_exists.cpp

@@ -9,8 +9,8 @@
 
 #include <sys/stat.h>
 
-IGL_INLINE bool igl::file_exists(const char* filename)
+IGL_INLINE bool igl::file_exists(const std::string filename)
 {
   struct stat status;
-  return (stat(filename,&status)==0);
+  return (stat(filename.c_str(),&status)==0);
 }

include/igl/file_exists.h

@@ -8,6 +8,7 @@
 #ifndef IGL_FILE_EXISTS_H
 #define IGL_FILE_EXISTS_H
 #include "igl_inline.h"
+#include <string>
 namespace igl
 {
   // Check if a file or directory exists like PHP's file_exists function:
@@ -16,7 +17,7 @@ namespace igl
   //   filename  path to file
   // Returns true if file exists and is readable and false if file doesn't
   // exist or *is not readable*
-  IGL_INLINE bool file_exists(const char * filename);
+  IGL_INLINE bool file_exists(const std::string filename);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/find_cross_field_singularities.cpp

@@ -36,7 +36,7 @@ IGL_INLINE void igl::find_cross_field_singularities(const Eigen::PlainObjectBase
   {
     ///check that is on border..
     if (V_border[vid])
-      return;
+      continue;
 
     int missmatch=0;
     for (unsigned int i=0;i<VF[vid].size();i++)

include/igl/fit_rigid.cpp

@@ -0,0 +1,25 @@
+// 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/.
+#include "fit_rigid.h"
+#include "polar_svd.h"
+
+IGL_INLINE void igl::fit_rigid(
+  const Eigen::MatrixXd & A,
+  const Eigen::MatrixXd & B,
+  Eigen::Rotation2Dd & R,
+  Eigen::RowVector2d & t)
+{
+  // build covariance matrix
+  const auto & Amean = A.colwise().mean();
+  const auto & Bmean = B.colwise().mean();
+  const auto & S = (B.rowwise() - Bmean).transpose() * (A.rowwise() - Amean);
+  Eigen::Matrix2d Rm,Tm;
+  polar_svd(S.eval(),Rm,Tm);
+  R.fromRotationMatrix(Rm);
+  t = Amean - Bmean*Rm;
+}

include/igl/fit_rigid.h

@@ -0,0 +1,30 @@
+// 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_FIT_RIGID_H
+#define IGL_FIT_RIGID_H
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <Eigen/Geometry>
+
+namespace igl
+{
+  // Fit a rigid 
+  IGL_INLINE void fit_rigid(
+    const Eigen::MatrixXd & A,
+    const Eigen::MatrixXd & B,
+    Eigen::Rotation2Dd & R,
+    Eigen::RowVector2d & t);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "fit_rigid.cpp"
+#endif
+
+#endif
+

include/igl/get_modifiers.cpp

@@ -0,0 +1,42 @@
+#include "get_modifiers.h"
+
+/* glutGetModifiers return mask. */
+#ifndef GLUT_ACTIVE_SHIFT
+#  define GLUT_ACTIVE_SHIFT 1
+#endif
+#ifndef GLUT_ACTIVE_CTRL
+#  define GLUT_ACTIVE_CTRL 2
+#endif
+#ifndef GLUT_ACTIVE_ALT
+#  define GLUT_ACTIVE_ALT 4
+#endif
+#ifndef GLUT_ACTIVE_COMMAND
+#  define GLUT_ACTIVE_COMMAND 8
+#endif
+
+#ifdef __APPLE__
+#include <Carbon/Carbon.h>
+#endif
+
+IGL_INLINE int igl::get_modifiers()
+{
+  int mod = 0;
+#ifdef __APPLE__
+  // http://stackoverflow.com/a/18082326/148668
+  KeyMap keyStates;
+  const auto & carbon_is_keydown = [&keyStates]( uint16_t vKey )->bool
+  {
+    uint8_t index = vKey / 32 ;
+    uint8_t shift = vKey % 32 ;
+    return keyStates[index].bigEndianValue & (1 << shift) ;
+  };
+  GetKeys(keyStates) ;
+  mod |= (carbon_is_keydown(kVK_Command)?GLUT_ACTIVE_COMMAND:0);
+  mod |= (carbon_is_keydown(kVK_Shift)?GLUT_ACTIVE_SHIFT:0);
+  mod |= (carbon_is_keydown(kVK_Option)?GLUT_ACTIVE_ALT:0);
+  mod |= (carbon_is_keydown(kVK_Control)?GLUT_ACTIVE_CTRL:0);
+#else
+#  error "Not supported.
+#endif
+  return mod;
+}

include/igl/get_modifiers.h

@@ -0,0 +1,22 @@
+#ifndef GET_MODIFIERS_H
+#define GET_MODIFIERS_H
+#include "igl_inline.h"
+namespace igl
+{
+  enum Modifier
+  {
+    MODIFIER_OPTION = 1,
+    MODIFIER_SHIFT = 3,
+    MODIFIER_CONTROL = 5,
+    MODIFIER_COMMAND = 9,
+    NUM_MODIFIERS = 4,
+  };
+  // Retrieve current modifier constellation. 
+  //
+  // Returns int that's an "or" of the active modifiers above.
+  IGL_INLINE int get_modifiers();
+}
+#ifndef IGL_STATIC_LIBRARY
+#include "get_modifiers.cpp"
+#endif
+#endif

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

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "jet.h"
 #include "colon.h"
@@ -81,7 +81,7 @@ IGL_INLINE void igl::jet(const T x_in, T & r, T & g, T & b)
   {
     r = 0;
     g = gone*(x-1./8.)/(3./8.-1./8.);
-    b = bone; 
+    b = bone;
   }else if(x<5./8.)
   {
     r = rone*(x-3./8.)/(5./8.-3./8.);
@@ -134,4 +134,5 @@ template void igl::jet<double>(double, double&, double&, double&);
 template void igl::jet<float>(float, float*);
 template void igl::jet<Eigen::Array<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Array<double, -1, 1, 0, -1, 1> > const&, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::jet<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::jet<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/line_field_missmatch.cpp

@@ -0,0 +1,143 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2014 Nico Pietroni <nico.pietroni@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 "line_field_missmatch.h"
+
+#include <vector>
+#include <deque>
+#include <igl/comb_line_field.h>
+#include <igl/rotate_vectors.h>
+#include <igl/comb_cross_field.h>
+#include <igl/comb_line_field.h>
+#include <igl/per_face_normals.h>
+#include <igl/is_border_vertex.h>
+#include <igl/vertex_triangle_adjacency.h>
+#include <igl/triangle_triangle_adjacency.h>
+#include <igl/rotation_matrix_from_directions.h>
+#include <igl/local_basis.h>
+
+namespace igl {
+template <typename DerivedV, typename DerivedF, typename DerivedO>
+class MissMatchCalculatorLine
+{
+public:
+
+    const Eigen::PlainObjectBase<DerivedV> &V;
+    const Eigen::PlainObjectBase<DerivedF> &F;
+    const Eigen::PlainObjectBase<DerivedV> &PD1;
+    const Eigen::PlainObjectBase<DerivedV> &PD2;
+    Eigen::PlainObjectBase<DerivedV> N;
+
+private:
+    // internal
+    std::vector<bool> V_border; // bool
+    std::vector<std::vector<int> > VF;
+    std::vector<std::vector<int> > VFi;
+    Eigen::PlainObjectBase<DerivedF> TT;
+    Eigen::PlainObjectBase<DerivedF> TTi;
+
+
+private:
+
+    //compute the mismatch between 2 faces
+    inline int MissMatchByLine(const int f0,
+                               const int f1)
+    {
+        Eigen::Matrix<typename DerivedV::Scalar, 3, 1> dir0 = PD1.row(f0);
+        Eigen::Matrix<typename DerivedV::Scalar, 3, 1> dir1 = PD1.row(f1);
+        Eigen::Matrix<typename DerivedV::Scalar, 3, 1> n0 = N.row(f0);
+        Eigen::Matrix<typename DerivedV::Scalar, 3, 1> n1 = N.row(f1);
+
+        Eigen::Matrix<typename DerivedV::Scalar, 3, 1> dir1Rot = igl::rotation_matrix_from_directions(n1,n0)*dir1;
+        dir1Rot.normalize();
+
+        // TODO: this should be equivalent to the other code below, to check!
+        // Compute the angle between the two vectors
+        //    double a0 = atan2(dir0.dot(B2.row(f0)),dir0.dot(B1.row(f0)));
+        //    double a1 = atan2(dir1Rot.dot(B2.row(f0)),dir1Rot.dot(B1.row(f0)));
+        //
+        //    double angle_diff = a1-a0;   //VectToAngle(f0,dir1Rot);
+
+        double angle_diff = atan2(dir1Rot.dot(PD2.row(f0)),dir1Rot.dot(PD1.row(f0)));
+
+        double step=M_PI;
+        int i=(int)floor((angle_diff/step)+0.5);
+        assert((i>=-2)&&(i<=2));
+        int k=0;
+        if (i>=0)
+            k=i%2;
+        else
+            k=(2+i)%2;
+
+        assert((k==0)||(k==1));
+        return (k*2);
+    }
+
+public:
+
+    inline MissMatchCalculatorLine(const Eigen::PlainObjectBase<DerivedV> &_V,
+                               const Eigen::PlainObjectBase<DerivedF> &_F,
+                               const Eigen::PlainObjectBase<DerivedV> &_PD1,
+                               const Eigen::PlainObjectBase<DerivedV> &_PD2
+                               ):
+        V(_V),
+        F(_F),
+        PD1(_PD1),
+        PD2(_PD2)
+    {
+        igl::per_face_normals(V,F,N);
+        V_border = igl::is_border_vertex(V,F);
+        igl::vertex_triangle_adjacency(V,F,VF,VFi);
+        igl::triangle_triangle_adjacency(V,F,TT,TTi);
+    }
+
+    inline void calculateMissmatchLine(Eigen::PlainObjectBase<DerivedO> &Handle_MMatch)
+    {
+        Handle_MMatch.setConstant(F.rows(),3,-1);
+        for (unsigned int i=0;i<F.rows();i++)
+        {
+            for (int j=0;j<3;j++)
+            {
+                if (i==TT(i,j) || TT(i,j) == -1)
+                    Handle_MMatch(i,j)=0;
+                else
+                    Handle_MMatch(i,j) = MissMatchByLine(i,TT(i,j));
+            }
+        }
+    }
+
+};
+}
+
+
+template <typename DerivedV, typename DerivedF, typename DerivedO>
+IGL_INLINE void igl::line_field_missmatch(const Eigen::PlainObjectBase<DerivedV> &V,
+                                const Eigen::PlainObjectBase<DerivedF> &F,
+                                const Eigen::PlainObjectBase<DerivedV> &PD1,
+                                const bool isCombed,
+                                Eigen::PlainObjectBase<DerivedO> &missmatch)
+{
+    Eigen::PlainObjectBase<DerivedV> PD1_combed;
+    Eigen::PlainObjectBase<DerivedV> PD2_combed;
+
+    if (!isCombed)
+        igl::comb_line_field(V,F,PD1,PD1_combed);
+    else
+    {
+        PD1_combed = PD1;
+    }
+    Eigen::MatrixXd B1,B2,B3;
+    igl::local_basis(V,F,B1,B2,B3);
+    PD2_combed = igl::rotate_vectors(PD1_combed, Eigen::VectorXd::Constant(1,M_PI/2), B1, B2);
+    igl::MissMatchCalculatorLine<DerivedV, DerivedF, DerivedO> sf(V, F, PD1_combed, PD2_combed);
+    sf.calculateMissmatchLine(missmatch);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+#endif

include/igl/line_field_missmatch.h

@@ -0,0 +1,42 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2014 Nico Pietroni <nico.pietroni@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_LINE_FIELD_MISSMATCH_H
+#define IGL_LINE_FIELD_MISSMATCH_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Calculates the missmatch (integer), at each face edge, of a cross field defined on the mesh faces.
+  // The integer missmatch is a multiple of pi/2 that transforms the cross on one side of the edge to
+  // the cross on the other side. It represents the deviation from a Lie connection across the edge.
+
+  // Inputs:
+  //   V         #V by 3 eigen Matrix of mesh vertex 3D positions
+  //   F         #F by 3 eigen Matrix of face (quad) indices
+  //   PD1       #F by 3 eigen Matrix of the first per face cross field vector
+  //   PD2       #F by 3 eigen Matrix of the second per face cross field vector
+  //   isCombed  boolean, specifying whether the field is combed (i.e. matching has been precomputed.
+  //             If not, the field is combed first.
+  // Output:
+  //   Handle_MMatch    #F by 3 eigen Matrix containing the integer missmatch of the cross field
+  //                    across all face edges
+  //
+
+    template <typename DerivedV, typename DerivedF, typename DerivedO>
+    IGL_INLINE void line_field_missmatch(const Eigen::PlainObjectBase<DerivedV> &V,
+                                         const Eigen::PlainObjectBase<DerivedF> &F,
+                                         const Eigen::PlainObjectBase<DerivedV> &PD1,
+                                         const bool isCombed,
+                                         Eigen::PlainObjectBase<DerivedO> &missmatch);
+}
+#ifndef IGL_STATIC_LIBRARY
+#include "line_field_missmatch.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_face_normals.cpp

@@ -52,4 +52,5 @@ IGL_INLINE void igl::per_face_normals(
 template void igl::per_face_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_face_normals<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, 3, 1, 0, 3, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 1, 0, 3, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::per_face_normals<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::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> >&);
+template void igl::per_face_normals<Eigen::Matrix<float, -1, -1, 1, -1, -1>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<float, -1, -1, 1, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 1, -1, -1> >&);
 #endif

include/igl/per_vertex_normals.cpp

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "per_vertex_normals.h"
 
@@ -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);
@@ -98,11 +98,7 @@ 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, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, igl::PerVertexNormalsWeightingType, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -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> >&);
-template void igl::per_vertex_normals<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, igl::PerVertexNormalsWeightingType, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+// generated by autoexplicit.sh
 #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/png/render_to_png_async.cpp

@@ -39,7 +39,9 @@ static IGL_INLINE bool render_to_png_async_helper(
     }
   }
 
-  return img->save(png_file.c_str(),fast);
+  bool ret = img->save(png_file.c_str(),fast);
+  delete img;
+  return ret;
 }
 
 IGL_INLINE std::thread igl::render_to_png_async(
@@ -61,5 +63,7 @@ IGL_INLINE std::thread igl::render_to_png_async(
     GL_UNSIGNED_BYTE,
     img->data());
   // Part that should be asynchronous  
-  return std::thread(render_to_png_async_helper,img,png_file,alpha,fast);
+  std::thread t(render_to_png_async_helper,img,png_file,alpha,fast);
+  t.detach();
+  return t;
 }

include/igl/polar_svd.cpp

@@ -69,4 +69,5 @@ template void igl::polar_svd<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Ma
 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::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> >&);
 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::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> >&);
 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::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> >&);
+template void igl::polar_svd<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, 2, 2, 0, 2, 2>, Eigen::Matrix<double, 2, 2, 0, 2, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 2, 2, 0, 2, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, 2, 2, 0, 2, 2> >&);
 #endif

include/igl/project.cpp

@@ -114,7 +114,7 @@ IGL_INLINE Eigen::PlainObjectBase<Derivedobj> igl::project(
 #endif
 #endif
 
-Eigen::Vector3f igl::project(const Eigen::Vector3f&  obj,
+IGL_INLINE Eigen::Vector3f igl::project(const Eigen::Vector3f&  obj,
                         const Eigen::Matrix4f& model,
                         const Eigen::Matrix4f& proj,
                         const Eigen::Vector4f&  viewport)
@@ -147,6 +147,8 @@ template Eigen::PlainObjectBase<Eigen::Matrix<float, 3, 1, 0, 3, 1> > igl::proje
 template Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> > igl::project<Eigen::Matrix<double, 1, -1, 1, 1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> > const&);
 template int igl::project<Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&);
 template Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > igl::project<Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&);
+template Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> > igl::project<Eigen::Matrix<double, 1, 2, 1, 1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> > const&);
+template Eigen::PlainObjectBase<Eigen::Matrix<double, 2, 1, 0, 2, 1> > igl::project<Eigen::Matrix<double, 2, 1, 0, 2, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 2, 1, 0, 2, 1> > const&);
 #endif
 #endif
 

include/igl/project.h

@@ -44,7 +44,7 @@ namespace igl
 // Returns:
 //   screen space x, y, and z coordinates respectively
 // Returns return value of gluProject call
-  Eigen::Vector3f project(const Eigen::Vector3f&  obj,
+  IGL_INLINE Eigen::Vector3f project(const Eigen::Vector3f&  obj,
                           const Eigen::Matrix4f& model,
                           const Eigen::Matrix4f& proj,
                           const Eigen::Vector4f&  viewport);

include/igl/project_to_line.cpp

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

include/igl/project_to_line_segment.cpp

@@ -44,4 +44,5 @@ IGL_INLINE void igl::project_to_line_segment(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 template void igl::project_to_line_segment<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, 1, 0, 1, 1>, Eigen::Matrix<double, 1, 1, 0, 1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 1, 0, 1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 1, 0, 1, 1> >&);
+template void igl::project_to_line_segment<Eigen::Matrix<double, 1, 2, 1, 1, 2>, Eigen::Matrix<double, 1, 2, 1, 1, 2>, Eigen::Matrix<double, 1, 2, 1, 1, 2>, Eigen::Matrix<double, 1, 1, 0, 1, 1>, Eigen::Matrix<double, 1, 1, 0, 1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 1, 0, 1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 1, 0, 1, 1> >&);
 #endif

include/igl/readDMAT.cpp

@@ -205,4 +205,5 @@ template bool igl::readDMAT<Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::basic_s
 template bool igl::readDMAT<Eigen::Matrix<double, 4, 1, 0, 4, 1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<double, 4, 1, 0, 4, 1> >&);
 template bool igl::readDMAT<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> >&);
 template bool igl::readDMAT<Eigen::Matrix<int, -1, 1, 0, -1, 1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template bool igl::readDMAT<Eigen::Matrix<int, -1, 2, 0, -1, 2> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&);
 #endif

include/igl/rotation_matrix_from_directions.cpp

@@ -10,18 +10,18 @@
 #include <Eigen/Geometry>
 
 template <typename Scalar>
-IGL_INLINE Eigen::Matrix<Scalar, 3, 3> igl::rotation_matrix_from_directions(Eigen::Matrix<Scalar, 3, 1> v0,
-                                                                        Eigen::Matrix<Scalar, 3, 1> v1,
+IGL_INLINE Eigen::Matrix<Scalar, 3, 3> igl::rotation_matrix_from_directions(const Eigen::Matrix<Scalar, 3, 1> v0,
+                                                                        const Eigen::Matrix<Scalar, 3, 1> v1,
                                                                         bool normalized)
 {
   Eigen::Matrix<Scalar, 3, 3> rotM;
   const double epsilon=0.00001;
-  if (!normalized)
-  {
-    v0.normalize();
-    v1.normalize();
-  }
-  Scalar dot=v0.dot(v1);
+  // if (!normalized)
+  // {
+    // v0.normalize();
+    // v1.normalize();
+  // }
+  Scalar dot=v0.normalized().dot(v1.normalized());
   ///control if there is no rotation
   if (dot>((double)1-epsilon))
   {
@@ -57,5 +57,5 @@ IGL_INLINE Eigen::Matrix<Scalar, 3, 3> igl::rotation_matrix_from_directions(Eige
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template Eigen::Matrix<double, 3, 3, 0, 3, 3> igl::rotation_matrix_from_directions<double>(Eigen::Matrix<double, 3, 1, 0, 3, 1>, Eigen::Matrix<double, 3, 1, 0, 3, 1>, bool);
+template Eigen::Matrix<double, 3, 3, 0, 3, 3> igl::rotation_matrix_from_directions<double>(const Eigen::Matrix<double, 3, 1, 0, 3, 1>, const Eigen::Matrix<double, 3, 1, 0, 3, 1>, const bool);
 #endif

include/igl/rotation_matrix_from_directions.h

@@ -14,7 +14,7 @@
 
 namespace igl {
   /// Given 2 vectors centered on origin calculate the rotation matrix from first to the second
-  
+
   // Inputs:
   //   v0, v1         the two #3 by 1 vectors
   //   normalized     boolean, if false, then the vectors are normalized prior to the calculation
@@ -22,9 +22,9 @@ namespace igl {
   //                  3 by 3 rotation matrix that takes v0 to v1
   //
   template <typename Scalar>
-  IGL_INLINE Eigen::Matrix<Scalar, 3, 3> rotation_matrix_from_directions(Eigen::Matrix<Scalar, 3, 1> v0,
-                                                                     Eigen::Matrix<Scalar, 3, 1> v1,
-                                                                     bool normalized=true);
+  IGL_INLINE Eigen::Matrix<Scalar, 3, 3> rotation_matrix_from_directions(const Eigen::Matrix<Scalar, 3, 1> v0,
+                                                                     const Eigen::Matrix<Scalar, 3, 1> v1,
+                                                                     const bool normalized=true);
 }
 
 

include/igl/sort.cpp

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "sort.h"
 
@@ -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,11 @@ 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> >&);
+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> >&);
 #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<Index> >(F.cols()));
+  TTi.resize(m,vector<vector<Index> >(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/triangles_from_strip.cpp

@@ -0,0 +1,29 @@
+#include "triangles_from_strip.h"
+#include <iostream>
+
+template <typename DerivedS, typename DerivedF>
+IGL_INLINE void igl::triangles_from_strip(
+  const Eigen::MatrixBase<DerivedS>& S,
+  Eigen::PlainObjectBase<DerivedF>& F)
+{
+  using namespace std;
+  F.resize(S.size()-2,3);
+  for(int s = 0;s < S.size()-2;s++)
+  {
+    if(s%2 == 0)
+    {
+      F(s,0) = S(s+2);
+      F(s,1) = S(s+1);
+      F(s,2) = S(s+0);
+    }else
+    {
+      F(s,0) = S(s+0);
+      F(s,1) = S(s+1);
+      F(s,2) = S(s+2);
+    }
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+#endif

include/igl/triangles_from_strip.h

@@ -0,0 +1,33 @@
+// 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_TRIANGLES_FROM_STRIP_H
+#define IGL_TRIANGLES_FROM_STRIP_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // TRIANGLES_FROM_STRIP Create a list of triangles from a stream of indices
+  // along a strip.
+  //
+  // Inputs:
+  //   S  #S list of indices
+  // Outputs:
+  //   F  #S-2 by 3 list of triangle indices
+  //
+  template <typename DerivedS, typename DerivedF>
+  IGL_INLINE void triangles_from_strip(
+    const Eigen::MatrixBase<DerivedS>& S,
+    Eigen::PlainObjectBase<DerivedF>& F);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "triangles_from_strip.cpp"
+#endif
+
+#endif
+

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

tutorial/506_FrameField/main.cpp

@@ -201,7 +201,7 @@ int main(int argc, char *argv[])
   igl::barycenter(V_deformed, F, B_deformed);
 
   // Find the closest crossfield to the deformed frame field
-  igl::frame_to_cross_field(V,F,FF1_deformed,FF2_deformed,X1_deformed);
+igl::frame_to_cross_field(V_deformed,F,FF1_deformed,FF2_deformed,X1_deformed);
 
   // Find a smooth crossfield that interpolates the deformed constraints
   MatrixXd bc_x(b.size(),3);