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

Conflicts:
	include/igl/list_to_matrix.cpp
	include/igl/peel_outer_hull_layers.cpp

Former-commit-id: a90aa948c5bd56342b0d82ebee5cb0e6b8ed90a2

README.md

@@ -1,6 +1,6 @@
 # libigl - A simple C++ geometry processing library
 
-![](tutorial/images/libigl-logo.jpg)
+![](libigl-teaser.png)
 
 <https://github.com/libigl/libigl/>
 
@@ -23,7 +23,7 @@ group prototypes a lot in MATLAB, and we have a useful [conversion
 table](matlab-to-eigen.html) from
 MATLAB to libigl/Eigen.
 
-# Tutorial
+## Tutorial
 
 As of version 1.0, libigl includes an introductory
 [tutorial](tutorial/tutorial.html) that covers
@@ -79,7 +79,19 @@ libigl depends only on the [Eigen](http://eigen.tuxfamily.org) library.
 
 For more information see our [tutorial](tutorial/tutorial.html).
 
-# Download
+### GCC and the optional CGAL dependency
+The `include/igl/cgal/*.h` headers depend on CGAL. It has come to our attention
+that CGAL does not work properly with GCC 4.8. To the best of our knowledge,
+GCC 4.7 and clang will work correctly.
+
+### OpenMP and Windows
+Some of our functions will take advantage of OpenMP if available. However, it
+has come to our attention that Visual Studio + Eigen does not work properly
+with OpenMP. Since OpenMP only improves performance without affecting
+functionality we recommend avoiding OpenMP on Windows or proceeding with
+caution.
+
+## Download
 You can keep up to date by cloning a read-only copy of our GitHub
 [repository](https://github.com/libigl).
 
@@ -105,11 +117,11 @@ BibTeX entry:
   title = {{libigl}: A simple {C++} geometry processing library},
   author = {Alec Jacobson and Daniele Panozzo and others},
   note = {http://libigl.github.io/libigl/},
-  year = {2014},
+  year = {2015},
 }
 ```
 
-# Contact
+## Contact
 
 Libigl is a group endeavor led by [Alec
 Jacobson](http://www.cs.columbia.edu/~jacobson/) and [Daniele
@@ -125,6 +137,8 @@ spending time maintaining this.
 If you find bugs or have problems please use our [github issue tracking
 page](https://github.com/libigl/libigl/issues).
 
-### Copyright
-2014 Alec Jacobson, Daniele Panozzo, Olga Diamanti, Kenshi
+## Copyright
+2015 Alec Jacobson, Daniele Panozzo, Olga Diamanti, Christian Schüller, Kenshi
 Takayama, Leo Sacht, Wenzel Jacob, Nico Pietroni, Amir Vaxman
+
+![](tutorial/images/libigl-logo.jpg)

RELEASE_HISTORY.txt

@@ -1,3 +1,4 @@
+1.1.3  Bug fixes in active set and boundary_conditions
 1.1.1  PLY file format support
 1.1.0  Mesh boolean operations using CGAL and cork, implementing [Attene 14]
 1.0.3  Bone heat method

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.1.1
+1.1.3

build/Makefile.conf

@@ -107,7 +107,7 @@ endif
 ifeq ($(IGL_USERNAME),daniele)
 	IGL_WITH_MATLAB=0
 	AFLAGS=-m64
-	GG=g++
+	GG=g++-4.8
 	EIGEN3_INC=-I/usr/local/include/eigen3
 endif
 

examples/skeleton-builder/example.cpp

@@ -1,4 +1,3 @@
-#include <igl/draw_skeleton_3d.h>
 #include <igl/draw_skeleton_vector_graphics.h>
 #include <igl/two_axis_valuator_fixed_up.h>
 #include <igl/read_triangle_mesh.h>
@@ -41,13 +40,13 @@
 #include <igl/writeTGF.h>
 #include <igl/file_exists.h>
 #include <igl/centroid.h>
+#include <igl/draw_skeleton_3d.h>
 
 #include <Eigen/Core>
 #include <Eigen/Geometry>
 
 #ifdef __APPLE__
 #include <GLUT/glut.h>
-#include <Carbon/Carbon.h>
 #else
 #include <GL/glut.h>
 #endif
@@ -310,7 +309,7 @@ void display()
             }
           }
         }
-        draw_skeleton_3d(s.C,s.BE,MatrixXd(),colors);
+        draw_skeleton_3d(s.C,s.BE,MatrixXd(),colors,bbd*0.5);
         break;
       }
       case SKEL_STYLE_TYPE_VECTOR_GRAPHICS:
@@ -707,6 +706,7 @@ void init_relative()
   Vmid = 0.5*(Vmax + Vmin);
   centroid(V,F,Vcen);
   bbd = (Vmax-Vmin).norm();
+  cout<<"bbd: "<<bbd<<endl;
   camera.push_away(2);
 }
 

examples/skeleton-posing/example.cpp

@@ -52,7 +52,6 @@
 #ifndef GLUT_ACTIVE_COMMAND
 #  define GLUT_ACTIVE_COMMAND 9
 #endif
-#include <Carbon/Carbon.h>
 #else
 #include <GL/glut.h>
 #endif
@@ -296,7 +295,7 @@ void display()
       default:
       case SKEL_STYLE_TYPE_3D:
       {
-        draw_skeleton_3d(C,BE,T,s.colors);
+        draw_skeleton_3d(C,BE,T,s.colors,bbd*0.5);
         break;
       }
       case SKEL_STYLE_TYPE_VECTOR_GRAPHICS:
@@ -631,6 +630,28 @@ bool save_weights()
   }
 }
 
+bool save_mesh()
+{
+  using namespace std;
+  using namespace igl;
+  using namespace Eigen;
+  MatrixXd T;
+  forward_kinematics(C,BE,P,s.mouse.rotations(),T);
+  MatrixXd U = M*T;
+  string output_filename;
+  next_filename(output_pose_prefix,4,".obj",output_filename);
+  if(writeOBJ(output_filename,U,F))
+  {
+    cout<<GREENGIN("Current mesh written to "+
+      output_filename+".")<<endl;
+    return true;
+  }else
+  {
+    cout<<REDRUM("Writing to "+output_filename+" failed.")<<endl;
+    return false;
+  }
+}
+
 void key(unsigned char key, int mouse_x, int mouse_y)
 {
   using namespace std;
@@ -675,6 +696,10 @@ void key(unsigned char key, int mouse_x, int mouse_y)
       break;
     }
     case 'S':
+    {
+      save_mesh();
+      break;
+    }
     case 's':
     {
       save_pose();
@@ -897,7 +922,8 @@ int main(int argc, char * argv[])
   cout<<"D,d                    Deselect all."<<endl;
   cout<<"R                      Reset selected rotation."<<endl;
   cout<<"r                      Reset all rotations."<<endl;
-  cout<<"S,s                    Save current pose."<<endl;
+  cout<<"S                      Save current posed mesh."<<endl;
+  cout<<"s                      Save current skeleton pose."<<endl;
   cout<<"W,w                    Save current weights."<<endl;
   cout<<"Z,z                    Snap to canonical view."<<endl;
   cout<<"⌘ Z                    Undo."<<endl;

examples/skeleton/example.cpp

@@ -293,7 +293,7 @@ void display()
   {
     default:
     case SKEL_STYLE_TYPE_3D:
-      draw_skeleton_3d(C,BE,T);
+      draw_skeleton_3d(C,BE,T,MAYA_VIOLET,bbd*0.5);
       break;
     case SKEL_STYLE_TYPE_VECTOR_GRAPHICS:
       draw_skeleton_vector_graphics(C,BE,T);

include/igl/active_set.cpp

@@ -8,6 +8,7 @@
 #include "active_set.h"
 #include "min_quad_with_fixed.h"
 #include "slice.h"
+#include "slice_into.h"
 #include "cat.h"
 #include "matlab_format.h"
 
@@ -44,6 +45,9 @@ IGL_INLINE igl::SolverStatus igl::active_set(
   )
 {
 //#define ACTIVE_SET_CPP_DEBUG
+#ifdef ACTIVE_SET_CPP_DEBUG
+#  warning "ACTIVE_SET_CPP_DEBUG"
+#endif
   using namespace igl;
   using namespace Eigen;
   using namespace std;
@@ -259,39 +263,52 @@ IGL_INLINE igl::SolverStatus igl::active_set(
     }
 #endif
     
+    Eigen::PlainObjectBase<DerivedZ> sol;
+    if(known_i.size() == A.rows())
+    {
+      // Everything's fixed?
 #ifdef ACTIVE_SET_CPP_DEBUG
-    cout<<"  min_quad_with_fixed_precompute"<<endl;
+      cout<<"  everything's fixed."<<endl;
 #endif
-    if(!min_quad_with_fixed_precompute(A,known_i,Aeq_i,params.Auu_pd,data))
+      Z.resize(A.rows(),Y_i.cols());
+      slice_into(Y_i,known_i,1,Z);
+      sol.resize(0,Y_i.cols());
+      assert(Aeq_i.rows() == 0 && "All fixed but linearly constrained");
+    }else
     {
-      cerr<<"Error: min_quad_with_fixed precomputation failed."<<endl;
-      if(iter > 0 && Aeq_i.rows() > Aeq.rows())
+#ifdef ACTIVE_SET_CPP_DEBUG
+      cout<<"  min_quad_with_fixed_precompute"<<endl;
+#endif
+      if(!min_quad_with_fixed_precompute(A,known_i,Aeq_i,params.Auu_pd,data))
       {
-        cerr<<"  *Are you sure rows of [Aeq;Aieq] are linearly independent?*"<<
-          endl;
+        cerr<<"Error: min_quad_with_fixed precomputation failed."<<endl;
+        if(iter > 0 && Aeq_i.rows() > Aeq.rows())
+        {
+          cerr<<"  *Are you sure rows of [Aeq;Aieq] are linearly independent?*"<<
+            endl;
+        }
+        ret = SOLVER_STATUS_ERROR;
+        break;
       }
-      ret = SOLVER_STATUS_ERROR;
-      break;
-    }
 #ifdef ACTIVE_SET_CPP_DEBUG
-    cout<<"  min_quad_with_fixed_solve"<<endl;
+      cout<<"  min_quad_with_fixed_solve"<<endl;
 #endif
-    Eigen::PlainObjectBase<DerivedZ> sol;
-    if(!min_quad_with_fixed_solve(data,B,Y_i,Beq_i,Z,sol))
-    {
-      cerr<<"Error: min_quad_with_fixed solve failed."<<endl;
-      ret = SOLVER_STATUS_ERROR;
-      break;
-    }
-    //cout<<matlab_format((Aeq*Z-Beq).eval(),"cr")<<endl;
-    //cout<<matlab_format(Z,"Z")<<endl;
+      if(!min_quad_with_fixed_solve(data,B,Y_i,Beq_i,Z,sol))
+      {
+        cerr<<"Error: min_quad_with_fixed solve failed."<<endl;
+        ret = SOLVER_STATUS_ERROR;
+        break;
+      }
+      //cout<<matlab_format((Aeq*Z-Beq).eval(),"cr")<<endl;
+      //cout<<matlab_format(Z,"Z")<<endl;
 #ifdef ACTIVE_SET_CPP_DEBUG
-    cout<<"  post"<<endl;
+      cout<<"  post"<<endl;
 #endif
+      // Computing Lagrange multipliers needs to be adjusted slightly if A is not symmetric
+      assert(data.Auu_sym);
+    }
 
     // Compute Lagrange multiplier values for known_i
-    // This needs to be adjusted slightly if A is not symmetric
-    assert(data.Auu_sym);
     SparseMatrix<AT> Ak;
     // Slow
     slice(A,known_i,1,Ak);

include/igl/angles.cpp

@@ -0,0 +1,49 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2015 Daniele Panozzo <daniele.panozzo@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "angles.h"
+#include <cassert>
+
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename Derivedtheta>
+void igl::angles(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<Derivedtheta>& theta)
+{
+  theta.resize(F.rows(),F.cols());
+
+  auto corner = [](const Eigen::PlainObjectBase<DerivedV>& x, const Eigen::PlainObjectBase<DerivedV>& y, const Eigen::PlainObjectBase<DerivedV>& z)
+  {
+    Eigen::RowVector3d v1 = (x-y).normalized();
+    Eigen::RowVector3d v2 = (z-y).normalized();
+
+    // http://stackoverflow.com/questions/10133957/signed-angle-between-two-vectors-without-a-reference-plane
+    double s = v1.cross(v2).norm();
+    double c = v1.dot(v2);
+
+    return atan2(s, c);
+  };
+
+  for(unsigned i=0; i<F.rows(); ++i)
+  {
+    for(unsigned j=0; j<F.cols(); ++j)
+    {
+      theta(i,j) = corner(
+        V.row(F(i,int(j-1+F.cols())%F.cols())),
+        V.row(F(i,j)),
+        V.row(F(i,(j+1+F.cols())%F.cols()))
+        );
+    }
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+#endif

include/igl/angles.h

@@ -0,0 +1,38 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2015 Daniele Panozzo <daniele.panozzo@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_ANGLES_H
+#define IGL_ANGLES_H
+#warning "Deprecated. Use igl/internal_angles.h instead"
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // ANGLES Compute angles for each corner of each triangle
+  //
+  // Inputs:
+  //   V  #V by dim list of vertex positions
+  //   F  #V by 3[4] list of triangle[quads] indices
+  // Outputs:
+  //   theta  #F by 3[4] list of angles for each corner (in radians)
+  //
+  template <
+  typename DerivedV,
+  typename DerivedF,
+  typename Derivedtheta>
+  IGL_INLINE void angles(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<Derivedtheta>& theta);
+
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "angles.cpp"
+#endif
+
+#endif

include/igl/any_of.h

@@ -15,6 +15,8 @@ namespace igl
   // Inputs:
   //   S  matrix
   // Returns whether any entries are true
+  //
+  // Seems that Eigen (now) implements this for `Eigen::Array` 
   template <typename Mat>
   IGL_INLINE bool any_of(const Mat & S);
 }

include/igl/bbw/bbw.cpp

@@ -104,11 +104,11 @@ IGL_INLINE bool igl::bbw(
     {
       case QP_SOLVER_IGL_ACTIVE_SET:
       {
-        //if(data.verbosity >= 1)
-        //{
+        if(data.verbosity >= 1)
+        {
           cout<<"BBW: max_iter: "<<data.active_set_params.max_iter<<endl;
-          cout<<"BBW: max_iter: "<<eff_params.max_iter<<endl;
-        //}
+          cout<<"BBW: eff_max_iter: "<<eff_params.max_iter<<endl;
+        }
         if(data.verbosity >= 1)
         {
           cout<<"BBW: Computing initial weights for "<<m<<" handle"<<

include/igl/boolean/mesh_boolean.cpp

@@ -1,10 +1,12 @@
 #include "mesh_boolean.h"
-#include <igl/peal_outer_hull_layers.h>
+#include <igl/peel_outer_hull_layers.h>
 #include <igl/cgal/remesh_self_intersections.h>
 #include <igl/remove_unreferenced.h>
 #include <igl/unique_simplices.h>
 #include <iostream>
 
+//#define IGL_MESH_BOOLEAN_DEBUG
+
 template <
   typename DerivedVA,
   typename DerivedFA,
@@ -97,10 +99,16 @@ IGL_INLINE void igl::mesh_boolean(
   typedef Matrix<Index,Dynamic,2> MatrixX2I;
   typedef Matrix<Index,Dynamic,1> VectorXI;
   typedef Matrix<typename DerivedJ::Scalar,Dynamic,1> VectorXJ;
+#ifdef IGL_MESH_BOOLEAN_DEBUG
+  cout<<"mesh boolean..."<<endl;
+#endif
   MatrixX3S V(VA.rows()+VB.rows(),3);
   MatrixX3I F(FA.rows()+FB.rows(),3);
   V.block(0,0,VA.rows(),VA.cols()) = VA;
   V.block(VA.rows(),0,VB.rows(),VB.cols()) = VB;
+#ifdef IGL_MESH_BOOLEAN_DEBUG
+  cout<<"prepare selfintersect input..."<<endl;
+#endif
   switch(type)
   {
     // Minus is implemented by flipping B and computing union
@@ -138,6 +146,9 @@ IGL_INLINE void igl::mesh_boolean(
     }
   };
 
+#ifdef IGL_MESH_BOOLEAN_DEBUG
+  cout<<"resolve..."<<endl;
+#endif
   MatrixX3S CV;
   MatrixX3I CF;
   VectorXJ CJ;
@@ -157,12 +168,18 @@ IGL_INLINE void igl::mesh_boolean(
     return;
   }
 
+#ifdef IGL_MESH_BOOLEAN_DEBUG
+  cout<<"peel..."<<endl;
+#endif
   Matrix<bool,Dynamic,1> from_A(CF.rows());
-  // Peal layers keeping track of odd and even flips
+  // peel layers keeping track of odd and even flips
   Matrix<bool,Dynamic,1> odd;
   Matrix<bool,Dynamic,1> flip;
-  peal_outer_hull_layers(CV,CF,odd,flip);
+  peel_outer_hull_layers(CV,CF,odd,flip);
 
+#ifdef IGL_MESH_BOOLEAN_DEBUG
+  cout<<"categorize..."<<endl;
+#endif
   const Index m = CF.rows();
   // Faces of output vG[i] = j means ith face of output should be jth face in F
   std::vector<Index> vG;
@@ -205,6 +222,9 @@ IGL_INLINE void igl::mesh_boolean(
     G.row(g) = Gflip[g] ? CF.row(vG[g]).reverse().eval() : CF.row(vG[g]);
     GJ(g) = CJ(vG[g]);
   }
+#ifdef IGL_MESH_BOOLEAN_DEBUG
+  cout<<"clean..."<<endl;
+#endif
   // remove duplicates: cancel out in all cases? assertion in others?
   {
     MatrixXi oldG = G;

include/igl/boundary_conditions.cpp

@@ -35,10 +35,9 @@ IGL_INLINE bool igl::boundary_conditions(
     return false;
   }
 
-
   vector<int> bci;
   vector<int> bcj;
-  vector<int> bcv;
+  vector<double> bcv;
 
   // loop over points
   for(int p = 0;p<P.size();p++)
@@ -95,7 +94,6 @@ IGL_INLINE bool igl::boundary_conditions(
     }
   }
 
-  // Cage edges are not considered yet
   // loop over cage edges
   for(int e = 0;e<CE.rows();e++)
   {
@@ -153,7 +151,7 @@ IGL_INLINE bool igl::boundary_conditions(
   for(i = 0;i<bc.rows();i++)
   {
     double sum = bc.row(i).sum();
-    assert(sum != 0);
+    assert(sum != 0 && "Some boundary vertex getting all zero BCs");
     bc.row(i).array() /= sum;
   }
 

include/igl/boundary_loop.cpp

@@ -1,105 +1,137 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-//
+// 
 // Copyright (C) 2014 Stefan Brugger <stefanbrugger@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 "boundary_loop.h"
-
+#include "slice.h"
 #include "triangle_triangle_adjacency.h"
 #include "vertex_triangle_adjacency.h"
+#include "is_border_vertex.h"
+#include <set>
 
+template <typename DerivedF, typename Index>
 IGL_INLINE void igl::boundary_loop(
-    const Eigen::MatrixXd& V,
-    const Eigen::MatrixXi& F,
-    Eigen::VectorXi& b)
+    const Eigen::PlainObjectBase<DerivedF> & F, 
+    std::vector<std::vector<Index> >& L)
 {
-  std::vector<int> bnd;
-  bnd.clear();
-  std::vector<bool> isVisited(V.rows(),false);
-
-  // Actually mesh only needs to be manifold near boundary, so this is
-  // over zealous (see gptoolbox's outline_loop for a more general
-  // (and probably faster) implementation)
-  assert(is_edge_manifold(V,F) && "Mesh must be manifold");
-  Eigen::MatrixXi TT,TTi;
-  std::vector<std::vector<int> > VF, VFi;
-  igl::triangle_triangle_adjacency(V,F,TT,TTi);
-  igl::vertex_triangle_adjacency(V,F,VF,VFi);
-
-  // Extract one boundary edge
-  bool done = false;
-  for (int i = 0; i < TT.rows() && !done; i++)
+  using namespace std;
+  using namespace Eigen;
+  using namespace igl;
+
+  MatrixXd Vdummy(F.maxCoeff(),1);
+  MatrixXi TT,TTi;
+  vector<std::vector<int> > VF, VFi;
+  triangle_triangle_adjacency(Vdummy,F,TT,TTi);
+  vertex_triangle_adjacency(Vdummy,F,VF,VFi);
+
+  vector<bool> unvisited = is_border_vertex(Vdummy,F);
+  set<int> unseen;
+  for (int i = 0; i < unvisited.size(); ++i)
   {
-    for (int j = 0; j < TT.cols(); j++)
-    {
-      if (TT(i,j) < 0)
-      {
-        int idx1, idx2;
-        idx1 = F(i,j);
-        idx2 = F(i,(j+1) % F.cols());
-        bnd.push_back(idx1);
-        bnd.push_back(idx2);
-        isVisited[idx1] = true;
-        isVisited[idx2] = true;
-        done = true;
-        break;
-      }
-    }
+    if (unvisited[i])
+      unseen.insert(unseen.end(),i);
   }
 
-  // Traverse boundary
-  while(1)
+  while (!unseen.empty())
   {
-    bool changed = false;
-    int lastV;
-    lastV = bnd[bnd.size()-1];
+    vector<Index> l;
 
-    for (int i = 0; i < (int)VF[lastV].size(); i++)
-    {
-      int curr_neighbor = VF[lastV][i];
+    // Get first vertex of loop
+    int start = *unseen.begin();
+    unseen.erase(unseen.begin());
+    unvisited[start] = false;
+    l.push_back(start);
 
-      if (TT.row(curr_neighbor).minCoeff() < 0.) // Face contains boundary edge
+    bool done = false;
+    while (!done)
+    {
+      // Find next vertex
+      bool newBndEdge = false;
+      int v = l[l.size()-1];
+      int next;
+      for (int i = 0; i < (int)VF[v].size() && !newBndEdge; i++)
       {
-        int idx_lastV_in_face;
-        if (F(curr_neighbor,0) == lastV) idx_lastV_in_face = 0;
-        if (F(curr_neighbor,1) == lastV) idx_lastV_in_face = 1;
-        if (F(curr_neighbor,2) == lastV) idx_lastV_in_face = 2;
-
-        int idx_prev = (idx_lastV_in_face + F.cols()-1) % F.cols();
-        int idx_next = (idx_lastV_in_face + 1) % F.cols();
-        bool isPrevVisited = isVisited[F(curr_neighbor,idx_prev)];
-        bool isNextVisited = isVisited[F(curr_neighbor,idx_next)];
-
-        bool gotBndEdge = false;
-        int next_bnd_vertex;
-        if (!isNextVisited && TT(curr_neighbor,idx_lastV_in_face) < 0)
-        {
-          next_bnd_vertex = idx_next;
-          gotBndEdge = true;
-        }
-        else if (!isPrevVisited && TT(curr_neighbor,(idx_lastV_in_face+2) % F.cols()) < 0)
-        {
-          next_bnd_vertex = idx_prev;
-          gotBndEdge = true;
-        }
+        int fid = VF[v][i];
 
-        if (gotBndEdge)
+        if (TT.row(fid).minCoeff() < 0.) // Face contains boundary edge
         {
-          changed = true;
-          bnd.push_back(F(curr_neighbor,next_bnd_vertex));
-          isVisited[F(curr_neighbor,next_bnd_vertex)] = true;
-          break;
+          int vLoc;
+          if (F(fid,0) == v) vLoc = 0;
+          if (F(fid,1) == v) vLoc = 1;
+          if (F(fid,2) == v) vLoc = 2;
+
+          int vPrev = F(fid,(vLoc + F.cols()-1) % F.cols());
+          int vNext = F(fid,(vLoc + 1) % F.cols());
+
+          bool newBndEdge = false;
+          if (unvisited[vPrev] && TT(fid,(vLoc+2) % F.cols()) < 0)
+          {
+            next = vPrev;
+            newBndEdge = true;
+          }
+          else if (unvisited[vNext] && TT(fid,vLoc) < 0)
+          {
+            next = vNext;
+            newBndEdge = true;
+          }
         }
       }
-    }
 
-    if (!changed)
-      break;
+      if (newBndEdge)
+      {
+        l.push_back(next);
+        unseen.erase(next);
+        unvisited[next] = false;
+      }
+      else
+        done = true;
+    }
+    L.push_back(l);
   }
+}
 
-  b.resize(bnd.size());
-  for(unsigned i=0;i<bnd.size();++i)
-    b(i) = bnd[i];
+template <typename DerivedF, typename Index>
+IGL_INLINE void igl::boundary_loop(
+  const Eigen::PlainObjectBase<DerivedF>& F, 
+  std::vector<Index>& L)
+{
+  using namespace Eigen;
+  using namespace std;
+
+  vector<vector<int> > Lall;
+  boundary_loop(F,Lall);
+
+  int idxMax = -1;
+  int maxLen = 0;
+  for (int i = 0; i < Lall.size(); ++i)
+  {
+    if (Lall[i].size() > maxLen)
+    {
+      maxLen = Lall[i].size();
+      idxMax = i;
+    }
+  }   
+
+  L.resize(Lall[idxMax].size());
+  for (int i = 0; i < Lall[idxMax].size(); ++i)
+    L[i] = Lall[idxMax][i];
 }
+
+template <typename DerivedF, typename DerivedL>
+IGL_INLINE void igl::boundary_loop(
+  const Eigen::PlainObjectBase<DerivedF>& F, 
+  Eigen::PlainObjectBase<DerivedL>& L)
+{
+  using namespace Eigen;
+  using namespace std;
+
+  vector<int> Lvec;
+  boundary_loop(F,Lvec);
+
+  L.resize(Lvec.size());
+  for (int i = 0; i < Lvec.size(); ++i)
+    L(i) = Lvec[i];
+}

include/igl/boundary_loop.h

@@ -7,30 +7,60 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #ifndef IGL_BOUNDARY_LOOP_H
 #define IGL_BOUNDARY_LOOP_H
-#include <igl/igl_inline.h>
+#include "igl_inline.h"
 
 #include <Eigen/Dense>
 #include <vector>
 
 namespace igl
 {
+  // Compute list of ordered boundary loops for a manifold mesh.
+  //
+  // Templates:
+  //  Index  index type
+  // Inputs:
+  //   F  #V by dim list of mesh faces
+  // Outputs:
+  //   L  list of loops where L[i] = ordered list of boundary vertices in loop i
+  //
+  template <typename DerivedF, typename Index>
+  IGL_INLINE void boundary_loop(
+    const Eigen::PlainObjectBase<DerivedF>& F, 
+    std::vector<std::vector<Index> >& L);
+
+
+  // Compute ordered boundary loops for a manifold mesh and return the 
+  // longest loop in terms of vertices.
+  //
+  // Templates:
+  //  Index  index type
+  // Inputs:
+  //   F  #V by dim list of mesh faces
+  // Outputs:
+  //   L  ordered list of boundary vertices of longest boundary loop
+  //
+  template <typename DerivedF, typename Index>
+  IGL_INLINE void boundary_loop(
+    const Eigen::PlainObjectBase<DerivedF>& F, 
+    std::vector<Index>& L);
 
-  // Compute sorted list of boundary vertices for a manifold mesh with single
-  // boundary.
+  // Compute ordered boundary loops for a manifold mesh and return the 
+  // longest loop in terms of vertices.
   //
+  // Templates:
+  //  Index  index type
   // Inputs:
-  //   V  #V by dim list of mesh vertex positions
   //   F  #V by dim list of mesh faces
   // Outputs:
-  //   bnd   sorted list of boundary vertex indices
+  //   L  ordered list of boundary vertices of longest boundary loop
+  //
+  template <typename DerivedF, typename DerivedL>
   IGL_INLINE void boundary_loop(
-  	const Eigen::MatrixXd& V, 
-  	const Eigen::MatrixXi& F, 
-    Eigen::VectorXi& bnd);
+    const Eigen::PlainObjectBase<DerivedF>& F, 
+    Eigen::PlainObjectBase<DerivedL>& L);
 }
 
 #ifndef IGL_STATIC_LIBRARY
 #  include "boundary_loop.cpp"
 #endif
-
 #endif

include/igl/cgal/CGAL_includes.hpp

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

include/igl/cgal/SelfIntersectMesh.h

@@ -212,10 +212,10 @@ namespace igl
       // Getters:
     public:
       //const IndexList& get_lIF() const{ return lIF;}
-      static inline void box_intersect(
+      static inline void box_intersect_static(
         SelfIntersectMesh * SIM, 
-        const SelfIntersectMesh::Box &a, 
-        const SelfIntersectMesh::Box &b);
+        const Box &a, 
+        const Box &b);
   };
 }
 
@@ -227,9 +227,8 @@ namespace igl
 #include <igl/get_seconds.h>
 #include <igl/C_STR.h>
 
-#include <boost/function.hpp>
-#include <boost/bind.hpp>
 
+#include <functional>
 #include <algorithm>
 #include <exception>
 #include <cassert>
@@ -277,7 +276,7 @@ inline void igl::SelfIntersectMesh<
   DerivedFF,
   DerivedIF,
   DerivedJ,
-  DerivedIM>::box_intersect(
+  DerivedIM>::box_intersect_static(
   Self * SIM, 
   const typename Self::Box &a, 
   const typename Self::Box &b)
@@ -350,8 +349,12 @@ inline igl::SelfIntersectMesh<
     boxes.push_back(Box(tit->bbox(), tit));
   }
   // Leapfrog callback
-  boost::function<void(const Box &a,const Box &b)> cb
-    = boost::bind(&box_intersect, this, _1,_2);
+  std::function<void(const Box &a,const Box &b)> cb = 
+    std::bind(&box_intersect_static, this, 
+      // Explicitly use std namespace to avoid confusion with boost (who puts
+      // _1 etc. in global namespace)
+      std::placeholders::_1,
+      std::placeholders::_2);
   //cout<<"boxes and bind: "<<tictoc()<<endl;
   // Run the self intersection algorithm with all defaults
   try{
@@ -911,7 +914,7 @@ inline bool igl::SelfIntersectMesh<
     try
     {
       CGAL::Object result = CGAL::intersection(A,B);
-      if(result)
+      if(!result.empty())
       {
         if(CGAL::object_cast<Segment_3 >(&result))
         {

include/igl/cgal/complex_to_mesh.cpp

@@ -11,6 +11,7 @@
 #include <igl/remove_unreferenced.h>
 
 #include <CGAL/Surface_mesh_default_triangulation_3.h>
+#include <CGAL/Triangulation_cell_base_with_circumcenter_3.h>
 #include <set>
 #include <stack>
 

include/igl/cgal/mesh_to_polyhedron.cpp

@@ -0,0 +1,45 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 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 "mesh_to_polyhedron.h"
+#include <CGAL/Polyhedron_3.h>
+#include <CGAL/Polyhedron_incremental_builder_3.h>
+template <typename Polyhedron>
+IGL_INLINE bool igl::mesh_to_polyhedron(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  Polyhedron & poly)
+{
+  typedef typename Polyhedron::HalfedgeDS HalfedgeDS;
+  // Postcondition: hds is a valid polyhedral surface.
+  CGAL::Polyhedron_incremental_builder_3<HalfedgeDS> B(poly.hds());
+  B.begin_surface(V.rows(),F.rows());
+  typedef typename HalfedgeDS::Vertex   Vertex;
+  typedef typename Vertex::Point Point;
+  assert(V.cols() == 3 && "V must be #V by 3");
+  for(int v = 0;v<V.rows();v++)
+  {
+    B.add_vertex(Point(V(v,0),V(v,1),V(v,2)));
+  }
+  assert(F.cols() == 3 && "F must be #F by 3");
+  for(int f=0;f<F.rows();f++)
+  {
+    B.begin_facet();
+    for(int c = 0;c<3;c++)
+    {
+      B.add_vertex_to_facet(F(f,c));
+    }
+    B.end_facet();
+  }
+  if(B.error())
+  {
+    B.rollback();
+    return false;
+  }
+  B.end_surface();
+  return poly.is_valid();
+}

include/igl/cgal/mesh_to_polyhedron.h

@@ -0,0 +1,36 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 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_MESH_TO_POLYHEDRON_H
+#define IGL_MESH_TO_POLYHEDRON_H
+#include <igl/igl_inline.h>
+#include <Eigen/Core>
+
+namespace igl
+{
+  // Convert a mesh (V,F) to a CGAL Polyhedron
+  //
+  // Templates:
+  //   Polyhedron  CGAL Polyhedron type (e.g. Polyhedron_3)
+  // Inputs:
+  //   V  #V by 3 list of vertex positions
+  //   F  #F by 3 list of triangle indices
+  // Outputs:
+  //   poly  cgal polyhedron
+  // Returns true only if (V,F) can be converted to a valid polyhedron (i.e. if
+  // (V,F) is vertex and edge manifold).
+  template <typename Polyhedron>
+  IGL_INLINE bool mesh_to_polyhedron(
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    Polyhedron & poly);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "mesh_to_polyhedron.cpp"
+#endif
+
+#endif

include/igl/cgal/point_mesh_squared_distance.cpp

@@ -55,6 +55,16 @@ IGL_INLINE void igl::point_mesh_squared_distance_precompute(
   assert(V.cols() == 3);
   // Must be triangles
   assert(F.cols() == 3);
+
+  // WTF ALERT!!!! 
+  //
+  // There's a bug in clang probably or at least in cgal. Without calling this
+  // line (I guess invoking some compilation from <vector>), clang will vomit
+  // errors inside CGAL.
+  //
+  // http://stackoverflow.com/questions/27748442/is-clangs-c11-support-reliable
+  T.reserve(0);
+
   // Make list of cgal triangles
   mesh_to_cgal_triangle_list(V,F,T);
   tree.clear();
@@ -110,7 +120,7 @@ IGL_INLINE void igl::point_mesh_squared_distance(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template void igl::point_mesh_squared_distance<CGAL::Epeck>( const Eigen::MatrixXd & P, const Eigen::MatrixXd & V, const Eigen::MatrixXi & F, Eigen::VectorXd & sqrD, Eigen::VectorXi & I, Eigen::MatrixXd & C);
 template void igl::point_mesh_squared_distance_precompute<CGAL::Epick>(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, CGAL::AABB_tree<CGAL::AABB_traits<CGAL::Epick, CGAL::AABB_triangle_primitive<CGAL::Epick, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >::iterator, CGAL::Boolean_tag<false> > > >&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >&);
-
+template void igl::point_mesh_squared_distance<CGAL::Epeck>( const Eigen::MatrixXd & P, const Eigen::MatrixXd & V, const Eigen::MatrixXi & F, Eigen::VectorXd & sqrD, Eigen::VectorXi & I, Eigen::MatrixXd & C);
+template void igl::point_mesh_squared_distance<CGAL::Epick>(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, 1, 0, -1, 1>&, Eigen::Matrix<int, -1, 1, 0, -1, 1>&, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
 #endif

include/igl/cgal/polyhedron_to_mesh.cpp

@@ -0,0 +1,57 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 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 "polyhedron_to_mesh.h"
+#include <CGAL/Polyhedron_3.h>
+
+template <typename Polyhedron>
+IGL_INLINE void igl::polyhedron_to_mesh(
+  const Polyhedron & poly,
+  Eigen::MatrixXd & V,
+  Eigen::MatrixXi & F)
+{
+  using namespace std;
+  V.resize(poly.size_of_vertices(),3);
+  F.resize(poly.size_of_facets(),3);
+  typedef typename Polyhedron::Vertex_const_iterator Vertex_iterator;
+  std::map<Vertex_iterator,size_t> vertex_to_index;
+  {
+    size_t v = 0;
+    for(
+      typename Polyhedron::Vertex_const_iterator p = poly.vertices_begin();
+      p != poly.vertices_end();
+      p++)
+    {
+      V(v,0) = p->point().x();
+      V(v,1) = p->point().y();
+      V(v,2) = p->point().z();
+      vertex_to_index[p] = v;
+      v++;
+    }
+  }
+  {
+    size_t f = 0;
+    for(
+      typename Polyhedron::Facet_const_iterator facet = poly.facets_begin();
+      facet != poly.facets_end();
+      ++facet)
+    {
+      typename Polyhedron::Halfedge_around_facet_const_circulator he = 
+        facet->facet_begin();
+      // Facets in polyhedral surfaces are at least triangles.
+      assert(CGAL::circulator_size(he) == 3 && "Facets should be triangles");
+      size_t c = 0;
+      do {
+        //// This is stooopidly slow
+        // F(f,c) = std::distance(poly.vertices_begin(), he->vertex());
+        F(f,c) = vertex_to_index[he->vertex()];
+        c++;
+      } while ( ++he != facet->facet_begin());
+      f++;
+    }
+  }
+}

include/igl/cgal/polyhedron_to_mesh.h

@@ -0,0 +1,34 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 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_POLYHEDRON_TO_MESH_H
+#define IGL_POLYHEDRON_TO_MESH_H
+#include <igl/igl_inline.h>
+#include <Eigen/Core>
+
+namespace igl
+{
+  // Convert a CGAL Polyhedron to a mesh (V,F)
+  //
+  // Templates:
+  //   Polyhedron  CGAL Polyhedron type (e.g. Polyhedron_3)
+  // Inputs:
+  //   poly  cgal polyhedron
+  // Outputs:
+  //   V  #V by 3 list of vertex positions
+  //   F  #F by 3 list of triangle indices
+  template <typename Polyhedron>
+  IGL_INLINE void polyhedron_to_mesh(
+    const Polyhedron & poly,
+    Eigen::MatrixXd & V,
+    Eigen::MatrixXi & F);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "polyhedron_to_mesh.cpp"
+#endif
+
+#endif

include/igl/cgal/signed_distance.cpp

@@ -282,6 +282,9 @@ IGL_INLINE void igl::signed_distance_winding_number(
 }
 
 #ifdef IGL_STATIC_LIBRARY
+// This template is necessary for the others to compile with clang
+// http://stackoverflow.com/questions/27748442/is-clangs-c11-support-reliable
+template void igl::signed_distance<CGAL::Epick>( const Eigen::MatrixXd & , const Eigen::MatrixXd & , const Eigen::MatrixXi & , const SignedDistanceType , Eigen::VectorXd & , Eigen::VectorXi &, Eigen::MatrixXd & , Eigen::MatrixXd & );
 template CGAL::Epick::FT igl::signed_distance_winding_number<CGAL::Epick>(CGAL::AABB_tree<CGAL::AABB_traits<CGAL::Epick, CGAL::AABB_triangle_primitive<CGAL::Epick, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >::iterator, CGAL::Boolean_tag<false> > > > const&, igl::WindingNumberAABB<Eigen::Matrix<double, 3, 1, 0, 3, 1> > const&, CGAL::Epick::Point_3 const&);
 template CGAL::Epick::FT igl::signed_distance_pseudonormal<CGAL::Epick>(CGAL::AABB_tree<CGAL::AABB_traits<CGAL::Epick, CGAL::AABB_triangle_primitive<CGAL::Epick, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >::iterator, CGAL::Boolean_tag<false> > > > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, CGAL::Epick::Point_3 const&);
 #endif

include/igl/collapse_small_triangles.cpp

@@ -0,0 +1,132 @@
+#include "collapse_small_triangles.h"
+
+#include "bounding_box_diagonal.h"
+#include "doublearea.h"
+#include "edge_lengths.h"
+#include "colon.h"
+#include "faces_first.h"
+
+#include <limits>
+
+#include <iostream>
+
+void igl::collapse_small_triangles(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const double eps,
+  Eigen::MatrixXi & FF)
+{
+  using namespace Eigen;
+  using namespace std;
+
+  // Compute bounding box diagonal length
+  double bbd = bounding_box_diagonal(V);
+  MatrixXd l;
+  edge_lengths(V,F,l);
+  VectorXd dblA;
+  doublearea(l,dblA);
+
+  // Minimum area tolerance
+  const double min_dblarea = 2.0*eps*bbd*bbd;
+
+  Eigen::VectorXi FIM = colon<int>(0,V.rows()-1);
+  int num_edge_collapses = 0;
+  // Loop over triangles
+  for(int f = 0;f<F.rows();f++)
+  {
+    if(dblA(f) < min_dblarea)
+    {
+      double minl = 0;
+      int minli = -1;
+      // Find shortest edge
+      for(int e = 0;e<3;e++)
+      {
+        if(minli==-1 || l(f,e)<minl)
+        {
+          minli = e;
+          minl = l(f,e);
+        }
+      }
+      double maxl = 0;
+      int maxli = -1;
+      // Find longest edge
+      for(int e = 0;e<3;e++)
+      {
+        if(maxli==-1 || l(f,e)>maxl)
+        {
+          maxli = e;
+          maxl = l(f,e);
+        }
+      }
+      // Be sure that min and max aren't the same
+      maxli = (minli==maxli?(minli+1)%3:maxli);
+
+      // Collapse min edge maintaining max edge: i-->j
+      // Q: Why this direction?
+      int i = maxli;
+      int j = ((minli+1)%3 == maxli ? (minli+2)%3: (minli+1)%3);
+      assert(i != minli);
+      assert(j != minli);
+      assert(i != j);
+      FIM(F(f,i)) = FIM(F(f,j));
+      num_edge_collapses++;
+    }
+  }
+
+  // Reindex faces
+  MatrixXi rF = F;
+  // Loop over triangles
+  for(int f = 0;f<rF.rows();f++)
+  {
+    for(int i = 0;i<rF.cols();i++)
+    {
+      rF(f,i) = FIM(rF(f,i));
+    }
+  }
+
+  FF.resize(rF.rows(),rF.cols());
+  int num_face_collapses=0;
+  // Only keep uncollapsed faces
+  {
+    int ff = 0;
+    // Loop over triangles
+    for(int f = 0;f<rF.rows();f++)
+    {
+      bool collapsed = false;
+      // Check if any indices are the same
+      for(int i = 0;i<rF.cols();i++)
+      {
+        for(int j = i+1;j<rF.cols();j++)
+        {
+          if(rF(f,i)==rF(f,j))
+          {
+            collapsed = true;
+            num_face_collapses++;
+            break;
+          }
+        }
+      }
+      if(!collapsed)
+      {
+        FF.row(ff++) = rF.row(f);
+      }
+    }
+    // Use conservative resize
+    FF.conservativeResize(ff,FF.cols());
+  }
+  //cout<<"num_edge_collapses: "<<num_edge_collapses<<endl;
+  //cout<<"num_face_collapses: "<<num_face_collapses<<endl;
+  if(num_edge_collapses == 0)
+  {
+    // There must have been a "collapsed edge" in the input
+    assert(num_face_collapses==0);
+    // Base case
+    return;
+  }
+
+  //// force base case
+  //return;
+
+  MatrixXi recFF = FF;
+  return collapse_small_triangles(V,recFF,eps,FF);
+}

include/igl/collapse_small_triangles.h

@@ -0,0 +1,40 @@
+// 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_COLLAPSE_SMALL_TRIANGLES_H
+#define IGL_COLLAPSE_SMALL_TRIANGLES_H
+#include <Eigen/Dense>
+namespace igl
+{
+  // Given a triangle mesh (V,F) compute a new mesh (VV,FF) which contains the
+  // original faces and vertices of (V,F) except any small triangles have been
+  // removed via collapse.
+  //
+  // We are *not* following the rules in "Mesh Optimization" [Hoppe et al]
+  // Section 4.2. But for our purposes we don't care about this criteria.
+  //
+  // Inputs:
+  //   V  #V by 3 list of vertex positions
+  //   F  #F by 3 list of triangle indices into V
+  //   eps  epsilon for smallest allowed area treated as fraction of squared bounding box
+  //     diagonal
+  // Outputs:
+  //   FF  #FF by 3 list of triangle indices into V
+  //
+  //
+  void collapse_small_triangles(
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    const double eps,
+    Eigen::MatrixXi & FF);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "collapse_small_triangles.cpp"
+#endif
+  
+#endif

include/igl/doublearea.cpp

@@ -17,10 +17,13 @@ IGL_INLINE void igl::doublearea(
   const Eigen::PlainObjectBase<DerivedF> & F,
   Eigen::PlainObjectBase<DeriveddblA> & dblA)
 {
+  if (F.cols() == 4) // quads are handled by a specialized function
+    return doublearea_quad(V,F,dblA);
+
   const int dim = V.cols();
   // Only support triangles
   assert(F.cols() == 3);
-  const int m = F.rows();
+  const size_t m = F.rows();
   // Compute edge lengths
   Eigen::PlainObjectBase<DerivedV> l;
   // "Lecture Notes on Geometric Robustness" Shewchuck 09, Section 3.1
@@ -42,7 +45,7 @@ IGL_INLINE void igl::doublearea(
     case 3:
     {
       dblA = Eigen::PlainObjectBase<DeriveddblA>::Zero(m,1);
-      for(int f = 0;f<F.rows();f++)
+      for(size_t f = 0;f<m;f++)
       {
         for(int d = 0;d<3;d++)
         {
@@ -56,7 +59,7 @@ IGL_INLINE void igl::doublearea(
     case 2:
     {
       dblA.resize(m,1);
-      for(int f = 0;f<m;f++)
+      for(size_t f = 0;f<m;f++)
       {
         dblA(f) = proj_doublearea(0,1,f);
       }
@@ -115,13 +118,13 @@ IGL_INLINE void igl::doublearea(
   // Only support triangles
   assert(ul.cols() == 3);
   // Number of triangles
-  const int m = ul.rows();
+  const size_t m = ul.rows();
   Eigen::PlainObjectBase<Derivedl> l;
   MatrixXi _;
   sort(ul,2,false,l,_);
   // semiperimeters
   Matrix<typename Derivedl::Scalar,Dynamic,1> s = l.rowwise().sum()*0.5;
-  assert(s.rows() == m);
+  assert((size_t)s.rows() == m);
   // resize output
   dblA.resize(l.rows(),1);
   // Minimum number of iterms per openmp thread
@@ -129,7 +132,7 @@ IGL_INLINE void igl::doublearea(
   #  define IGL_OMP_MIN_VALUE 1000
   #endif
   #pragma omp parallel for if (m>IGL_OMP_MIN_VALUE)
-  for(int i = 0;i<m;i++)
+  for(size_t i = 0;i<m;i++)
   {
     //// Heron's formula for area
     //const typename Derivedl::Scalar arg =
@@ -148,6 +151,35 @@ IGL_INLINE void igl::doublearea(
   }
 }
 
+template <typename DerivedV, typename DerivedF, typename DeriveddblA>
+IGL_INLINE void igl::doublearea_quad(
+const Eigen::PlainObjectBase<DerivedV> & V,
+const Eigen::PlainObjectBase<DerivedF> & F,
+Eigen::PlainObjectBase<DeriveddblA> & dblA)
+{
+  assert(V.cols() == 3); // Only supports points in 3D
+  assert(F.cols() == 4); // Only support quads
+  const size_t m = F.rows();
+
+  // Split the quads into triangles
+  Eigen::MatrixXi Ft(F.rows()*2,3);
+
+  for(size_t i=0; i<m;++i)
+  {
+    Ft.row(i*2    ) << F(i,0), F(i,1), F(i,2);
+    Ft.row(i*2 + 1) << F(i,2), F(i,3), F(i,0);
+  }
+
+  // Compute areas
+  Eigen::VectorXd doublearea_tri;
+  igl::doublearea(V,Ft,doublearea_tri);
+
+  dblA.resize(F.rows(),1);
+  for(unsigned i=0; i<F.rows();++i)
+    dblA(i) = doublearea_tri(i*2) + doublearea_tri(i*2 + 1);
+
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 // generated by autoexplicit.sh

include/igl/doublearea.h

@@ -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/.
 #ifndef IGL_DOUBLEAREA_H
 #define IGL_DOUBLEAREA_H
@@ -11,7 +11,7 @@
 #include <Eigen/Dense>
 namespace igl
 {
-  // DOUBLEAREA computes twice the area for each input triangle
+  // DOUBLEAREA computes twice the area for each input triangle[quad]
   //
   // Templates:
   //   DerivedV  derived type of eigen matrix for V (e.g. derived from
@@ -22,16 +22,16 @@ namespace igl
   //     MatrixXd)
   // Inputs:
   //   V  #V by dim list of mesh vertex positions
-  //   F  #F by simplex_size list of mesh faces (must be triangles)
+  //   F  #F by simplex_size list of mesh faces (must be triangles or quads)
   // Outputs:
-  //   dblA  #F list of triangle double areas (SIGNED only for 2D input)
+  //   dblA  #F list of triangle[quad] double areas (SIGNED only for 2D input)
   //
   // Known bug: For dim==3 complexity is O(#V + #F)!! Not just O(#F). This is a big deal
   // if you have 1million unreferenced vertices and 1 face
   template <typename DerivedV, typename DerivedF, typename DeriveddblA>
-  IGL_INLINE void doublearea( 
-    const Eigen::PlainObjectBase<DerivedV> & V, 
-    const Eigen::PlainObjectBase<DerivedF> & F, 
+  IGL_INLINE void doublearea(
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedF> & F,
     Eigen::PlainObjectBase<DeriveddblA> & dblA);
   // Stream of triangles
   template <
@@ -39,9 +39,9 @@ namespace igl
     typename DerivedB,
     typename DerivedC,
     typename DerivedD>
-  IGL_INLINE void doublearea( 
-    const Eigen::PlainObjectBase<DerivedA> & A, 
-    const Eigen::PlainObjectBase<DerivedB> & B, 
+  IGL_INLINE void doublearea(
+    const Eigen::PlainObjectBase<DerivedA> & A,
+    const Eigen::PlainObjectBase<DerivedB> & B,
     const Eigen::PlainObjectBase<DerivedC> & C,
     Eigen::PlainObjectBase<DerivedD> & D);
   // Single triangle in 2D!
@@ -51,9 +51,9 @@ namespace igl
     typename DerivedA,
     typename DerivedB,
     typename DerivedC>
-  IGL_INLINE typename DerivedA::Scalar doublearea_single( 
-    const Eigen::PlainObjectBase<DerivedA> & A, 
-    const Eigen::PlainObjectBase<DerivedB> & B, 
+  IGL_INLINE typename DerivedA::Scalar doublearea_single(
+    const Eigen::PlainObjectBase<DerivedA> & A,
+    const Eigen::PlainObjectBase<DerivedB> & B,
     const Eigen::PlainObjectBase<DerivedC> & C);
   // Same as above but use instrinsic edge lengths rather than (V,F) mesh
   // Templates:
@@ -64,14 +64,29 @@ namespace igl
   //   DeriveddblA  derived type of eigen matrix for dblA (e.g. derived from
   //     MatrixXd)
   // Inputs:
-  //   l  #F by dim list of edge lengths using 
+  //   l  #F by dim list of edge lengths using
   //     for triangles, columns correspond to edges 23,31,12
   // Outputs:
   //   dblA  #F list of triangle double areas
   template <typename Derivedl, typename DeriveddblA>
-  IGL_INLINE void doublearea( 
-    const Eigen::PlainObjectBase<Derivedl> & l, 
+  IGL_INLINE void doublearea(
+    const Eigen::PlainObjectBase<Derivedl> & l,
     Eigen::PlainObjectBase<DeriveddblA> & dblA);
+
+  // DOUBLEAREA_QUAD computes twice the area for each input quadrilateral
+  //
+  // Inputs:
+  //   V  #V by dim list of mesh vertex positions
+  //   F  #F by simplex_size list of mesh faces (must be quadrilaterals)
+  // Outputs:
+  //   dblA  #F list of quadrilateral double areas
+  //
+  template <typename DerivedV, typename DerivedF, typename DeriveddblA>
+  IGL_INLINE void doublearea_quad(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DeriveddblA> & dblA);
+
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/draw_skeleton_3d.cpp

@@ -21,7 +21,8 @@ IGL_INLINE void igl::draw_skeleton_3d(
   const Eigen::PlainObjectBase<DerivedC> & C,
   const Eigen::PlainObjectBase<DerivedBE> & BE,
   const Eigen::PlainObjectBase<DerivedT> & T,
-  const Eigen::PlainObjectBase<Derivedcolor> & color)
+  const Eigen::PlainObjectBase<Derivedcolor> & color,
+  const double half_bbd)
 {
   // Note: Maya's skeleton *does* scale with the mesh suggesting a scale
   // parameter. Further, its joint balls are not rotated with the bones.
@@ -29,7 +30,7 @@ IGL_INLINE void igl::draw_skeleton_3d(
   using namespace std;
   if(color.size() == 0)
   {
-    return draw_skeleton_3d(C,BE,T,MAYA_SEA_GREEN);
+    return draw_skeleton_3d(C,BE,T,MAYA_SEA_GREEN,half_bbd);
   }
   assert(color.cols() == 4 || color.size() == 4);
   if(T.size() == 0)
@@ -42,15 +43,13 @@ IGL_INLINE void igl::draw_skeleton_3d(
     {
       return;
     }
-    return draw_skeleton_3d(C,BE,T,color);
+    return draw_skeleton_3d(C,BE,T,color,half_bbd);
   }
   assert(T.rows() == BE.rows()*(C.cols()+1));
   assert(T.cols() == C.cols());
-  // old settings
-  int old_lighting=0;
-  double old_line_width=1;
-  glGetIntegerv(GL_LIGHTING,&old_lighting);
-  glGetDoublev(GL_LINE_WIDTH,&old_line_width);
+  // push old settings
+  glPushAttrib(GL_LIGHTING_BIT);
+  glPushAttrib(GL_LINE_BIT);
   glDisable(GL_LIGHTING);
   glLineWidth(1.0);
 
@@ -103,7 +102,7 @@ IGL_INLINE void igl::draw_skeleton_3d(
     auto s = C.row(BE(e,0));
     auto d = C.row(BE(e,1));
     auto b = (d-s).transpose().eval();
-    double r = 0.02;
+    double r = 0.02*half_bbd;
     Matrix4d Te = Matrix4d::Identity();
     Te.block(0,0,3,4) = T.block(e*4,0,4,3).transpose();
     Quaterniond q;
@@ -113,8 +112,7 @@ IGL_INLINE void igl::draw_skeleton_3d(
     glTranslated(s(0),s(1),s(2));
     if(color.size() == 4)
     {
-      Vector4d d = color.template cast<double>();
-      glColor4dv(d.data());
+      glColor4d( color(0), color(1), color(2), color(3));
     }else
     {
       glColor4d( color(e,0), color(e,1), color(e,2), color(e,3));
@@ -136,8 +134,8 @@ IGL_INLINE void igl::draw_skeleton_3d(
     glPopMatrix();
   }
   // Reset settings
-  (old_lighting ? glEnable(GL_LIGHTING) : glDisable(GL_LIGHTING));
-  glLineWidth(old_line_width);
+  glPopAttrib();
+  glPopAttrib();
 }
 
 template <typename DerivedC, typename DerivedBE, typename DerivedT>
@@ -161,5 +159,6 @@ IGL_INLINE void igl::draw_skeleton_3d(
 // Explicit template instanciation
 template void igl::draw_skeleton_3d<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&);
 template void igl::draw_skeleton_3d<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&);
-template void igl::draw_skeleton_3d<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -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<float, -1, -1, 0, -1, -1> > const&);
+template void igl::draw_skeleton_3d<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<float, 4, 1, 0, 4, 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<float, 4, 1, 0, 4, 1> > const&, double);
+template void igl::draw_skeleton_3d<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -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<float, -1, -1, 0, -1, -1> > const&, double);
 #endif

include/igl/draw_skeleton_3d.h

@@ -19,6 +19,7 @@ namespace igl
   //   BE  #BE by 2 list of bone edge indices into C
   //   T  #BE*(dim+1) by dim  matrix of stacked transposed bone transformations
   //   color  #BE|1 by 4 list of color
+  //   half_bbd  half bounding box diagonal to determine scaling {1.0}
   template <
     typename DerivedC, 
     typename DerivedBE, 
@@ -28,7 +29,8 @@ namespace igl
     const Eigen::PlainObjectBase<DerivedC> & C,
     const Eigen::PlainObjectBase<DerivedBE> & BE,
     const Eigen::PlainObjectBase<DerivedT> & T,
-    const Eigen::PlainObjectBase<Derivedcolor> & color);
+    const Eigen::PlainObjectBase<Derivedcolor> & color,
+    const double half_bbd=0.5);
   // Default color
   template <typename DerivedC, typename DerivedBE, typename DerivedT>
   IGL_INLINE void draw_skeleton_3d(

include/igl/embree/unproject_onto_mesh.cpp

@@ -19,7 +19,7 @@ IGL_INLINE bool igl::unproject_onto_mesh(
   const Eigen::Vector4f& viewport,
   const igl::EmbreeIntersector & ei,
   int& fid,
-  int& vid)
+  Eigen::Vector3f& bc)
 {
   using namespace std;
   using namespace Eigen;
@@ -31,18 +31,32 @@ IGL_INLINE bool igl::unproject_onto_mesh(
 
   if (hits.size()> 0)
   {
-    Vector3d bc(1.0-hits[0].u-hits[0].v, hits[0].u, hits[0].v);
-    int i;
-    bc.maxCoeff(&i);
-
+    bc << 1.0-hits[0].u-hits[0].v, hits[0].u, hits[0].v;
     fid = hits[0].id;
-    vid = F(fid,i);
     return true;
   }
 
   return false;
 }
 
+IGL_INLINE bool igl::unproject_onto_mesh(
+  const Eigen::Vector2f& pos,
+  const Eigen::MatrixXi& F,
+  const Eigen::Matrix4f& model,
+  const Eigen::Matrix4f& proj,
+  const Eigen::Vector4f& viewport,
+  const igl::EmbreeIntersector & ei,
+  int& fid,
+  int& vid)
+{
+  Eigen::Vector3f bc;
+  bool hit = unproject_onto_mesh(pos,F,model,proj,viewport,ei,fid,bc);
+  int i;
+  bc.maxCoeff(&i);
+  vid = F(fid,i);
+  return hit;
+}
+
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instanciation

include/igl/embree/unproject_onto_mesh.h

@@ -17,6 +17,30 @@ namespace igl
 {
   // Forward define
   class EmbreeIntersector;
+  // Unproject a screen location (using the given model, proj and viewport) to find
+  // the first hit on a mesh.
+  //
+  // Inputs:
+  //    pos        screen space coordinates
+  //    F          #F by 3 face matrix
+  //    model      model matrix
+  //    proj       projection matrix
+  //    viewport   vieweport vector
+  //    ei         EmbreeIntersector containing (V,F)
+  // Outputs:
+  //    fid        id of the first face hit
+  //    bc         barycentric coordinates of hit
+  // Returns true if there is a hit
+  IGL_INLINE bool unproject_onto_mesh(
+    const Eigen::Vector2f& pos,
+    const Eigen::MatrixXi& F,
+    const Eigen::Matrix4f& model,
+    const Eigen::Matrix4f& proj,
+    const Eigen::Vector4f& viewport,
+    const igl::EmbreeIntersector & ei,
+    int& fid,
+    Eigen::Vector3f& bc);
+  
   // Unproject a screen location (using the given model, proj and viewport) to find
   // the first face on the mesh and the closest vertex
   //
@@ -31,9 +55,6 @@ namespace igl
   //    fid        id of the first face hit
   //    vid        vertex id of the closest vertex hit
   // Returns true if there is a hit
-  //
-  // Known bugs: This should return the barycentric coordinates in F(fid,:)
-  // rather than vid.
   IGL_INLINE bool unproject_onto_mesh(
     const Eigen::Vector2f& pos,
     const Eigen::MatrixXi& F,
@@ -43,6 +64,7 @@ namespace igl
     const igl::EmbreeIntersector & ei,
     int& fid,
     int& vid);
+
 }
 #ifndef IGL_STATIC_LIBRARY
 #  include "unproject_onto_mesh.cpp"

include/igl/facet_components.cpp

@@ -13,7 +13,8 @@ IGL_INLINE void igl::facet_components(
   vector<vector<vector<Index > > > TT;
   vector<vector<vector<Index > > > TTi;
   triangle_triangle_adjacency(F,TT,TTi);
-  return facet_components(TT,C);
+  Eigen::VectorXi counts;
+  return facet_components(TT,C,counts);
 }
 
 template <
@@ -68,9 +69,9 @@ IGL_INLINE void igl::facet_components(
     }
     id++;
   }
-  assert(id == vcounts.size());
+  assert((size_t) id == vcounts.size());
   const size_t ncc = vcounts.size();
-  assert(C.maxCoeff()+1 == ncc);
+  assert((size_t)C.maxCoeff()+1 == ncc);
   counts.resize(ncc,1);
   for(size_t i = 0;i<ncc;i++)
   {

include/igl/fit_rigid.cpp

@@ -6,7 +6,7 @@
 // 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"
+#include "procrustes.h"
 
 IGL_INLINE void igl::fit_rigid(
   const Eigen::MatrixXd & A,
@@ -14,12 +14,8 @@ IGL_INLINE void igl::fit_rigid(
   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;
+  using namespace Eigen;
+  Matrix2d Rmat;
+  procrustes(A,B,Rmat,t);
+  R.fromRotationMatrix(Rmat);
 }

include/igl/fit_rigid.h

@@ -8,12 +8,14 @@
 #ifndef IGL_FIT_RIGID_H
 #define IGL_FIT_RIGID_H
 
+#warning "Deprecated. Use igl/procrustes.h instead"
 #include "igl_inline.h"
 #include <Eigen/Core>
 #include <Eigen/Geometry>
 
 namespace igl
 {
+  // Deprecated: please use procrustes(...) instead.
   // Fit a rigid 
   IGL_INLINE void fit_rigid(
     const Eigen::MatrixXd & A,

include/igl/frame_field_deformer.h

@@ -13,8 +13,8 @@
 
 namespace igl
 {
-  // Deform a mesh to transform the given per-face frame field to be as close as possible
-  // to a cross field, in the least square sense.
+  // Deform a mesh to transform the given per-face frame field to be as close
+  // as possible to a cross field, in the least square sense.
   //
   // Inputs:
   //   V       #V by 3 coordinates of the vertices

include/igl/get_modifiers.cpp

@@ -1,42 +0,0 @@
-#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

@@ -1,22 +1,69 @@
 #ifndef GET_MODIFIERS_H
 #define GET_MODIFIERS_H
-#include "igl_inline.h"
+//#include "igl_inline.h"
 namespace igl
 {
   enum Modifier
   {
     MODIFIER_OPTION = 1,
-    MODIFIER_SHIFT = 3,
-    MODIFIER_CONTROL = 5,
-    MODIFIER_COMMAND = 9,
+    MODIFIER_SHIFT = 2,
+    MODIFIER_CONTROL = 4,
+    MODIFIER_COMMAND = 8,
     NUM_MODIFIERS = 4,
   };
   // Retrieve current modifier constellation. 
   //
   // Returns int that's an "or" of the active modifiers above.
-  IGL_INLINE int get_modifiers();
+  //
+  // FORCED INLINE
+  inline int get_modifiers();
 }
-#ifndef IGL_STATIC_LIBRARY
-#include "get_modifiers.cpp"
+
+// Implementation 
+
+/* 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/HIToolbox/Events.h>
+#include <Carbon/Carbon.h>
 #endif
+
+#warning "igl::get_modifiers is deprecated. If using GLUT, try Alec's glut patch www.alecjacobson.com/weblog/?p=3659 and use glutGetModifiers"
+
+// FORCED INLINE
+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
+#  warning "igl::get_modifiers not supported on your OS, some demos may not work correctly."
+#endif
+  return mod;
+}
+
 #endif

include/igl/is_irregular_vertex.cpp

@@ -0,0 +1,43 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2015 Daniele Panozzo <daniele.panozzo@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "is_irregular_vertex.h"
+#include <vector>
+
+#include "is_border_vertex.h"
+
+template <typename DerivedV, typename DerivedF>
+IGL_INLINE std::vector<bool> igl::is_irregular_vertex(const Eigen::PlainObjectBase<DerivedV> &V, const Eigen::PlainObjectBase<DerivedF> &F)
+{
+  Eigen::VectorXi count = Eigen::VectorXi::Zero(F.maxCoeff());
+
+  for(unsigned i=0; i<F.rows();++i)
+  {
+    for(unsigned j=0; j<F.cols();++j)
+    {
+      if (F(i,j) < F(i,(j+1)%F.cols())) // avoid duplicate edges
+      {
+        count(F(i,j  )) += 1;
+        count(F(i,(j+1)%F.cols())) += 1;
+      }
+    }
+  }
+
+  std::vector<bool> border = is_border_vertex(V,F);
+
+  std::vector<bool> res(count.size());
+
+  for (unsigned i=0; i<res.size(); ++i)
+    res[i] = !border[i] && count[i] != (F.cols() == 3 ? 6 : 4 );
+
+  return res;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+template std::vector<bool, std::allocator<bool> > igl::is_irregular_vertex<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&);
+#endif

include/igl/is_irregular_vertex.h

@@ -0,0 +1,33 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2015 Daniele Panozzo <daniele.panozzo@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_IS_IRREGULAR_VERTEX_H
+#define IGL_IS_IRREGULAR_VERTEX_H
+#include "igl_inline.h"
+
+#include <Eigen/Core>
+#include <vector>
+
+namespace igl
+{
+  // Determine if a vertex is irregular, i.e. it has more than 6 (triangles)
+  // or 4 (quads) incident edges. Vertices on the boundary are ignored.
+  //
+  // Inputs:
+  //   V  #V by dim list of vertex positions
+  //   F  #F by 3[4] list of triangle[quads] indices
+  // Returns #V vector of bools revealing whether vertices are singular
+  //
+  template <typename DerivedV, typename DerivedF>
+  IGL_INLINE std::vector<bool> is_irregular_vertex(const Eigen::PlainObjectBase<DerivedV> &V, const Eigen::PlainObjectBase<DerivedF> &F);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "is_irregular_vertex.cpp"
+#endif
+
+#endif

include/igl/jet.cpp

@@ -56,7 +56,7 @@
 template <typename T>
 IGL_INLINE void igl::jet(const T x, T * rgb)
 {
-  igl::jet(x,rgb[0],rgb[1],rgb[2]);
+  return jet(x,rgb[0],rgb[1],rgb[2]);
 }
 
 template <typename T>

include/igl/list_to_matrix.cpp

@@ -152,4 +152,5 @@ template bool igl::list_to_matrix<int, Eigen::PlainObjectBase<Eigen::Matrix<int,
 template bool igl::list_to_matrix<double, Eigen::Matrix<double, 1, -1, 1, 1, -1> >(std::vector<double, std::allocator<double> > const&, Eigen::Matrix<double, 1, -1, 1, 1, -1>&);
 template bool igl::list_to_matrix<unsigned long, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > >(std::vector<unsigned long, std::allocator<unsigned long> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template bool igl::list_to_matrix<double, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(std::__1::vector<double, std::__1::allocator<double> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1>&);
+template bool igl::list_to_matrix<float, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > >(std::vector<std::vector<float, std::allocator<float> >, std::allocator<std::vector<float, std::allocator<float> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/map_vertices_to_circle.cpp

@@ -8,8 +8,6 @@
 
 #include "map_vertices_to_circle.h"
 
-#include "boundary_loop.h"
-
 IGL_INLINE void igl::map_vertices_to_circle(
   const Eigen::MatrixXd& V,
   const Eigen::VectorXi& bnd,

include/igl/massmatrix.cpp

@@ -8,6 +8,7 @@
 #include "massmatrix.h"
 #include "normalize_row_sums.h"
 #include "sparse.h"
+#include "doublearea.h"
 #include "repmat.h"
 #include <Eigen/Geometry>
 #include <iostream>
@@ -51,15 +52,8 @@ IGL_INLINE void igl::massmatrix(
       l(i,1) = (V.row(F(i,2))-V.row(F(i,0))).norm();
       l(i,2) = (V.row(F(i,0))-V.row(F(i,1))).norm();
     }
-    // semiperimeters
-    Matrix<Scalar,Dynamic,1> s = l.rowwise().sum()*0.5;
-    assert(s.rows() == m);
-    // Heron's forumal for area
-    Matrix<Scalar,Dynamic,1> dblA(m);
-    for(int i = 0;i<m;i++)
-    {
-      dblA(i) = 2.0*sqrt(s(i)*(s(i)-l(i,0))*(s(i)-l(i,1))*(s(i)-l(i,2)));
-    }
+    Matrix<Scalar,Dynamic,1> dblA;
+    doublearea(l,dblA);
 
     switch(eff_type)
     {

include/igl/max_size.cpp

@@ -31,4 +31,5 @@ template int igl::max_size<std::vector<int, std::allocator<int> > >(std::vector<
 // generated by autoexplicit.sh
 template int igl::max_size<std::vector<double, std::allocator<double> > >(std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > > const&);
 template int igl::max_size<std::vector<bool, std::allocator<bool> > >(std::vector<std::vector<bool, std::allocator<bool> >, std::allocator<std::vector<bool, std::allocator<bool> > > > const&);
+template int igl::max_size<std::vector<float, std::allocator<float> > >(std::vector<std::vector<float, std::allocator<float> >, std::allocator<std::vector<float, std::allocator<float> > > > const&);
 #endif

include/igl/min_quad_with_fixed.cpp

@@ -96,6 +96,7 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
 
   SparseMatrix<T> Auu;
   slice(A,data.unknown,data.unknown,Auu);
+  assert(Auu.size() > 0 && "All DOFs seem to be fixed.");
 
   // Positive definiteness is *not* determined, rather it is given as a
   // parameter

include/igl/min_size.cpp

@@ -37,4 +37,5 @@ template int igl::min_size<std::vector<int, std::allocator<int> > >(std::vector<
 // generated by autoexplicit.sh
 template int igl::min_size<std::vector<double, std::allocator<double> > >(std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > > const&);
 template int igl::min_size<std::vector<bool, std::allocator<bool> > >(std::vector<std::vector<bool, std::allocator<bool> >, std::allocator<std::vector<bool, std::allocator<bool> > > > const&);
+template int igl::min_size<std::vector<float, std::allocator<float> > >(std::vector<std::vector<float, std::allocator<float> >, std::allocator<std::vector<float, std::allocator<float> > > > const&);
 #endif

include/igl/outer_hull.cpp

@@ -15,6 +15,7 @@
 #include <map>
 #include <queue>
 #include <iostream>
+//#define IGL_OUTER_HULL_DEBUG
 
 template <
   typename DerivedV,
@@ -40,6 +41,13 @@ IGL_INLINE void igl::outer_hull(
   typedef Matrix<typename DerivedJ::Scalar,Dynamic,DerivedJ::ColsAtCompileTime> MatrixXJ;
   const Index m = F.rows();
 
+#ifdef IGL_OUTER_HULL_DEBUG
+  cout<<"outer hull..."<<endl;
+#endif
+
+#ifdef IGL_OUTER_HULL_DEBUG
+  cout<<"edge map..."<<endl;
+#endif
   typedef Matrix<typename DerivedF::Scalar,Dynamic,2> MatrixX2I;
   typedef Matrix<typename DerivedF::Index,Dynamic,1> VectorXI;
   MatrixX2I E,uE;
@@ -50,6 +58,9 @@ IGL_INLINE void igl::outer_hull(
   vector<vector<vector<Index > > > TT,_1;
   triangle_triangle_adjacency(E,EMAP,uE2E,false,TT,_1);
   VectorXI counts;
+#ifdef IGL_OUTER_HULL_DEBUG
+  cout<<"facet components..."<<endl;
+#endif
   facet_components(TT,C,counts);
   assert(C.maxCoeff()+1 == counts.rows());
   const size_t ncc = counts.rows();
@@ -61,7 +72,14 @@ IGL_INLINE void igl::outer_hull(
     typename DerivedV::Scalar,
     DerivedF::RowsAtCompileTime,
     3> N;
+#ifdef IGL_OUTER_HULL_DEBUG
+  cout<<"normals..."<<endl;
+#endif
   per_face_normals(V,F,N);
+
+#ifdef IGL_OUTER_HULL_DEBUG
+  cout<<"reindex..."<<endl;
+#endif
   // H contains list of faces on outer hull;
   vector<bool> FH(m,false);
   vector<bool> EH(3*m,false);
@@ -80,11 +98,17 @@ IGL_INLINE void igl::outer_hull(
     vIM[C(f)](g[C(f)]++) = f;
   }
 
+#ifdef IGL_OUTER_HULL_DEBUG
+  cout<<"barycenters..."<<endl;
+#endif
   // assumes that "resolve" has handled any coplanar cases correctly and nearly
   // coplanar cases can be sorted based on barycenter.
   MatrixXV BC;
   barycenter(V,F,BC);
 
+#ifdef IGL_OUTER_HULL_DEBUG
+  cout<<"loop over CCs..."<<endl;
+#endif
   for(Index id = 0;id<(Index)ncc;id++)
   {
     auto & IM = vIM[id];
@@ -92,6 +116,9 @@ IGL_INLINE void igl::outer_hull(
     // component
     int f;
     bool f_flip;
+#ifdef IGL_OUTER_HULL_DEBUG
+  cout<<"outer facet..."<<endl;
+#endif
     outer_facet(V,F,N,IM,f,f_flip);
     int FHcount = 0;
     // Q contains list of face edges to continue traversing upong
@@ -99,7 +126,10 @@ IGL_INLINE void igl::outer_hull(
     Q.push(f+0*m);
     Q.push(f+1*m);
     Q.push(f+2*m);
-    flip[f] = f_flip;
+    flip(f) = f_flip;
+#ifdef IGL_OUTER_HULL_DEBUG
+  cout<<"BFS..."<<endl;
+#endif
     while(!Q.empty())
     {
       // face-edge
@@ -123,11 +153,11 @@ IGL_INLINE void igl::outer_hull(
       }
       // find overlapping face-edges
       const auto & neighbors = uE2E[EMAP(e)];
-      const auto & fN = (flip[f]?-1.:1.)*N.row(f);
+      const auto & fN = (flip(f)?-1.:1.)*N.row(f);
       // source of edge according to f
-      const int fs = flip[f]?F(f,(c+2)%3):F(f,(c+1)%3);
+      const int fs = flip(f)?F(f,(c+2)%3):F(f,(c+1)%3);
       // destination of edge according to f
-      const int fd = flip[f]?F(f,(c+1)%3):F(f,(c+2)%3);
+      const int fd = flip(f)?F(f,(c+1)%3):F(f,(c+2)%3);
       const auto & eV = (V.row(fd)-V.row(fs)).normalized();
       // Loop over and find max dihedral angle
       typename DerivedV::Scalar max_di = -1;
@@ -144,8 +174,8 @@ IGL_INLINE void igl::outer_hull(
         // are faces consistently oriented
         //const int ns = F(nf,(nc+1)%3);
         const int nd = F(nf,(nc+2)%3);
-        const bool cons = (flip[f]?fd:fs) == nd;
-        const auto & nN = (cons? (flip[f]?-1:1.) : (flip[f]?1.:-1.) )*N.row(nf);
+        const bool cons = (flip(f)?fd:fs) == nd;
+        const auto & nN = (cons? (flip(f)?-1:1.) : (flip(f)?1.:-1.) )*N.row(nf);
         const auto & ndi = M_PI - atan2( fN.cross(nN).dot(eV), fN.dot(nN));
         if(ndi>=max_di)
         {
@@ -159,8 +189,8 @@ IGL_INLINE void igl::outer_hull(
         const int nf = max_ne%m;
         const int nc = max_ne/m;
         const int nd = F(nf,(nc+2)%3);
-        const bool cons = (flip[f]?fd:fs) == nd;
-        flip[nf] = (cons ? flip[f] : !flip[f]);
+        const bool cons = (flip(f)?fd:fs) == nd;
+        flip(nf) = (cons ? flip(f) : !flip(f));
         const int ne1 = nf+((nc+1)%3)*m;
         const int ne2 = nf+((nc+2)%3)*m;
         if(!EH[ne1])
@@ -185,16 +215,16 @@ IGL_INLINE void igl::outer_hull(
         const size_t f = IM(i);
         //if(f_flip)
         //{
-        //  flip[f] = !flip[f];
+        //  flip(f) = !flip(f);
         //}
         if(FH[f])
         {
-          vG[id].row(h) = (flip[f]?F.row(f).reverse().eval():F.row(f));
+          vG[id].row(h) = (flip(f)?F.row(f).reverse().eval():F.row(f));
           vJ[id](h,0) = f;
           h++;
         }
       }
-      assert(h == FHcount);
+      assert((int)h == FHcount);
     }
   }
 
@@ -245,7 +275,7 @@ IGL_INLINE void igl::outer_hull(
     // q could be so close (<~1e-16) to B that the winding number is not a robust way to
     // determine inside/outsideness. We could try to find a _better_ q which is
     // farther away, but couldn't they all be bad?
-    MatrixXV q = BC.row(AJ(1));
+    MatrixXV q = BC.row(AJ(0));
     // In a perfect world, it's enough to test a single point.
     double w;
     winding_number_3(

include/igl/outline_ordered.cpp

@@ -1,62 +0,0 @@
-// This file is part of libigl, a simple c++ geometry processing library.
-// 
-// Copyright (C) 2014 Stefan Brugger <stefanbrugger@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 "outline_ordered.h"
-#include "exterior_edges.h"
-#include <set>
-
-template <typename DerivedF, typename Index>
-IGL_INLINE void igl::outline_ordered(
-    const Eigen::PlainObjectBase<DerivedF> & F, 
-    std::vector<std::vector<Index> >& L)
-{
-  using namespace std;
-  using namespace Eigen;
-  // Exterior edges include some non-manifold edges (poor function name). I
-  // suppose `outline_ordered` is not well defined for non-manifold meshes, but
-  // perhaps this should just call `boundary_facets`
-  MatrixXi E = exterior_edges(F);
-
-  set<int> unseen;
-  for (int i = 0; i < E.rows(); ++i)
-  {
-    unseen.insert(unseen.end(),i);
-  }
-
-  while (!unseen.empty())
-  {
-      vector<Index> l;
-
-      // Get first vertex of loop
-      int startEdge = *unseen.begin();
-      unseen.erase(unseen.begin());
-
-      int start = E(startEdge,0);
-      int next = E(startEdge,1);
-      l.push_back(start);
-
-      while (start != next)
-      {
-          l.push_back(next);
-
-          // Find next edge
-          int nextEdge;
-          set<int>::iterator it;
-          for (it=unseen.begin(); it != unseen.end() ; ++it)
-          {
-              if (E(*it,0) == next || E(*it,1) == next)
-              {
-                  nextEdge = *it;
-                  break;
-              }                  
-          }
-          unseen.erase(nextEdge);
-          next = (E(nextEdge,0) == next) ? E(nextEdge,1) : E(nextEdge,0);
-      }
-      L.push_back(l);
-  }
-}

include/igl/outline_ordered.h

@@ -1,35 +0,0 @@
-// This file is part of libigl, a simple c++ geometry processing library.
-// 
-// Copyright (C) 2014 Stefan Brugger <stefanbrugger@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_OUTLINE_ORDERED_H
-#define IGL_OUTLINE_ORDERED_H
-#include "igl_inline.h"
-
-#include <Eigen/Dense>
-#include <vector>
-
-namespace igl
-{
-  // Compute list of ordered boundary loops for a manifold mesh.
-  //
-  // Templates:
-  //  Index  index type
-  // Inputs:
-  //   F  #V by dim list of mesh faces
-  // Outputs:
-  //   L  list of loops where L[i] = ordered list of boundary vertices in loop i
-  //
-  template <typename DerivedF, typename Index>
-  IGL_INLINE void outline_ordered(
-    const Eigen::PlainObjectBase<DerivedF> & F, 
-    std::vector<std::vector<Index> >& L);
-}
-
-#ifndef IGL_STATIC_LIBRARY
-#  include "outline_ordered.cpp"
-#endif
-#endif

include/igl/parula.cpp

@@ -0,0 +1,59 @@
+#include "parula.h"
+
+
+template <typename T>
+IGL_INLINE void igl::parula(const T x, T * rgb)
+{
+  return parula(x,rgb[0],rgb[1],rgb[2]);
+}
+
+template <typename T>
+IGL_INLINE void igl::parula(const T f, T & r, T & g, T & b)
+{
+  // clamp to (0,1)
+  const float eff_f = (f>=1.?1.:(f<=0.?0.:f));
+  // continuous index into array
+  const float ff = eff_f*(PARULA_COLOR_MAP.rows()-1);
+  size_t s = floor(ff);
+  size_t d = ceil(ff);
+  const float t = (s==d ? 0. : (ff-s)/float(d-s));
+
+  assert(t>=0 && t<=1);
+  const auto rgb_s = PARULA_COLOR_MAP.row(s);
+  const auto rgb_d = PARULA_COLOR_MAP.row(d);
+  const auto rgb = rgb_s + t*(rgb_d-rgb_s);
+  r = rgb(0);
+  g = rgb(1);
+  b = rgb(2);
+}
+
+
+template <typename DerivedZ, typename DerivedC>
+IGL_INLINE void igl::parula(
+  const Eigen::PlainObjectBase<DerivedZ> & Z,
+  const bool normalize,
+  Eigen::PlainObjectBase<DerivedC> & C)
+{
+  const double min_z = (normalize?Z.minCoeff():0);
+  const double max_z = (normalize?Z.maxCoeff():1);
+  return parula(Z,min_z,max_z,C);
+}
+template <typename DerivedZ, typename DerivedC>
+IGL_INLINE void igl::parula(
+  const Eigen::PlainObjectBase<DerivedZ> & Z,
+  const double min_z,
+  const double max_z,
+  Eigen::PlainObjectBase<DerivedC> & C)
+{
+  C.resize(Z.rows(),3);
+  for(int r = 0;r<Z.rows();r++)
+  {
+    parula((-min_z+Z(r,0))/(max_z-min_z),C(r,0),C(r,1),C(r,2));
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instanciation
+template void igl::parula<double>(double, double*);
+template void igl::parula<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+#endif

include/igl/parula.h

@@ -0,0 +1,320 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 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_PARULA_H
+#define IGL_PARULA_H
+#include "igl_inline.h"
+//#ifndef IGL_NO_EIGEN
+#  include <Eigen/Dense>
+//#endif
+namespace igl
+{
+  // PARULA like MATLAB's parula
+  //
+  // Inputs:
+  //   m  number of colors 
+  // Outputs:
+  //   J  m by list of RGB colors between 0 and 1
+  //
+  // Wrapper for directly computing [r,g,b] values for a given factor f between
+  // 0 and 1
+  //
+  // Inputs:
+  //   f  factor determining color value as if 0 was min and 1 was max
+  // Outputs:
+  //   r  red value
+  //   g  green value
+  //   b  blue value
+  template <typename T>
+  IGL_INLINE void parula(const T f, T * rgb);
+  template <typename T>
+  IGL_INLINE void parula(const T f, T & r, T & g, T & b);
+  // Inputs:
+  //   Z  #Z list of factors 
+  //   normalize  whether to normalize Z to be tightly between [0,1]
+  // Outputs:
+  //   C  #C by 3 list of rgb colors
+  template <typename DerivedZ, typename DerivedC>
+  IGL_INLINE void parula(
+    const Eigen::PlainObjectBase<DerivedZ> & Z,
+    const bool normalize,
+    Eigen::PlainObjectBase<DerivedC> & C);
+  // Inputs:
+  //   min_z  value at blue
+  //   max_z  value at red
+  template <typename DerivedZ, typename DerivedC>
+  IGL_INLINE void parula(
+    const Eigen::PlainObjectBase<DerivedZ> & Z,
+    const double min_Z,
+    const double max_Z,
+    Eigen::PlainObjectBase<DerivedC> & C);
+  const Eigen::Matrix<float,256,4> PARULA_COLOR_MAP = 
+    (Eigen::Matrix<float,256,4>()<<
+      0.2081,0.1663,0.5292,1,
+      0.2091,0.1721,0.5411,1,
+      0.2101,0.1779,0.553,1,
+      0.2109,0.1837,0.565,1,
+      0.2116,0.1895,0.5771,1,
+      0.2121,0.1954,0.5892,1,
+      0.2124,0.2013,0.6013,1,
+      0.2125,0.2072,0.6135,1,
+      0.2123,0.2132,0.6258,1,
+      0.2118,0.2192,0.6381,1,
+      0.2111,0.2253,0.6505,1,
+      0.2099,0.2315,0.6629,1,
+      0.2084,0.2377,0.6753,1,
+      0.2063,0.244,0.6878,1,
+      0.2038,0.2503,0.7003,1,
+      0.2006,0.2568,0.7129,1,
+      0.1968,0.2632,0.7255,1,
+      0.1921,0.2698,0.7381,1,
+      0.1867,0.2764,0.7507,1,
+      0.1802,0.2832,0.7634,1,
+      0.1728,0.2902,0.7762,1,
+      0.1641,0.2975,0.789,1,
+      0.1541,0.3052,0.8017,1,
+      0.1427,0.3132,0.8145,1,
+      0.1295,0.3217,0.8269,1,
+      0.1147,0.3306,0.8387,1,
+      0.0986,0.3397,0.8495,1,
+      0.0816,0.3486,0.8588,1,
+      0.0646,0.3572,0.8664,1,
+      0.0482,0.3651,0.8722,1,
+      0.0329,0.3724,0.8765,1,
+      0.0213,0.3792,0.8796,1,
+      0.0136,0.3853,0.8815,1,
+      0.0086,0.3911,0.8827,1,
+      0.006,0.3965,0.8833,1,
+      0.0051,0.4017,0.8834,1,
+      0.0054,0.4066,0.8831,1,
+      0.0067,0.4113,0.8825,1,
+      0.0089,0.4159,0.8816,1,
+      0.0116,0.4203,0.8805,1,
+      0.0148,0.4246,0.8793,1,
+      0.0184,0.4288,0.8779,1,
+      0.0223,0.4329,0.8763,1,
+      0.0264,0.437,0.8747,1,
+      0.0306,0.441,0.8729,1,
+      0.0349,0.4449,0.8711,1,
+      0.0394,0.4488,0.8692,1,
+      0.0437,0.4526,0.8672,1,
+      0.0477,0.4564,0.8652,1,
+      0.0514,0.4602,0.8632,1,
+      0.0549,0.464,0.8611,1,
+      0.0582,0.4677,0.8589,1,
+      0.0612,0.4714,0.8568,1,
+      0.064,0.4751,0.8546,1,
+      0.0666,0.4788,0.8525,1,
+      0.0689,0.4825,0.8503,1,
+      0.071,0.4862,0.8481,1,
+      0.0729,0.4899,0.846,1,
+      0.0746,0.4937,0.8439,1,
+      0.0761,0.4974,0.8418,1,
+      0.0773,0.5012,0.8398,1,
+      0.0782,0.5051,0.8378,1,
+      0.0789,0.5089,0.8359,1,
+      0.0794,0.5129,0.8341,1,
+      0.0795,0.5169,0.8324,1,
+      0.0793,0.521,0.8308,1,
+      0.0788,0.5251,0.8293,1,
+      0.0778,0.5295,0.828,1,
+      0.0764,0.5339,0.827,1,
+      0.0746,0.5384,0.8261,1,
+      0.0724,0.5431,0.8253,1,
+      0.0698,0.5479,0.8247,1,
+      0.0668,0.5527,0.8243,1,
+      0.0636,0.5577,0.8239,1,
+      0.06,0.5627,0.8237,1,
+      0.0562,0.5677,0.8234,1,
+      0.0523,0.5727,0.8231,1,
+      0.0484,0.5777,0.8228,1,
+      0.0445,0.5826,0.8223,1,
+      0.0408,0.5874,0.8217,1,
+      0.0372,0.5922,0.8209,1,
+      0.0342,0.5968,0.8198,1,
+      0.0317,0.6012,0.8186,1,
+      0.0296,0.6055,0.8171,1,
+      0.0279,0.6097,0.8154,1,
+      0.0265,0.6137,0.8135,1,
+      0.0255,0.6176,0.8114,1,
+      0.0248,0.6214,0.8091,1,
+      0.0243,0.625,0.8066,1,
+      0.0239,0.6285,0.8039,1,
+      0.0237,0.6319,0.801,1,
+      0.0235,0.6352,0.798,1,
+      0.0233,0.6384,0.7948,1,
+      0.0231,0.6415,0.7916,1,
+      0.023,0.6445,0.7881,1,
+      0.0229,0.6474,0.7846,1,
+      0.0227,0.6503,0.781,1,
+      0.0227,0.6531,0.7773,1,
+      0.0232,0.6558,0.7735,1,
+      0.0238,0.6585,0.7696,1,
+      0.0246,0.6611,0.7656,1,
+      0.0263,0.6637,0.7615,1,
+      0.0282,0.6663,0.7574,1,
+      0.0306,0.6688,0.7532,1,
+      0.0338,0.6712,0.749,1,
+      0.0373,0.6737,0.7446,1,
+      0.0418,0.6761,0.7402,1,
+      0.0467,0.6784,0.7358,1,
+      0.0516,0.6808,0.7313,1,
+      0.0574,0.6831,0.7267,1,
+      0.0629,0.6854,0.7221,1,
+      0.0692,0.6877,0.7173,1,
+      0.0755,0.6899,0.7126,1,
+      0.082,0.6921,0.7078,1,
+      0.0889,0.6943,0.7029,1,
+      0.0956,0.6965,0.6979,1,
+      0.1031,0.6986,0.6929,1,
+      0.1104,0.7007,0.6878,1,
+      0.118,0.7028,0.6827,1,
+      0.1258,0.7049,0.6775,1,
+      0.1335,0.7069,0.6723,1,
+      0.1418,0.7089,0.6669,1,
+      0.1499,0.7109,0.6616,1,
+      0.1585,0.7129,0.6561,1,
+      0.1671,0.7148,0.6507,1,
+      0.1758,0.7168,0.6451,1,
+      0.1849,0.7186,0.6395,1,
+      0.1938,0.7205,0.6338,1,
+      0.2033,0.7223,0.6281,1,
+      0.2128,0.7241,0.6223,1,
+      0.2224,0.7259,0.6165,1,
+      0.2324,0.7275,0.6107,1,
+      0.2423,0.7292,0.6048,1,
+      0.2527,0.7308,0.5988,1,
+      0.2631,0.7324,0.5929,1,
+      0.2735,0.7339,0.5869,1,
+      0.2845,0.7354,0.5809,1,
+      0.2953,0.7368,0.5749,1,
+      0.3064,0.7381,0.5689,1,
+      0.3177,0.7394,0.563,1,
+      0.3289,0.7406,0.557,1,
+      0.3405,0.7417,0.5512,1,
+      0.352,0.7428,0.5453,1,
+      0.3635,0.7438,0.5396,1,
+      0.3753,0.7446,0.5339,1,
+      0.3869,0.7454,0.5283,1,
+      0.3986,0.7461,0.5229,1,
+      0.4103,0.7467,0.5175,1,
+      0.4218,0.7473,0.5123,1,
+      0.4334,0.7477,0.5072,1,
+      0.4447,0.7482,0.5021,1,
+      0.4561,0.7485,0.4972,1,
+      0.4672,0.7487,0.4924,1,
+      0.4783,0.7489,0.4877,1,
+      0.4892,0.7491,0.4831,1,
+      0.5,0.7491,0.4786,1,
+      0.5106,0.7492,0.4741,1,
+      0.5212,0.7492,0.4698,1,
+      0.5315,0.7491,0.4655,1,
+      0.5418,0.749,0.4613,1,
+      0.5519,0.7489,0.4571,1,
+      0.5619,0.7487,0.4531,1,
+      0.5718,0.7485,0.449,1,
+      0.5816,0.7482,0.4451,1,
+      0.5913,0.7479,0.4412,1,
+      0.6009,0.7476,0.4374,1,
+      0.6103,0.7473,0.4335,1,
+      0.6197,0.7469,0.4298,1,
+      0.629,0.7465,0.4261,1,
+      0.6382,0.746,0.4224,1,
+      0.6473,0.7456,0.4188,1,
+      0.6564,0.7451,0.4152,1,
+      0.6653,0.7446,0.4116,1,
+      0.6742,0.7441,0.4081,1,
+      0.683,0.7435,0.4046,1,
+      0.6918,0.743,0.4011,1,
+      0.7004,0.7424,0.3976,1,
+      0.7091,0.7418,0.3942,1,
+      0.7176,0.7412,0.3908,1,
+      0.7261,0.7405,0.3874,1,
+      0.7346,0.7399,0.384,1,
+      0.743,0.7392,0.3806,1,
+      0.7513,0.7385,0.3773,1,
+      0.7596,0.7378,0.3739,1,
+      0.7679,0.7372,0.3706,1,
+      0.7761,0.7364,0.3673,1,
+      0.7843,0.7357,0.3639,1,
+      0.7924,0.735,0.3606,1,
+      0.8005,0.7343,0.3573,1,
+      0.8085,0.7336,0.3539,1,
+      0.8166,0.7329,0.3506,1,
+      0.8246,0.7322,0.3472,1,
+      0.8325,0.7315,0.3438,1,
+      0.8405,0.7308,0.3404,1,
+      0.8484,0.7301,0.337,1,
+      0.8563,0.7294,0.3336,1,
+      0.8642,0.7288,0.33,1,
+      0.872,0.7282,0.3265,1,
+      0.8798,0.7276,0.3229,1,
+      0.8877,0.7271,0.3193,1,
+      0.8954,0.7266,0.3156,1,
+      0.9032,0.7262,0.3117,1,
+      0.911,0.7259,0.3078,1,
+      0.9187,0.7256,0.3038,1,
+      0.9264,0.7256,0.2996,1,
+      0.9341,0.7256,0.2953,1,
+      0.9417,0.7259,0.2907,1,
+      0.9493,0.7264,0.2859,1,
+      0.9567,0.7273,0.2808,1,
+      0.9639,0.7285,0.2754,1,
+      0.9708,0.7303,0.2696,1,
+      0.9773,0.7326,0.2634,1,
+      0.9831,0.7355,0.257,1,
+      0.9882,0.739,0.2504,1,
+      0.9922,0.7431,0.2437,1,
+      0.9952,0.7476,0.2373,1,
+      0.9973,0.7524,0.231,1,
+      0.9986,0.7573,0.2251,1,
+      0.9991,0.7624,0.2195,1,
+      0.999,0.7675,0.2141,1,
+      0.9985,0.7726,0.209,1,
+      0.9976,0.7778,0.2042,1,
+      0.9964,0.7829,0.1995,1,
+      0.995,0.788,0.1949,1,
+      0.9933,0.7931,0.1905,1,
+      0.9914,0.7981,0.1863,1,
+      0.9894,0.8032,0.1821,1,
+      0.9873,0.8083,0.178,1,
+      0.9851,0.8133,0.174,1,
+      0.9828,0.8184,0.17,1,
+      0.9805,0.8235,0.1661,1,
+      0.9782,0.8286,0.1622,1,
+      0.9759,0.8337,0.1583,1,
+      0.9736,0.8389,0.1544,1,
+      0.9713,0.8441,0.1505,1,
+      0.9692,0.8494,0.1465,1,
+      0.9672,0.8548,0.1425,1,
+      0.9654,0.8603,0.1385,1,
+      0.9638,0.8659,0.1343,1,
+      0.9623,0.8716,0.1301,1,
+      0.9611,0.8774,0.1258,1,
+      0.96,0.8834,0.1215,1,
+      0.9593,0.8895,0.1171,1,
+      0.9588,0.8958,0.1126,1,
+      0.9586,0.9022,0.1082,1,
+      0.9587,0.9088,0.1036,1,
+      0.9591,0.9155,0.099,1,
+      0.9599,0.9225,0.0944,1,
+      0.961,0.9296,0.0897,1,
+      0.9624,0.9368,0.085,1,
+      0.9641,0.9443,0.0802,1,
+      0.9662,0.9518,0.0753,1,
+      0.9685,0.9595,0.0703,1,
+      0.971,0.9673,0.0651,1,
+      0.9736,0.9752,0.0597,1,
+      0.9763,0.9831,0.0538,1).finished();
+};
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "parula.cpp"
+#endif
+
+#endif
+

include/igl/peal_outer_hull_layers.cpp → include/igl/peel_outer_hull_layers.cpp

@@ -1,7 +1,12 @@
-#include <vector>
-#include "peal_outer_hull_layers.h"
+#include "peel_outer_hull_layers.h"
 #include "outer_hull.h"
+#include "writePLY.h"
 #include <vector>
+#include <iostream>
+//#define IGL_PEEL_OUTER_HULL_LAYERS_DEBUG
+#ifdef IGL_PEEL_OUTER_HULL_LAYERS_DEBUG
+#include "STR.h"
+#endif
 
 using namespace std;
 template <
@@ -9,7 +14,7 @@ template <
   typename DerivedF,
   typename Derivedodd,
   typename Derivedflip>
-IGL_INLINE void igl::peal_outer_hull_layers(
+IGL_INLINE size_t igl::peel_outer_hull_layers(
   const Eigen::PlainObjectBase<DerivedV > & V,
   const Eigen::PlainObjectBase<DerivedF > & F,
   Eigen::PlainObjectBase<Derivedodd > & odd,
@@ -22,7 +27,13 @@ IGL_INLINE void igl::peal_outer_hull_layers(
   typedef Matrix<Index,Dynamic,1> MatrixXI;
   typedef Matrix<typename Derivedflip::Scalar,Dynamic,Derivedflip::ColsAtCompileTime> MatrixXflip;
   const Index m = F.rows();
+#ifdef IGL_PEEL_OUTER_HULL_LAYERS_DEBUG
+  cout<<"peel outer hull layers..."<<endl;
+#endif
 
+#ifdef IGL_PEEL_OUTER_HULL_LAYERS_DEBUG
+  cout<<"resize output ..."<<endl;
+#endif
   // keep track of iteration parity and whether flipped in hull
   MatrixXF Fr = F;
   odd.resize(m,1);
@@ -40,7 +51,15 @@ IGL_INLINE void igl::peal_outer_hull_layers(
     MatrixXF Fo;
     MatrixXI Jo;
     MatrixXflip flipr;
+#ifdef IGL_PEEL_OUTER_HULL_LAYERS_DEBUG
+  cout<<"calling outer hull..."<<endl;
+  writePLY(STR("outer-hull-input-"<<iter<<".ply"),V,Fr);
+#endif
     outer_hull(V,Fr,Fo,Jo,flipr);
+#ifdef IGL_PEEL_OUTER_HULL_LAYERS_DEBUG
+  writePLY(STR("outer-hull-output-"<<iter<<".ply"),V,Fo);
+  cout<<"reindex, flip..."<<endl;
+#endif
     assert(Fo.rows() == Jo.rows());
     // all faces in Fo of Fr
     vector<bool> in_outer(Fr.rows(),false);
@@ -72,10 +91,11 @@ IGL_INLINE void igl::peal_outer_hull_layers(
     odd_iter = !odd_iter;
     iter++;
   }
+  return iter;
 }
 
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template void igl::peal_outer_hull_layers<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<bool, -1, 1, 0, -1, 1>, Eigen::Matrix<bool, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&);
+template size_t igl::peel_outer_hull_layers<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<bool, -1, 1, 0, -1, 1>, Eigen::Matrix<bool, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/peal_outer_hull_layers.h → include/igl/peel_outer_hull_layers.h

@@ -1,27 +1,28 @@
-#ifndef PEAL_OUTER_HULL_LAYERS_H
-#define PEAL_OUTER_HULL_LAYERS_H
+#ifndef PEEL_OUTER_HULL_LAYERS_H
+#define PEEL_OUTER_HULL_LAYERS_H
 #include "igl_inline.h"
 #include <Eigen/Core>
 namespace igl
 {
   // Computes necessary generic information for boolean operations by
-  // successively "pealing" off the "outer hull" of a mesh (V,F) resulting from
+  // successively "peeling" off the "outer hull" of a mesh (V,F) resulting from
   // "resolving" all (self-)intersections.
   //
   // Inputs:
   //   V  #V by 3 list of vertex positions
   //   F  #F by 3 list of triangle indices into V
   // Outputs:
-  //   odd  #F list of whether facet belongs to an odd iteration peal (otherwise
-  //     an even iteration peal)
+  //   odd  #F list of whether facet belongs to an odd iteration peel (otherwise
+  //     an even iteration peel)
   //   flip  #F list of whether a facet's orientation was flipped when facet
-  //     "pealed" into its associated outer hull layer.
+  //     "peeled" into its associated outer hull layer.
+  // Returns number of peels
   template <
     typename DerivedV,
     typename DerivedF,
     typename Derivedodd,
     typename Derivedflip>
-  IGL_INLINE void peal_outer_hull_layers(
+  IGL_INLINE size_t peel_outer_hull_layers(
     const Eigen::PlainObjectBase<DerivedV > & V,
     const Eigen::PlainObjectBase<DerivedF > & F,
     Eigen::PlainObjectBase<Derivedodd > & odd,
@@ -29,6 +30,6 @@ namespace igl
 }
 
 #ifndef IGL_STATIC_LIBRARY
-#  include "peal_outer_hull_layers.cpp"
+#  include "peel_outer_hull_layers.cpp"
 #endif
 #endif

include/igl/ply.h.REMOVED.git-id

@@ -1 +1 @@
-05e282fe767ea0bd8987f68ae4b79d0ab4c26a09
+6693ed81cd72749c90548f5078b27623d108a17f

include/igl/procrustes.cpp

@@ -0,0 +1,136 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 Stefan Brugger <stefanbrugger@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 "procrustes.h"
+#include "polar_svd.h"
+#include "polar_dec.h"
+ 
+template <
+  typename DerivedX, 
+  typename DerivedY, 
+  typename Scalar, 
+  typename DerivedR, 
+  typename DerivedT>
+IGL_INLINE void igl::procrustes(
+    const Eigen::PlainObjectBase<DerivedX>& X,
+    const Eigen::PlainObjectBase<DerivedY>& Y,
+    bool includeScaling,
+    bool includeReflections,
+    Scalar& scale,
+    Eigen::PlainObjectBase<DerivedR>& R,
+    Eigen::PlainObjectBase<DerivedT>& t)
+{
+  using namespace Eigen;
+  assert (X.rows() == Y.rows() && "Same number of points");
+  assert(X.cols() == Y.cols() && "Points have same dimensions");
+
+  // Center data
+  const VectorXd Xmean = X.colwise().mean();      
+  const VectorXd Ymean = Y.colwise().mean();      
+  MatrixXd XC = X.rowwise() - Xmean.transpose();
+  MatrixXd YC = Y.rowwise() - Ymean.transpose();
+
+  // Scale
+  scale = 1.;
+  if (includeScaling)
+  {
+     double scaleX = XC.norm() / XC.rows();
+     double scaleY = YC.norm() / YC.rows();
+     scale = scaleY/scaleX;
+     XC *= scale;
+     assert (abs(XC.norm() / XC.rows() - scaleY) < 1e-8);
+  }
+
+  // Rotation 
+  MatrixXd S = XC.transpose() * YC; 
+  MatrixXd T;
+  if (includeReflections)
+  {
+    polar_dec(S,R,T);
+  }else
+  {
+    polar_svd(S,R,T);
+  }
+  R.transposeInPlace();
+
+  // Translation
+  t = Ymean - scale*R.transpose()*Xmean;
+}
+
+
+template <
+  typename DerivedX, 
+  typename DerivedY, 
+  typename Scalar, 
+  int DIM, 
+  int TType>
+IGL_INLINE void igl::procrustes(
+    const Eigen::PlainObjectBase<DerivedX>& X,
+    const Eigen::PlainObjectBase<DerivedY>& Y,
+    bool includeScaling,
+    bool includeReflections,
+    Eigen::Transform<Scalar,DIM,TType>& T)
+{
+  using namespace Eigen;
+  double scale;
+  MatrixXd R;
+  VectorXd t;
+  procrustes(X,Y,includeScaling,includeReflections,scale,R,t);
+
+  // Combine
+  T = Translation<Scalar,DIM>(t) * R * Scaling(scale);
+}
+
+template <
+  typename DerivedX, 
+  typename DerivedY, 
+  typename DerivedR, 
+  typename DerivedT>
+IGL_INLINE void igl::procrustes(
+    const Eigen::PlainObjectBase<DerivedX>& X,
+    const Eigen::PlainObjectBase<DerivedY>& Y,
+    bool includeScaling,
+    bool includeReflections,
+    Eigen::PlainObjectBase<DerivedR>& S,
+    Eigen::PlainObjectBase<DerivedT>& t)
+{
+  double scale;
+  procrustes(X,Y,includeScaling,includeReflections,scale,S,t);
+  S *= scale;
+}
+
+template <
+  typename DerivedX, 
+  typename DerivedY, 
+  typename DerivedR, 
+  typename DerivedT>
+IGL_INLINE void igl::procrustes(
+    const Eigen::PlainObjectBase<DerivedX>& X,
+    const Eigen::PlainObjectBase<DerivedY>& Y,
+    Eigen::PlainObjectBase<DerivedR>& R,
+    Eigen::PlainObjectBase<DerivedT>& t)
+{
+  procrustes(X,Y,false,false,R,t);
+}
+
+template <
+  typename DerivedX, 
+  typename DerivedY, 
+  typename Scalar, 
+  typename DerivedT>
+IGL_INLINE void igl::procrustes(
+    const Eigen::PlainObjectBase<DerivedX>& X,
+    const Eigen::PlainObjectBase<DerivedY>& Y,
+    Eigen::Rotation2D<Scalar>& R,
+    Eigen::PlainObjectBase<DerivedT>& t)
+{
+  using namespace Eigen;
+  assert (X.cols() == 2 && Y.cols() == 2 && "Points must have dimension 2");
+  Matrix2d Rmat;
+  procrustes(X,Y,false,false,Rmat,t);
+  R.fromRotationMatrix(Rmat);
+}

include/igl/procrustes.h

@@ -0,0 +1,137 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 Stefan Brugger <stefanbrugger@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_PROCRUSTES_H
+#define IGL_PROCRUSTES_H
+#include "igl_inline.h"
+
+#include <Eigen/Dense>
+#include <Eigen/Geometry>
+
+namespace igl
+{
+  // Solve Procrustes problem in d dimensions.  Given two point sets X,Y in R^d
+  // find best scale s, orthogonal R  and translation t s.t. |s*X*R + t - Y|^2
+  // is minimized.
+  //
+  // Templates:
+  //    DerivedV point type
+  //    Scalar   scalar type
+  //    DerivedR type of R
+  //    DerivedT type of t
+  // Inputs:
+  //    X  #V by DIM first list of points
+  //    Y  #V by DIM second list of points
+  //    includeScaling  if scaling should be allowed
+  //    includeReflections  if R is allowed to be a reflection
+  // Outputs:
+  //    scale  scaling
+  //    R      orthogonal matrix
+  //    t      translation
+  //
+  // Example:
+  //   MatrixXd X, Y; (containing 3d points as rows)
+  //   double scale;
+  //   MatrixXd R;
+  //   VectorXd t;
+  //   igl::procrustes(X,Y,true,false,scale,R,t);
+  //   R *= scale;
+  //   MatrixXd Xprime = (X * R).rowwise() + t.transpose();
+  //
+  template <
+    typename DerivedX, 
+    typename DerivedY, 
+    typename Scalar, 
+    typename DerivedR, 
+    typename DerivedT>
+  IGL_INLINE void procrustes(
+    const Eigen::PlainObjectBase<DerivedX>& X,
+    const Eigen::PlainObjectBase<DerivedY>& Y,
+    bool includeScaling,
+    bool includeReflections,
+    Scalar& scale,
+    Eigen::PlainObjectBase<DerivedR>& R,
+    Eigen::PlainObjectBase<DerivedT>& t);
+  // Same as above but returns Eigen transformation object.
+  //
+  // Templates:
+  //    DerivedV point type
+  //    Scalar   scalar type
+  //    DIM      point dimension
+  //    TType    type of transformation
+  //             (Isometry,Affine,AffineCompact,Projective)
+  // Inputs:
+  //    X  #V by DIM first list of points
+  //    Y  #V by DIM second list of points
+  //    includeScaling  if scaling should be allowed
+  //    includeReflections  if R is allowed to be a reflection
+  // Outputs:
+  //    T  transformation that minimizes error    
+  //
+  // Example:
+  //   MatrixXd X, Y; (containing 3d points as rows)
+  //   AffineCompact3d T;
+  //   igl::procrustes(X,Y,true,false,T);
+  //   MatrixXd Xprime = (X * T.linear()).rowwise() + T.translation().transpose();
+  template <
+    typename DerivedX, 
+    typename DerivedY, 
+    typename Scalar, 
+    int DIM, 
+    int TType>
+  IGL_INLINE void procrustes(
+    const Eigen::PlainObjectBase<DerivedX>& X,
+    const Eigen::PlainObjectBase<DerivedY>& Y,
+    bool includeScaling,
+    bool includeReflections,
+    Eigen::Transform<Scalar,DIM,TType>& T);
+
+
+  // Convenient wrapper that returns S=scale*R instead of scale and R separately
+  template <
+    typename DerivedX, 
+    typename DerivedY, 
+    typename DerivedR, 
+    typename DerivedT>
+  IGL_INLINE void procrustes(
+    const Eigen::PlainObjectBase<DerivedX>& X,
+    const Eigen::PlainObjectBase<DerivedY>& Y,
+    bool includeScaling,
+    bool includeReflections,
+    Eigen::PlainObjectBase<DerivedR>& S,
+    Eigen::PlainObjectBase<DerivedT>& t);
+
+  // Convenient wrapper for rigid case (no scaling, no reflections)
+  template <
+    typename DerivedX, 
+    typename DerivedY, 
+    typename DerivedR, 
+    typename DerivedT>
+  IGL_INLINE void procrustes(
+    const Eigen::PlainObjectBase<DerivedX>& X,
+    const Eigen::PlainObjectBase<DerivedY>& Y,
+    Eigen::PlainObjectBase<DerivedR>& R,
+    Eigen::PlainObjectBase<DerivedT>& t);
+
+  // Convenient wrapper for 2D case.
+  template <
+    typename DerivedX, 
+    typename DerivedY, 
+    typename Scalar, 
+    typename DerivedT>
+  IGL_INLINE void procrustes(
+    const Eigen::PlainObjectBase<DerivedX>& X,
+    const Eigen::PlainObjectBase<DerivedY>& Y,
+    Eigen::Rotation2D<Scalar>& R,
+    Eigen::PlainObjectBase<DerivedT>& t);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+  #include "procrustes.cpp"
+#endif
+
+#endif

include/igl/project_to_line.cpp

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

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

include/igl/readPLY.cpp

@@ -184,6 +184,18 @@ IGL_INLINE bool igl::readPLY(
     list_to_matrix(vUV,UV);
 }
 
+template <
+  typename DerivedV,
+  typename DerivedF>
+IGL_INLINE bool igl::readPLY(
+  const std::string & filename,
+  Eigen::PlainObjectBase<DerivedV> & V,
+  Eigen::PlainObjectBase<DerivedF> & F)
+{
+  Eigen::MatrixXd N,UV;
+  return readPLY(filename,V,F,N,UV);
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 template bool igl::readPLY<double, int, double, double>(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&, std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > >&, std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > >&);

include/igl/readPLY.h

@@ -46,6 +46,13 @@ namespace igl
     Eigen::PlainObjectBase<DerivedF> & F,
     Eigen::PlainObjectBase<DerivedN> & N,
     Eigen::PlainObjectBase<DerivedUV> & UV);
+  template <
+    typename DerivedV,
+    typename DerivedF>
+  IGL_INLINE bool readPLY(
+    const std::string & filename,
+    Eigen::PlainObjectBase<DerivedV> & V,
+    Eigen::PlainObjectBase<DerivedF> & F);
 }
 #ifndef IGL_STATIC_LIBRARY
 #  include "readPLY.cpp"

include/igl/readSTL.cpp

@@ -88,7 +88,8 @@ IGL_INLINE bool igl::readSTL(
       int ret;
       char facet[IGL_LINE_MAX],normal[IGL_LINE_MAX];
       vector<TypeN > n(3);
-      ret = fscanf(stl_file,"%s %s %lg %lg %lg",facet,normal,&n[0],&n[1],&n[2]);
+      double nd[3];
+      ret = fscanf(stl_file,"%s %s %lg %lg %lg",facet,normal,nd,nd+1,nd+2);
       if(string("endsolid") == facet)
       {
         break;
@@ -98,6 +99,8 @@ IGL_INLINE bool igl::readSTL(
         cerr<<"IOError: "<<filename<<" bad format (1)."<<endl;
         goto close_false;
       }
+      // copy casts to Type
+      n[0] = nd[0]; n[1] = nd[1]; n[2] = nd[2];
       N.push_back(n);
       char outer[IGL_LINE_MAX], loop[IGL_LINE_MAX];
       ret = fscanf(stl_file,"%s %s",outer,loop);
@@ -117,13 +120,16 @@ IGL_INLINE bool igl::readSTL(
         }else if(ret == 1 && string("vertex") == word)
         {
           vector<TypeV> v(3);
-          int ret = fscanf(stl_file,"%lg %lg %lg",&v[0],&v[1],&v[2]);
+          double vd[3];
+          int ret = fscanf(stl_file,"%lg %lg %lg",vd,vd+1,vd+2);
           if(ret != 3)
           {
             cerr<<"IOError: "<<filename<<" bad format (3)."<<endl;
             goto close_false;
           }
           f.push_back(V.size());
+          // copy casts to Type
+          v[0] = vd[0]; v[1] = vd[1]; v[2] = vd[2];
           V.push_back(v);
         }else
         {
@@ -219,4 +225,5 @@ close_true:
 template bool igl::readSTL<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 // generated by autoexplicit.sh
 template bool igl::readSTL<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template bool igl::readSTL<Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/slice_tets.cpp

@@ -0,0 +1,221 @@
+#include "slice_tets.h"
+#include <igl/sort.h>
+#include <igl/cat.h>
+#include <igl/per_face_normals.h>
+#include <cassert>
+#include <algorithm>
+#include <vector>
+
+template <
+  typename DerivedV, 
+  typename DerivedT, 
+  typename Derivedplane,
+  typename DerivedU,
+  typename DerivedG,
+  typename DerivedJ,
+  typename BCType>
+IGL_INLINE void igl::slice_tets(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedT>& T,
+  const Eigen::PlainObjectBase<Derivedplane> & plane,
+  Eigen::PlainObjectBase<DerivedU>& U,
+  Eigen::PlainObjectBase<DerivedG>& G,
+  Eigen::PlainObjectBase<DerivedJ>& J,
+  Eigen::SparseMatrix<BCType> & BC)
+{
+  using namespace Eigen;
+  using namespace std;
+  assert(V.cols() == 3 && "V should be #V by 3");
+  assert(T.cols() == 4 && "T should be #T by 4");
+  assert(plane.size() == 4 && "Plane equation should be 4 coefficients");
+
+  // number of tets
+  const size_t m = T.rows();
+
+  typedef typename DerivedV::Scalar Scalar;
+  typedef typename DerivedT::Scalar Index;
+  typedef Matrix<Scalar,Dynamic,1> VectorXS;
+  typedef Matrix<Scalar,Dynamic,4> MatrixX4S;
+  typedef Matrix<Scalar,Dynamic,3> MatrixX3S;
+  typedef Matrix<Scalar,Dynamic,2> MatrixX2S;
+  typedef Matrix<Index,Dynamic,4> MatrixX4I;
+  typedef Matrix<Index,Dynamic,3> MatrixX3I;
+  typedef Matrix<Index,Dynamic,1> VectorXI;
+  typedef Matrix<bool,Dynamic,1> VectorXb;
+  
+  // Value of plane's implicit function at all vertices
+  VectorXS IV = 
+    (V.col(0)*plane(0) + 
+     V.col(1)*plane(1) + 
+     V.col(2)*plane(2)).array()
+    + plane(3);
+  MatrixX4S IT(m,4);
+  for(size_t t = 0;t<m;t++)
+  {
+    for(size_t c = 0;c<4;c++)
+    {
+      IT(t,c) = IV(T(t,c));
+    }
+  }
+
+  const auto & extract_rows = [](
+    const PlainObjectBase<DerivedT> & T,
+    const MatrixX4S & IT,
+    const VectorXb & I,
+    MatrixX4I  & TI,
+    MatrixX4S & ITI,
+    VectorXI & JI)
+  {
+    const Index num_I = std::count(I.data(),I.data()+I.size(),true);
+    TI.resize(num_I,4);
+    ITI.resize(num_I,4);
+    JI.resize(num_I,1);
+    {
+      size_t k = 0;
+      for(size_t t = 0;t<(size_t)T.rows();t++)
+      {
+        if(I(t))
+        {
+          TI.row(k) = T.row(t);
+          ITI.row(k) = IT.row(t);
+          JI(k) = t;
+          k++;
+        }
+      }
+      assert(k == num_I);
+    }
+  };
+
+  VectorXb I13 = (IT.array()<0).rowwise().count()==1;
+  VectorXb I31 = (IT.array()>0).rowwise().count()==1;
+  VectorXb I22 = (IT.array()<0).rowwise().count()==2;
+  MatrixX4I T13,T31,T22;
+  MatrixX4S IT13,IT31,IT22;
+  VectorXI J13,J31,J22;
+  extract_rows(T,IT,I13,T13,IT13,J13);
+  extract_rows(T,IT,I31,T31,IT31,J31);
+  extract_rows(T,IT,I22,T22,IT22,J22);
+
+  const auto & apply_sort = [] (
+     const MatrixX4I & T, 
+     const MatrixX4I & sJ, 
+     MatrixX4I & sT)
+  {
+    sT.resize(T.rows(),4);
+    for(size_t t = 0;t<(size_t)T.rows();t++)
+    {
+      for(size_t c = 0;c<4;c++)
+      {
+        sT(t,c) = T(t,sJ(t,c));
+      }
+    }
+  };
+
+  const auto & one_below = [&V,&apply_sort](
+    const MatrixX4I & T,
+    const MatrixX4S & IT,
+    MatrixX3I & G,
+    SparseMatrix<BCType> & BC)
+  {
+    // Number of tets
+    const size_t m = T.rows();
+    MatrixX4S sIT;
+    MatrixX4I sJ;
+    sort(IT,2,true,sIT,sJ);
+    MatrixX4I sT;
+    apply_sort(T,sJ,sT);
+    MatrixX3S lambda = 
+      sIT.rightCols(3).array() /
+      (sIT.rightCols(3).colwise()-sIT.col(0)).array();
+    vector<Triplet<BCType> > IJV;
+    IJV.reserve(m*3*2);
+    for(size_t c = 0;c<3;c++)
+    {
+      for(size_t t = 0;t<(size_t)m;t++)
+      {
+        IJV.push_back(Triplet<BCType>(c*m+t,  sT(t,0),  lambda(t,c)));
+        IJV.push_back(Triplet<BCType>(c*m+t,sT(t,c+1),1-lambda(t,c)));
+      }
+    }
+    BC.resize(m*3,V.rows());
+    BC.reserve(m*3*2);
+    BC.setFromTriplets(IJV.begin(),IJV.end());
+    G.resize(m,3);
+    for(size_t c = 0;c<3;c++)
+    {
+      G.col(c).setLinSpaced(m,0+c*m,(m-1)+c*m);
+    }
+  };
+
+  const auto & two_below = [&V,&apply_sort](
+    const MatrixX4I & T,
+    const MatrixX4S & IT,
+    MatrixX3I & G,
+    SparseMatrix<BCType> & BC)
+  {
+    // Number of tets
+    const size_t m = T.rows();
+    MatrixX4S sIT;
+    MatrixX4I sJ;
+    sort(IT,2,true,sIT,sJ);
+    MatrixX4I sT;
+    apply_sort(T,sJ,sT);
+    MatrixX2S lambda = 
+      sIT.rightCols(2).array() /
+      (sIT.rightCols(2).colwise()-sIT.col(0)).array();
+    MatrixX2S gamma = 
+      sIT.rightCols(2).array() /
+      (sIT.rightCols(2).colwise()-sIT.col(1)).array();
+    vector<Triplet<BCType> > IJV;
+    IJV.reserve(m*4*2);
+    for(size_t c = 0;c<2;c++)
+    {
+      for(size_t t = 0;t<(size_t)m;t++)
+      {
+        IJV.push_back(Triplet<BCType>(0*2*m+c*m+t,  sT(t,0),  lambda(t,c)));
+        IJV.push_back(Triplet<BCType>(0*2*m+c*m+t,sT(t,c+2),1-lambda(t,c)));
+        IJV.push_back(Triplet<BCType>(1*2*m+c*m+t,  sT(t,1),   gamma(t,c)));
+        IJV.push_back(Triplet<BCType>(1*2*m+c*m+t,sT(t,c+2),1- gamma(t,c)));
+      }
+    }
+    BC.resize(m*4,V.rows());
+    BC.reserve(m*4*2);
+    BC.setFromTriplets(IJV.begin(),IJV.end());
+    G.resize(2*m,3);
+    G.block(0,0,m,1) = VectorXI::LinSpaced(m,0+0*m,(m-1)+0*m);
+    G.block(0,1,m,1) = VectorXI::LinSpaced(m,0+1*m,(m-1)+1*m);
+    G.block(0,2,m,1) = VectorXI::LinSpaced(m,0+3*m,(m-1)+3*m);
+    G.block(m,0,m,1) = VectorXI::LinSpaced(m,0+0*m,(m-1)+0*m);
+    G.block(m,1,m,1) = VectorXI::LinSpaced(m,0+3*m,(m-1)+3*m);
+    G.block(m,2,m,1) = VectorXI::LinSpaced(m,0+2*m,(m-1)+2*m);
+  };
+
+  MatrixX3I G13,G31,G22;
+  SparseMatrix<BCType> BC13,BC31,BC22;
+  one_below(T13,IT13,G13,BC13);
+  one_below(T31,-IT31,G31,BC31);
+  two_below(T22,IT22,G22,BC22);
+
+  BC = cat(1,cat(1,BC13,BC31),BC22);
+  U = BC*V;
+  G.resize(G13.rows()+G31.rows()+G22.rows(),3);
+  G<<G13,(G31.array()+BC13.rows()),(G22.array()+BC13.rows()+BC31.rows());
+  MatrixX3S N;
+  per_face_normals(U,G,N);
+  Matrix<Scalar,1,3> planeN(plane(0),plane(1),plane(2));
+  VectorXb flip = (N.array().rowwise() * planeN.array()).rowwise().sum()<0;
+  for(size_t g = 0;g<(size_t)G.rows();g++)
+  {
+    if(flip(g))
+    {
+      G.row(g) = G.row(g).reverse().eval();
+    }
+  }
+
+  J.resize(G.rows());
+  J<<J13,J31,J22,J22;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+#endif

include/igl/slice_tets.h

@@ -0,0 +1,60 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 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_SLICE_TETS_H
+#define IGL_SLICE_TETS_H
+#include "igl_inline.h"
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+#include <vector>
+
+namespace igl
+{
+  // SLICE_TETS Slice through a tet mesh (V,T) along a given plane (via its
+  // implicit equation).
+  //
+  // Inputs:
+  //   V  #V by 3 list of tet mesh vertices
+  //   T  #T by 4 list of tet indices into V 
+  //   plane  list of 4 coefficients in the plane equation: [x y z 1]'*plane = 0
+  //   Optional:
+  //     'Manifold' followed by whether to stitch together triangles into a
+  //       manifold mesh {true}: results in more compact U but slightly slower.
+  // Outputs:
+  //   U  #U by 3 list of triangle mesh vertices along slice
+  //   G  #G by 3 list of triangles indices into U
+  //   J  #G list of indices into T revealing from which tet each faces comes
+  //   BC  #U by #V list of barycentric coordinates (or more generally: linear
+  //     interpolation coordinates) so that U = BC*V
+  // 
+  template <
+    typename DerivedV, 
+    typename DerivedT, 
+    typename Derivedplane,
+    typename DerivedU,
+    typename DerivedG,
+    typename DerivedJ,
+    typename BCType>
+  IGL_INLINE void slice_tets(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedT>& T,
+    const Eigen::PlainObjectBase<Derivedplane> & plane,
+    Eigen::PlainObjectBase<DerivedU>& U,
+    Eigen::PlainObjectBase<DerivedG>& G,
+    Eigen::PlainObjectBase<DerivedJ>& J,
+    Eigen::SparseMatrix<BCType> & BC);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "slice_tets.cpp"
+#endif
+
+#endif
+
+

include/igl/sort.h

@@ -65,10 +65,10 @@ namespace igl
   //   index_map  an index map such that sorted[i] = unsorted[index_map[i]]
   template <class T>
   IGL_INLINE void sort(
-      const std::vector<T> &unsorted,
-      const bool ascending,
-      std::vector<T> &sorted,
-      std::vector<size_t> &index_map);
+    const std::vector<T> &unsorted,
+    const bool ascending,
+    std::vector<T> &sorted,
+    std::vector<size_t> &index_map);
 
 }
 

include/igl/viewer/TODOs.txt

@@ -16,6 +16,7 @@
 - zoom with pan rather than scaling
 - refresh draw while resizing
 - use constructor initializer list rather than complicated constructor
+- support per-element alpha values
 + snap to canonical view key shortcut is not working
 + resize TwBar with window
 + trackball should be able to drag over TwBar

include/igl/viewer/Viewer.cpp

@@ -493,8 +493,12 @@ namespace igl
 
     // first try to load it with a plugin
     for (unsigned int i = 0; i<plugins.size(); ++i)
+    {
       if (plugins[i]->load(mesh_file_name_string))
+      {
         return true;
+      }
+    }
 
     data.clear();
 
@@ -510,7 +514,9 @@ namespace igl
     if (extension == "off" || extension =="OFF")
     {
       if (!igl::readOFF(mesh_file_name_string, data.V, data.F))
+      {
         return false;
+      }
     }
     else if (extension == "obj" || extension =="OBJ")
     {
@@ -521,7 +527,9 @@ namespace igl
       Eigen::MatrixXi UV_F;
 
       if (!(igl::readOBJ(mesh_file_name_string, data.V, data.F, corner_normals, fNormIndices, UV_V, UV_F)))
+      {
         return false;
+      }
     }
     else
     {
@@ -989,9 +997,15 @@ namespace igl
 
     opengl.init();
 
+    // Alec: It seems silly to overload launch to take a filename as an
+    // argument. load_mesh_from_file has many side effects so it makes
+    // debugging launch difficult.
+
     // Load the mesh passed as input
     if (filename.size() > 0)
+    {
       load_mesh_from_file(filename.c_str());
+    }
 
     core.align_camera_center(data.V,data.F);
 

include/igl/viewer/ViewerCore.cpp

@@ -10,6 +10,7 @@
 #include <igl/quat_to_mat.h>
 #include <igl/massmatrix.h>
 #include <Eigen/Geometry>
+#include <iostream>
 
 
 Eigen::Matrix4f lookAt (
@@ -143,7 +144,11 @@ IGL_INLINE void igl::ViewerCore::align_camera_center(
   const Eigen::MatrixXi& F)
 {
   get_scale_and_shift_to_fit_mesh(V,F,model_zoom,model_translation);
-  object_scale = (V.colwise().maxCoeff() - V.colwise().minCoeff()).norm();
+  // Rather than crash on empty mesh...
+  if(V.size() > 0)
+  {
+    object_scale = (V.colwise().maxCoeff() - V.colwise().minCoeff()).norm();
+  }
 }
 
 IGL_INLINE void igl::ViewerCore::get_scale_and_shift_to_fit_mesh(

include/igl/writeOBJ.cpp

@@ -97,7 +97,7 @@ IGL_INLINE bool igl::writeOBJ(
   {
     for(int i = 0;i<(int)TC.rows();i++)
     {
-      fprintf(obj_file, "vt %0.15g %0.15g\n",TC(i,0),1-TC(i,1));
+      fprintf(obj_file, "vt %0.15g %0.15g\n",TC(i,0),TC(i,1));
     }
     fprintf(obj_file,"\n");
   }

include/igl/writePLY.cpp

@@ -60,7 +60,7 @@ IGL_INLINE bool igl::writePLY(
   const bool has_texture_coords = UV.rows() > 0;
   std::vector<Vertex> vlist(V.rows());
   std::vector<Face> flist(F.rows());
-  for(size_t i = 0;i<V.rows();i++)
+  for(size_t i = 0;i<(size_t)V.rows();i++)
   {
     vlist[i].x = V(i,0);
     vlist[i].y = V(i,1);
@@ -77,11 +77,11 @@ IGL_INLINE bool igl::writePLY(
       vlist[i].t = UV(i,1);
     }
   }
-  for(size_t i = 0;i<F.rows();i++)
+  for(size_t i = 0;i<(size_t)F.rows();i++)
   {
     flist[i].nverts = F.cols();
     flist[i].verts = new int[F.cols()];
-    for(size_t c = 0;c<F.cols();c++)
+    for(size_t c = 0;c<(size_t)F.cols();c++)
     {
       flist[i].verts[c] = F(i,c);
     }
@@ -133,13 +133,28 @@ IGL_INLINE bool igl::writePLY(
 
   ply_close(ply);
   fclose(fp);
-  for(size_t i = 0;i<F.rows();i++)
+  for(size_t i = 0;i<(size_t)F.rows();i++)
   {
     delete[] flist[i].verts;
   }
   return true;
 }
 
+template <
+  typename DerivedV,
+  typename DerivedF>
+IGL_INLINE bool igl::writePLY(
+  const std::string & filename,
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  const bool ascii)
+{
+  Eigen::MatrixXd N,UV;
+  return writePLY(filename,V,F,N,UV,ascii);
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
+template bool igl::writePLY<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, bool);
+template bool igl::writePLY<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, bool);
 #endif

include/igl/writePLY.h

@@ -34,6 +34,14 @@ namespace igl
     const Eigen::PlainObjectBase<DerivedN> & N,
     const Eigen::PlainObjectBase<DerivedUV> & UV,
     const bool ascii = true);
+  template <
+    typename DerivedV,
+    typename DerivedF>
+  IGL_INLINE bool writePLY(
+    const std::string & filename,
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    const bool ascii = true);
 }
 #ifndef IGL_STATIC_LIBRARY
 #  include "writePLY.cpp"

index.html

@@ -15,7 +15,7 @@
 <h1 id="libigl-asimplecgeometryprocessinglibrary">libigl - A simple C++ geometry processing library</h1>
 
 <figure>
-<img src="tutorial/images/libigl-logo.jpg" alt="" />
+<img src="libigl-teaser.png" alt="" />
 <figcaption></figcaption></figure>
 
 <p><a href="https://github.com/libigl/libigl/">https://github.com/libigl/libigl/</a></p>
@@ -39,7 +39,7 @@ group prototypes a lot in MATLAB, and we have a useful <a href="matlab-to-eigen.
 table</a> from
 MATLAB to libigl/Eigen.</p>
 
-<h1 id="tutorial">Tutorial</h1>
+<h2 id="tutorial">Tutorial</h2>
 
 <p>As of version 1.0, libigl includes an introductory
 <a href="tutorial/tutorial.html">tutorial</a> that covers
@@ -94,7 +94,21 @@ libigl depends only on the <a href="http://eigen.tuxfamily.org">Eigen</a> librar
 
 <p>For more information see our <a href="tutorial/tutorial.html">tutorial</a>.</p>
 
-<h1 id="download">Download</h1>
+<h3 id="gccandtheoptionalcgaldependency">GCC and the optional CGAL dependency</h3>
+
+<p>The <code>include/igl/cgal/*.h</code> headers depend on CGAL. It has come to our attention
+that CGAL does not work properly with GCC 4.8. To the best of our knowledge,
+GCC 4.7 and clang will work correctly.</p>
+
+<h3 id="openmpandwindows">OpenMP and Windows</h3>
+
+<p>Some of our functions will take advantage of OpenMP if available. However, it
+has come to our attention that Visual Studio + Eigen does not work properly
+with OpenMP. Since OpenMP only improves performance without affecting
+functionality we recommend avoiding OpenMP on Windows or proceeding with
+caution.</p>
+
+<h2 id="download">Download</h2>
 
 <p>You can keep up to date by cloning a read-only copy of our GitHub
 <a href="https://github.com/libigl">repository</a>.</p>
@@ -122,11 +136,11 @@ BibTeX entry:</p>
   title = {{libigl}: A simple {C++} geometry processing library},
   author = {Alec Jacobson and Daniele Panozzo and others},
   note = {http://libigl.github.io/libigl/},
-  year = {2014},
+  year = {2015},
 }
 </code></pre>
 
-<h1 id="contact">Contact</h1>
+<h2 id="contact">Contact</h2>
 
 <p>Libigl is a group endeavor led by <a href="http://www.cs.columbia.edu/~jacobson/">Alec
 Jacobson</a> and <a href="http://www.inf.ethz.ch/personal/dpanozzo/">Daniele
@@ -142,10 +156,14 @@ spending time maintaining this.</p>
 <p>If you find bugs or have problems please use our <a href="https://github.com/libigl/libigl/issues">github issue tracking
 page</a>.</p>
 
-<h3 id="copyright">Copyright</h3>
+<h2 id="copyright">Copyright</h2>
 
-<p>2014 Alec Jacobson, Daniele Panozzo, Olga Diamanti, Kenshi
+<p>2015 Alec Jacobson, Daniele Panozzo, Olga Diamanti, Christian Schüller, Kenshi
 Takayama, Leo Sacht, Wenzel Jacob, Nico Pietroni, Amir Vaxman</p>
 
+<figure>
+<img src="tutorial/images/libigl-logo.jpg" alt="" />
+<figcaption></figcaption></figure>
+
 </body>
 </html>

libigl-teaser.png.REMOVED.git-id

@@ -0,0 +1 @@
+a946fc13d18614f2696a1ed94978ab61e76a722b

tutorial/403_BoundedBiharmonicWeights/CMakeLists.txt

@@ -3,18 +3,23 @@ project(403_BoundedBiharmonicWeights)
 
 include("../CMakeLists.shared")
 
+#find_package(EMBREE REQUIRED)
+#include_directories(${EMBREE_INCLUDE_DIRS})
+
 if(NOT MOSEK_FOUND)
   add_definitions(-DIGL_NO_MOSEK)
   if(LIBIGL_USE_STATIC_LIBRARY)
     set(LIBIGLMOSEK_LIBRARY "")
   endif(LIBIGL_USE_STATIC_LIBRARY)
+  set(MOSEK_LIBRARIES "")
+else(NOT MOSEK_FOUND)
+  include_directories( ${MOSEK_INCLUDE_DIR} )
 endif(NOT MOSEK_FOUND)
 
-include_directories( ${MOSEK_INCLUDE_DIR} )
-
 set(SOURCES
 ${PROJECT_SOURCE_DIR}/main.cpp
 )
 
 add_executable(${PROJECT_NAME}_bin ${SOURCES} ${SHARED_SOURCES})
 target_link_libraries(${PROJECT_NAME}_bin ${SHARED_LIBRARIES} ${MOSEK_LIBRARIES})
+#target_link_libraries(${PROJECT_NAME}_bin ${SHARED_LIBRARIES} ${MOSEK_LIBRARIES} ${EMBREE_LIBRARIES})

tutorial/403_BoundedBiharmonicWeights/main.cpp

@@ -22,6 +22,7 @@
 #include <igl/readTGF.h>
 #include <igl/viewer/Viewer.h>
 #include <igl/bbw/bbw.h>
+//#include <igl/embree/bone_heat.h>
 
 #include <Eigen/Geometry>
 #include <Eigen/StdVector>
@@ -143,6 +144,16 @@ int main(int argc, char *argv[])
   {
     return false;
   }
+
+  //MatrixXd Vsurf = V.topLeftCorner(F.maxCoeff()+1,V.cols());
+  //MatrixXd Wsurf;
+  //if(!igl::bone_heat(Vsurf,F,C,VectorXi(),BE,MatrixXi(),Wsurf))
+  //{
+  //  return false;
+  //}
+  //W.setConstant(V.rows(),Wsurf.cols(),1);
+  //W.topLeftCorner(Wsurf.rows(),Wsurf.cols()) = Wsurf = Wsurf = Wsurf = Wsurf;
+
   // Normalize weights to sum to one
   igl::normalize_row_sums(W,W);
   // precompute linear blend skinning matrix

tutorial/501_HarmonicParam/main.cpp

@@ -35,7 +35,7 @@ int main(int argc, char *argv[])
 
   // Find the open boundary
   Eigen::VectorXi bnd;
-  igl::boundary_loop(V,F,bnd);
+  igl::boundary_loop(F,bnd);
 
   // Map the boundary to a circle, preserving edge proportions
   Eigen::MatrixXd bnd_uv;

tutorial/502_LSCMParam/main.cpp

@@ -39,7 +39,7 @@ int main(int argc, char *argv[])
 
   // Fix two points on the boundary
   VectorXi bnd,b(2,1);
-  igl::boundary_loop(V,F,bnd);
+  igl::boundary_loop(F,bnd);
   b(0) = bnd(0);
   b(1) = bnd(round(bnd.size()/2));
   MatrixXd bc(2,2);

tutorial/503_ARAPParam/main.cpp

@@ -46,7 +46,7 @@ int main(int argc, char *argv[])
 
   // Compute the initial solution for ARAP (harmonic parametrization)
   Eigen::VectorXi bnd;
-  igl::boundary_loop(V,F,bnd);
+  igl::boundary_loop(F,bnd);
   Eigen::MatrixXd bnd_uv;
   igl::map_vertices_to_circle(V,bnd,bnd_uv);
 

tutorial/605_Tetgen/main.cpp

@@ -65,7 +65,7 @@ int main(int argc, char *argv[])
   igl::readOFF("../shared/fertility.off",V,F);
 
   // Tetrahedralize the interior
-  igl::tetrahedralize(V,F,"pq1.414", TV,TT,TF);
+  igl::tetrahedralize(V,F,"pq1.414Y", TV,TT,TF);
 
   // Compute barycenters
   igl::barycenter(TV,TT,B);

tutorial/701_Statistics/CMakeLists.txt

@@ -0,0 +1,11 @@
+cmake_minimum_required(VERSION 2.6)
+project(701_Statistics)
+
+include("../CMakeLists.shared")
+
+set(SOURCES
+${PROJECT_SOURCE_DIR}/main.cpp
+)
+
+add_executable(${PROJECT_NAME}_bin ${SOURCES} ${SHARED_SOURCES})
+target_link_libraries(${PROJECT_NAME}_bin ${SHARED_LIBRARIES} ${CGAL_LIBRARIES})

tutorial/701_Statistics/main.cpp

@@ -0,0 +1,53 @@
+#include <igl/readOBJ.h>
+
+#include <Eigen/Core>
+#include <iostream>
+
+#include <igl/is_irregular_vertex.h>
+#include <igl/doublearea.h>
+#include <igl/internal_angles.h>
+
+int main(int argc, char *argv[])
+{
+  using namespace Eigen;
+  using namespace std;
+
+  MatrixXd V;
+  MatrixXi F;
+
+  igl::readOBJ("../shared/horse_quad.obj",V,F);
+
+  // Count the number of irregular vertices, the border is ignored
+  vector<bool> irregular = igl::is_irregular_vertex(V,F);
+
+  int vertex_count = V.rows();
+  int irregular_vertex_count = std::count(irregular.begin(),irregular.end(),true);
+  double irregular_ratio = double(irregular_vertex_count)/vertex_count;
+
+  printf("Irregular vertices: \n%d/%d (%.2f%%)\n",irregular_vertex_count,vertex_count, irregular_ratio*100);
+
+  // Compute areas, min, max and standard deviation
+  VectorXd area;
+  igl::doublearea(V,F,area);
+  area = area.array() / 2;
+  
+  double area_avg   = area.mean();
+  double area_min   = area.minCoeff() / area_avg;
+  double area_max   = area.maxCoeff() / area_avg;
+  double area_sigma = sqrt( ((area.array()-area_avg)/area_avg).square().mean() );
+
+  printf("Areas (Min/Max)/Avg_Area Sigma: \n%.2f/%.2f (%.2f)\n",area_min,area_max,area_sigma);
+
+  // Compute per face angles, min, max and standard deviation
+  MatrixXd angles;
+  igl::angles(V,F,angles);
+  angles = 360.0 * (angles/(2*M_PI)); // Convert to degrees
+  
+  double angle_avg   = angles.mean();
+  double angle_min   = angles.minCoeff();
+  double angle_max   = angles.maxCoeff();
+  double angle_sigma = sqrt( (angles.array()-angle_avg).square().mean() );
+  
+  printf("Angles in degrees (Min/Max) Sigma: \n%.2f/%.2f (%.2f)\n",angle_min,angle_max,angle_sigma);
+
+}

tutorial/702_WindingNumber/CMakeLists.txt

@@ -0,0 +1,11 @@
+cmake_minimum_required(VERSION 2.6)
+project(702_WindingNumber)
+
+include("../CMakeLists.shared")
+
+set(SOURCES
+${PROJECT_SOURCE_DIR}/main.cpp
+)
+
+add_executable(${PROJECT_NAME}_bin ${SOURCES} ${SHARED_SOURCES})
+target_link_libraries(${PROJECT_NAME}_bin ${SHARED_LIBRARIES})

tutorial/702_WindingNumber/main.cpp

@@ -0,0 +1,139 @@
+#include <igl/readMESH.h>
+#include <igl/winding_number.h>
+#include <igl/barycenter.h>
+#include <igl/boundary_facets.h>
+#include <igl/parula.h>
+#include <igl/slice_tets.h>
+#include <igl/slice.h>
+#include <igl/viewer/Viewer.h>
+#include <Eigen/Sparse>
+#include <iostream>
+
+Eigen::MatrixXd V,BC;
+Eigen::VectorXd W;
+Eigen::MatrixXi T,F,G;
+double slice_z = 0.5;
+enum OverLayType
+{
+  OVERLAY_NONE = 0,
+  OVERLAY_INPUT = 1,
+  OVERLAY_OUTPUT = 2,
+  NUM_OVERLAY = 3,
+} overlay = OVERLAY_NONE;
+
+void update_visualization(igl::Viewer & viewer)
+{
+  using namespace Eigen;
+  using namespace std;
+  Eigen::Vector4d plane(
+    0,0,1,-((1-slice_z)*V.col(2).minCoeff()+slice_z*V.col(2).maxCoeff()));
+  MatrixXd V_vis;
+  MatrixXi F_vis;
+  VectorXi J;
+  SparseMatrix<double> bary;
+  igl::slice_tets(V,T,plane,V_vis,F_vis,J,bary);
+  VectorXd W_vis;
+  igl::slice(W,J,W_vis);
+  MatrixXd C_vis;
+  // color without normalizing
+  igl::parula(W_vis,false,C_vis);
+
+
+  const auto & append_mesh = [&C_vis,&F_vis,&V_vis](
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    const RowVector3d & color)
+  {
+    F_vis.conservativeResize(F_vis.rows()+F.rows(),3);
+    F_vis.bottomRows(F.rows()) = F.array()+V_vis.rows();
+    V_vis.conservativeResize(V_vis.rows()+V.rows(),3);
+    V_vis.bottomRows(V.rows()) = V;
+    C_vis.conservativeResize(C_vis.rows()+F.rows(),3);
+    C_vis.bottomRows(F.rows()).rowwise() = color;
+  };
+  switch(overlay)
+  {
+    case OVERLAY_INPUT:
+      append_mesh(V,F,RowVector3d(1.,0.894,0.227));
+      break;
+    case OVERLAY_OUTPUT:
+      append_mesh(V,G,RowVector3d(0.8,0.8,0.8));
+      break;
+    default:
+      break;
+  }
+  viewer.data.clear();
+  viewer.data.set_mesh(V_vis,F_vis);
+  viewer.data.set_colors(C_vis);
+  viewer.data.set_face_based(true);
+}
+
+bool key_down(igl::Viewer& viewer, unsigned char key, int mod)
+{
+  switch(key)
+  {
+    default:
+      return false;
+    case ' ':
+      overlay = (OverLayType)((1+(int)overlay)%NUM_OVERLAY);
+      break;
+    case '.':
+      slice_z = std::min(slice_z+0.01,0.99);
+      break;
+    case ',':
+      slice_z = std::max(slice_z-0.01,0.01);
+      break;
+  }
+  update_visualization(viewer);
+  return true;
+}
+
+int main(int argc, char *argv[])
+{
+  using namespace Eigen;
+  using namespace std;
+
+  cout<<"Usage:"<<endl;
+  cout<<"[space]  toggle showing input mesh, output mesh or slice "<<endl;
+  cout<<"         through tet-mesh of convex hull."<<endl;
+  cout<<"'.'/','  push back/pull forward slicing plane."<<endl;
+  cout<<endl;
+
+  // Load mesh: (V,T) tet-mesh of convex hull, F contains facets of input
+  // surface mesh _after_ self-intersection resolution
+  igl::readMESH("../shared/big-sigcat.mesh",V,T,F);
+
+  // Compute barycenters of all tets
+  igl::barycenter(V,T,BC);
+
+  // Compute generalized winding number at all barycenters
+  cout<<"Computing winding number over all "<<T.rows()<<" tets..."<<endl;
+  igl::winding_number(V,F,BC,W);
+
+  // Extract interior tets
+  MatrixXi CT((W.array()>0.5).count(),4);
+  {
+    size_t k = 0;
+    for(size_t t = 0;t<T.rows();t++)
+    {
+      if(W(t)>0.5)
+      {
+        CT.row(k) = T.row(t);
+        k++;
+      }
+    }
+  }
+  // find bounary facets of interior tets
+  igl::boundary_facets(CT,G);
+  // boundary_facets seems to be reversed...
+  G = G.rowwise().reverse().eval();
+
+  // normalize
+  W = (W.array() - W.minCoeff())/(W.maxCoeff()-W.minCoeff());
+
+  // Plot the generated mesh
+  igl::Viewer viewer;
+  update_visualization(viewer);
+  viewer.callback_key_down = &key_down;
+  viewer.launch();
+}

tutorial/CMakeLists.shared

@@ -2,7 +2,7 @@ if(NOT CMAKELISTS_SHARED_INCLUDED)
 set(CMAKELISTS_SHARED_INCLUDED TRUE)
 SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11")
 find_package(OpenMP)
-if (OPENMP_FOUND)
+if (OPENMP_FOUND AND NOT WIN32)
   set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")
   set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
 endif()

tutorial/CMakeLists.txt

@@ -77,3 +77,5 @@ add_subdirectory("608_LIM")
 if(CGAL_FOUND)
 add_subdirectory("609_Boolean")
 endif()
+add_subdirectory("701_Statistics")
+add_subdirectory("702_WindingNumber")

tutorial/cmake/FindGLEW.cmake

@@ -7,6 +7,7 @@
 
 FIND_PATH(GLEW_INCLUDE_DIR GL/glew.h
    ${PROJECT_SOURCE_DIR}/../../external/glew/include
+   ${PROJECT_SOURCE_DIR}/../external/glew/include
    /usr/include
    /usr/local/include
    $ENV{GLEWROOT}/include

tutorial/cmake/FindTETGEN.cmake

@@ -6,10 +6,10 @@
 #  TETGEN_SOURCES - the TETGEN source files
 
 FIND_PATH(TETGEN_INCLUDE_DIR tetgen.h
-   /usr/include
-   /usr/local/include
    ${PROJECT_SOURCE_DIR}/../libigl/external/tetgen/
    ${PROJECT_SOURCE_DIR}/../../external/tetgen/
+   /usr/include
+   /usr/local/include
    NO_DEFAULT_PATH
 )
 

tutorial/images/big-sigcat-winding-number.gif.REMOVED.git-id

@@ -0,0 +1 @@
+848269f54baa921cd2b5799d8cc77e3024bae9d5

tutorial/shared/big-sigcat.mesh.REMOVED.git-id

@@ -0,0 +1 @@
+0ca1d1e69f9b544b53246a4645a4a711c40de98c

tutorial/shared/horse_quad.obj.REMOVED.git-id

@@ -0,0 +1 @@
+7ab1d930fa0c28df9118e924b7e30c9af6875474

tutorial/tutorial.html.REMOVED.git-id

@@ -1 +1 @@
-a30fd5f786682eb951f8170826b3b3212240db31
+2e5cd62765c38c507481d1880609faab4f97a9a9

tutorial/tutorial.md.REMOVED.git-id

@@ -1 +1 @@
-227bf977c70fc9d17fb2b2823c0818e6205893d2
+fa82ad61f8b29ace7195ce6be520e1735648430e