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

Conflicts:
	include/igl/embree/reorient_facets_raycast.cpp
	include/igl/local_basis.cpp
	include/igl/min_quad_with_fixed.h
	include/igl/xml/XMLSerializer.h

Former-commit-id: 1f42dca56021ecca8d3b3ea99ebee40d5dea3cbd

.gitignore

@@ -57,3 +57,9 @@ external/tetgen/tetgen
 external/tinyxml2/test
 external/tinyxml2/tinyxml2.pc
 external/yimg/showpng
+README.html
+tutorial/readme.html
+tutorial/*/build/*
+tutorial/*/Makefile
+external/glew/build
+*buildXcode*

.mailmap

@@ -11,3 +11,9 @@ Alec Jacobson <alecjacobson@gmail.com> ajx <devnull@localhost>
 Alec Jacobson <alecjacobson@gmail.com> mangledorf <alecjacobson@gmail.com>
 Daniele Panozzo <panozzo@inf.ethz.ch> Daniele Panozzo <daniele.panozzo@gmail.com>
 Daniele Panozzo <panozzo@inf.ethz.ch> dpanozzo <devnull@localhost>
+Wenzel Jakob <wenzel@inf.ethz.ch> Wenzel Jakob <wenzel@cs.cornell.edu>
+Olga Diamanti <olga.diamanti@inf.ethz.ch> dolga <devnull@localhost>
+schuellc <schuellchr@gmail.com> schuellc <devnull@localhost>
+Kenshi Takayama <kenshi84@gmail.com> Kenshi Takayama (kenshi <Kenshi Takayama (kenshi@gmail.com)>
+Kenshi Takayama <kenshi84@gmail.com> kenshi <kenshi@jackal.ethz.ch>
+Kenshi Takayama <kenshi84@gmail.com> kenshi84 <kenshi84@gmail.com>

Makefile

@@ -54,6 +54,10 @@ ifeq ($(IGL_WITH_TETGEN),1)
 	EXTRA_DIRS+=include/igl/tetgen
 	EXTRAS += tetgen
 endif
+ifeq ($(IGL_WITH_VIEWER),1)
+	EXTRA_DIRS+=include/igl/viewer
+	EXTRAS += viewer
+endif
 ifeq ($(IGL_WITH_XML),1)
 	EXTRA_DIRS+=include/igl/xml
 	EXTRAS += xml

Makefile.conf

@@ -47,6 +47,7 @@ endif
 ifeq ($(IGL_USERNAME),ajx)
 	MOSEKPLATFORM=osx64x86
 	MOSEKVERSION=7
+	IGL_WITH_VIEWER=1
 	IGL_WITH_TETGEN=1
 	IGL_WITH_EMBREE=1
 	IGL_WITH_MATLAB=1

README.md

@@ -25,6 +25,7 @@ listed below.
 - Mosek  libiglmosek extra only
 - Matlab  libiglmatlab extra only
 - boost  libiglboost, libiglcgal extra only
+- SSE/AVX  libiglsvd3x3 extra only
 - CGAL  libiglcgal extra only
     * boost 
     * gmp
@@ -154,6 +155,9 @@ To build the a static YImg library on Mac OS X issue:
     cd external/yimg
     make
 
+You may need to install libpng. Systems with X11 might find this already
+installed at `/usr/X11/lib`.
+
 
 ### Windows (Experimental) ###
 To build a static library (.lib) on windows, open Visual Studio 2010.

RELEASE_HISTORY.txt

@@ -1,3 +1,4 @@
+0.4.6  Generalized Winding Numbers
 0.4.5  CGAL extra: mesh selfintersection
 0.4.4  STL file format support
 0.4.3  ARAP implementation

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
-0.4.5
+0.4.6

documentation/implemented-papers.tex

@@ -40,6 +40,9 @@ choosing skinning transformations \cite{Jacobson:FAST:2012}.
 ``skeletal subspace deformation'', or ``enveloping''. This technique is often
 attributed to \cite{Magnenat-Thalmann:1988:JLD}.
 
+\paragraph{\texttt{winding\_number}} implements ``Generalized Winding Numbers''
+\cite{Jacobson:WN:2013}
+
 \bibliographystyle{acmsiggraph}
 \bibliography{references} 
 

documentation/references.bib.REMOVED.git-id

@@ -1 +1 @@
-92b6b20756fd71dd62d7b52066a7bcf07e4acd28
+346cd2c8206ee2c83e7aef50adb832e483320d09

examples/ambient-occlusion-mex/compile.m

@@ -4,6 +4,6 @@ mex -v -o ambient_occlusion -DMEX -largeArrayDims ...
   CXXFLAGS="\$CXXFLAGS -m64 -msse4.2 -fopenmp" ...
   -I/usr/local/igl/libigl/external/embree/embree/ ...
   -I/usr/local/igl/libigl/external/embree/ ...
-  -L/usr/local/igl/libigl/external/embree/bin -lembree -lsys ... 
+  -L/usr/local/igl/libigl/external/embree/build -lembree -lsys ... 
   mexFunction.cpp parse_rhs.cpp;
 

examples/arap/Makefile

@@ -13,6 +13,8 @@ LIBIGL_INC=-I$(LIBIGL)/include
 LIBIGL_LIB=-L$(LIBIGL)/lib -ligl -liglpng -liglsvd3x3
 
 EIGEN3_INC=-I/opt/local/include/eigen3 -I/opt/local/include/eigen3/unsupported
+#CFLAGS+=-DIGL_HEADER_ONLY
+
 
 #EMBREE=$(LIBIGL)/external/embree
 #EMBREE_INC=-I$(EMBREE)/rtcore -I$(EMBREE)/common
@@ -23,14 +25,25 @@ EMBREE_LIB=-L$(EMBREE)/build -lembree -lsys
 
 # YIMAGE Library
 YIMG=$(LIBIGL)/external/yimg/
-YIMG_LIB=-L$(YIMG) -lyimg -lz -L/usr/X11/lib -lpng -bind_at_load
+YIMG_LIB=-L$(YIMG) -lyimg -lz -L/opt/local/lib -lpng
 YIMG_INC=-I/usr/X11/include -I$(YIMG)
 
 ANTTWEAKBAR_INC=-I$(LIBIGL)/external/AntTweakBar/include
-ANTTWEAKBAR_LIB=-L$(LIBIGL)/external/AntTweakBar/lib -lAntTweakBar -framework AppKit
-
-INC=$(LIBIGL_INC) $(ANTTWEAKBAR_INC) $(EIGEN3_INC) $(YIMG_INC)
-LIB=$(OPENGL_LIB) $(GLUT_LIB) $(ANTTWEAKBAR_LIB) $(LIBIGL_LIB) $(YIMG_LIB)
+ANTTWEAKBAR_LIB=-L$(LIBIGL)/external/AntTweakBar/lib -lAntTweakBar
+UNAME := $(shell uname)
+# Apple needs to load the AppKit framework for anttweakbar and maybe bind at
+# load for png
+ifeq ($(UNAME), Darwin)
+	YIMG_LIB+=-bind_at_load
+	ANTTWEAKBAR_LIB+=-framework AppKit
+endif
+
+# SVD 
+SINGULAR_VALUE_DECOMPOSITION_INC=\
+	-I$(LIBIGL)/external/Singular_Value_Decomposition/
+
+INC=$(LIBIGL_INC) $(ANTTWEAKBAR_INC) $(EIGEN3_INC) $(YIMG_INC) $(SINGULAR_VALUE_DECOMPOSITION_INC)
+LIB=$(OPENGL_LIB) $(GLUT_LIB) $(ANTTWEAKBAR_LIB) $(LIBIGL_LIB) $(YIMG_LIB) 
 
 example: example.o
 	g++ $(OPENMP) $(AFLAGS) $(CFLAGS) $(LIB) -o example $^

examples/arap/example.cpp

@@ -28,7 +28,6 @@
 #include <igl/material_colors.h>
 #include <igl/barycenter.h>
 #include <igl/matlab_format.h>
-#include <igl/material_colors.h>
 #include <igl/ReAntTweakBar.h>
 #include <igl/pathinfo.h>
 #include <igl/Camera.h>
@@ -70,6 +69,9 @@ Eigen::Quaterniond animation_to_quat;
 // Use vector for range-based `for`
 std::vector<State> undo_stack;
 std::vector<State> redo_stack;
+bool paused = false;
+double t = 0;
+double pause_time = 0.0;
 
 void push_undo()
 {
@@ -133,7 +135,7 @@ float light_pos[4] = {0.1,0.1,-0.9,0};
 Eigen::MatrixXd V,U,N,C,mid;
 Eigen::VectorXi S;
 igl::ARAPData arap_data;
-Eigen::MatrixXi F;
+Eigen::MatrixXi F,T;
 int selected_col = 0;
 // Faces
 // Bounding box diagonal length
@@ -141,8 +143,6 @@ double bbd;
 int tot_num_samples = 0;
 #define REBAR_NAME "temp.rbr"
 igl::ReTwBar rebar; // Pointer to the tweak bar
-bool flip_y = false;
-bool rotate_xy = false;
 
 int num_in_selection(const Eigen::VectorXi & S)
 {
@@ -166,6 +166,14 @@ bool init_arap()
   assert(S.rows() == V.rows());
   C.resize(S.rows(),3);
   MatrixXd bc = MatrixXd::Zero(b.size(),S.maxCoeff()+1);
+  MatrixXi * Ele;
+  if(T.rows()>0)
+  {
+    Ele = &T;
+  }else
+  {
+    Ele = &F;
+  }
   // get b from S
   {
     int bi = 0;
@@ -176,12 +184,17 @@ bool init_arap()
         b(bi) = v;
         bc(bi,S(v)) = 1;
         bi++;
-        if(S(v) == 0)
-        {
-          C.row(v) = RowVector3d(0.039,0.31,1);
-        }else
+        switch(S(v))
         {
-          C.row(v) = RowVector3d(1,0.41,0.70);
+          case 0:
+            C.row(v) = RowVector3d(0.039,0.31,1);
+            break;
+          case 1:
+            C.row(v) = RowVector3d(1,0.41,0.70);
+            break;
+          default:
+            C.row(v) = RowVector3d(0.4,0.8,0.3);
+            break;
         }
       }else
       {
@@ -197,12 +210,15 @@ bool init_arap()
   VectorXi _S;
   VectorXd _D;
   MatrixXd W;
-  if(!harmonic(V,F,b,bc,1,W))
+  if(!harmonic(V,*Ele,b,bc,1,W))
   {
     return false;
   }
-  partition(W,100,arap_data.G,_S,_D);
-  return arap_precomputation(V,F,b,arap_data);
+  //arap_data.with_dynamics = true;
+  //arap_data.h = 0.5;
+  //arap_data.max_iter = 100;
+  //partition(W,100,arap_data.G,_S,_D);
+  return arap_precomputation(V,*Ele,V.cols(),b,arap_data);
 }
 
 bool update_arap()
@@ -213,6 +229,10 @@ bool update_arap()
   MatrixXd bc(num_in_selection(S),V.cols());
   // get b from S
   {
+    if(!paused)
+    {
+      t = get_seconds()-pause_time;
+    }
     int bi = 0;
     for(int v = 0;v<S.rows(); v++)
     {
@@ -223,25 +243,25 @@ bool update_arap()
         {
           case 0:
           {
-            //const double r = mid(0)*0.25;
-            //bc(bi,0) += r*cos(0.5*get_seconds()*2.*PI);
-            //bc(bi,1) -= r+r*sin(0.5*get_seconds()*2.*PI);
+            const double r = mid(0)*0.25;
+            bc(bi,0) += r*cos(0.5*t*2.*PI);
+            bc(bi,1) -= r+r*sin(0.5*t*2.*PI);
             break;
           }
           case 1:
           {
-            //const double r = mid(1)*0.15;
-            //bc(bi,1) += r+r*cos(0.15*get_seconds()*2.*PI);
-            //bc(bi,2) -= r*sin(0.15*get_seconds()*2.*PI);
+            const double r = mid(1)*0.15;
+            bc(bi,1) += r+r*cos(0.15*t*2.*PI);
+            bc(bi,2) -= r*sin(0.15*t*2.*PI);
 
             //// Pull-up
             //bc(bi,0) += 0.42;//mid(0)*0.5;
             //bc(bi,1) += 0.55;//mid(0)*0.5;
-            // Bend
-            Vector3d t(-1,0,0);
-            Quaterniond q(AngleAxisd(PI/1.5,Vector3d(0,1.0,0.1).normalized()));
-            const Vector3d a = bc.row(bi);
-            bc.row(bi) = (q*(a-t) + t) + Vector3d(1.5,0.1,0.9);
+            //// Bend
+            //Vector3d t(-1,0,0);
+            //Quaterniond q(AngleAxisd(PI/1.5,Vector3d(0,1.0,0.1).normalized()));
+            //const Vector3d a = bc.row(bi);
+            //bc.row(bi) = (q*(a-t) + t) + Vector3d(1.5,0.1,0.9);
                 
 
             break;
@@ -576,6 +596,7 @@ void mouse_drag(int mouse_x, int mouse_y)
 void key(unsigned char key, int mouse_x, int mouse_y)
 {
   using namespace std;
+  using namespace igl;
   switch(key)
   {
     // ESC
@@ -584,6 +605,20 @@ void key(unsigned char key, int mouse_x, int mouse_y)
     // ^C
     case char(3):
       exit(0);
+    case ' ':
+      {
+        static double pause_start,pause_stop;
+        paused = !paused;
+        if(paused)
+        {
+          pause_start = get_seconds();
+        }else
+        {
+          pause_stop = get_seconds();
+          pause_time += (pause_stop-pause_start);
+        }
+        break;
+      }
     default:
       if(!TwEventKeyboardGLUT(key,mouse_x,mouse_y))
       {
@@ -603,6 +638,8 @@ int main(int argc, char * argv[])
   // init mesh
   string filename = "../shared/decimated-knight.obj";
   string sfilename = "../shared/decimated-knight-selection.dmat";
+  //string filename = "../shared/decimated-knight.mesh";
+  //string sfilename = "../shared/decimated-knight-1-selection.dmat";
   if(argc < 3)
   {
     cerr<<"Usage:"<<endl<<"    ./example input.obj selection.dmat"<<endl;
@@ -613,11 +650,27 @@ int main(int argc, char * argv[])
     filename = argv[1];
     sfilename = argv[2];
   }
+  string d,b,ext,f;
+  pathinfo(filename,d,b,ext,f);
+  // Convert extension to lower case
+  transform(ext.begin(), ext.end(), ext.begin(), ::tolower);
 
   vector<vector<double > > vV,vN,vTC;
   vector<vector<int > > vF,vTF,vFN;
   // Convert extension to lower case
-  if(!igl::readOBJ(filename,vV,vTC,vN,vF,vTF,vFN))
+  if(ext == "obj")
+  {
+    if(!igl::readOBJ(filename,vV,vTC,vN,vF,vTF,vFN))
+    {
+      return 1;
+    }
+  }else if(ext == "mesh")
+  {
+    if(!igl::readMESH(filename,V,T,F))
+    {
+      return 1;
+    }
+  }else
   {
     return 1;
   }
@@ -660,8 +713,6 @@ int main(int argc, char * argv[])
     "igl_trackball,two-a...-fixed-up");
   rebar.TwAddVarCB( "rotation_type", RotationTypeTW,
     set_rotation_type,get_rotation_type,NULL,"keyIncr=] keyDecr=[");
-  rebar.TwAddVarRW("flip_y", TW_TYPE_BOOLCPP, &flip_y,"key=f");
-  rebar.TwAddVarRW("rotate_xy", TW_TYPE_BOOLCPP, &rotate_xy,"key=r");
   rebar.load(REBAR_NAME);
 
   glutInitDisplayString( "rgba depth double samples>=8 ");

include/igl/C_STR.h

@@ -5,8 +5,8 @@
 // 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 C_STR_H
-#define C_STR_H
+#ifndef IGL_C_STR_H
+#define IGL_C_STR_H
 // http://stackoverflow.com/a/2433143/148668
 // Suppose you have a function:
 //   void func(const char * c);

include/igl/InElementAABB.h

@@ -0,0 +1,397 @@
+#ifndef IGL_IN_ELEMENT_AABB_H
+#define IGL_IN_ELEMENT_AABB_H
+
+#include <Eigen/Core>
+#include <memory>
+
+namespace igl
+{
+  class InElementAABB
+  {
+    public:
+      std::shared_ptr<InElementAABB> m_left, m_right;
+      Eigen::RowVectorXd m_bb_min,m_bb_max;
+      // -1 non-leaf
+      int m_element;
+      InElementAABB():
+        m_left(NULL), m_right(NULL),
+        m_bb_min(), m_bb_max(), m_element(-1)
+      {}
+      // Build an Axis-Aligned Bounding Box tree for a given mesh and given
+      // serialization of a previous AABB tree.
+      //
+      // Inputs:
+      //   V  #V by dim list of mesh vertex positions. 
+      //   Ele  #Ele by dim+1 list of mesh indices into #V. 
+      //   bb_mins  max_tree by dim list of bounding box min corner positions
+      //   bb_maxs  max_tree by dim list of bounding box max corner positions
+      //   elements  max_tree list of element or (not leaf id) indices into Ele
+      //   i  recursive call index {0}
+      inline void init(
+          const Eigen::MatrixXd & V,
+          const Eigen::MatrixXi & Ele, 
+          const Eigen::MatrixXd & bb_mins,
+          const Eigen::MatrixXd & bb_maxs,
+          const Eigen::VectorXi & elements,
+          const int i = 0);
+      // Wrapper for root with empty serialization
+      inline void init(
+        const Eigen::MatrixXd & V,
+        const Eigen::MatrixXi & Ele);
+      // Build an Axis-Aligned Bounding Box tree for a given mesh.
+      //
+      // Inputs:
+      //   V  #V by dim list of mesh vertex positions. 
+      //   Ele  #Ele by dim+1 list of mesh indices into #V. 
+      //   SI  #Ele by dim list revealing for each coordinate where Ele's
+      //     barycenters would be sorted: SI(e,d) = i --> the dth coordinate of
+      //     the barycenter of the eth element would be placed at position i in a
+      //     sorted list.
+      //   I  #I list of indices into Ele of elements to include (for recursive
+      //     calls)
+      // 
+      inline void init(
+          const Eigen::MatrixXd & V,
+          const Eigen::MatrixXi & Ele, 
+          const Eigen::MatrixXi & SI,
+          const Eigen::VectorXi & I);
+      // Find the indices of elements containing given point.
+      //
+      // Inputs:
+      //   V  #V by dim list of mesh vertex positions. **Should be same as used to
+      //     construct mesh.**
+      //   Ele  #Ele by dim+1 list of mesh indices into #V. **Should be same as used to
+      //     construct mesh.**
+      //   q  dim row-vector query position
+      //   first  whether to only return first element containing q
+      // Returns:
+      //   list of indices of elements containing q
+      inline std::vector<int> find(
+        const Eigen::MatrixXd & V,
+        const Eigen::MatrixXi & Ele, 
+        const Eigen::RowVectorXd & q, 
+        const bool first=false) const;
+  
+      // If number of elements m then total tree size should be 2*h where h is
+      // the deepest depth 2^ceil(log(#Ele*2-1))
+      inline int subtree_size();
+  
+      // Serialize this class into 3 arrays (so we can pass it pack to matlab)
+      //
+      // Outputs:
+      //   bb_mins  max_tree by dim list of bounding box min corner positions
+      //   bb_maxs  max_tree by dim list of bounding box max corner positions
+      //   elements  max_tree list of element or (not leaf id) indices into Ele
+      //   i  recursive call index into these arrays {0}
+      inline void serialize(
+        Eigen::MatrixXd & bb_mins,
+        Eigen::MatrixXd & bb_maxs,
+        Eigen::VectorXi & elements,
+        const int i = 0);
+  };
+}
+
+// Implementation
+#include <igl/volume.h>
+#include <igl/colon.h>
+#include <igl/doublearea.h>
+#include <igl/matlab_format.h>
+#include <igl/colon.h>
+#include <igl/sort.h>
+#include <igl/barycenter.h>
+#include <iostream>
+
+inline void igl::InElementAABB::init(
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & Ele, 
+    const Eigen::MatrixXd & bb_mins,
+    const Eigen::MatrixXd & bb_maxs,
+    const Eigen::VectorXi & elements,
+    const int i)
+{
+  using namespace std;
+  using namespace Eigen;
+  using namespace igl;
+  if(bb_mins.size() > 0)
+  {
+    assert(bb_mins.rows() == bb_maxs.rows() && "Serial tree arrays must match");
+    assert(bb_mins.cols() == V.cols() && "Serial tree array dim must match V");
+    assert(bb_mins.cols() == bb_maxs.cols() && "Serial tree arrays must match");
+    assert(bb_mins.rows() == elements.rows() &&
+      "Serial tree arrays must match");
+    // construct from serialization
+    m_bb_min = bb_mins.row(i);
+    m_bb_max = bb_maxs.row(i);
+    m_element = elements(i);
+    // Not leaf then recurse
+    if(m_element == -1)
+    {
+      m_left = make_shared<InElementAABB>();
+      m_left->init( V,Ele,bb_mins,bb_maxs,elements,2*i+1);
+      m_right = make_shared<InElementAABB>();
+      m_right->init( V,Ele,bb_mins,bb_maxs,elements,2*i+2);
+    }
+  }else
+  {
+    VectorXi allI = colon<int>(0,Ele.rows()-1);
+    MatrixXd BC;
+    barycenter(V,Ele,BC);
+    MatrixXi SI(BC.rows(),BC.cols());
+    {
+      MatrixXd _;
+      MatrixXi IS;
+      igl::sort(BC,1,true,_,IS);
+      // Need SI(i) to tell which place i would be sorted into
+      const int dim = IS.cols();
+      for(int i = 0;i<IS.rows();i++)
+      {
+        for(int d = 0;d<dim;d++)
+        {
+          SI(IS(i,d),d) = i;
+        }
+      }
+    }
+    init(V,Ele,SI,allI);
+  }
+}
+
+inline void igl::InElementAABB::init(
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & Ele)
+{
+  using namespace Eigen;
+  return init(V,Ele,MatrixXd(),MatrixXd(),VectorXi(),0);
+}
+
+inline void igl::InElementAABB::init(
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & Ele, 
+    const Eigen::MatrixXi & SI,
+    const Eigen::VectorXi & I)
+{
+  using namespace Eigen;
+  using namespace std;
+  using namespace igl;
+  const int dim = V.cols();
+  const double inf = numeric_limits<double>::infinity();
+  m_bb_min.setConstant(1,dim, inf);
+  m_bb_max.setConstant(1,dim,-inf);
+  // Compute bounding box
+  for(int i = 0;i<I.rows();i++)
+  {
+    for(int c = 0;c<Ele.cols();c++)
+    {
+      for(int d = 0;d<dim;d++)
+      {
+        m_bb_min(d) = min(m_bb_min(d),V(Ele(I(i),c),d));
+        m_bb_max(d) = max(m_bb_max(d),V(Ele(I(i),c),d));
+      }
+    }
+  }
+  switch(I.size())
+  {
+    case 0:
+      {
+        assert(false);
+      }
+    case 1:
+      {
+        m_element = I(0);
+        break;
+      }
+    default:
+      {
+        // Compute longest direction
+        int max_d = -1;
+        double max_len = -inf;
+        for(int d = 0;d<dim;d++)
+        {
+          const auto diff = (m_bb_max[d] - m_bb_min[d]);
+          if( diff > max_len )
+          {
+            max_len = diff;
+            max_d = d;
+          }
+        }
+        // Can't use median on BC directly because many may have same value,
+        // but can use median on sorted BC indices
+        VectorXi SIdI(I.rows());
+        for(int i = 0;i<I.rows();i++)
+        {
+          SIdI(i) = SI(I(i),max_d);
+        }
+        // Since later I use <= I think I don't need to worry about odd/even
+        // Pass by copy to avoid changing input
+        const auto median = [](VectorXi A)->double
+        {
+          size_t n = A.size()/2;
+          nth_element(A.data(),A.data()+n,A.data()+A.size());
+          if(A.rows() % 2 == 1)
+          {
+            return A(n);
+          }else
+          {
+            nth_element(A.data(),A.data()+n-1,A.data()+A.size());
+            return 0.5*(A(n)+A(n-1));
+          }
+        };
+        const double med = median(SIdI);
+        VectorXi LI((I.rows()+1)/2),RI(I.rows()/2);
+        assert(LI.rows()+RI.rows() == I.rows());
+        // Distribute left and right
+        {
+          int li = 0;
+          int ri = 0;
+          for(int i = 0;i<I.rows();i++)
+          {
+            if(SIdI(i)<=med)
+            {
+              LI(li++) = I(i);
+            }else
+            {
+              RI(ri++) = I(i);
+            }
+          }
+        }
+        if(LI.rows()>0)
+        {
+          m_left = make_shared<InElementAABB>();
+          m_left->init(V,Ele,SI,LI);
+        }
+        if(RI.rows()>0)
+        {
+          m_right = make_shared<InElementAABB>();
+          m_right->init(V,Ele,SI,RI);
+        }
+      }
+  }
+}
+
+inline std::vector<int> igl::InElementAABB::find(
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & Ele, 
+    const Eigen::RowVectorXd & q, 
+    const bool first) const
+{
+  using namespace std;
+  using namespace igl;
+  using namespace Eigen;
+  bool inside = true;
+  const int dim = m_bb_max.size();
+  assert(q.size() == m_bb_max.size());
+  const double epsilon = 1e-14;
+  for(int d = 0;d<q.size()&&inside;d++)
+  {
+    inside &= (q(d)-m_bb_min(d))>=epsilon;
+    inside &= (m_bb_max(d)-q(d))>=epsilon;
+  }
+  if(!inside)
+  {
+    return std::vector<int>();
+  }
+  if(m_element != -1)
+  {
+    // Initialize to some value > -epsilon
+    double a1=1,a2=1,a3=1,a4=1;
+    switch(dim)
+    {
+      case 3:
+        {
+          // Barycentric coordinates
+          const RowVector3d V1 = V.row(Ele(m_element,0));
+          const RowVector3d V2 = V.row(Ele(m_element,1));
+          const RowVector3d V3 = V.row(Ele(m_element,2));
+          const RowVector3d V4 = V.row(Ele(m_element,3));
+          a1 = volume_single(V2,V4,V3,(RowVector3d)q);
+          a2 = volume_single(V1,V3,V4,(RowVector3d)q);
+          a3 = volume_single(V1,V4,V2,(RowVector3d)q);
+          a4 = volume_single(V1,V2,V3,(RowVector3d)q);
+          break;
+        }
+      case 2:
+        {
+          // Barycentric coordinates
+          const Vector2d V1 = V.row(Ele(m_element,0));
+          const Vector2d V2 = V.row(Ele(m_element,1));
+          const Vector2d V3 = V.row(Ele(m_element,2));
+          a1 = doublearea_single(V1,V2,(Vector2d)q);
+          a2 = doublearea_single(V2,V3,(Vector2d)q);
+          a3 = doublearea_single(V3,V1,(Vector2d)q);
+          break;
+        }
+      default:assert(false);
+    }
+    if(
+        a1>=-epsilon && 
+        a2>=-epsilon && 
+        a3>=-epsilon && 
+        a4>=-epsilon)
+    {
+      return std::vector<int>(1,m_element);
+    }else
+    {
+      return std::vector<int>();
+    }
+  }
+  std::vector<int> left = m_left->find(V,Ele,q,first);
+  if(first && !left.empty())
+  {
+    return left;
+  }
+  std::vector<int> right = m_right->find(V,Ele,q,first);
+  if(first)
+  {
+    return right;
+  }
+  left.insert(left.end(),right.begin(),right.end());
+  return left;
+}
+
+inline int igl::InElementAABB::subtree_size()
+{
+  // 1 for self
+  int n = 1;
+  int n_left = 0,n_right = 0;
+  if(m_left != NULL)
+  {
+    n_left = m_left->subtree_size();
+  }
+  if(m_right != NULL)
+  {
+    n_right = m_right->subtree_size();
+  }
+  n += 2*std::max(n_left,n_right);
+  return n;
+}
+
+
+inline void igl::InElementAABB::serialize(
+    Eigen::MatrixXd & bb_mins,
+    Eigen::MatrixXd & bb_maxs,
+    Eigen::VectorXi & elements,
+    const int i)
+{
+  using namespace std;
+  // Calling for root then resize output
+  if(i==0)
+  {
+    const int m = subtree_size();
+    //cout<<"m: "<<m<<endl;
+    bb_mins.resize(m,m_bb_min.size());
+    bb_maxs.resize(m,m_bb_max.size());
+    elements.resize(m,1);
+  }
+  //cout<<i<<" ";
+  bb_mins.row(i) = m_bb_min;
+  bb_maxs.row(i) = m_bb_max;
+  elements(i) = m_element;
+  if(m_left != NULL)
+  {
+    m_left->serialize(bb_mins,bb_maxs,elements,2*i+1);
+  }
+  if(m_right != NULL)
+  {
+    m_right->serialize(bb_mins,bb_maxs,elements,2*i+2);
+  }
+}
+#endif

include/igl/MCTables.hh

@@ -24,8 +24,8 @@
  \*===========================================================================*/
 
 //=============================================================================
-#ifndef ISOEX_MC_TABLES_HH
-#define ISOEX_MC_TABLES_HH
+#ifndef IGL_ISOEX_MC_TABLES_HH
+#define IGL_ISOEX_MC_TABLES_HH
 //=============================================================================
 
 

include/igl/MouseController.h

@@ -5,8 +5,8 @@
 // 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 MOUSECONTROLLER_H
-#define MOUSECONTROLLER_H
+#ifndef IGL_MOUSECONTROLLER_H
+#define IGL_MOUSECONTROLLER_H
 // Needs to be included before others
 #include <Eigen/StdVector>
 #include <igl/RotateWidget.h>

include/igl/OpenGL_convenience.h

@@ -5,8 +5,8 @@
 // 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 OPENGL_CONVENIENCE_H
-#define OPENGL_CONVENIENCE_H
+#ifndef IGL_OPENGL_CONVENIENCE_H
+#define IGL_OPENGL_CONVENIENCE_H
 #ifndef IGL_NO_OPENGL
 
 // Always use this:

include/igl/STR.h

@@ -5,8 +5,8 @@
 // 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 STR_H
-#define STR_H
+#ifndef IGL_STR_H
+#define IGL_STR_H
 // http://stackoverflow.com/a/2433143/148668
 #include <string>
 #include <sstream>

include/igl/SortableRow.h

@@ -12,56 +12,59 @@
 // reordering
 #include <Eigen/Core>
 
-// Templates:
-//   T  should be a matrix that implments .size(), and operator(int i)
-template <typename T>
-class SortableRow
+namespace igl
 {
-  public:
-    T data;
-  public:
-    SortableRow():data(){};
-    SortableRow(const T & data):data(data){};
-    bool operator<(const SortableRow & that) const
-    {
-      // Get reference so that I can use parenthesis
-      const SortableRow<T> & THIS = *this;
-      // Lexicographical
-      int minc = (THIS.data.size() < that.data.size()? 
-          THIS.data.size() : that.data.size());
-      // loop over columns
-      for(int i = 0;i<minc;i++)
+  // Templates:
+  //   T  should be a matrix that implments .size(), and operator(int i)
+  template <typename T>
+  class SortableRow
+  {
+    public:
+      T data;
+    public:
+      SortableRow():data(){};
+      SortableRow(const T & data):data(data){};
+      bool operator<(const SortableRow & that) const
       {
-        if(THIS.data(i) == that.data(i))
+        // Get reference so that I can use parenthesis
+        const SortableRow<T> & THIS = *this;
+        // Lexicographical
+        int minc = (THIS.data.size() < that.data.size()? 
+            THIS.data.size() : that.data.size());
+        // loop over columns
+        for(int i = 0;i<minc;i++)
         {
-          continue;
+          if(THIS.data(i) == that.data(i))
+          {
+            continue;
+          }
+          return THIS.data(i) < that.data(i);
         }
-        return THIS.data(i) < that.data(i);
-      }
-      // All characters the same, comes done to length
-      return THIS.data.size()<that.data.size();
-    };
-    bool operator==(const SortableRow & that) const
-    {
-      // Get reference so that I can use parenthesis
-      const SortableRow<T> & THIS = *this;
-      if(THIS.data.size() != that.data.size())
+        // All characters the same, comes done to length
+        return THIS.data.size()<that.data.size();
+      };
+      bool operator==(const SortableRow & that) const
       {
-        return false;
-      }
-      for(int i = 0;i<THIS.data.size();i++)
-      {
-        if(THIS.data(i) != that.data(i))
+        // Get reference so that I can use parenthesis
+        const SortableRow<T> & THIS = *this;
+        if(THIS.data.size() != that.data.size())
         {
           return false;
         }
-      }
-      return true;
-    };
-    bool operator!=(const SortableRow & that) const
-    {
-      return !(*this == that);
-    };
-};
+        for(int i = 0;i<THIS.data.size();i++)
+        {
+          if(THIS.data(i) != that.data(i))
+          {
+            return false;
+          }
+        }
+        return true;
+      };
+      bool operator!=(const SortableRow & that) const
+      {
+        return !(*this == that);
+      };
+  };
+}
 
 #endif

include/igl/Timer.h

@@ -5,177 +5,174 @@
 // 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/.
-// High Resolution Timer.
-//
-// Resolution on Mac (clock tick)
-// Resolution on Linux (1 us not tested)
-// Resolution on Windows (clock tick not tested)
-
-#ifndef IGL_TIMER_H
-#define IGL_TIMER_H
-
-#ifdef WIN32   // Windows system specific
-#include <windows.h>
-#elif __APPLE__ // Unix based system specific
-#include <mach/mach_time.h> // for mach_absolute_time
-#else
-#include <sys/time.h>
-#endif
-
-namespace igl
-{
-  class Timer
-  {
-  public:
-    // default constructor
-    Timer():
-      stopped(0),
-#ifdef WIN32
-      frequency(),
-      startCount(),
-      endCount()
-#elif __APPLE__
-      startCount(0),
-      endCount(0)
-#else
-      startCount(),
-      endCount()
-#endif
-    {
-#ifdef WIN32
-      QueryPerformanceFrequency(&frequency);
-      startCount.QuadPart = 0;
-      endCount.QuadPart = 0;
-#elif __APPLE__
-      startCount = 0;
-      endCount = 0;
-#else
-      startCount.tv_sec = startCount.tv_usec = 0;
-      endCount.tv_sec = endCount.tv_usec = 0;
-#endif
-      
-      stopped = 0;
-    }
-    // default destructor
-    ~Timer()                     
-    {
-      
-    }
-    
-#ifdef __APPLE__
-    //Raw mach_absolute_times going in, difference in seconds out
-    double subtractTimes( uint64_t endTime, uint64_t startTime )
-    {
-      uint64_t difference = endTime - startTime;
-      static double conversion = 0.0;
-      
-      if( conversion == 0.0 )
-      {
-        mach_timebase_info_data_t info;
-        kern_return_t err = mach_timebase_info( &info );
-        
-        //Convert the timebase into seconds
-        if( err == 0  )
-          conversion = 1e-9 * (double) info.numer / (double) info.denom;
-      }
-      
-      return conversion * (double) difference;
-    }
-#endif
-    
-    // start timer
-    void   start()               
-    {
-      stopped = 0; // reset stop flag
-#ifdef WIN32
-      QueryPerformanceCounter(&startCount);
-#elif __APPLE__
-      startCount = mach_absolute_time();
-#else
-      gettimeofday(&startCount, NULL);
-#endif
-      
-    }
-    
-    // stop the timer
-    void   stop()                
-    {
-      stopped = 1; // set timer stopped flag
-      
-#ifdef WIN32
-      QueryPerformanceCounter(&endCount);
-#elif __APPLE__
-      endCount = mach_absolute_time();
-#else
-      gettimeofday(&endCount, NULL);
-#endif
-      
-    }
-    // get elapsed time in second
-    double getElapsedTime()      
-    {
-      return this->getElapsedTimeInSec();
-    }
-    // get elapsed time in second (same as getElapsedTime)
-    double getElapsedTimeInSec() 
-    {
-      return this->getElapsedTimeInMicroSec() * 0.000001;
-    }
-    
-    // get elapsed time in milli-second
-    double getElapsedTimeInMilliSec()
-    {
-      return this->getElapsedTimeInMicroSec() * 0.001;
-    }
-    // get elapsed time in micro-second
-    double getElapsedTimeInMicroSec()          
-    {
-      double startTimeInMicroSec = 0;
-      double endTimeInMicroSec = 0;
-
-#ifdef WIN32
-      if(!stopped)
-        QueryPerformanceCounter(&endCount);
-      
-      startTimeInMicroSec = 
-        startCount.QuadPart * (1000000.0 / frequency.QuadPart);
-      endTimeInMicroSec = endCount.QuadPart * (1000000.0 / frequency.QuadPart);
-#elif __APPLE__
-      if (!stopped)
-        endCount = mach_absolute_time();
-      
-      return subtractTimes(endCount,startCount)/1e-6;
-#else
-      if(!stopped)
-        gettimeofday(&endCount, NULL);
-      
-      startTimeInMicroSec = 
-        (startCount.tv_sec * 1000000.0) + startCount.tv_usec;
-      endTimeInMicroSec = (endCount.tv_sec * 1000000.0) + endCount.tv_usec;
-#endif
-      
-      return endTimeInMicroSec - startTimeInMicroSec;
-    }
-    
-    
-  protected:
-    
-    
-  private:
-    // stop flag 
-    int    stopped;               
-#ifdef WIN32
-    // ticks per second
-    LARGE_INTEGER frequency;      
-    LARGE_INTEGER startCount;     
-    LARGE_INTEGER endCount;       
-#elif __APPLE__
-    uint64_t startCount;           
-    uint64_t endCount;             
-#else
-    timeval startCount;           
-    timeval endCount;             
-#endif
-  };
-}
-#endif // TIMER_H_DEF
+// High Resolution Timer.
+//
+// Resolution on Mac (clock tick)
+// Resolution on Linux (1 us not tested)
+// Resolution on Windows (clock tick not tested)
+
+#ifndef IGL_TIMER_H
+#define IGL_TIMER_H
+
+#ifdef WIN32   // Windows system specific
+#include <windows.h>
+#elif __APPLE__ // Unix based system specific
+#include <mach/mach_time.h> // for mach_absolute_time
+#else
+#include <sys/time.h>
+#endif
+
+namespace igl
+{
+  class Timer
+  {
+  public:
+    // default constructor
+    Timer():
+      stopped(0),
+#ifdef WIN32
+      frequency(),
+      startCount(),
+      endCount()
+#elif __APPLE__
+      startCount(0),
+      endCount(0)
+#else
+      startCount(),
+      endCount()
+#endif
+    {
+#ifdef WIN32
+      QueryPerformanceFrequency(&frequency);
+      startCount.QuadPart = 0;
+      endCount.QuadPart = 0;
+#elif __APPLE__
+      startCount = 0;
+      endCount = 0;
+#else
+      startCount.tv_sec = startCount.tv_usec = 0;
+      endCount.tv_sec = endCount.tv_usec = 0;
+#endif
+
+      stopped = 0;
+    }
+    // default destructor
+    ~Timer()                     
+    {
+
+    }
+
+#ifdef __APPLE__
+    //Raw mach_absolute_times going in, difference in seconds out
+    double subtractTimes( uint64_t endTime, uint64_t startTime )
+    {
+      uint64_t difference = endTime - startTime;
+      static double conversion = 0.0;
+
+      if( conversion == 0.0 )
+      {
+        mach_timebase_info_data_t info;
+        kern_return_t err = mach_timebase_info( &info );
+
+        //Convert the timebase into seconds
+        if( err == 0  )
+          conversion = 1e-9 * (double) info.numer / (double) info.denom;
+      }
+
+      return conversion * (double) difference;
+    }
+#endif
+
+    // start timer
+    void   start()               
+    {
+      stopped = 0; // reset stop flag
+#ifdef WIN32
+      QueryPerformanceCounter(&startCount);
+#elif __APPLE__
+      startCount = mach_absolute_time();
+#else
+      gettimeofday(&startCount, NULL);
+#endif
+
+    }
+
+    // stop the timer
+    void   stop()                
+    {
+      stopped = 1; // set timer stopped flag
+
+#ifdef WIN32
+      QueryPerformanceCounter(&endCount);
+#elif __APPLE__
+      endCount = mach_absolute_time();
+#else
+      gettimeofday(&endCount, NULL);
+#endif
+
+    }
+    // get elapsed time in second
+    double getElapsedTime()      
+    {
+      return this->getElapsedTimeInSec();
+    }
+    // get elapsed time in second (same as getElapsedTime)
+    double getElapsedTimeInSec() 
+    {
+      return this->getElapsedTimeInMicroSec() * 0.000001;
+    }
+
+    // get elapsed time in milli-second
+    double getElapsedTimeInMilliSec()
+    {
+      return this->getElapsedTimeInMicroSec() * 0.001;
+    }
+    // get elapsed time in micro-second
+    double getElapsedTimeInMicroSec()          
+    {
+      double startTimeInMicroSec = 0;
+      double endTimeInMicroSec = 0;
+
+#ifdef WIN32
+      if(!stopped)
+        QueryPerformanceCounter(&endCount);
+
+      startTimeInMicroSec = 
+        startCount.QuadPart * (1000000.0 / frequency.QuadPart);
+      endTimeInMicroSec = endCount.QuadPart * (1000000.0 / frequency.QuadPart);
+#elif __APPLE__
+      if (!stopped)
+        endCount = mach_absolute_time();
+
+      return subtractTimes(endCount,startCount)/1e-6;
+#else
+      if(!stopped)
+        gettimeofday(&endCount, NULL);
+
+      startTimeInMicroSec = 
+        (startCount.tv_sec * 1000000.0) + startCount.tv_usec;
+      endTimeInMicroSec = (endCount.tv_sec * 1000000.0) + endCount.tv_usec;
+#endif
+
+      return endTimeInMicroSec - startTimeInMicroSec;
+    }
+
+  private:
+    // stop flag 
+    int    stopped;               
+#ifdef WIN32
+    // ticks per second
+    LARGE_INTEGER frequency;      
+    LARGE_INTEGER startCount;     
+    LARGE_INTEGER endCount;       
+#elif __APPLE__
+    uint64_t startCount;           
+    uint64_t endCount;             
+#else
+    timeval startCount;           
+    timeval endCount;             
+#endif
+  };
+}
+#endif // TIMER_H_DEF
+

include/igl/WindingNumberAABB.h

@@ -0,0 +1,340 @@
+#ifndef IGL_WINDINGNUMBERAABB_H
+#define IGL_WINDINGNUMBERAABB_H
+#include "WindingNumberTree.h"
+
+namespace igl
+{
+  template <typename Point>
+  class WindingNumberAABB : public WindingNumberTree<Point>
+  {
+    protected:
+      Point min_corner;
+      Point max_corner;
+      double total_positive_area;
+    public: 
+      enum SplitMethod
+      {
+        CENTER_ON_LONGEST_AXIS = 0,
+        MEDIAN_ON_LONGEST_AXIS = 1,
+        NUM_SPLIT_METHODS = 3
+      } split_method;
+    public:
+      inline WindingNumberAABB(
+        const Eigen::MatrixXd & V,
+        const Eigen::MatrixXi & F);
+      inline WindingNumberAABB(
+        const WindingNumberTree<Point> & parent,
+        const Eigen::MatrixXi & F);
+      // Initialize some things
+      inline void init();
+      inline bool inside(const Point & p) const;
+      inline virtual void grow();
+      // Compute min and max corners
+      inline void compute_min_max_corners();
+      inline double max_abs_winding_number(const Point & p) const;
+      inline double max_simple_abs_winding_number(const Point & p) const;
+  };
+}
+
+// Implementation
+
+#include "winding_number.h"
+
+#include "barycenter.h"
+#include "median.h"
+#include "doublearea.h"
+#include "per_face_normals.h"
+
+#include <limits>
+#include <vector>
+#include <iostream>
+
+// Minimum number of faces in a hierarchy element (this is probably dependent
+// on speed of machine and compiler optimization)
+#ifndef WindingNumberAABB_MIN_F
+#  define WindingNumberAABB_MIN_F 100
+#endif
+
+template <typename Point>
+inline void igl::WindingNumberAABB<Point>::init()
+{
+  using namespace igl;
+  using namespace Eigen;
+  assert(max_corner.size() == 3);
+  assert(min_corner.size() == 3);
+  compute_min_max_corners();
+  VectorXd dblA;
+  doublearea(this->getV(),this->getF(),dblA);
+  total_positive_area = dblA.sum()/2.0;
+}
+
+template <typename Point>
+inline igl::WindingNumberAABB<Point>::WindingNumberAABB(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F):
+  WindingNumberTree<Point>(V,F),
+  min_corner(),
+  max_corner(),
+  total_positive_area(std::numeric_limits<double>::infinity()),
+  split_method(MEDIAN_ON_LONGEST_AXIS)
+{
+  init();
+}
+
+template <typename Point>
+inline igl::WindingNumberAABB<Point>::WindingNumberAABB(
+  const WindingNumberTree<Point> & parent,
+  const Eigen::MatrixXi & F):
+  WindingNumberTree<Point>(parent,F),
+  min_corner(),
+  max_corner(),
+  total_positive_area(std::numeric_limits<double>::infinity()),
+  split_method(MEDIAN_ON_LONGEST_AXIS)
+{
+  init();
+}
+
+template <typename Point>
+inline void igl::WindingNumberAABB<Point>::grow()
+{
+  using namespace std;
+  using namespace Eigen;
+  using namespace igl;
+  //cout<<"cap.rows(): "<<this->getcap().rows()<<endl;
+  //cout<<"F.rows(): "<<this->getF().rows()<<endl;
+
+  // Base cases
+  if(
+    this->getF().rows() <= (WindingNumberAABB_MIN_F>0?WindingNumberAABB_MIN_F:0) ||
+    (this->getcap().rows() - 2) >= this->getF().rows())
+  {
+    // Don't grow
+    return;
+  }
+
+  // Compute longest direction
+  int max_d = -1;
+  double max_len = -numeric_limits<double>::infinity();
+  for(int d = 0;d<min_corner.size();d++)
+  {
+    if( (max_corner[d] - min_corner[d]) > max_len )
+    {
+      max_len = (max_corner[d] - min_corner[d]);
+      max_d = d;
+    }
+  }
+  // Compute facet barycenters
+  MatrixXd BC;
+  barycenter(this->getV(),this->getF(),BC);
+
+
+  // Blerg, why is selecting rows so difficult
+
+  vector<int> id( this->getF().rows());
+  double split_value;
+  // Split in longest direction
+  switch(split_method)
+  {
+    case MEDIAN_ON_LONGEST_AXIS:
+      // Determine median
+      median(BC.col(max_d),split_value);
+      break;
+    default:
+      assert(false);
+    case CENTER_ON_LONGEST_AXIS:
+      split_value = 0.5*(max_corner[max_d] + min_corner[max_d]);
+      break;
+  }
+  //cout<<"c: "<<0.5*(max_corner[max_d] + min_corner[max_d])<<" "<<
+  //  "m: "<<split_value<<endl;;
+
+  for(int i = 0;i<this->getF().rows();i++)
+  {
+    if(BC(i,max_d) <= split_value)
+    {
+      id[i] = 0; //left
+    }else
+    {
+      id[i] = 1; //right
+    }
+  }
+
+  const int lefts = (int) count(id.begin(),id.end(),0);
+  const int rights = (int) count(id.begin(),id.end(),1);
+  if(lefts == 0 || rights == 0)
+  {
+    // badly balanced base case (could try to recut)
+    return;
+  }
+  assert(lefts+rights == this->getF().rows());
+  MatrixXi leftF(lefts,  this->getF().cols());
+  MatrixXi rightF(rights,this->getF().cols());
+  int left_i = 0;
+  int right_i = 0;
+  for(int i = 0;i<this->getF().rows();i++)
+  {
+    if(id[i] == 0)
+    {
+      leftF.row(left_i++) = this->getF().row(i);
+    }else if(id[i] == 1)
+    {
+      rightF.row(right_i++) = this->getF().row(i);
+    }else
+    {
+      assert(false);
+    }
+  }
+  assert(right_i == rightF.rows());
+  assert(left_i == leftF.rows());
+  // Finally actually grow children and Recursively grow
+  WindingNumberAABB<Point> * leftWindingNumberAABB = new WindingNumberAABB<Point>(*this,leftF);
+  leftWindingNumberAABB->grow();
+  this->children.push_back(leftWindingNumberAABB);
+  WindingNumberAABB<Point> * rightWindingNumberAABB = new WindingNumberAABB<Point>(*this,rightF);
+  rightWindingNumberAABB->grow();
+  this->children.push_back(rightWindingNumberAABB);
+}
+
+template <typename Point>
+inline bool igl::WindingNumberAABB<Point>::inside(const Point & p) const
+{
+  assert(p.size() == max_corner.size());
+  assert(p.size() == min_corner.size());
+  for(int i = 0;i<p.size();i++)
+  {
+    if( p(i) < min_corner(i) || p(i) >= max_corner(i))
+    {
+      return false;
+    }
+  }
+  return true;
+}
+
+template <typename Point>
+inline void igl::WindingNumberAABB<Point>::compute_min_max_corners()
+{
+  using namespace std;
+  // initialize corners
+  for(int d = 0;d<min_corner.size();d++)
+  {
+    min_corner[d] =  numeric_limits<double>::infinity();
+    max_corner[d] = -numeric_limits<double>::infinity();
+  }
+
+  this->center = Point(0,0,0);
+  // Loop over facets
+  for(int i = 0;i<this->getF().rows();i++)
+  {
+    for(int j = 0;j<this->getF().cols();j++)
+    {
+      for(int d = 0;d<min_corner.size();d++)
+      {
+        min_corner[d] = 
+          this->getV()(this->getF()(i,j),d) < min_corner[d] ?  
+            this->getV()(this->getF()(i,j),d) : min_corner[d];
+        max_corner[d] = 
+          this->getV()(this->getF()(i,j),d) > max_corner[d] ?  
+            this->getV()(this->getF()(i,j),d) : max_corner[d];
+      }
+      // This is biased toward vertices incident on more than one face, but
+      // perhaps that's good
+      this->center += this->getV().row(this->getF()(i,j));
+    }
+  }
+  // Average
+  this->center.array() /= this->getF().size();
+
+  //cout<<"min_corner: "<<this->min_corner.transpose()<<endl;
+  //cout<<"Center: "<<this->center.transpose()<<endl;
+  //cout<<"max_corner: "<<this->max_corner.transpose()<<endl;
+  //cout<<"Diag center: "<<((this->max_corner + this->min_corner)*0.5).transpose()<<endl;
+  //cout<<endl;
+
+  this->radius = (max_corner-min_corner).norm()/2.0;
+}
+
+template <typename Point>
+inline double igl::WindingNumberAABB<Point>::max_abs_winding_number(const Point & p) const
+{
+  using namespace std;
+  // Only valid if not inside
+  if(inside(p))
+  {
+    return numeric_limits<double>::infinity();
+  }
+  // Q: we know the total positive area so what's the most this could project
+  // to? Remember it could be layered in the same direction.
+  return numeric_limits<double>::infinity();
+}
+
+template <typename Point>
+inline double igl::WindingNumberAABB<Point>::max_simple_abs_winding_number(const Point & p) const
+{
+  using namespace std;
+  using namespace Eigen;
+  using namespace igl;
+  // Only valid if not inside
+  if(inside(p))
+  {
+    return numeric_limits<double>::infinity();
+  }
+  // Max simple is the same as sum of positive winding number contributions of
+  // bounding box
+
+  // begin precomputation
+  //MatrixXd BV((int)pow(2,3),3);
+  MatrixXd BV((int)(1<<3),3);
+  BV <<
+    min_corner[0],min_corner[1],min_corner[2],
+    min_corner[0],min_corner[1],max_corner[2],
+    min_corner[0],max_corner[1],min_corner[2],
+    min_corner[0],max_corner[1],max_corner[2],
+    max_corner[0],min_corner[1],min_corner[2],
+    max_corner[0],min_corner[1],max_corner[2],
+    max_corner[0],max_corner[1],min_corner[2],
+    max_corner[0],max_corner[1],max_corner[2];
+  MatrixXi BF(2*2*3,3);
+  BF <<
+    0,6,4,
+    0,2,6,
+    0,3,2,
+    0,1,3,
+    2,7,6,
+    2,3,7,
+    4,6,7,
+    4,7,5,
+    0,4,5,
+    0,5,1,
+    1,5,7,
+    1,7,3;
+  MatrixXd BFN;
+  per_face_normals(BV,BF,BFN);
+  // end of precomputation
+
+  // Only keep those with positive dot products
+  MatrixXi PBF(BF.rows(),BF.cols());
+  int pbfi = 0;
+  Point p2c = 0.5*(min_corner+max_corner)-p;
+  for(int i = 0;i<BFN.rows();i++)
+  {
+    if(p2c.dot(BFN.row(i)) > 0)
+    {
+      PBF.row(pbfi++) = BF.row(i);
+    }
+  }
+  PBF.conservativeResize(pbfi,PBF.cols());
+  double w = numeric_limits<double>::infinity();
+  igl::winding_number_3(
+    BV.data(),
+    BV.rows(),
+    PBF.data(),
+    PBF.rows(),
+    p.data(),
+    1,
+    &w);
+  return w;
+}
+
+//// Explicit instanciation
+//template class igl::WindingNumberAABB<Eigen::Vector3d >;
+#endif

include/igl/WindingNumberMethod.h

@@ -0,0 +1,13 @@
+#ifndef IGL_WINDINGNUMBERMETHOD_H
+#define IGL_WINDINGNUMBERMETHOD_H
+namespace igl
+{
+  enum WindingNumberMethod
+  {
+    EXACT_WINDING_NUMBER_METHOD = 0, // Exact hierarchical evaluation
+    APPROX_SIMPLE_WINDING_NUMBER_METHOD = 1,
+    APPROX_CACHE_WINDING_NUMBER_METHOD = 2, // Approximate hierarchical evaluation
+    NUM_WINDING_NUMBER_METHODS = 3
+  };
+}
+#endif

include/igl/WindingNumberTree.h

@@ -0,0 +1,455 @@
+#ifndef IGL_BOUNDINGTREE_H
+#define IGL_BOUNDINGTREE_H
+#include <list>
+#include <map>
+#include <Eigen/Dense>
+#include "WindingNumberMethod.h"
+
+namespace igl
+{
+  // Templates:
+  //   Point  type for points in space, e.g. Eigen::Vector3d
+  template <typename Point>
+  class WindingNumberTree
+  {
+    public:
+      // Method to use (see enum above)
+      //static double min_max_w;
+      static std::map< 
+        std::pair<const WindingNumberTree*,const WindingNumberTree*>, double>
+          cached;
+    protected:
+      WindingNumberMethod method;
+      const WindingNumberTree * parent;
+      std::list<WindingNumberTree * > children;
+      //// List of boundary edges (recall edges are vertices in 2d)
+      //const Eigen::MatrixXi boundary;
+      // Base mesh vertices
+      Eigen::MatrixXd & V;
+      // Base mesh vertices with duplicates removed
+      Eigen::MatrixXd SV;
+      // Facets in this bounding volume
+      Eigen::MatrixXi F;
+      // Tesselated boundary curve
+      Eigen::MatrixXi cap;
+      // Upper Bound on radius of enclosing ball
+      double radius;
+      // (Approximate) center (of mass)
+      Point center;
+    public:
+      // For root
+      inline WindingNumberTree(
+        const Eigen::MatrixXd & V,
+        const Eigen::MatrixXi & F);
+      // For chilluns 
+      inline WindingNumberTree(
+        const WindingNumberTree<Point> & parent,
+        const Eigen::MatrixXi & F);
+      inline virtual ~WindingNumberTree();
+      // Set method
+      inline void set_method( const WindingNumberMethod & m);
+    public:
+      inline const Eigen::MatrixXd & getV() const;
+      inline const Eigen::MatrixXi & getF() const;
+      inline const Eigen::MatrixXi & getcap() const;
+      // Grow the Tree recursively
+      inline virtual void grow();
+      // Determine whether a given point is inside the bounding 
+      //
+      // Inputs:
+      //   p  query point 
+      // Returns true if the point p is inside this bounding volume
+      inline virtual bool inside(const Point & p) const;
+      // Compute the (partial) winding number of a given point p
+      // According to method
+      //  
+      // Inputs:
+      //   p  query point 
+      // Returns winding number 
+      inline double winding_number(const Point & p) const;
+      // Same as above, but always computes winding number using exact method
+      // (sum over every facet)
+      inline double winding_number_all(const Point & p) const;
+      // Same as above, but always computes using sum over tesslated boundary
+      inline double winding_number_boundary(const Point & p) const;
+      //// Same as winding_number above, but if max_simple_abs_winding_number is
+      //// less than some threshold min_max_w just return 0 (colloquially the "fast
+      //// multipole method)
+      ////
+      ////
+      //// Inputs:
+      ////   p  query point 
+      ////   min_max_w  minimum max simple w to be processed
+      //// Returns approximate winding number
+      //double winding_number_approx_simple(
+      //  const Point & p, 
+      //  const double min_max_w);
+      // Print contents of Tree
+      //
+      // Optional input:
+      //   tab  tab to show depth
+      inline void print(const char * tab="");
+      // Determine max absolute winding number
+      //
+      // Inputs:
+      //   p  query point 
+      // Returns max winding number of 
+      inline virtual double max_abs_winding_number(const Point & p) const; 
+      // Same as above, but stronger assumptions on (V,F). Assumes (V,F) is a
+      // simple polyhedron
+      inline virtual double max_simple_abs_winding_number(const Point & p) const;
+      // Compute or read cached winding number for point p with respect to mesh
+      // in bounding box, recursing according to approximation criteria
+      //
+      // Inputs:
+      //   p  query point 
+      //   that  WindingNumberTree containing mesh w.r.t. which we're computing w.n.
+      // Returns cached winding number
+      inline virtual double cached_winding_number(const WindingNumberTree & that, const Point & p) const;
+  };
+}
+
+// Implementation
+
+#include "WindingNumberTree.h"
+#include "winding_number.h"
+#include "triangle_fan.h"
+#include "exterior_edges.h"
+
+#include <igl/PI.h>
+#include <igl/remove_duplicate_vertices.h>
+
+#include <iostream>
+#include <limits>
+
+//template <typename Point>
+//WindingNumberMethod WindingNumberTree<Point>::method = EXACT_WINDING_NUMBER_METHOD;
+//template <typename Point>
+//double WindingNumberTree<Point>::min_max_w = 0;
+template <typename Point>
+std::map< std::pair<const igl::WindingNumberTree<Point>*,const igl::WindingNumberTree<Point>*>, double>
+  igl::WindingNumberTree<Point>::cached;
+
+static Eigen::MatrixXd dummyV;
+template <typename Point>
+inline igl::WindingNumberTree<Point>::WindingNumberTree(
+  const Eigen::MatrixXd & _V,
+  const Eigen::MatrixXi & _F):
+  method(EXACT_WINDING_NUMBER_METHOD),
+  parent(NULL),
+  V(dummyV),
+  SV(),
+  F(),
+  //boundary(igl::boundary_faces<Eigen::MatrixXi,Eigen::MatrixXi>(F))
+  cap(),
+  radius(std::numeric_limits<double>::infinity()),
+  center(0,0,0)
+{
+  // Remove any exactly duplicate vertices
+  // Q: Can this ever increase the complexity of the boundary?
+  // Q: Would we gain even more by remove almost exactly duplicate vertices?
+  Eigen::MatrixXi SF,SVI,SVJ;
+  igl::remove_duplicate_vertices(_V,_F,0.0,SV,SVI,SVJ,F);
+  triangle_fan(exterior_edges(F),cap);
+  V = SV;
+}
+
+template <typename Point>
+inline igl::WindingNumberTree<Point>::WindingNumberTree(
+  const igl::WindingNumberTree<Point> & parent,
+  const Eigen::MatrixXi & _F):
+  method(parent.method),
+  parent(&parent),
+  V(parent.V),
+  SV(),
+  F(_F),
+  cap(triangle_fan(exterior_edges(_F)))
+{
+}
+
+template <typename Point>
+inline igl::WindingNumberTree<Point>::~WindingNumberTree()
+{
+  using namespace std;
+  // Delete children
+  typename list<WindingNumberTree<Point>* >::iterator cit = children.begin();
+  while(cit != children.end())
+  {
+    // clear the memory of this item
+    delete (* cit);
+    // erase from list, returns next element in iterator
+    cit = children.erase(cit);
+  }
+}
+      
+template <typename Point>
+inline void igl::WindingNumberTree<Point>::set_method(const WindingNumberMethod & m)
+{
+  this->method = m;
+  for(auto child : children)
+  {
+    child->set_method(m);
+  }
+}
+
+template <typename Point>
+inline const Eigen::MatrixXd & igl::WindingNumberTree<Point>::getV() const
+{
+  return V;
+}
+
+template <typename Point>
+inline const Eigen::MatrixXi & igl::WindingNumberTree<Point>::getF() const
+{
+  return F;
+}
+
+template <typename Point>
+inline const Eigen::MatrixXi & igl::WindingNumberTree<Point>::getcap() const
+{
+  return cap;
+}
+
+template <typename Point>
+inline void igl::WindingNumberTree<Point>::grow()
+{
+  // Don't grow
+  return;
+}
+
+template <typename Point>
+inline bool igl::WindingNumberTree<Point>::inside(const Point & /*p*/) const
+{
+  return true;
+}
+
+template <typename Point>
+inline double igl::WindingNumberTree<Point>::winding_number(const Point & p) const
+{
+  using namespace std;
+  //cout<<"+"<<boundary.rows();
+  // If inside then we need to be careful
+  if(inside(p))
+  {
+    // If not a leaf then recurse
+    if(children.size()>0)
+    {
+      // Recurse on each child and accumulate
+      double sum = 0;
+      for(
+        typename list<WindingNumberTree<Point>* >::const_iterator cit = children.begin();
+        cit != children.end();
+        cit++)
+      {
+        switch(method)
+        {
+          case EXACT_WINDING_NUMBER_METHOD:
+            sum += (*cit)->winding_number(p);
+            break;
+          case APPROX_SIMPLE_WINDING_NUMBER_METHOD:
+          case APPROX_CACHE_WINDING_NUMBER_METHOD:
+            //if((*cit)->max_simple_abs_winding_number(p) > min_max_w)
+            //{
+              sum += (*cit)->winding_number(p);
+            //}
+            break;
+          default:
+            assert(false);
+            break;
+        }
+      }
+      return sum;
+    }else
+    {
+      return winding_number_all(p);
+    }
+  }else{
+    // Otherwise we can just consider boundary
+    // Q: If we using the "multipole" method should we also subdivide the
+    // boundary case?
+    if((cap.rows() - 2) < F.rows())
+    {
+      switch(method)
+      {
+        case EXACT_WINDING_NUMBER_METHOD:
+          return winding_number_boundary(p);
+        case APPROX_SIMPLE_WINDING_NUMBER_METHOD:
+        {
+          double dist = (p-center).norm();
+          // Radius is already an overestimate of inside
+          if(dist>1.0*radius)
+          {
+            return 0;
+          }else
+          {
+            return winding_number_boundary(p);
+          }
+        }
+        case APPROX_CACHE_WINDING_NUMBER_METHOD:
+        {
+          return parent->cached_winding_number(*this,p);
+        }
+        default: assert(false);break;
+      }
+    }else
+    {
+      // doesn't pay off to use boundary
+      return winding_number_all(p);
+    }
+  }
+  return 0;
+}
+
+template <typename Point>
+inline double igl::WindingNumberTree<Point>::winding_number_all(const Point & p) const
+{
+  double w = 0;
+  igl::winding_number_3(
+    V.data(),
+    V.rows(),
+    F.data(),
+    F.rows(),
+    p.data(),
+    1,
+    &w);
+  return w;
+}
+
+template <typename Point>
+inline double igl::WindingNumberTree<Point>::winding_number_boundary(const Point & p) const
+{
+  using namespace Eigen;
+  using namespace std;
+
+  double w = 0;
+  // `cap` is already flipped inside out, so we don't need to flip sign of w
+  igl::winding_number_3(
+    V.data(),
+    V.rows(),
+    cap.data(),
+    cap.rows(),
+    &p[0],
+    1,
+    &w);
+  return w;
+}
+
+//template <typename Point>
+//inline double igl::WindingNumberTree<Point>::winding_number_approx_simple(
+//  const Point & p, 
+//  const double min_max_w)
+//{
+//  using namespace std;
+//  if(max_simple_abs_winding_number(p) > min_max_w)
+//  {
+//    return winding_number(p);
+//  }else
+//  {
+//    cout<<"Skipped! "<<max_simple_abs_winding_number(p)<<"<"<<min_max_w<<endl;
+//    return 0;
+//  }
+//}
+
+template <typename Point>
+inline void igl::WindingNumberTree<Point>::print(const char * tab)
+{
+  using namespace std;
+  // Print all facets
+  cout<<tab<<"["<<endl<<F<<endl<<"]";
+  // Print children
+  for(
+      typename list<WindingNumberTree<Point>* >::iterator cit = children.begin();
+      cit != children.end();
+      cit++)
+  {
+    cout<<","<<endl;
+    (*cit)->print((string(tab)+"").c_str());
+  }
+}
+
+template <typename Point>
+inline double igl::WindingNumberTree<Point>::max_abs_winding_number(const Point & p) const
+{
+  return std::numeric_limits<double>::infinity();
+}
+
+template <typename Point>
+inline double igl::WindingNumberTree<Point>::max_simple_abs_winding_number(const Point & p) const
+{
+  using namespace std;
+  return numeric_limits<double>::infinity();
+}
+
+template <typename Point>
+inline double igl::WindingNumberTree<Point>::cached_winding_number(
+  const igl::WindingNumberTree<Point> & that,
+  const Point & p) const
+{
+  using namespace std;
+  using namespace igl;
+  // Simple metric for "far".
+  //   this             that
+  //                   --------
+  //   -----          /   |    \ .
+  //  /  r  \        /    R     \ .
+  // | p !   |      |     !      |
+  //  \_____/        \          /
+  //                  \________/
+  //
+  // 
+  // a = angle formed by trapazoid formed by raising sides with lengths r and R
+  // at respective centers.
+  //
+  // a = atan2(R-r,d), where d is the distance between centers
+
+  // That should be bigger (what about parent? what about sister?)
+  bool far = this->radius<that.radius;
+  if(far)
+  {
+    double a = atan2(
+      that.radius - this->radius,
+      (that.center - this->center).norm());
+    assert(a>0);
+    far = (a<PI/8.0);
+  }
+
+  if(far)
+  {
+    // Not implemented yet
+    pair<const WindingNumberTree*,const WindingNumberTree*> this_that(this,&that);
+    // Need to compute it for first time?
+    if(cached.count(this_that)==0)
+    {
+      cached[this_that] = 
+        that.winding_number_boundary(this->center);
+    }
+    return cached[this_that];
+  }else if(children.size() == 0)
+  {
+    // not far and hierarchy ended too soon: can't use cache
+    return that.winding_number_boundary(p);
+  }else
+  {
+    for(
+      typename list<WindingNumberTree<Point>* >::const_iterator cit = children.begin();
+      cit != children.end();
+      cit++)
+    {
+      if((*cit)->inside(p))
+      {
+        return (*cit)->cached_winding_number(that,p);
+      }
+    }
+    // Not inside any children? This can totally happen because bounding boxes
+    // are set to bound contained facets. So sibilings may overlap and their
+    // union may not contain their parent (though, their union is certainly a
+    // subset of their parent).
+    assert(false);
+  }
+  return 0;
+}
+
+// Explicit instanciation
+//template class igl::WindingNumberTree<Eigen::Vector3d >;
+
+#endif

include/igl/active_set.h

@@ -33,8 +33,8 @@ namespace igl
   //   Y  list of fixed values corresponding to known rows in Z
   //   Aeq  meq by n list of linear equality constraint coefficients
   //   Beq  meq by 1 list of linear equality constraint constant values
-  //   Aieq  mieq by n list of linear equality constraint coefficients
-  //   Bieq  mieq by 1 list of linear equality constraint constant values
+  //   Aieq  mieq by n list of linear inequality constraint coefficients
+  //   Bieq  mieq by 1 list of linear inequality constraint constant values
   //   lx  n by 1 list of lower bounds [] implies -Inf
   //   ux  n by 1 list of upper bounds [] implies Inf
   //   params  struct of additional parameters (see below)

include/igl/add_barycenter.cpp

@@ -19,34 +19,34 @@ IGL_INLINE void igl::add_barycenter(
     Eigen::PlainObjectBase<Index> & FD)
 {
   using namespace Eigen;
-	// Compute face barycenter
-	Eigen::MatrixXd BC;
-	igl::barycenter(V,F,BC);
-
-	// Add the barycenters to the vertices
-	VD.resize(V.rows()+F.rows(),3);
-	VD.block(0,0,V.rows(),3) = V;
-	VD.block(V.rows(),0,F.rows(),3) = BC;
-
-	// Each face is split four ways
-	FD.resize(F.rows()*3,3);
-
-	for (unsigned i=0; i<F.rows(); ++i)
-	{
-		int i0 = F(i,0);
-		int i1 = F(i,1);
-		int i2 = F(i,2);
-		int i3 = V.rows() + i;
-
-		Vector3i F0,F1,F2;
-		F0 << i0,i1,i3;
-		F1 << i1,i2,i3;
-		F2 << i2,i0,i3;
-
-		FD.row(i*3 + 0) = F0;
-		FD.row(i*3 + 1) = F1;
-		FD.row(i*3 + 2) = F2;
-	}
+  // Compute face barycenter
+  Eigen::MatrixXd BC;
+  igl::barycenter(V,F,BC);
+
+  // Add the barycenters to the vertices
+  VD.resize(V.rows()+F.rows(),3);
+  VD.block(0,0,V.rows(),3) = V;
+  VD.block(V.rows(),0,F.rows(),3) = BC;
+
+  // Each face is split four ways
+  FD.resize(F.rows()*3,3);
+
+  for (unsigned i=0; i<F.rows(); ++i)
+  {
+    int i0 = F(i,0);
+    int i1 = F(i,1);
+    int i2 = F(i,2);
+    int i3 = V.rows() + i;
+
+    Vector3i F0,F1,F2;
+    F0 << i0,i1,i3;
+    F1 << i1,i2,i3;
+    F2 << i2,i0,i3;
+
+    FD.row(i*3 + 0) = F0;
+    FD.row(i*3 + 1) = F1;
+    FD.row(i*3 + 2) = F2;
+  }
 
 
 }

include/igl/angular_distance.cpp

@@ -8,7 +8,7 @@
 #include "angular_distance.h"
 #include <igl/EPS.h>
 #include <igl/PI.h>
-double igl::angular_distance(
+IGL_INLINE double igl::angular_distance(
   const Eigen::Quaterniond & A,
   const Eigen::Quaterniond & B)
 {

include/igl/angular_distance.h

@@ -18,7 +18,7 @@ namespace igl
   //   A  unit quaternion
   //   B  unit quaternion
   // Returns angular distance
-  double angular_distance(
+  IGL_INLINE double angular_distance(
     const Eigen::Quaterniond & A,
     const Eigen::Quaterniond & B);
 }

include/igl/arap_rhs.cpp

@@ -14,13 +14,14 @@
 IGL_INLINE void igl::arap_rhs(
   const Eigen::MatrixXd & V, 
   const Eigen::MatrixXi & F,
+  const int dim,
   const igl::ARAPEnergyType energy,
   Eigen::SparseMatrix<double>& K)
 {
   using namespace igl;
   using namespace Eigen;
   // Number of dimensions
-  int dim = V.cols();
+  int Vdim = V.cols();
   //// Number of mesh vertices
   //int n = V.rows();
   //// Number of mesh elements
@@ -41,7 +42,7 @@ IGL_INLINE void igl::arap_rhs(
     default:
       fprintf(
         stderr,
-        "covariance_scatter_matrix.h: Error: Unsupported arap energy %d\n",
+        "arap_rhs.h: Error: Unsupported arap energy %d\n",
         energy);
       return;
   }
@@ -49,19 +50,37 @@ IGL_INLINE void igl::arap_rhs(
   SparseMatrix<double> KX,KY,KZ;
   arap_linear_block(V,F,0,energy,KX);
   arap_linear_block(V,F,1,energy,KY);
-  if(dim == 2)
+  if(Vdim == 2)
   {
     K = cat(2,repdiag(KX,dim),repdiag(KY,dim));
-  }else if(dim == 3)
+  }else if(Vdim == 3)
   {
     arap_linear_block(V,F,2,energy,KZ);
-    K = cat(2,cat(2,repdiag(KX,dim),repdiag(KY,dim)),repdiag(KZ,dim));
+    if(dim == 3)
+    {
+      K = cat(2,cat(2,repdiag(KX,dim),repdiag(KY,dim)),repdiag(KZ,dim));
+    }else if(dim ==2)
+    {
+      SparseMatrix<double> ZZ(KX.rows()*2,KX.cols());
+      K = cat(2,cat(2,
+            cat(2,repdiag(KX,dim),ZZ),
+            cat(2,repdiag(KY,dim),ZZ)),
+            cat(2,repdiag(KZ,dim),ZZ));
+    }else
+    {
+      assert(false);
+      fprintf(
+      stderr,
+      "arap_rhs.h: Error: Unsupported dimension %d\n",
+      dim);
+    }
   }else
   {
+    assert(false);
     fprintf(
      stderr,
-     "covariance_scatter_matrix.h: Error: Unsupported dimension %d\n",
-     dim);
+     "arap_rhs.h: Error: Unsupported dimension %d\n",
+     Vdim);
     return;
   }
   

include/igl/arap_rhs.h

@@ -18,8 +18,10 @@ namespace igl
   // ARAP_RHS build right-hand side constructor of global poisson solve for
   // various Arap energies
   // Inputs:
-  //   V  #V by dim list of initial domain positions
+  //   V  #V by Vdim list of initial domain positions
   //   F  #F by 3 list of triangle indices into V
+  //   dim  dimension being used at solve time. For deformation usually dim =
+  //     V.cols(), for surface parameterization V.cols() = 3 and dim = 2
   //   energy  igl::ARAPEnergyType enum value defining which energy is being
   //     used. See igl::ARAPEnergyType.h for valid options and explanations.
   // Outputs:
@@ -27,9 +29,10 @@ namespace igl
   //     b = K * reshape(permute(R,[3 1 2]),size(V|F,1)*size(V,2)*size(V,2),1);
   //   
   // See also: arap_linear_block
-  void arap_rhs(
+  IGL_INLINE void arap_rhs(
     const Eigen::MatrixXd & V, 
     const Eigen::MatrixXi & F,
+    const int dim,
     const igl::ARAPEnergyType energy,
     Eigen::SparseMatrix<double>& K);
 }

include/igl/areamatrix.cpp

@@ -0,0 +1,56 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2014 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 "areamatrix.h"
+#include <vector>
+
+// Bug in unsupported/Eigen/SparseExtra needs iostream first
+#include <iostream>
+#include <unsupported/Eigen/SparseExtra>
+
+//#include <igl/boundary_vertices_sorted.h>
+#include <igl/boundary_faces.h>
+
+template <typename DerivedF, typename Scalar>
+IGL_INLINE void igl::areamatrix(
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  Eigen::SparseMatrix<Scalar>& A)
+{
+  using namespace Eigen;
+  using namespace std;
+
+  // number of vertices
+  const int n = F.maxCoeff()+1;
+
+	SparseMatrix<Scalar> aux (n * 2, n * 2);
+	SparseMatrix<Scalar> auxT(n * 2, n * 2);
+
+	vector<Triplet<Scalar> > auxTripletList;
+	vector<Triplet<Scalar> > auxTTripletList;
+
+  MatrixXi E;
+  boundary_faces(F,E);
+
+	for(int k = 0; k < E.rows(); k++)
+  {
+		int i = E(k,0);
+		int j = E(k,1);
+		auxTripletList.push_back(Triplet<Scalar>(i+n, j, -0.5));
+		auxTripletList.push_back(Triplet<Scalar>(i, j+n, 0.5));
+		auxTTripletList.push_back(Triplet<Scalar>(j, i+n, -0.5));
+		auxTTripletList.push_back(Triplet<Scalar>(j+n, i, 0.5));
+	}
+
+	aux.setFromTriplets(auxTripletList.begin(), auxTripletList.end());
+	auxT.setFromTriplets(auxTTripletList.begin(), auxTTripletList.end());
+	A = (aux + auxT)/2;
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+template void igl::areamatrix<Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&);
+#endif

include/igl/areamatrix.h

@@ -0,0 +1,41 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2014 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_AREAMATRIX_H
+#define IGL_AREAMATRIX_H
+#include "igl_inline.h"
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+namespace igl
+{
+  // Constructs the symmetric area matrix A, s.t.  [V.col(0)' V.col(1)'] * A *
+  // [V.col(0); V.col(1)] is the **vector area** of the mesh (V,F).
+  //
+  // Templates:
+  //   DerivedV  derived type of eigen matrix for V (e.g. derived from
+  //     MatrixXd)
+  //   DerivedF  derived type of eigen matrix for F (e.g. derived from
+  //     MatrixXi)
+  //   Scalar  scalar type for eigen sparse matrix (e.g. double)
+  // Inputs:
+  //   F  #F by 3 list of mesh faces (must be triangles)
+  // Outputs:
+  //   A  #Vx2 by #Vx2 area matrix
+  //
+  template <typename DerivedF, typename Scalar>
+  IGL_INLINE void areamatrix(
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    Eigen::SparseMatrix<Scalar>& A);
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "areamatrix.cpp"
+#endif
+
+#endif

include/igl/barycentric2global.h

@@ -14,24 +14,25 @@
 
 namespace igl 
 {
-	// Converts barycentric coordinates in the embree form to 3D coordinates
-	// Embree stores barycentric coordinates as triples: fid, bc1, bc2
-	// fid is the id of a face, bc1 is the displacement of the point wrt the 
-	// first vertex v0 and the edge v1-v0. Similarly, bc2 is the displacement
-	// wrt v2-v0.
-  	// 
-  	// Input:
-  	// V:  #Vx3 Vertices of the mesh
-  	// F:  #Fxe Faces of the mesh
-  	// bc: #Xx3 Barycentric coordinates, one row per point
-  	//
-  	// Output:
-  	// #X: #Xx3 3D coordinates of all points in bc
+  // Converts barycentric coordinates in the embree form to 3D coordinates
+  // Embree stores barycentric coordinates as triples: fid, bc1, bc2
+  // fid is the id of a face, bc1 is the displacement of the point wrt the 
+  // first vertex v0 and the edge v1-v0. Similarly, bc2 is the displacement
+  // wrt v2-v0.
+  // 
+  // Input:
+  // V:  #Vx3 Vertices of the mesh
+  // F:  #Fxe Faces of the mesh
+  // bc: #Xx3 Barycentric coordinates, one row per point
+  //
+  // Output:
+  // #X: #Xx3 3D coordinates of all points in bc
   template <typename Scalar, typename Index>
-  IGL_INLINE Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic> barycentric2global(
-    const Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic> & V, 
-	  const Eigen::Matrix<Index,Eigen::Dynamic,Eigen::Dynamic>   & F, 
-    const Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic>  & bc);
+  IGL_INLINE Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic> 
+    barycentric2global(
+      const Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic> & V, 
+      const Eigen::Matrix<Index,Eigen::Dynamic,Eigen::Dynamic>   & F, 
+      const Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic>  & bc);
 }
 
 #ifdef IGL_HEADER_ONLY

include/igl/barycentric_coordinates.cpp

@@ -0,0 +1,87 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "barycentric_coordinates.h"
+#include "volume.h"
+#include "doublearea.h"
+
+template <
+  typename DerivedP,
+  typename DerivedA,
+  typename DerivedB,
+  typename DerivedC,
+  typename DerivedD,
+  typename DerivedL>
+IGL_INLINE void igl::barycentric_coordinates(
+  const Eigen::PlainObjectBase<DerivedP> & P,
+  const Eigen::PlainObjectBase<DerivedA> & A,
+  const Eigen::PlainObjectBase<DerivedB> & B,
+  const Eigen::PlainObjectBase<DerivedC> & C,
+  const Eigen::PlainObjectBase<DerivedD> & D,
+  Eigen::PlainObjectBase<DerivedL> & L)
+{
+  using namespace Eigen;
+  assert(P.cols() == 3 && "query must be in 3d");
+  assert(A.cols() == 3 && "corners must be in 3d");
+  assert(B.cols() == 3 && "corners must be in 3d");
+  assert(C.cols() == 3 && "corners must be in 3d");
+  assert(D.cols() == 3 && "corners must be in 3d");
+  assert(P.rows() == A.rows() && "Must have same number of queries as corners");
+  assert(A.rows() == B.rows() && "Corners must be same size");
+  assert(A.rows() == C.rows() && "Corners must be same size");
+  assert(A.rows() == D.rows() && "Corners must be same size");
+  typedef Matrix<typename DerivedL::Scalar,DerivedL::RowsAtCompileTime,1> 
+    VectorXS;
+  // Total volume
+  VectorXS vol,LA,LB,LC,LD;
+  volume(B,D,C,P,LA);
+  volume(A,C,D,P,LB);
+  volume(A,D,B,P,LC);
+  volume(A,B,C,P,LD);
+  volume(A,B,C,D,vol);
+  L.resize(P.rows(),4);
+  L<<LA,LB,LC,LD;
+  L.array().colwise() /= vol.array();
+}
+
+template <
+  typename DerivedP,
+  typename DerivedA,
+  typename DerivedB,
+  typename DerivedC,
+  typename DerivedL>
+IGL_INLINE void igl::barycentric_coordinates(
+  const Eigen::PlainObjectBase<DerivedP> & P,
+  const Eigen::PlainObjectBase<DerivedA> & A,
+  const Eigen::PlainObjectBase<DerivedB> & B,
+  const Eigen::PlainObjectBase<DerivedC> & C,
+  Eigen::PlainObjectBase<DerivedL> & L)
+{
+  using namespace Eigen;
+  assert(P.cols() == 2 && "query must be in 2d");
+  assert(A.cols() == 2 && "corners must be in 2d");
+  assert(B.cols() == 2 && "corners must be in 2d");
+  assert(C.cols() == 2 && "corners must be in 2d");
+  assert(P.rows() == A.rows() && "Must have same number of queries as corners");
+  assert(A.rows() == B.rows() && "Corners must be same size");
+  assert(A.rows() == C.rows() && "Corners must be same size");
+  typedef Matrix<typename DerivedL::Scalar,DerivedL::RowsAtCompileTime,1> 
+    VectorXS;
+  // Total area
+  VectorXS dblA,LA,LB,LC;
+  doublearea(P,B,C,LA);
+  doublearea(A,P,C,LB);
+  doublearea(A,B,P,LC);
+  doublearea(A,B,C,dblA);
+  L.resize(P.rows(),3);
+  L<<LA,LB,LC;
+  L.array().colwise() /= dblA.array();
+}
+
+#ifndef IGL_HEADER_ONLY
+template void igl::barycentric_coordinates<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+#endif

include/igl/barycentric_coordinates.h

@@ -0,0 +1,68 @@
+// 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_BARYCENTRIC_COORDINATES_H
+#define IGL_BARYCENTRIC_COORDINATES_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Compute barycentric coordinates in a tet
+  //
+  // Inputs:
+  //   P  #P by 3 Query points in 3d
+  //   A  #P by 3 Tet corners in 3d
+  //   B  #P by 3 Tet corners in 3d
+  //   C  #P by 3 Tet corners in 3d
+  //   D  #P by 3 Tet corners in 3d
+  // Outputs:
+  //   L  #P by 4 list of barycentric coordinates
+  //   
+  template <
+    typename DerivedP,
+    typename DerivedA,
+    typename DerivedB,
+    typename DerivedC,
+    typename DerivedD,
+    typename DerivedL>
+  IGL_INLINE void barycentric_coordinates(
+    const Eigen::PlainObjectBase<DerivedP> & P,
+    const Eigen::PlainObjectBase<DerivedA> & A,
+    const Eigen::PlainObjectBase<DerivedB> & B,
+    const Eigen::PlainObjectBase<DerivedC> & C,
+    const Eigen::PlainObjectBase<DerivedD> & D,
+    Eigen::PlainObjectBase<DerivedL> & L);
+  // Compute barycentric coordinates in a triangle
+  //
+  // Inputs:
+  //   P  #P by 2 Query points in 2d
+  //   A  #P by 2 Triangle corners in 2d
+  //   B  #P by 2 Triangle corners in 2d
+  //   C  #P by 2 Triangle corners in 2d
+  // Outputs:
+  //   L  #P by e list of barycentric coordinates
+  //   
+  template <
+    typename DerivedP,
+    typename DerivedA,
+    typename DerivedB,
+    typename DerivedC,
+    typename DerivedL>
+  IGL_INLINE void barycentric_coordinates(
+    const Eigen::PlainObjectBase<DerivedP> & P,
+    const Eigen::PlainObjectBase<DerivedA> & A,
+    const Eigen::PlainObjectBase<DerivedB> & B,
+    const Eigen::PlainObjectBase<DerivedC> & C,
+    Eigen::PlainObjectBase<DerivedL> & L);
+
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "barycentric_coordinates.cpp"
+#endif
+
+#endif

include/igl/boost/bfs_orient.cpp

@@ -11,7 +11,7 @@
 #include <queue>
 
 template <typename DerivedF, typename DerivedFF, typename DerivedC>
-void igl::bfs_orient(
+IGL_INLINE void igl::bfs_orient(
   const Eigen::PlainObjectBase<DerivedF> & F,
   Eigen::PlainObjectBase<DerivedFF> & FF,
   Eigen::PlainObjectBase<DerivedC> & C)

include/igl/boost/bfs_orient.h

@@ -24,7 +24,7 @@ namespace igl
   //
   //
   template <typename DerivedF, typename DerivedFF, typename DerivedC>
-  void bfs_orient(
+  IGL_INLINE void bfs_orient(
     const Eigen::PlainObjectBase<DerivedF> & F,
     Eigen::PlainObjectBase<DerivedFF> & FF,
     Eigen::PlainObjectBase<DerivedC> & C);

include/igl/boundary_vertices_sorted.cpp

@@ -0,0 +1,105 @@
+// 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 "boundary_vertices_sorted.h"
+
+#include "tt.h"
+#include "vf.h"
+
+IGL_INLINE void igl::boundary_vertices_sorted(
+    const Eigen::MatrixXd& V,
+    const Eigen::MatrixXi& F,
+    Eigen::VectorXi& b)
+{
+  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_manifold(V,F) && "Mesh must be manifold");
+  Eigen::MatrixXi TT,TTi;
+  std::vector<std::vector<int> > VF, VFi;
+  igl::tt(V,F,TT,TTi);
+  igl::vf(V,F,VF,VFi);
+
+  // Extract one boundary edge
+  bool done = false;
+  for (int i = 0; i < TT.rows() && !done; 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;
+      }
+    }
+  }
+
+  // Traverse boundary
+  while(1)
+  {
+    bool changed = false;
+    int lastV;
+    lastV = bnd[bnd.size()-1];
+
+    for (int i = 0; i < VF[lastV].size(); i++)
+    {
+      int curr_neighbor = VF[lastV][i];
+
+      if (TT.row(curr_neighbor).minCoeff() < 0.) // Face contains boundary edge
+      {
+        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;
+        }
+
+        if (gotBndEdge)
+        {
+          changed = true;
+          bnd.push_back(F(curr_neighbor,next_bnd_vertex));
+          isVisited[F(curr_neighbor,next_bnd_vertex)] = true;
+          break;
+        }
+      }
+    }
+
+    if (!changed)
+      break;
+  }
+
+  b.resize(bnd.size());
+  for(unsigned i=0;i<bnd.size();++i)
+    b(i) = bnd[i];
+}

include/igl/boundary_vertices_sorted.h

@@ -0,0 +1,36 @@
+// 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_BOUNDARY_VERTICES_SORTED_H
+#define IGL_BOUNDARY_VERTICES_SORTED_H
+#include <igl/igl_inline.h>
+
+#include <Eigen/Dense>
+#include <vector>
+
+namespace igl
+{
+
+  // Compute sorted list of boundary vertices for a manifold mesh with single
+  // boundary.
+  //
+  // 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
+  IGL_INLINE void boundary_vertices_sorted(
+  	const Eigen::MatrixXd& V, 
+  	const Eigen::MatrixXi& F, 
+    Eigen::VectorXi& bnd);
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "boundary_vertices_sorted.cpp"
+#endif
+
+#endif

include/igl/cgal/CGAL_includes.hpp

@@ -1,3 +1,10 @@
+// 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_CGAL_INCLUDES_H
 #define IGL_CGAL_INCLUDES_H
 

include/igl/cgal/SelfIntersectMesh.h

@@ -1,3 +1,10 @@
+// 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_SELFINTERSECTMESH_H
 #define IGL_SELFINTERSECTMESH_H
 
@@ -78,13 +85,18 @@ namespace igl
     public:
       // Constructs (VV,FF) a new mesh with self-intersections of (V,F)
       // subdivided
+      //
+      // See also: selfintersect.h
       inline SelfIntersectMesh(
         const Eigen::MatrixXd & V,
         const Eigen::MatrixXi & F,
         const SelfintersectParam & params,
         Eigen::MatrixXd & VV,
         Eigen::MatrixXi & FF,
-        Eigen::MatrixXi & IF);
+        Eigen::MatrixXi & IF,
+        Eigen::VectorXi & J,
+        Eigen::VectorXi & IM
+        );
     private:
       // Helper function to mark a face as offensive
       //
@@ -243,7 +255,9 @@ inline igl::SelfIntersectMesh<Kernel>::SelfIntersectMesh(
   const SelfintersectParam & params,
   Eigen::MatrixXd & VV,
   Eigen::MatrixXi & FF,
-  Eigen::MatrixXi & IF):
+  Eigen::MatrixXi & IF,
+  Eigen::VectorXi & J,
+  Eigen::VectorXi & IM):
   V(V),
   F(F),
   count(0),
@@ -348,8 +362,8 @@ inline igl::SelfIntersectMesh<Kernel>::SelfIntersectMesh(
           //  P[o].to_3d(vit->point())<<endl;
 #ifndef NDEBUG
           // I want to be sure that the original corners really show up as the
-          // original corners of the CDT. I.e. I don't trust CGAL to maintain the
-          // order
+          // original corners of the CDT. I.e. I don't trust CGAL to maintain
+          // the order
           assert(T[f].vertex(i) == P[o].to_3d(vit->point()));
 #endif
           // For first three, use original index in F
@@ -442,6 +456,7 @@ inline igl::SelfIntersectMesh<Kernel>::SelfIntersectMesh(
 #endif
   // Append faces
   FF.resize(F.rows()-offending.size()+NF_count,3);
+  J.resize(FF.rows());
   // First append non-offending original faces
   // There's an Eigen way to do this in one line but I forget
   int off = 0;
@@ -449,7 +464,9 @@ inline igl::SelfIntersectMesh<Kernel>::SelfIntersectMesh(
   {
     if(!offensive[f])
     {
-      FF.row(off++) = F.row(f);
+      FF.row(off) = F.row(f);
+      J(off) = f;
+      off++;
     }
   }
   assert(off == (int)(F.rows()-offending.size()));
@@ -457,6 +474,7 @@ inline igl::SelfIntersectMesh<Kernel>::SelfIntersectMesh(
   for(int o = 0;o<(int)offending.size();o++)
   {
     FF.block(off,0,NF[o].rows(),3) = NF[o];
+    J.block(off,0,NF[o].rows(),1).setConstant(offending[o]);
     off += NF[o].rows();
   }
   // Append vertices
@@ -483,6 +501,38 @@ inline igl::SelfIntersectMesh<Kernel>::SelfIntersectMesh(
       i++;
     }
   }
+  IM.resize(VV.rows(),1);
+  map<Point_3,int> vv2i;
+  // Safe to check for duplicates using double for original vertices: if
+  // incoming reps are different then the points are unique.
+  for(int v = 0;v<V.rows();v++)
+  {
+    const Point_3 p(V(v,0),V(v,1),V(v,2));
+    if(vv2i.count(p)==0)
+    {
+      vv2i[p] = v;
+    }
+    assert(vv2i.count(p) == 1);
+    IM(v) = vv2i[p];
+  }
+  // Must check for duplicates of new vertices using exact.
+  {
+    int v = V.rows();
+    for(
+      typename list<Point_3>::const_iterator nvit = NV.begin();
+      nvit != NV.end();
+      nvit++)
+    {
+      const Point_3 & p = *nvit;
+      if(vv2i.count(p)==0)
+      {
+        vv2i[p] = v;
+      }
+      assert(vv2i.count(p) == 1);
+      IM(v) = vv2i[p];
+      v++;
+    }
+  }
 
   // Q: Does this give the same result as TETGEN?
   // A: For the cow and beast, yes.

include/igl/cgal/intersect_other.cpp

@@ -1,3 +1,10 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
 #include "intersect_other.h"
 #include "CGAL_includes.hpp"
 #include "mesh_to_cgal_triangle_list.h"
@@ -121,3 +128,7 @@ IGL_INLINE void igl::intersect_other(
   }
 
 }
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

include/igl/cgal/intersect_other.h

@@ -1,3 +1,10 @@
+// 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_INTERSECT_OTHER_H
 #define IGL_INTERSECT_OTHER_H
 #include <igl/igl_inline.h>

include/igl/cgal/mesh_to_cgal_triangle_list.cpp

@@ -1,3 +1,10 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
 #include "mesh_to_cgal_triangle_list.h"
 
 #include <cassert>

include/igl/cgal/mesh_to_cgal_triangle_list.h

@@ -1,3 +1,10 @@
+// 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_MESH_TO_CGAL_TRIANGLE_LIST_H
 #define IGL_MESH_TO_CGAL_TRIANGLE_LIST_H
 #include <igl/igl_inline.h>
@@ -7,6 +14,9 @@ namespace igl
 {
   // Convert a mesh (V,F) to a list of CGAL triangles
   //
+  // Templates:
+  //   Kernal  CGAL computation and construction kernel (e.g.
+  //     CGAL::Exact_predicates_exact_constructions_kernel)
   // Inputs:
   //   V  #V by 3 list of vertex positions
   //   F  #F by 3 list of triangle indices

include/igl/cgal/point_mesh_squared_distance.cpp

@@ -0,0 +1,63 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "point_mesh_squared_distance.h"
+#include "mesh_to_cgal_triangle_list.h"
+
+template <typename Kernel>
+IGL_INLINE void igl::point_mesh_squared_distance(
+  const Eigen::MatrixXd & P,
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  Eigen::VectorXd & sqrD,
+  Eigen::VectorXi & I,
+  Eigen::MatrixXd & C)
+{
+  using namespace std;
+  typedef CGAL::Point_3<Kernel>    Point_3;
+  typedef CGAL::Triangle_3<Kernel> Triangle_3; 
+  typedef typename std::vector<Triangle_3>::iterator Iterator;
+  typedef CGAL::AABB_triangle_primitive<Kernel, Iterator> Primitive;
+  typedef CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
+  typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
+  typedef typename Tree::Point_and_primitive_id Point_and_primitive_id;
+
+  // Must be 3D
+  assert(V.cols() == 3);
+  assert(P.cols() == 3);
+  // Must be triangles
+  assert(F.cols() == 3);
+  // Make list of cgal triangles
+  Tree tree;
+  vector<Triangle_3> T;
+  mesh_to_cgal_triangle_list(V,F,T);
+  tree.insert(T.begin(),T.end());
+
+  tree.accelerate_distance_queries();
+  const int n = P.rows();
+  sqrD.resize(n,1);
+  I.resize(n,1);
+  C.resize(n,P.cols());
+  for(int p = 0;p < n;p++)
+  {
+    Point_3 query(P(p,0),P(p,1),P(p,2));
+    // Find closest point and primitive id
+    Point_and_primitive_id pp = tree.closest_point_and_primitive(query);
+    Point_3 closest_point = pp.first;
+    C(p,0) = CGAL::to_double(closest_point[0]);
+    C(p,1) = CGAL::to_double(closest_point[1]);
+    C(p,2) = CGAL::to_double(closest_point[2]);
+    sqrD(p) = CGAL::to_double((closest_point-query).squared_length());
+    I(p) = pp.second - T.begin();
+  }
+}
+
+#ifndef IGL_HEADER_ONLY
+// 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);
+
+#endif

include/igl/cgal/point_mesh_squared_distance.h

@@ -0,0 +1,44 @@
+// 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_POINT_MESH_SQUARED_DISTANCE_H
+#define IGL_POINT_MESH_SQUARED_DISTANCE_H
+#include <igl/igl_inline.h>
+#include <Eigen/Core>
+#include "CGAL_includes.hpp"
+namespace igl
+{
+  // Compute distances from a set of points P to a triangle mesh (V,F)
+  //
+  // Templates:
+  //   Kernal  CGAL computation and construction kernel (e.g.
+  //     CGAL::Simple_cartesian<double>)
+  // Inputs:
+  //   P  #P by 3 list of query point positions
+  //   V  #V by 3 list of vertex positions
+  //   F  #F by 3 list of triangle indices
+  // Outputs:
+  //   sqrD  #P list of smallest squared distances
+  //   I  #P list of facet indices corresponding to smallest distances
+  //   C  #P by 3 list of closest points
+  //
+  // Known bugs: This only computes distances to triangles. So unreferenced
+  // vertices are ignored.
+  template <typename Kernel>
+  IGL_INLINE void point_mesh_squared_distance(
+    const Eigen::MatrixXd & P,
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    Eigen::VectorXd & sqrD,
+    Eigen::VectorXi & I,
+    Eigen::MatrixXd & C);
+}
+#ifdef IGL_HEADER_ONLY
+#  include "point_mesh_squared_distance.cpp"
+#endif
+
+#endif

include/igl/cgal/selfintersect.cpp

@@ -1,3 +1,10 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
 #include "selfintersect.h"
 #include "SelfIntersectMesh.h"
 #include <igl/C_STR.h>
@@ -9,7 +16,9 @@ IGL_INLINE void igl::selfintersect(
   const SelfintersectParam & params,
   Eigen::MatrixXd & VV,
   Eigen::MatrixXi & FF,
-  Eigen::MatrixXi & IF)
+  Eigen::MatrixXi & IF,
+  Eigen::VectorXi & J,
+  Eigen::VectorXi & IM)
 {
   using namespace std;
   if(params.detect_only)
@@ -27,7 +36,7 @@ IGL_INLINE void igl::selfintersect(
 //#endif
 
     typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
-    SelfIntersectMesh<Kernel> SIM = SelfIntersectMesh<Kernel>(V,F,params,VV,FF,IF);
+    SelfIntersectMesh<Kernel> SIM = SelfIntersectMesh<Kernel>(V,F,params,VV,FF,IF,J,IM);
 
 //#ifdef __APPLE__
 //    signal(SIGFPE,SIG_DFL);
@@ -36,6 +45,6 @@ IGL_INLINE void igl::selfintersect(
   }else
   {
     typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
-    SelfIntersectMesh<Kernel> SIM = SelfIntersectMesh<Kernel>(V,F,params,VV,FF,IF);
+    SelfIntersectMesh<Kernel> SIM = SelfIntersectMesh<Kernel>(V,F,params,VV,FF,IF,J,IM);
   }
 }

include/igl/cgal/selfintersect.h

@@ -1,3 +1,10 @@
+// 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_SELFINTERSECT_H
 #define IGL_SELFINTERSECT_H
 #include <igl/igl_inline.h>
@@ -39,13 +46,22 @@ namespace igl
   //   FF  #FF by 3 list of triangle indices into V
   //   IF  #intersecting face pairs by 2  list of intersecting face pairs,
   //     indexing F
+  //   J  #FF list of indices into F denoting birth triangle
+  //   IM  #VV list of indices into VV of unique vertices.
+  //
+  // Known bugs: If an existing edge in (V,F) lies exactly on another face then
+  // any resulting additional vertices along that edge may not get properly
+  // connected so that the output mesh has the same global topology. This is
+  // because 
   IGL_INLINE void selfintersect(
     const Eigen::MatrixXd & V,
     const Eigen::MatrixXi & F,
     const SelfintersectParam & params,
     Eigen::MatrixXd & VV,
     Eigen::MatrixXi & FF,
-    Eigen::MatrixXi & IF);
+    Eigen::MatrixXi & IF,
+    Eigen::VectorXi & J,
+    Eigen::VectorXi & IM);
 }
 
 #ifdef IGL_HEADER_ONLY

include/igl/cotmatrix.cpp

@@ -6,6 +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 "cotmatrix.h"
+#include <vector>
 
 // For error printing
 #include <cstdio>
@@ -23,9 +24,9 @@ IGL_INLINE void igl::cotmatrix(
 {
   using namespace igl;
   using namespace Eigen;
-  Eigen::DynamicSparseMatrix<double> foo;
+  using namespace std;
 
-  DynamicSparseMatrix<Scalar, RowMajor> dyn_L (V.rows(), V.rows());
+  L.resize(V.rows(),V.rows());
   Matrix<int,Dynamic,2> edges;
   int simplex_size = F.cols();
   // 3 for triangles, 4 for tets
@@ -35,7 +36,7 @@ IGL_INLINE void igl::cotmatrix(
     // This is important! it could decrease the comptuation time by a factor of 2
     // Laplacian for a closed 2d manifold mesh will have on average 7 entries per
     // row
-    dyn_L.reserve(7*V.rows());
+    L.reserve(10*V.rows());
     edges.resize(3,2);
     edges << 
       1,2,
@@ -43,7 +44,7 @@ IGL_INLINE void igl::cotmatrix(
       0,1;
   }else if(simplex_size == 4)
   {
-    dyn_L.reserve(17*V.rows());
+    L.reserve(17*V.rows());
     edges.resize(6,2);
     edges << 
       1,2,
@@ -60,6 +61,8 @@ IGL_INLINE void igl::cotmatrix(
   Matrix<Scalar,Dynamic,Dynamic> C;
   cotangent(V,F,C);
   
+  vector<Triplet<Scalar> > IJV;
+  IJV.reserve(F.rows()*edges.rows()*4);
   // Loop over triangles
   for(int i = 0; i < F.rows(); i++)
   {
@@ -68,14 +71,13 @@ IGL_INLINE void igl::cotmatrix(
     {
       int source = F(i,edges(e,0));
       int dest = F(i,edges(e,1));
-      dyn_L.coeffRef(source,dest) += C(i,e);
-      dyn_L.coeffRef(dest,source) += C(i,e);
-      dyn_L.coeffRef(source,source) += -C(i,e);
-      dyn_L.coeffRef(dest,dest) += -C(i,e);
+      IJV.push_back(Triplet<Scalar>(source,dest,C(i,e)));
+      IJV.push_back(Triplet<Scalar>(dest,source,C(i,e)));
+      IJV.push_back(Triplet<Scalar>(source,source,-C(i,e)));
+      IJV.push_back(Triplet<Scalar>(dest,dest,-C(i,e)));
     }
   }
-    // Corner indices of this triangle
-  L = SparseMatrix<Scalar>(dyn_L);
+  L.setFromTriplets(IJV.begin(),IJV.end());
 }
 
 #ifndef IGL_HEADER_ONLY

include/igl/cotmatrix.h

@@ -9,7 +9,6 @@
 #define IGL_COTMATRIX_H
 #include "igl_inline.h"
 
-#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
 #include <Eigen/Dense>
 #include <Eigen/Sparse>
 

include/igl/covariance_scatter_matrix.h

@@ -17,12 +17,12 @@ namespace igl
 {
   // Construct the covariance scatter matrix for a given arap energy
   // Inputs:
-  //   V  #V by dim list of initial domain positions
+  //   V  #V by Vdim list of initial domain positions
   //   F  #F by 3 list of triangle indices into V
   //   energy  ARAPEnergyType enum value defining which energy is being used.
   //     See ARAPEnergyType.h for valid options and explanations.
   // Outputs:
-  //   CSM dim*#V by dim*#V sparse matrix containing special laplacians along
+  //   CSM dim*#V/#F by dim*#V sparse matrix containing special laplacians along
   //     the diagonal so that when multiplied by V gives covariance matrix
   //     elements, can be used to speed up covariance matrix computation
   IGL_INLINE void covariance_scatter_matrix(

include/igl/cross.cpp

@@ -17,3 +17,29 @@ IGL_INLINE void igl::cross(
   out[1] = a[2]*b[0]-a[0]*b[2];
   out[2] = a[0]*b[1]-a[1]*b[0];
 }
+
+template <
+  typename DerivedA,
+  typename DerivedB,
+  typename DerivedC>
+IGL_INLINE void igl::cross(
+  const Eigen::PlainObjectBase<DerivedA> & A,
+  const Eigen::PlainObjectBase<DerivedB> & B,
+  Eigen::PlainObjectBase<DerivedC> & C)
+{
+  assert(A.cols() == 3 && "#cols should be 3");
+  assert(B.cols() == 3 && "#cols should be 3");
+  assert(A.rows() == B.rows() && "#rows in A and B should be equal");
+  C.resize(A.rows(),3);
+  for(int d = 0;d<3;d++)
+  {
+    C.col(d) = 
+      A.col((d+1)%3).array() * B.col((d+2)%3).array() -
+      A.col((d+2)%3).array() * B.col((d+1)%3).array();
+  }
+}
+
+#ifndef IGL_HEADER_ONLY
+template void igl::cross<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::cross<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::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> >&);
+#endif

include/igl/cross.h

@@ -8,6 +8,7 @@
 #ifndef IGL_CROSS_H
 #define IGL_CROSS_H
 #include "igl_inline.h"
+#include <Eigen/Core>
 namespace igl
 {
   // Computes out = cross(a,b)
@@ -16,10 +17,22 @@ namespace igl
   //   b  right 3d vector
   // Outputs:
   //   out  result 3d vector
+  IGL_INLINE void cross( const double *a, const double *b, double *out);
+  // Computes C = cross(A,B,2);
+  //
+  // Inputs:
+  //   A  #A by 3 list of row-vectors
+  //   B  #A by 3 list of row-vectors
+  // Outputs:
+  //   C  #A by 3 list of row-vectors
+  template <
+    typename DerivedA,
+    typename DerivedB,
+    typename DerivedC>
   IGL_INLINE void cross(
-    const double *a, 
-    const double *b,
-    double *out);
+    const Eigen::PlainObjectBase<DerivedA> & A,
+    const Eigen::PlainObjectBase<DerivedB> & B,
+    Eigen::PlainObjectBase<DerivedC> & C);
 }
 
 #ifdef IGL_HEADER_ONLY

include/igl/dated_copy.cpp

@@ -1,3 +1,10 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
 #include "dated_copy.h"
 #include "dirname.h"
 #include "basename.h"
@@ -9,7 +16,7 @@
 #include <unistd.h>
 #include <iostream>
 
-bool igl::dated_copy(const std::string & src_path, const std::string & dir)
+IGL_INLINE bool igl::dated_copy(const std::string & src_path, const std::string & dir)
 {
   using namespace std;
   using namespace igl;
@@ -75,7 +82,7 @@ bool igl::dated_copy(const std::string & src_path, const std::string & dir)
   return true;
 }
 
-bool igl::dated_copy(const std::string & src_path)
+IGL_INLINE bool igl::dated_copy(const std::string & src_path)
 {
   return dated_copy(src_path,dirname(src_path));
 }

include/igl/dated_copy.h

@@ -20,9 +20,9 @@ namespace igl
   // Example:
   //   dated_copy("/path/to/foo","/bar/");
   //   // copies /path/to/foo to /bar/foo-2013-12-12T18-10-56
-  bool dated_copy(const std::string & src_path, const std::string & dir);
+  IGL_INLINE bool dated_copy(const std::string & src_path, const std::string & dir);
   // Wrapper using directory of source file
-  bool dated_copy(const std::string & src_path);
+  IGL_INLINE bool dated_copy(const std::string & src_path);
 }
 #ifdef IGL_HEADER_ONLY
 #  include "dated_copy.cpp"

include/igl/doublearea.cpp

@@ -60,6 +60,41 @@ IGL_INLINE void igl::doublearea(
   }
 }
 
+template <
+  typename DerivedA,
+  typename DerivedB,
+  typename DerivedC,
+  typename DerivedD>
+IGL_INLINE void doublearea( 
+  const Eigen::PlainObjectBase<DerivedA> & A, 
+  const Eigen::PlainObjectBase<DerivedB> & B, 
+  const Eigen::PlainObjectBase<DerivedC> & C,
+  Eigen::PlainObjectBase<DerivedD> & D)
+{
+  assert(A.cols() == 2 && "corners should be 2d");
+  assert(B.cols() == 2 && "corners should be 2d");
+  assert(C.cols() == 2 && "corners should be 2d");
+  assert(A.rows() == B.rows() && "corners should have same length");
+  assert(A.rows() == C.rows() && "corners should have same length");
+  const auto & R = A-C;
+  const auto & S = B-C;
+  D = R.col(0).array()*S.col(1).array() - R.col(1).array()*S.col(0).array();
+}
+
+template <
+  typename DerivedA,
+  typename DerivedB,
+  typename DerivedC>
+IGL_INLINE typename DerivedA::Scalar igl::doublearea_single(
+  const Eigen::PlainObjectBase<DerivedA> & A, 
+  const Eigen::PlainObjectBase<DerivedB> & B, 
+  const Eigen::PlainObjectBase<DerivedC> & C)
+{
+  auto r = A-C;
+  auto s = B-C;
+  return r(0)*s(1) - r(1)*s(0);
+}
+
 template <typename Derivedl, typename DeriveddblA>
 IGL_INLINE void igl::doublearea( 
   const Eigen::PlainObjectBase<Derivedl> & ul,
@@ -108,4 +143,5 @@ template void igl::doublearea<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::M
 template void igl::doublearea<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template void igl::doublearea<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::doublearea<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -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<double, -1, 1, 0, -1, 1> >&);
+template Eigen::Matrix<double, 2, 1, 0, 2, 1>::Scalar igl::doublearea_single<Eigen::Matrix<double, 2, 1, 0, 2, 1>, Eigen::Matrix<double, 2, 1, 0, 2, 1>, Eigen::Matrix<double, 2, 1, 0, 2, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 2, 1, 0, 2, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 2, 1, 0, 2, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 2, 1, 0, 2, 1> > const&);
 #endif

include/igl/doublearea.h

@@ -24,14 +24,35 @@ namespace igl
   //   V  #V by dim list of mesh vertex positions
   //   F  #F by simplex_size list of mesh faces (must be triangles)
   // Outputs:
-  //   dblA  #F list of triangle double areas
+  //   dblA  #F list of triangle double areas (SIGNED only for 2D input)
   //
-  // Note: THESE ARE *NOT* SIGNED. In matlab doublearea is signed in 2d
   template <typename DerivedV, typename DerivedF, typename DeriveddblA>
   IGL_INLINE void doublearea( 
     const Eigen::PlainObjectBase<DerivedV> & V, 
     const Eigen::PlainObjectBase<DerivedF> & F, 
     Eigen::PlainObjectBase<DeriveddblA> & dblA);
+  // Stream of triangles
+  template <
+    typename DerivedA,
+    typename DerivedB,
+    typename DerivedC,
+    typename DerivedD>
+  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!
+  //
+  // This should handle streams of corners not just single corners
+  template <
+    typename DerivedA,
+    typename DerivedB,
+    typename DerivedC>
+  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:
   //   DerivedV  derived type of eigen matrix for V (e.g. derived from

include/igl/dqs.cpp

@@ -1,3 +1,10 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
 #include "dqs.h"
 #include <Eigen/Geometry>
 template <

include/igl/dqs.h

@@ -1,3 +1,10 @@
+// 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_DQS_H
 #define IGL_DQS_H
 #include "igl_inline.h"

include/igl/draw_rectangular_marquee.cpp

@@ -9,7 +9,7 @@
 #include "OpenGL_convenience.h"
 #include "material_colors.h"
 
-void igl::draw_rectangular_marquee(
+IGL_INLINE void igl::draw_rectangular_marquee(
   const int from_x,
   const int from_y,
   const int to_x,

include/igl/draw_rectangular_marquee.h

@@ -7,6 +7,7 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #ifndef IGL_DRAW_RECTANGULAR_MARQUEE_H
 #define IGL_DRAW_RECTANGULAR_MARQUEE_H
+#include "igl_inline.h"
 namespace igl
 {
   // Draw a rectangular marquee (selection box) in screen space. This sets up
@@ -17,7 +18,7 @@ namespace igl
   //   from_y  y coordinate of from point
   //   to_x  x coordinate of to point
   //   to_y  y coordinate of to point
-  void draw_rectangular_marquee(
+  IGL_INLINE void draw_rectangular_marquee(
     const int from_x,
     const int from_y,
     const int to_x,

include/igl/embree/ambient_occlusion.cpp

@@ -14,7 +14,7 @@ template <
   typename DerivedP,
   typename DerivedN,
   typename DerivedS >
-void igl::ambient_occlusion(
+IGL_INLINE void igl::ambient_occlusion(
   const igl::EmbreeIntersector & ei,
   const Eigen::PlainObjectBase<DerivedP> & P,
   const Eigen::PlainObjectBase<DerivedN> & N,
@@ -61,7 +61,7 @@ template <
   typename DerivedP,
   typename DerivedN,
   typename DerivedS >
-void igl::ambient_occlusion(
+IGL_INLINE void igl::ambient_occlusion(
   const Eigen::PlainObjectBase<DerivedV> & V,
   const Eigen::PlainObjectBase<DerivedF> & F,
   const Eigen::PlainObjectBase<DerivedP> & P,

include/igl/embree/ambient_occlusion.h

@@ -27,7 +27,7 @@ namespace igl
     typename DerivedP,
     typename DerivedN,
     typename DerivedS >
-  void ambient_occlusion(
+  IGL_INLINE void ambient_occlusion(
     const igl::EmbreeIntersector & ei,
     const Eigen::PlainObjectBase<DerivedP> & P,
     const Eigen::PlainObjectBase<DerivedN> & N,
@@ -41,7 +41,7 @@ namespace igl
     typename DerivedP,
     typename DerivedN,
     typename DerivedS >
-  void ambient_occlusion(
+  IGL_INLINE void ambient_occlusion(
     const Eigen::PlainObjectBase<DerivedV> & V,
     const Eigen::PlainObjectBase<DerivedF> & F,
     const Eigen::PlainObjectBase<DerivedP> & P,

include/igl/embree/project_mesh.cpp

@@ -35,11 +35,11 @@ namespace igl
   template <typename ScalarMatrix, typename IndexMatrix>
   IGL_INLINE ScalarMatrix project_mesh
   (
-		const ScalarMatrix & V_source,
-  	const IndexMatrix  & F_source,
+    const ScalarMatrix & V_source,
+    const IndexMatrix  & F_source,
     const ScalarMatrix & V_target,
     const IndexMatrix  & F_target
-	);
+  );
 
   // Project the point cloud V_source onto the triangle mesh
   // V_target,F_target. 

include/igl/embree/reorient_facets_raycast.cpp

@@ -10,7 +10,6 @@
 #include "../doublearea.h"
 #include "../random_dir.h"
 #include "../boost/bfs_orient.h"
-//#include "orient_outward_ao.h"
 #include "EmbreeIntersector.h"
 #include <iostream>
 #include <random>

include/igl/embree/reorient_facets_raycast.h

@@ -5,8 +5,8 @@
 // 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_ORIENT_OUTWARD_AO_H
-#define IGL_ORIENT_OUTWARD_AO_H
+#ifndef IGL_REORIENT_FACETS_RAYCAST_H
+#define IGL_REORIENT_FACETS_RAYCAST_H
 #include "../igl_inline.h"
 #include <Eigen/Core>
 namespace igl

include/igl/embree/unproject_in_mesh.cpp

@@ -12,7 +12,7 @@
 
 template <
   typename Derivedobj>
-int igl::unproject_in_mesh(
+IGL_INLINE int igl::unproject_in_mesh(
   const int x,
   const int y,
   const igl::EmbreeIntersector & ei,
@@ -24,7 +24,7 @@ int igl::unproject_in_mesh(
 
 template <
   typename Derivedobj>
-int igl::unproject_in_mesh(
+IGL_INLINE int igl::unproject_in_mesh(
   const int x,
   const int y,
   const igl::EmbreeIntersector & ei,

include/igl/embree/unproject_in_mesh.h

@@ -30,7 +30,7 @@ namespace igl
   //
   template <
     typename Derivedobj>
-  int unproject_in_mesh(
+  IGL_INLINE int unproject_in_mesh(
     const int x,
     const int y,
     const igl::EmbreeIntersector & ei,
@@ -38,7 +38,7 @@ namespace igl
 
   template <
     typename Derivedobj>
-  int unproject_in_mesh(
+  IGL_INLINE int unproject_in_mesh(
     const int x,
     const int y,
     const igl::EmbreeIntersector & ei,

include/igl/exterior_edges.cpp

@@ -0,0 +1,104 @@
+#include "exterior_edges.h"
+#include "all_edges.h"
+
+#include <cassert>
+#include <map>
+#include <utility>
+
+//template <typename T> inline int sgn(T val) {
+//      return (T(0) < val) - (val < T(0));
+//}
+
+//static void mod2(std::pair<const std::pair<const int, const int>, int>& p)
+//{
+//  using namespace std;
+//  // Be sure that sign of mod matches sign of argument
+//  p.second = p.second%2 ? sgn(p.second) : 0;
+//}
+
+//// http://stackoverflow.com/a/5517869/148668
+//struct Compare
+//{
+//   int i;
+//   Compare(const int& i) : i(i) {}
+//};
+//bool operator==(const std::pair<std::pair<const int, const int>,int>&p, const Compare& c)
+//{
+//  return c.i == p.second;
+//}
+//bool operator==(const Compare& c, const std::pair<std::pair<const int, const int>, int> &p)
+//{
+//  return c.i == p.second;
+//}
+
+IGL_INLINE void igl::exterior_edges(
+  const Eigen::MatrixXi & F,
+  Eigen::MatrixXi & E)
+{
+  using namespace Eigen;
+  using namespace std;
+  assert(F.cols() == 3);
+  MatrixXi all_E;
+  all_edges(F,all_E);
+  // Count occurances looking only at pairs (i,j) where i<j, so we count and
+  // edge i-->j as +1 and j-->i as -1
+  map<pair<const int,const int>,int> C;
+  // Loop over edges
+  for(int e = 0;e<all_E.rows();e++)
+  {
+    int i = all_E(e,0);
+    int j = all_E(e,1);
+    if(i<j)
+    {
+      // Forward direction --> +1
+      C[pair<const int,const int>(i,j)]++;
+    }else
+    {
+      // Backward direction --> -1
+      C[pair<const int,const int>(j,i)]--;
+    }
+  }
+  // Q: Why mod off factors of 2? +1/-1 already takes care of interior edges?
+  //// Mod off any factors of 2
+  //for_each(C.begin(),C.end(),mod2);
+  //int zeros = (int) count(C.begin(),C.end(),Compare(0));
+  //E.resize(C.size() - zeros,2);
+  E.resize(all_E.rows(),all_E.cols());
+  int e = 0;
+  // Find all edges with -1 or 1 occurances, flipping direction for -1
+  for(
+    map<pair<const int, const int>,int>::const_iterator cit=C.begin();
+    cit!=C.end();
+    cit++)
+  {
+    int i,j;
+    if(cit->second > 0)
+    {
+      i = cit->first.first;
+      j = cit->first.second;
+    } else if(cit->second < 0)
+    {
+      i = cit->first.second;
+      j = cit->first.first;
+    } else if(cit->second == 0)
+    {
+      continue;
+    }
+    for(int k = 0;k<abs(cit->second);k++)
+    {
+      E(e,0) = i;
+      E(e,1) = j;
+      e++;
+    }
+  }
+  E.conservativeResize(e,2);
+  assert(e == E.rows());
+}
+
+IGL_INLINE Eigen::MatrixXi igl::exterior_edges( const Eigen::MatrixXi & F)
+{
+  using namespace Eigen;
+  MatrixXi E;
+  exterior_edges(F,E);
+  return E;
+}

include/igl/exterior_edges.h

@@ -0,0 +1,26 @@
+#ifndef IGL_EXTERIOR_EDGES_H
+#define IGL_EXTERIOR_EDGES_H
+#include "igl_inline.h"
+#include <Eigen/Dense>
+namespace igl
+{
+  // EXTERIOR_EDGES Determines boundary "edges" and also edges with an
+  // odd number of occurances where seeing edge (i,j) counts as +1 and seeing
+  // the opposite edge (j,i) counts as -1
+  //
+  // Inputs:
+  //   F  #F by simplex_size list of "faces"
+  // Outputs:
+  //   E  #E by simplex_size-1  list of exterior edges
+  //
+  IGL_INLINE void exterior_edges(
+    const Eigen::MatrixXi & F,
+    Eigen::MatrixXi & E);
+  // Inline version
+  Eigen::MatrixXi exterior_edges( const Eigen::MatrixXi & F);
+}
+#ifdef IGL_HEADER_ONLY
+#  include "exterior_edges.h"
+#endif
+
+#endif

include/igl/face_areas.cpp

@@ -1,3 +1,10 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
 #include "face_areas.h"
 #include "edge_lengths.h"
 #include "doublearea.h"

include/igl/file_dialog_open.cpp

@@ -0,0 +1,87 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 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 "file_dialog_open.h"
+#include <cstdio>
+#include <cstring>
+
+
+#ifdef _WIN32
+ #include <Commdlg.h>
+#endif
+
+IGL_INLINE std::string igl::file_dialog_open()
+{
+  const int FILE_DIALOG_MAX_BUFFER = 1024;
+  char buffer[FILE_DIALOG_MAX_BUFFER];
+  
+#ifdef __APPLE__
+  // For apple use applescript hack
+  FILE * output = popen(
+    "osascript -e \""
+    "   tell application \\\"System Events\\\"\n"
+    "           activate\n"
+    "           set existing_file to choose file\n"
+    "   end tell\n"
+    "   set existing_file_path to (POSIX path of (existing_file))\n"
+    "\" 2>/dev/null | tr -d '\n' ","r");
+  while ( fgets(buffer, FILE_DIALOG_MAX_BUFFER, output) != NULL )
+  {
+  }
+#elif _WIN32
+  
+  // Use native windows file dialog box
+  // (code contributed by Tino Weinkauf)
+
+  OPENFILENAME ofn;       // common dialog box structure
+  char szFile[260];       // buffer for file name
+  HWND hwnd;              // owner window
+  HANDLE hf;              // file handle
+
+  // Initialize OPENFILENAME
+  ZeroMemory(&ofn, sizeof(ofn));
+  ofn.lStructSize = sizeof(ofn);
+  ofn.hwndOwner = NULL;//hwnd;
+  ofn.lpstrFile = new char[100];
+  // Set lpstrFile[0] to '\0' so that GetOpenFileName does not 
+  // use the contents of szFile to initialize itself.
+  ofn.lpstrFile[0] = '\0';
+  ofn.nMaxFile = sizeof(szFile);
+  ofn.lpstrFilter = "*.*\0";//off\0*.off\0obj\0*.obj\0mp\0*.mp\0";
+  ofn.nFilterIndex = 1;
+  ofn.lpstrFileTitle = NULL;
+  ofn.nMaxFileTitle = 0;
+  ofn.lpstrInitialDir = NULL;
+  ofn.Flags = OFN_PATHMUSTEXIST | OFN_FILEMUSTEXIST;
+
+  // Display the Open dialog box. 
+  int pos = 0;
+  if (GetOpenFileName(&ofn)==TRUE)
+  {
+    while(ofn.lpstrFile[pos] != '\0')
+    {
+      buffer[pos] = (char)ofn.lpstrFile[pos];
+      pos++;
+    }
+    buffer[pos] = 0;
+  } 
+  
+#else
+  
+  // For linux use zenity
+  FILE * output = popen("/usr/bin/zenity --file-selection","r");
+  while ( fgets(buffer, FILE_DIALOG_MAX_BUFFER, output) != NULL )
+  {
+  }
+  
+  if (strlen(buffer) > 0)
+  {
+    buffer[strlen(buffer)-1] = 0;
+  }
+#endif
+  return std::string(buffer);
+}

include/igl/file_dialog_open.h

@@ -0,0 +1,30 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 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_FILE_DIALOG_OPEN_H
+#define IGL_FILE_DIALOG_OPEN_H
+#include "igl_inline.h"
+
+#include <string>
+
+namespace igl
+{
+  // Returns a string with a path to an existing file
+  // The string is returned empty if no file is selected
+  // (on Linux machines, it assumes that Zenity is installed)
+  //
+  // Usage:
+  //   std::string str = get_open_file_path();
+  IGL_INLINE std::string file_dialog_open();
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "file_dialog_open.cpp"
+#endif
+
+#endif
+

include/igl/file_dialog_save.cpp

@@ -0,0 +1,88 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 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 "file_dialog_save.h"
+#include <cstdio>
+#include <cstring>
+
+
+#ifdef _WIN32
+ #include <Commdlg.h>
+#endif
+
+IGL_INLINE std::string igl::file_dialog_save()
+{
+  const int FILE_DIALOG_MAX_BUFFER = 1024;
+  char buffer[FILE_DIALOG_MAX_BUFFER];
+#ifdef __APPLE__
+  // For apple use applescript hack
+  // There is currently a bug in Applescript that strips extensions off
+  // of chosen existing files in the "choose file name" dialog
+  // I'm assuming that will be fixed soon
+  FILE * output = popen(
+    "osascript -e \""
+    "   tell application \\\"System Events\\\"\n"
+    "           activate\n"
+    "           set existing_file to choose file name\n"
+    "   end tell\n"
+    "   set existing_file_path to (POSIX path of (existing_file))\n"
+    "\" 2>/dev/null | tr -d '\n' ","r");
+  while ( fgets(buffer, FILE_DIALOG_MAX_BUFFER, output) != NULL )
+  {
+  }
+#elif _WIN32
+
+  // Use native windows file dialog box
+  // (code contributed by Tino Weinkauf)
+
+  OPENFILENAME ofn;       // common dialog box structure
+  char szFile[260];       // buffer for file name
+  HWND hwnd;              // owner window
+  HANDLE hf;              // file handle
+
+  // Initialize OPENFILENAME
+  ZeroMemory(&ofn, sizeof(ofn));
+  ofn.lStructSize = sizeof(ofn);
+  ofn.hwndOwner = NULL;//hwnd;
+  ofn.lpstrFile = new char[100];
+  // Set lpstrFile[0] to '\0' so that GetOpenFileName does not 
+  // use the contents of szFile to initialize itself.
+  ofn.lpstrFile[0] = '\0';
+  ofn.nMaxFile = sizeof(szFile);
+  ofn.lpstrFilter = "";
+  ofn.nFilterIndex = 1;
+  ofn.lpstrFileTitle = NULL;
+  ofn.nMaxFileTitle = 0;
+  ofn.lpstrInitialDir = NULL;
+  ofn.Flags = OFN_PATHMUSTEXIST | OFN_FILEMUSTEXIST;
+
+  // Display the Open dialog box. 
+  int pos = 0;
+  if (GetSaveFileName(&ofn)==TRUE)
+  {
+    while(ofn.lpstrFile[pos] != '\0')
+    {
+      buffer[pos] = (char)ofn.lpstrFile[pos];
+      pos++;
+    }
+    buffer[pos] = 0;
+  }
+
+#else
+  // For every other machine type use zenity
+  FILE * output = popen("/usr/bin/zenity --file-selection --save","r");
+  while ( fgets(buffer, FILE_DIALOG_MAX_BUFFER, output) != NULL )
+  {
+  }
+  
+  if (strlen(buffer) > 0)
+  {
+    buffer[strlen(buffer)-1] = 0;
+  }
+#endif
+  return std::string(buffer);
+}

include/igl/file_dialog_save.h

@@ -0,0 +1,31 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 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_FILE_DIALOG_SAVE_H
+#define IGL_FILE_DIALOG_SAVE_H
+#include "igl_inline.h"
+
+#include <string>
+
+namespace igl
+{
+  // Returns a string with a path to a new/existing file
+  // The string is returned empty if no file is selected
+  // (on Linux machines, it assumes that Zenity is installed)
+  //
+  // Usage:
+  //   char buffer[FILE_DIALOG_MAX_BUFFER];
+  //   get_save_file_path(buffer);
+  IGL_INLINE std::string file_dialog_save();
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "file_dialog_save.cpp"
+#endif
+
+#endif
+

include/igl/fit_plane.cpp

@@ -14,36 +14,37 @@ void igl::fit_plane(
     Eigen::RowVector3d & N,
     Eigen::RowVector3d & C)
 {
-    assert(V.rows()>0);
+  assert(V.rows()>0);
 
-    Eigen::Vector3d sum = V.colwise().sum();
+  Eigen::Vector3d sum = V.colwise().sum();
 
-    Eigen::Vector3d center = sum.array()/(double(V.rows()));
+  Eigen::Vector3d center = sum.array()/(double(V.rows()));
 
-    C = center;
+  C = center;
 
-    double sumXX=0.0f,sumXY=0.0f,sumXZ=0.0f;
-    double sumYY=0.0f,sumYZ=0.0f;
-    double sumZZ=0.0f;
+  double sumXX=0.0f,sumXY=0.0f,sumXZ=0.0f;
+  double sumYY=0.0f,sumYZ=0.0f;
+  double sumZZ=0.0f;
 
-    for(int i=0;i<V.rows();i++){
-        double diffX=V(i,0)-center(0);
-        double diffY=V(i,1)-center(1);
-        double diffZ=V(i,2)-center(2);
-        sumXX+=diffX*diffX;
-        sumXY+=diffX*diffY;
-        sumXZ+=diffX*diffZ;
-        sumYY+=diffY*diffY;
-        sumYZ+=diffY*diffZ;
-        sumZZ+=diffZ*diffZ;
-    }
+  for(int i=0;i<V.rows();i++)
+  {
+    double diffX=V(i,0)-center(0);
+    double diffY=V(i,1)-center(1);
+    double diffZ=V(i,2)-center(2);
+    sumXX+=diffX*diffX;
+    sumXY+=diffX*diffY;
+    sumXZ+=diffX*diffZ;
+    sumYY+=diffY*diffY;
+    sumYZ+=diffY*diffZ;
+    sumZZ+=diffZ*diffZ;
+  }
 
-    Eigen::MatrixXd m(3,3);
-    m << sumXX,sumXY,sumXZ,
-        sumXY,sumYY,sumYZ,
-        sumXZ,sumYZ,sumZZ;
+  Eigen::MatrixXd m(3,3);
+  m << sumXX,sumXY,sumXZ,
+    sumXY,sumYY,sumYZ,
+    sumXZ,sumYZ,sumZZ;
 
-	Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es(m);
+  Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es(m);
   
   N = es.eigenvectors().col(0);
 }

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

@@ -1 +1 @@
-da428b7e4052aa0f654f311c46b07a27bc50d653
+d7410116de836a2f96c3de72373ca1d089065867

include/igl/gaussian_curvature.cpp

@@ -0,0 +1,55 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "gaussian_curvature.h"
+#include "internal_angles.h"
+#include "PI.h"
+#include <iostream>
+template <typename DerivedV, typename DerivedF, typename DerivedK>
+IGL_INLINE void igl::gaussian_curvature(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedK> & K)
+{
+  using namespace Eigen;
+  using namespace std;
+  // internal corner angles
+  Matrix<
+    typename DerivedV::Scalar,
+    DerivedF::RowsAtCompileTime,
+    DerivedF::ColsAtCompileTime> A;
+  internal_angles(V,F,A);
+  K.resize(V.rows(),1);
+  K.setConstant(V.rows(),1,2.*PI);
+  assert(A.rows() == F.rows());
+  assert(A.cols() == F.cols());
+  assert(K.rows() == V.rows());
+  assert(F.maxCoeff() < V.rows());
+  assert(K.cols() == 1);
+  const int Frows = F.rows();
+  //K_G(x_i) = (2π - ∑θj)
+//#ifndef IGL_GAUSSIAN_CURVATURE_OMP_MIN_VALUE
+//#  define IGL_GAUSSIAN_CURVATURE_OMP_MIN_VALUE 1000
+//#endif
+//#pragma omp parallel for if (Frows>IGL_GAUSSIAN_CURVATURE_OMP_MIN_VALUE)
+  for(int f = 0;f<Frows;f++)
+  {
+    // throw normal at each corner
+    for(int j = 0; j < 3;j++)
+    {
+      // Q: Does this need to be critical?
+      // H: I think so, sadly. Maybe there's a way to use reduction
+//#pragma omp critical
+      K(F(f,j),0) -=  A(f,j);
+    }
+  }
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+template void igl::gaussian_curvature<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> >&);
+#endif

include/igl/gaussian_curvature.h

@@ -0,0 +1,32 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_GAUSSIAN_CURVATURE_H
+#define IGL_GAUSSIAN_CURVATURE_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Compute discrete gaussian curvature (angle deficit)
+  // Inputs:
+  //   V  #V by 3 eigen Matrix of mesh vertex 3D positions
+  //   F  #F by 3 eigen Matrix of face (triangle) indices
+  // Output:
+  //   K  #V by 1 eigen Matrix of discrete gaussian curvature values
+  template <typename DerivedV, typename DerivedF, typename DerivedK>
+  IGL_INLINE void gaussian_curvature(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    Eigen::PlainObjectBase<DerivedK> & K);
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "gaussian_curvature.cpp"
+#endif
+
+#endif
+

include/igl/harmonic.cpp

@@ -1,3 +1,10 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
 #include "harmonic.h"
 #include "cotmatrix.h"
 #include "massmatrix.h"
@@ -17,7 +24,16 @@ IGL_INLINE bool igl::harmonic(
   using namespace Eigen;
   SparseMatrix<double> L,M,Mi;
   cotmatrix(V,F,L);
-  massmatrix(V,F,MASSMATRIX_VORONOI,M);
+  switch(F.cols())
+  {
+    case 3:
+      massmatrix(V,F,MASSMATRIX_VORONOI,M);
+      break;
+    case 4:
+    default:
+      massmatrix(V,F,MASSMATRIX_BARYCENTRIC,M);
+      break;
+  }
   invert_diag(M,Mi);
   SparseMatrix<double> Q = -L;
   for(int p = 1;p<k;p++)

include/igl/histc.cpp

@@ -1,3 +1,10 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
 #include "histc.h"
 #include "matlab_format.h"
 #include <cassert>

include/igl/hsv_to_rgb.cpp

@@ -10,7 +10,7 @@
 
 
 template <typename T>
-void igl::hsv_to_rgb(const T * hsv, T * rgb)
+IGL_INLINE void igl::hsv_to_rgb(const T * hsv, T * rgb)
 {
   igl::hsv_to_rgb( 
       hsv[0],hsv[1],hsv[2],
@@ -18,7 +18,7 @@ void igl::hsv_to_rgb(const T * hsv, T * rgb)
 }
 
 template <typename T>
-void igl::hsv_to_rgb( 
+IGL_INLINE void igl::hsv_to_rgb( 
   const T & h, const T & s, const T & v, 
   T & r, T & g, T & b)
 {

include/igl/hsv_to_rgb.h

@@ -21,9 +21,9 @@ namespace igl
   //   g  green value ([0,1])
   //   b  blue value ([0,1])
   template <typename T>
-  void hsv_to_rgb(const T * hsv, T * rgb);
+  IGL_INLINE void hsv_to_rgb(const T * hsv, T * rgb);
   template <typename T>
-  void hsv_to_rgb( 
+  IGL_INLINE void hsv_to_rgb( 
     const T & h, const T & s, const T & v, 
     T & r, T & g, T & b);
 };

include/igl/in_element.cpp

@@ -0,0 +1,53 @@
+#include "in_element.h"
+
+IGL_INLINE void igl::in_element(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & Ele,
+  const Eigen::MatrixXd & Q,
+  const InElementAABB & aabb,
+  Eigen::VectorXi & I)
+{
+  using namespace std;
+  using namespace Eigen;
+  using namespace igl;
+  const int Qr = Q.rows();
+  I.setConstant(Qr,1,-1);
+#pragma omp parallel for
+  for(int e = 0;e<Qr;e++)
+  {
+    // find all
+    const auto R = aabb.find(V,Ele,Q.row(e),true);
+    if(!R.empty())
+    {
+      I(e) = R[0];
+    }
+  }
+}
+
+IGL_INLINE void igl::in_element(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & Ele,
+  const Eigen::MatrixXd & Q,
+  const InElementAABB & aabb,
+  Eigen::SparseMatrix<double> & I)
+{
+  using namespace std;
+  using namespace Eigen;
+  using namespace igl;
+  const int Qr = Q.rows();
+  std::vector<Triplet<double> > IJV;
+  IJV.reserve(Qr);
+#pragma omp parallel for
+  for(int e = 0;e<Qr;e++)
+  {
+    // find all
+    const auto R = aabb.find(V,Ele,Q.row(e),false);
+    for(const auto r : R)
+    {
+#pragma omp critical
+      IJV.push_back(Triplet<double>(e,r,1));
+    }
+  }
+  I.resize(Qr,Ele.rows());
+  I.setFromTriplets(IJV.begin(),IJV.end());
+}

include/igl/in_element.h

@@ -0,0 +1,46 @@
+#ifndef IGL_IN_ELEMENT_H
+#define IGL_IN_ELEMENT_H
+
+#include <igl/igl_inline.h>
+#include "InElementAABB.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  // Determine whether each point in a list of points is in the elements of a
+  // mesh.
+  //
+  // Inputs:
+  //   V  #V by dim list of mesh vertex positions. 
+  //   Ele  #Ele by dim+1 list of mesh indices into #V. 
+  //   Q  #Q by dim list of query point positions
+  //   aabb  axis-aligned bounding box tree object (see InElementAABB.h)
+  // Outputs:
+  //   I  #Q list of indices into Ele of first containing element (-1 means no
+  //     containing element)
+  IGL_INLINE void in_element(
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & Ele,
+    const Eigen::MatrixXd & Q,
+    const InElementAABB & aabb,
+    Eigen::VectorXi & I);
+  // Outputs:
+  //   I  #Q by #Ele sparse matrix revealing whether each element contains each
+  //     point: I(q,e) means point q is in element e
+  IGL_INLINE void in_element(
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & Ele,
+    const Eigen::MatrixXd & Q,
+    const InElementAABB & aabb,
+    Eigen::SparseMatrix<double> & I);
+  //
+  // Example:
+  //   InElementAABB aabb;
+  //   aabb.init(V,Ele);
+};
+
+#ifdef IGL_HEADER_ONLY
+#include "in_element.cpp"
+#endif
+
+#endif

include/igl/internal_angles.cpp

@@ -0,0 +1,48 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "internal_angles.h"
+#include "edge_lengths.h"
+
+template <typename DerivedV, typename DerivedF, typename DerivedK>
+IGL_INLINE void igl::internal_angles(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedK> & K)
+{
+  using namespace Eigen;
+  // Edge lengths
+  Matrix<
+    typename DerivedV::Scalar,
+    DerivedF::RowsAtCompileTime,
+    DerivedF::ColsAtCompileTime> L;
+  edge_lengths(V,F,L);
+
+  assert(F.cols() == 3 && "F should contain triangles");
+  return internal_angles(L,K);
+}
+
+template <typename DerivedL, typename DerivedK>
+IGL_INLINE void igl::internal_angles(
+  const Eigen::PlainObjectBase<DerivedL>& L,
+  Eigen::PlainObjectBase<DerivedK> & K)
+{
+  assert(L.cols() == 3 && "Edge-lengths should come from triangles");
+  K.resize(L.rows(),L.cols());
+  for(int d = 0;d<3;d++)
+  {
+    const auto & s1 = L.col(d).array();
+    const auto & s2 = L.col((d+1)%3).array();
+    const auto & s3 = L.col((d+2)%3).array();
+    K.col(d) = ((s3.square() + s2.square() - s1.square())/(2.*s3*s2)).acos();
+  }
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+template void igl::internal_angles<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> >&);
+#endif

include/igl/internal_angles.h

@@ -0,0 +1,40 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_INTERNAL_ANGLES_H
+#define IGL_INTERNAL_ANGLES_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Compute internal angles for a triangle mesh
+  // Inputs:
+  //   V  #V by 3 eigen Matrix of mesh vertex 3D positions
+  //   F  #F by 3 eigen Matrix of face (triangle) indices
+  // Output:
+  //   K  #F by 3 eigen Matrix of internal angles
+  //     for triangles, columns correspond to edges [1,2],[2,0],[0,1]
+  template <typename DerivedV, typename DerivedF, typename DerivedK>
+  IGL_INLINE void internal_angles(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    Eigen::PlainObjectBase<DerivedK> & K);
+  // Inputs:
+  //   L  #F by 3 list of edge lengths
+  template <typename DerivedL, typename DerivedK>
+  IGL_INLINE void internal_angles(
+    const Eigen::PlainObjectBase<DerivedL>& L,
+    Eigen::PlainObjectBase<DerivedK> & K);
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "internal_angles.cpp"
+#endif
+
+#endif
+
+

include/igl/is_planar.cpp

@@ -1,3 +1,10 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
 #include "is_planar.h"
 IGL_INLINE bool igl::is_planar(const Eigen::MatrixXd & V)
 {

include/igl/is_symmetric.cpp

@@ -15,6 +15,7 @@ IGL_INLINE bool igl::is_symmetric(const Eigen::SparseMatrix<T>& A)
   {
     return false;
   }
+  assert(A.size() != 0);
   Eigen::SparseMatrix<T> AT = A.transpose();
   Eigen::SparseMatrix<T> AmAT = A-AT;
   //// Eigen screws up something with LLT if you try to do
@@ -32,6 +33,7 @@ IGL_INLINE bool igl::is_symmetric(
   {
     return false;
   }
+  assert(A.size() != 0);
   const typename Eigen::PlainObjectBase<DerivedA>& AT = A.transpose();
   const typename Eigen::PlainObjectBase<DerivedA>& AmAT = A-AT;
   //// Eigen screws up something with LLT if you try to do
@@ -52,6 +54,7 @@ IGL_INLINE bool igl::is_symmetric(
   {
     return false;
   }
+  assert(A.size() != 0);
   SparseMatrix<AType> AT = A.transpose();
   SparseMatrix<AType> AmAT = A-AT;
   VectorXi AmATI,AmATJ;

include/igl/jet.cpp

@@ -54,13 +54,13 @@
 
 
 template <typename T>
-void igl::jet(const T x, T * rgb)
+IGL_INLINE void igl::jet(const T x, T * rgb)
 {
   igl::jet(x,rgb[0],rgb[1],rgb[2]);
 }
 
 template <typename T>
-void igl::jet(const T x, T & r, T & g, T & b)
+IGL_INLINE void igl::jet(const T x_in, T & r, T & g, T & b)
 {
   // Only important if the number of colors is small. In which case the rest is
   // still wrong anyway
@@ -69,6 +69,8 @@ void igl::jet(const T x, T & r, T & g, T & b)
   const double rone = 0.8;
   const double gone = 1.0;
   const double bone = 1.0;
+  T x = x_in;
+  x = (x_in<0?0:(x>1?1:x));
 
   if(x<1./8.)
   {
@@ -98,6 +100,21 @@ void igl::jet(const T x, T & r, T & g, T & b)
   }
 }
 
+template <typename DerivedZ, typename DerivedC>
+IGL_INLINE void igl::jet(
+  const Eigen::PlainObjectBase<DerivedZ> & Z,
+  const bool normalize,
+  Eigen::PlainObjectBase<DerivedC> & C)
+{
+  C.resize(Z.rows(),3);
+  const double min_z = (normalize?Z.minCoeff():0);
+  const double max_z = (normalize?Z.maxCoeff():-1);
+  for(int r = 0;r<Z.rows();r++)
+  {
+    jet((Z(r,0)-min_z)/(max_z-min_z),C(r,0),C(r,1),C(r,2));
+  }
+}
+
 #ifndef IGL_NO_HEADER
 template void igl::jet<double>(double, double*);
 template void igl::jet<double>(double, double&, double&, double&);

include/igl/jet.h

@@ -9,7 +9,7 @@
 #define IGL_JET_H
 #include "igl_inline.h"
 //#ifndef IGL_NO_EIGEN
-//#  include <Eigen/Dense>
+#  include <Eigen/Dense>
 //#endif
 namespace igl
 {
@@ -33,9 +33,19 @@ namespace igl
   //   g  green value
   //   b  blue value
   template <typename T>
-  void jet(const T f, T * rgb);
+  IGL_INLINE void jet(const T f, T * rgb);
   template <typename T>
-  void jet(const T f, T & r, T & g, T & b);
+  IGL_INLINE void jet(const T f, T & r, T & g, T & b);
+  // Inputs:
+  //   Z  #Z list of factos 
+  //   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 jet(
+    const Eigen::PlainObjectBase<DerivedZ> & Z,
+    const bool normalize,
+    Eigen::PlainObjectBase<DerivedC> & C);
 };
 
 #ifdef IGL_HEADER_ONLY

include/igl/kronecker_product.cpp

@@ -0,0 +1,65 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2014 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 "kronecker_product.h"
+
+// Bug in unsupported/Eigen/SparseExtra needs iostream first
+#include <iostream>
+#include <unsupported/Eigen/SparseExtra>
+
+template <typename Scalar>
+IGL_INLINE Eigen::SparseMatrix<Scalar> igl::kronecker_product(
+  const Eigen::SparseMatrix<Scalar> & A,
+  const Eigen::SparseMatrix<Scalar> & B)
+{
+  using namespace Eigen;
+  using namespace std;
+
+  // Convert B in triplets format
+  MatrixXd B_triplets(B.nonZeros(),3);
+  int count = 0;
+  for (int k=0; k<B.outerSize(); ++k)
+    for (SparseMatrix<double>::InnerIterator it(B,k); it; ++it)
+      B_triplets.row(count++) << it.row(), it.col(), it.value();
+
+  MatrixXd C_triplets(B_triplets.rows()*A.nonZeros(),3);
+  count = 0;
+  for (int k=0; k<A.outerSize(); ++k)
+    for (SparseMatrix<double>::InnerIterator it(A,k); it; ++it)
+    {
+      int i = it.row();
+      int j = it.col();
+      double v = it.value();
+
+      MatrixXd B_triplets_copy = B_triplets;
+      B_triplets_copy.col(0) = B_triplets_copy.col(0).array() + double(B.rows()*i);
+      B_triplets_copy.col(1) = B_triplets_copy.col(1).array() + double(B.cols()*j);
+      B_triplets_copy.col(2) = B_triplets_copy.col(2).array() * v;
+
+      C_triplets.block(count*B_triplets.rows(),0,
+                       B_triplets.rows(), B_triplets.cols()) = B_triplets_copy;
+
+      count++;
+    }
+
+  typedef Eigen::Triplet<double> T;
+  std::vector<T> triplets;
+  triplets.reserve(C_triplets.rows());
+  
+  for(unsigned i=0; i<C_triplets.rows(); ++i)
+    triplets.push_back(T(C_triplets(i,0),C_triplets(i,1),C_triplets(i,2)));
+  SparseMatrix<Scalar> C(A.rows()*B.rows(),A.cols()*B.cols());
+  C.setFromTriplets(triplets.begin(),triplets.end());
+  
+  return C;
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+// generated by autoexplicit.sh
+template Eigen::SparseMatrix<double, 0, int> igl::kronecker_product<double>(Eigen::SparseMatrix<double, 0, int> const&, Eigen::SparseMatrix<double, 0, int> const&);
+#endif

include/igl/kronecker_product.h

@@ -0,0 +1,34 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2014 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_KRONECKERPRODUCT_H
+#define IGL_KRONECKERPRODUCT_H
+#include "igl_inline.h"
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+namespace igl
+{
+  // Computes the Kronecker product between sparse matrices A and B.
+  //
+  // Inputs:
+  //   A  #M by #N sparse matrix
+  //   B  #P by #Q sparse matrix
+  // Returns #M*#P by #N*#Q sparse matrix
+  //
+  template <typename Scalar>
+  IGL_INLINE Eigen::SparseMatrix<Scalar> kronecker_product(
+    const Eigen::SparseMatrix<Scalar> & A,
+    const Eigen::SparseMatrix<Scalar> & B);
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "kronecker_product.cpp"
+#endif
+
+#endif

include/igl/lbs_matrix.cpp

@@ -1,3 +1,10 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
 #include "lbs_matrix.h"
 
 IGL_INLINE void igl::lbs_matrix(

include/igl/lens_flare.cpp

@@ -18,25 +18,26 @@
 
 // http://www.opengl.org/archives/resources/features/KilgardTechniques/LensFlare/glflare.c
 
-static void setup_texture(
-  const uint8_t * texture, 
-  const int width,
-  const int height,
-  GLuint texobj,
-  GLenum minFilter, GLenum maxFilter)
-{
-  glBindTexture(GL_TEXTURE_2D, texobj);
-  glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_MODULATE);
-  glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, minFilter);
-  glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, maxFilter);
-  glTexImage2D(GL_TEXTURE_2D, 0, 1, width, height, 0,
-    GL_LUMINANCE, GL_UNSIGNED_BYTE, texture);
-}
-
-void igl::lens_flare_load_textures(
+IGL_INLINE void igl::lens_flare_load_textures(
   std::vector<GLuint> & shine_id,
   std::vector<GLuint> & flare_id)
 {
+
+  const auto setup_texture =[](
+    const uint8_t * texture, 
+    const int width,
+    const int height,
+    GLuint texobj,
+    GLenum minFilter, GLenum maxFilter)
+  {
+    glBindTexture(GL_TEXTURE_2D, texobj);
+    glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_MODULATE);
+    glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, minFilter);
+    glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, maxFilter);
+    glTexImage2D(GL_TEXTURE_2D, 0, 1, width, height, 0,
+      GL_LUMINANCE, GL_UNSIGNED_BYTE, texture);
+  };
+
   glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
   shine_id.resize(10);
   glGenTextures(10,&shine_id[0]);
@@ -58,7 +59,7 @@ void igl::lens_flare_load_textures(
   }
 }
 
-void igl::lens_flare_create(
+IGL_INLINE void igl::lens_flare_create(
   const float * A,
   const float * B,
   const float * C,
@@ -84,7 +85,7 @@ void igl::lens_flare_create(
   flares[11] = Flare(5, -1.0f, 0.03f, A, 0.2);
 }
 
-void igl::lens_flare_draw(
+IGL_INLINE void igl::lens_flare_draw(
   const std::vector<igl::Flare> & flares,
   const std::vector<GLuint> & shine_ids,
   const std::vector<GLuint> & flare_ids,

include/igl/lens_flare.h

@@ -11,6 +11,7 @@
 #ifndef IGL_NO_OPENGL
 #include <igl/OpenGL_convenience.h>
 #include <Eigen/Core>
+#include "igl_inline.h"
 
 #include <vector>
 
@@ -45,7 +46,7 @@ namespace igl
   // Outputs:
   //   shine  list of texture ids for shines
   //   flare  list of texture ids for flares
-  void lens_flare_load_textures(
+  IGL_INLINE void lens_flare_load_textures(
     std::vector<GLuint> & shine_ids,
     std::vector<GLuint> & flare_ids);
   
@@ -57,7 +58,7 @@ namespace igl
   //   C  secondary color
   // Outputs:
   //   flares  list of flare objects
-  void lens_flare_create(
+  IGL_INLINE void lens_flare_create(
     const float * A,
     const float * B,
     const float * C,
@@ -74,7 +75,7 @@ namespace igl
   //   shine_tic  current "tic" in shine textures
   // Outputs:
   //   shine_tic  current "tic" in shine textures
-  void lens_flare_draw(
+  IGL_INLINE void lens_flare_draw(
     const std::vector<Flare> & flares,
     const std::vector<GLuint> & shine_ids,
     const std::vector<GLuint> & flare_ids,

include/igl/local_basis.cpp

@@ -1,3 +1,10 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
 #include "local_basis.h"
 
 #include <sstream>
@@ -7,27 +14,24 @@
 #include <vector>
 #include <Eigen/Geometry>
 
-using namespace Eigen;
-using namespace std;
 
-namespace igl
+template <typename DerivedV, typename DerivedF>
+IGL_INLINE void igl::local_basis(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedV>& B1,
+  Eigen::PlainObjectBase<DerivedV>& B2,
+  Eigen::PlainObjectBase<DerivedV>& B3
+  )
 {
+  using namespace Eigen;
+  using namespace std;
+  B1.resize(F.rows(),3);
+  B2.resize(F.rows(),3);
+  B3.resize(F.rows(),3);
 
-  template <typename DerivedV, typename DerivedF>
-  IGL_INLINE void local_basis(
-    const Eigen::PlainObjectBase<DerivedV>& V,
-    const Eigen::PlainObjectBase<DerivedF>& F,
-    Eigen::PlainObjectBase<DerivedV>& B1,
-    Eigen::PlainObjectBase<DerivedV>& B2,
-    Eigen::PlainObjectBase<DerivedV>& B3
-    )
+  for (unsigned i=0;i<F.rows();++i)
   {
-    B1.resize(F.rows(),3);
-    B2.resize(F.rows(),3);
-    B3.resize(F.rows(),3);
-
-    for (unsigned i=0;i<F.rows();++i)
-    {
       Eigen::Matrix<typename DerivedV::Scalar, 1, 3> v1 = (V.row(F(i,1)) - V.row(F(i,0))).normalized();
       Eigen::Matrix<typename DerivedV::Scalar, 1, 3> t = V.row(F(i,2)) - V.row(F(i,0));
       Eigen::Matrix<typename DerivedV::Scalar, 1, 3> v3 = v1.cross(t).normalized();
@@ -36,10 +40,7 @@ namespace igl
       B1.row(i) = v1;
       B2.row(i) = v2;
       B3.row(i) = v3;
-    }
-
   }
-
 }
 
 #ifndef IGL_HEADER_ONLY

include/igl/local_basis.h

@@ -1,3 +1,10 @@
+// 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_LOCALBASIS_H
 #define IGL_LOCALBASIS_H
 

include/igl/lscm.cpp

@@ -0,0 +1,69 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2014 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 "lscm.h"
+
+#include <igl/areamatrix.h>
+#include <igl/cotmatrix.h>
+//#include <igl/kronecker_product.h>
+#include <igl/repdiag.h>
+#include <igl/min_quad_with_fixed.h>
+#include <igl/matlab_format.h>
+#include <iostream>
+
+IGL_INLINE void igl::lscm(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  const Eigen::VectorXi& b,
+  const Eigen::MatrixXd& bc,
+  Eigen::MatrixXd & V_uv)
+{
+  using namespace Eigen;
+  using namespace std;
+  
+  // Assemble the area matrix (note that A is #Vx2 by #Vx2)
+  SparseMatrix<double> A;
+  igl::areamatrix(F,A);
+
+  // Assemble the cotan laplacian matrix
+  SparseMatrix<double> L;
+  igl::cotmatrix(V,F,L);
+
+  SparseMatrix<double> L_flat;
+  repdiag(L,2,L_flat);
+
+  VectorXi b_flat(b.size()*bc.cols(),1);
+  VectorXd bc_flat(bc.size(),1);
+  for(int c = 0;c<bc.cols();c++)
+  {
+    b_flat.block(c*b.size(),0,b.rows(),1) = c*V.rows() + b.array();
+    bc_flat.block(c*bc.rows(),0,bc.rows(),1) = bc.col(c);
+  }
+  
+  // Minimize the LSCM energy
+  SparseMatrix<double> Q = -0.5*L_flat + A;
+  const VectorXd B_flat = VectorXd::Zero(V.rows()*2);
+  igl::min_quad_with_fixed_data<double> data;
+  igl::min_quad_with_fixed_precompute(Q,b_flat,SparseMatrix<double>(),true,data);
+
+  MatrixXd W_flat;
+  if(!min_quad_with_fixed_solve(data,B_flat,bc_flat,VectorXd(),W_flat))
+    assert(false);
+
+
+  assert(W_flat.rows() == V.rows()*2);
+  V_uv.resize(V.rows(),2);
+  for (unsigned i=0;i<V_uv.cols();++i)
+  {
+    V_uv.col(i) = W_flat.block(V_uv.rows()*i,0,V_uv.rows(),1);
+  }
+
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

include/igl/lscm.h

@@ -0,0 +1,48 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2014 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_LSCM_H
+#define IGL_LSCM_H
+#include "igl_inline.h"
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+namespace igl
+{
+  // Compute a Least-squares conformal map parametrization following the algorithm
+  // presented in: Spectral Conformal Parameterization,
+  //               Patrick Mullen, Yiying Tong, Pierre Alliez and Mathieu Desbrun
+  //
+  // Inputs:
+  //   V  #V by 3 list of mesh vertex positions
+  //   F  #F by 3 list of mesh faces (must be triangles)
+  //   b  #b boundary indices into V
+  //   bc #b by 3 list of boundary values
+  // Outputs:
+  //   UV #V by 2 list of 2D mesh vertex positions in UV space
+  //
+  // X  Note: if b and bc are empty, lscm automatically removes the null space
+  // X        by fixing two farthest points on the boundary
+  // Boundary conditions are part of the api. It would be strange to secretly
+  // use other boundary conditions. This is also weird if there's more than one
+  // boundary loop.
+  // X  Note: (V,F) must be a genus-0 mesh, with a single boundary
+  // No longer the case. Should probably just be manifold.
+  IGL_INLINE void lscm(
+    const Eigen::MatrixXd& V,
+    const Eigen::MatrixXi& F,
+    const Eigen::VectorXi& b,
+    const Eigen::MatrixXd& bc,
+    Eigen::MatrixXd& V_uv);
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "lscm.cpp"
+#endif
+
+#endif