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

Conflicts:
	tutorial/cmake/FindMATLAB.cmake

Former-commit-id: 9e3660f49eb8ad27bf896362aa6c2ac10aa73728

README.md

@@ -1,16 +1,19 @@
 libigl - A simple c++ geometry processing library
 =================================================
 
-<http://igl.ethz.ch/projects/libigl/>
-<https://github.com/alecjacobson/libigl/>
+![](tutorial/images/libigl-logo.jpg)
+
+<http://libigl.github.io/libigl/>
+<https://github.com/libigl/libigl/>
 
 Copyright 2014 - Alec Jacobson, Daniele Panozzo, Olga Diamanti, Kenshi
 Takayama, Leo Sacht, Wenzel Jacob, etc.
 
-This is first and foremost a *header* library. Each header file should contain
-a single function.  The function may have multiple prototypes. All functions
-should use the igl namespace and should adhere to the conventions and styles
-listed below.
+Libigl is first and foremost a *header* library. Each header file should
+contain a single function.  This function may have multiple overloads and
+prototypes. All functions should use the `igl::` namespace and should adhere to
+the conventions and styles listed in the [style
+guidelines](style_guidelines.html).
 
 > **New:** As of 1 July 2014, we have released libigl as beta version 1.0.
 > There are a number of changes we collected for this release to minimize
@@ -63,8 +66,9 @@ Hello, mesh:
 
 ## Tutorial ##
 
-As of version 1.0, libigl includes an introductory tutorial that covers its basic
-functionalities. See [tutorial/tutorial.md](./tutorial/tutorial.md) to get started.
+As of version 1.0, libigl includes an introductory
+[tutorial](http://libigl.github.io/libigl/tutorial/tutorial.html) that covers
+its basic functionalities.
 
 ## Dependencies ##
 - Eigen3  Last tested with Eigen Version 3.2
@@ -405,7 +409,7 @@ BibTeX entry:
 @misc{libigl,
   title = {{libigl}: A simple {C++} geometry processing library},
   author = {Alec Jacobson and Daniele Panozzo and others},
-  note = {http://igl.ethz.ch/projects/libigl/},
+  note = {http://libigl.github.io/libigl/},
   year = {2014},
 }
 ```

RELEASE_HISTORY.txt

@@ -1,3 +1,4 @@
+1.1.0  Mesh boolean operations using CGAL and cork, implementing [Attene 14]
 1.0.3  Bone heat method
 1.0.2  Bug fix in winding number code
 1.0.1  Bug fixes and more CGAL support

VERSION.txt

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

build/Makefile

@@ -20,6 +20,9 @@ EXTRAS=
 ifeq ($(IGL_WITH_BBW),1)
 	EXTRAS += bbw
 endif
+ifeq ($(IGL_WITH_BOOLEAN),1)
+	EXTRAS += boolean
+endif
 ifeq ($(IGL_WITH_BOOST),1)
 	EXTRAS += boost
 endif

build/Makefile.conf

@@ -50,6 +50,7 @@ ifeq ($(IGL_USERNAME),ajx)
 	MOSEKVERSION=7
 	IGL_WITH_VIEWER=1
 	IGL_WITH_TETGEN=1
+	IGL_WITH_BOOLEAN=1
 	IGL_WITH_EMBREE=1
 	IGL_WITH_MATLAB=1
 	IGL_WITH_MOSEK=1

build/Makefile_boolean

@@ -0,0 +1,53 @@
+include Makefile.conf
+
+.PHONY: all
+all: libiglboolean
+debug: libiglboolean
+
+include Makefile.conf
+all: CFLAGS += -O3 -DNDEBUG -std=c++11
+debug: CFLAGS += -g -Wall -std=c++11
+
+.PHONY: libiglboolean
+libiglboolean: obj ../lib/libiglboolean.a
+
+SRC_DIR=../include/igl/boolean/
+CPP_FILES=$(wildcard $(SRC_DIR)*.cpp)
+OBJ_FILES=$(addprefix obj/,$(notdir $(CPP_FILES:.cpp=.o)))
+
+# include igl headers
+INC+=-I../include/
+
+# EXPECTS THAT CFLAGS IS ALREADY SET APPROPRIATELY 
+
+# Eigen dependency
+EIGEN3_INC=-I$(DEFAULT_PREFIX)/include/eigen3 -I$(DEFAULT_PREFIX)/include/eigen3/unsupported
+INC+=$(EIGEN3_INC)
+
+# cork dependency
+#CORK=../external/cork
+#CORK_INC=-I$(cork)/common -I$(cork)/rtcore/
+CORK=../external/cork/
+CORK_INC=-I$(CORK)/include/
+INC+=$(CORK_INC)
+
+# CGAL dependency
+CGAL=$(DEFAULT_PREFIX)
+CGAL_INC=-I$(CGAL)/include
+CGAL_FLAGS=-frounding-math -fsignaling-nans 
+CFLAGS+=$(CGAL_FLAGS)
+INC+=$(CGAL_INC)
+
+obj: 
+	mkdir -p obj
+
+../lib/libiglboolean.a: $(OBJ_FILES)
+	rm -f $@
+	ar cqs $@ $(OBJ_FILES)
+
+obj/%.o: $(SRC_DIR)/%.cpp $(SRC_DIR)/%.h
+	g++ $(AFLAGS) $(OPENMP) $(CFLAGS) -c -o $@ $< $(INC)
+
+clean:
+	rm -f $(OBJ_FILES)
+	rm -f ../lib/libiglboolean.a

include/igl/all_edges.cpp

@@ -46,4 +46,5 @@ IGL_INLINE void igl::all_edges(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 template void igl::all_edges<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::all_edges<Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 2, 0, -1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&);
 #endif

include/igl/boolean/MeshBooleanType.h

@@ -0,0 +1,17 @@
+#ifndef MESH_BOOLEAN_TYPE_H
+#define MESH_BOOLEAN_TYPE_H
+
+namespace igl
+{
+  enum MeshBooleanType
+  {
+    MESH_BOOLEAN_TYPE_UNION = 0,
+    MESH_BOOLEAN_TYPE_INTERSECT = 1,
+    MESH_BOOLEAN_TYPE_MINUS = 2,
+    MESH_BOOLEAN_TYPE_XOR = 3,
+    MESH_BOOLEAN_TYPE_RESOLVE = 4,
+    NUM_MESH_BOOLEAN_TYPES = 5
+  };
+};
+
+#endif

include/igl/boolean/from_cork_mesh.cpp

@@ -0,0 +1,37 @@
+#ifndef IGL_NO_CORK
+#include "from_cork_mesh.h"
+
+template <
+  typename DerivedV,
+  typename DerivedF>
+IGL_INLINE void igl::from_cork_mesh(
+  const CorkTriMesh & mesh,
+  Eigen::PlainObjectBase<DerivedV > & V,
+  Eigen::PlainObjectBase<DerivedF > & F)
+{
+  using namespace std;
+  F.resize(mesh.n_triangles,3);
+  V.resize(mesh.n_vertices,3);
+  for(size_t v = 0;v<mesh.n_vertices;v++)
+  {
+    for(size_t c = 0;c<3;c++)
+    {
+      V(v,c) = mesh.vertices[v*3+c];
+    }
+  }
+  for(size_t f = 0;f<mesh.n_triangles;f++)
+  {
+    for(size_t c = 0;c<3;c++)
+    {
+      F(f,c) = mesh.triangles[f*3+c];
+    }
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+template void igl::from_cork_mesh<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(CorkTriMesh const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::from_cork_mesh<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(CorkTriMesh const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
+#endif
+
+#endif

include/igl/boolean/from_cork_mesh.h

@@ -0,0 +1,28 @@
+#ifndef IGL_FROM_CORK_MESH_H
+#define IGL_FROM_CORK_MESH_H
+#ifndef IGL_NO_CORK
+#include <igl/igl_inline.h>
+#include <cork.h>
+#include <Eigen/Core>
+namespace igl
+{
+  // Convert cork's triangle mesh representation to a (V,F) mesh.
+  //
+  // Inputs:
+  //   mesh  cork representation of mesh
+  // Outputs:
+  //   V  #V by 3 list of vertex positions
+  //   F  #F by 3 list of triangle indices into V
+  template <
+    typename DerivedV,
+    typename DerivedF>
+  IGL_INLINE void from_cork_mesh(
+    const CorkTriMesh & mesh,
+    Eigen::PlainObjectBase<DerivedV > & V,
+    Eigen::PlainObjectBase<DerivedF > & F);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "from_cork_mesh.cpp"
+#endif
+#endif
+#endif

include/igl/boolean/mesh_boolean.cpp

@@ -0,0 +1,232 @@
+#include "mesh_boolean.h"
+#include <igl/peal_outer_hull_layers.h>
+#include <igl/cgal/remesh_self_intersections.h>
+#include <igl/remove_unreferenced.h>
+#include <igl/unique_simplices.h>
+#include <iostream>
+
+template <
+  typename DerivedVA,
+  typename DerivedFA,
+  typename DerivedVB,
+  typename DerivedFB,
+  typename DerivedVC,
+  typename DerivedFC,
+  typename DerivedJ>
+IGL_INLINE void igl::mesh_boolean(
+  const Eigen::PlainObjectBase<DerivedVA > & VA,
+  const Eigen::PlainObjectBase<DerivedFA > & FA,
+  const Eigen::PlainObjectBase<DerivedVB > & VB,
+  const Eigen::PlainObjectBase<DerivedFB > & FB,
+  const MeshBooleanType & type,
+  Eigen::PlainObjectBase<DerivedVC > & VC,
+  Eigen::PlainObjectBase<DerivedFC > & FC,
+  Eigen::PlainObjectBase<DerivedJ > & J)
+{
+  const std::function<void(
+    const Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3> &,
+    const Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&,
+          Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3> &,
+          Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&, 
+          Eigen::Matrix<typename DerivedJ::Scalar, Eigen::Dynamic,1>&)> 
+    empty_fun;
+  return mesh_boolean(VA,FA,VB,FB,type,empty_fun,VC,FC,J);
+}
+
+template <
+  typename DerivedVA,
+  typename DerivedFA,
+  typename DerivedVB,
+  typename DerivedFB,
+  typename DerivedVC,
+  typename DerivedFC>
+IGL_INLINE void igl::mesh_boolean(
+  const Eigen::PlainObjectBase<DerivedVA > & VA,
+  const Eigen::PlainObjectBase<DerivedFA > & FA,
+  const Eigen::PlainObjectBase<DerivedVB > & VB,
+  const Eigen::PlainObjectBase<DerivedFB > & FB,
+  const MeshBooleanType & type,
+  Eigen::PlainObjectBase<DerivedVC > & VC,
+  Eigen::PlainObjectBase<DerivedFC > & FC)
+{
+  Eigen::Matrix<typename DerivedFC::Index, Eigen::Dynamic,1> J;
+  const std::function<void(
+    const Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3> &,
+    const Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&,
+          Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3> &,
+          Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&, 
+          Eigen::Matrix<typename DerivedFC::Index, Eigen::Dynamic,1>&)> 
+    empty_fun;
+  return mesh_boolean(VA,FA,VB,FB,type,empty_fun,VC,FC,J);
+}
+
+template <
+  typename DerivedVA,
+  typename DerivedFA,
+  typename DerivedVB,
+  typename DerivedFB,
+  typename DerivedVC,
+  typename DerivedFC,
+  typename DerivedJ>
+IGL_INLINE void igl::mesh_boolean(
+  const Eigen::PlainObjectBase<DerivedVA > & VA,
+  const Eigen::PlainObjectBase<DerivedFA > & FA,
+  const Eigen::PlainObjectBase<DerivedVB > & VB,
+  const Eigen::PlainObjectBase<DerivedFB > & FB,
+  const MeshBooleanType & type,
+  const std::function<void(
+    const Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3>&,
+    const Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&,
+          Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3>&,
+          Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&,
+          Eigen::Matrix<typename DerivedJ::Scalar, Eigen::Dynamic,1>&)> 
+    & resolve_fun,
+  Eigen::PlainObjectBase<DerivedVC > & VC,
+  Eigen::PlainObjectBase<DerivedFC > & FC,
+  Eigen::PlainObjectBase<DerivedJ > & J)
+{
+  using namespace Eigen;
+  using namespace std;
+  using namespace igl;
+  MeshBooleanType eff_type = type;
+  // Concatenate A and B into a single mesh
+  typedef typename DerivedVC::Scalar Scalar;
+  typedef typename DerivedFC::Scalar Index;
+  typedef Matrix<Scalar,Dynamic,3> MatrixX3S;
+  typedef Matrix<Index,Dynamic,3> MatrixX3I;
+  typedef Matrix<Index,Dynamic,2> MatrixX2I;
+  typedef Matrix<Index,Dynamic,1> VectorXI;
+  typedef Matrix<typename DerivedJ::Scalar,Dynamic,1> VectorXJ;
+  MatrixX3S V(VA.rows()+VB.rows(),3);
+  MatrixX3I F(FA.rows()+FB.rows(),3);
+  V.block(0,0,VA.rows(),VA.cols()) = VA;
+  V.block(VA.rows(),0,VB.rows(),VB.cols()) = VB;
+  switch(type)
+  {
+    // Minus is implemented by flipping B and computing union
+    case MESH_BOOLEAN_TYPE_MINUS:
+      F.block(0,0,FA.rows(),FA.cols()) = FA.rowwise().reverse();
+      F.block(FA.rows(),0,FB.rows(),FB.cols()) = FB.array()+VA.rows();
+      //F.block(0,0,FA.rows(),3) = FA;
+      //F.block(FA.rows(),0,FB.rows(),3) = 
+      //  FB.rowwise().reverse().array()+VA.rows();
+      eff_type = MESH_BOOLEAN_TYPE_INTERSECT;
+      break;
+    default:
+      F.block(0,0,FA.rows(),FA.cols()) = FA;
+      F.block(FA.rows(),0,FB.rows(),FB.cols()) = FB.array()+VA.rows();
+      break;
+  }
+
+  // Resolve intersections (assumes A and B are solid)
+  const auto & libigl_resolve = [](
+    const MatrixX3S & V,
+    const MatrixX3I & F,
+    MatrixX3S & CV,
+    MatrixX3I & CF,
+    VectorXJ & J)
+  {
+    MatrixX3S SV;
+    MatrixX3I SF;
+    MatrixX2I SIF;
+    VectorXI SIM,UIM;
+    RemeshSelfIntersectionsParam params;
+    remesh_self_intersections(V,F,params,SV,SF,SIF,J,SIM);
+    for_each(SF.data(),SF.data()+SF.size(),[&SIM](int & a){a=SIM(a);});
+    {
+      remove_unreferenced(SV,SF,CV,CF,UIM);
+    }
+  };
+
+  MatrixX3S CV;
+  MatrixX3I CF;
+  VectorXJ CJ;
+  if(resolve_fun)
+  {
+    resolve_fun(V,F,CV,CF,CJ);
+  }else
+  {
+    libigl_resolve(V,F,CV,CF,CJ);
+  }
+
+  if(type == MESH_BOOLEAN_TYPE_RESOLVE)
+  {
+    FC = CF;
+    VC = CV;
+    J = CJ;
+    return;
+  }
+
+  Matrix<bool,Dynamic,1> from_A(CF.rows());
+  // Peal layers keeping track of odd and even flips
+  Matrix<bool,Dynamic,1> odd;
+  Matrix<bool,Dynamic,1> flip;
+  peal_outer_hull_layers(CV,CF,odd,flip);
+
+  const Index m = CF.rows();
+  // Faces of output vG[i] = j means ith face of output should be jth face in F
+  std::vector<Index> vG;
+  // Whether faces of output should be flipped, Gflip[i] = true means ith face
+  // of output should be F.row(vG[i]).reverse() rather than F.row(vG[i])
+  std::vector<bool> Gflip;
+  for(Index f = 0;f<m;f++)
+  {
+    switch(eff_type)
+    {
+      case MESH_BOOLEAN_TYPE_XOR:
+      case MESH_BOOLEAN_TYPE_UNION:
+        if((odd(f)&&!flip(f))||(!odd(f)&&flip(f)))
+        {
+          vG.push_back(f);
+          Gflip.push_back(false);
+        }else if(eff_type == MESH_BOOLEAN_TYPE_XOR)
+        {
+          vG.push_back(f);
+          Gflip.push_back(true);
+        }
+        break;
+      case MESH_BOOLEAN_TYPE_INTERSECT:
+        if((!odd(f) && !flip(f)) || (odd(f) && flip(f)))
+        {
+          vG.push_back(f);
+          Gflip.push_back(type == MESH_BOOLEAN_TYPE_MINUS);
+        }
+        break;
+      default:
+        assert(false && "Unknown type");
+        return;
+    }
+  }
+  const Index gm = vG.size();
+  MatrixX3I G(gm,3);
+  VectorXi GJ(gm,1);
+  for(Index g = 0;g<gm;g++)
+  {
+    G.row(g) = Gflip[g] ? CF.row(vG[g]).reverse().eval() : CF.row(vG[g]);
+    GJ(g) = CJ(vG[g]);
+  }
+  // remove duplicates: cancel out in all cases? assertion in others?
+  {
+    MatrixXi oldG = G;
+    VectorXi IA,_1;
+    unique_simplices(oldG,G,IA,_1);
+    assert(IA.rows() == G.rows());
+    J.resize(IA.rows(),1);
+    for(size_t j = 0;j<(size_t)J.size();j++)
+    {
+      J(j) = GJ(IA(j));
+    }
+  }
+  // remove unreferenced vertices
+  VectorXi newIM;
+  remove_unreferenced(CV,G,VC,FC,newIM);
+  //cerr<<"warning not removing unref"<<endl;
+  //VC = CV;
+  //FC = G;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+template void igl::mesh_boolean<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::MeshBooleanType const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::mesh_boolean<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::MeshBooleanType const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/boolean/mesh_boolean.h

@@ -0,0 +1,91 @@
+#ifndef MESH_BOOLEAN_H
+#define MESH_BOOLEAN_H
+
+#include <igl/igl_inline.h>
+#include "MeshBooleanType.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  //  MESH_BOOLEAN Compute boolean csg operations on "solid", consistently oriented
+  //  meshes.
+  // 
+  //  Inputs:
+  //    V  #V by 3 list of vertex positions of first mesh
+  //    F  #F by 3 list of triangle indices into V
+  //    U  #U by 3 list of vertex positions of second mesh
+  //    G  #G by 3 list of triangle indices into U
+  //    type  type of boolean operation
+  //  Outputs:
+  //    W  #W by 3 list of vertex positions of boolean result mesh
+  //    H  #H by 3 list of triangle indices into W
+  //    J  #H list of indices into [FA;FB] revealing "birth" facet
+  //  
+  //  See also: self_intersect
+  //     
+  template <
+    typename DerivedVA,
+    typename DerivedFA,
+    typename DerivedVB,
+    typename DerivedFB,
+    typename DerivedVC,
+    typename DerivedFC,
+    typename DerivedJ>
+  IGL_INLINE void mesh_boolean(
+    const Eigen::PlainObjectBase<DerivedVA > & VA,
+    const Eigen::PlainObjectBase<DerivedFA > & FA,
+    const Eigen::PlainObjectBase<DerivedVB > & VB,
+    const Eigen::PlainObjectBase<DerivedFB > & FB,
+    const MeshBooleanType & type,
+    Eigen::PlainObjectBase<DerivedVC > & VC,
+    Eigen::PlainObjectBase<DerivedFC > & FC,
+    Eigen::PlainObjectBase<DerivedJ > & J);
+  template <
+    typename DerivedVA,
+    typename DerivedFA,
+    typename DerivedVB,
+    typename DerivedFB,
+    typename DerivedVC,
+    typename DerivedFC>
+  IGL_INLINE void mesh_boolean(
+    const Eigen::PlainObjectBase<DerivedVA > & VA,
+    const Eigen::PlainObjectBase<DerivedFA > & FA,
+    const Eigen::PlainObjectBase<DerivedVB > & VB,
+    const Eigen::PlainObjectBase<DerivedFB > & FB,
+    const MeshBooleanType & type,
+    Eigen::PlainObjectBase<DerivedVC > & VC,
+    Eigen::PlainObjectBase<DerivedFC > & FC);
+  // Inputs:
+  //   resolve_fun  function handle for computing resolve of a
+  //     self-intersections of a mesh and outputting the new mesh.
+  template <
+    typename DerivedVA,
+    typename DerivedFA,
+    typename DerivedVB,
+    typename DerivedFB,
+    typename DerivedVC,
+    typename DerivedFC,
+    typename DerivedJ>
+  IGL_INLINE void mesh_boolean(
+    const Eigen::PlainObjectBase<DerivedVA > & VA,
+    const Eigen::PlainObjectBase<DerivedFA > & FA,
+    const Eigen::PlainObjectBase<DerivedVB > & VB,
+    const Eigen::PlainObjectBase<DerivedFB > & FB,
+    const MeshBooleanType & type,
+    const std::function<void(
+      const Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3> &,
+      const Eigen::Matrix<typename DerivedFC::Scalar,Eigen::Dynamic,3> &,
+            Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3> &,
+            Eigen::Matrix<typename DerivedFC::Scalar,Eigen::Dynamic,3> &,
+            Eigen::Matrix<typename  DerivedJ::Scalar,Eigen::Dynamic,1>&)> 
+      & resolve_fun,
+    Eigen::PlainObjectBase<DerivedVC > & VC,
+    Eigen::PlainObjectBase<DerivedFC > & FC,
+    Eigen::PlainObjectBase<DerivedJ > & J);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "mesh_boolean.cpp"
+#endif
+
+#endif

include/igl/boolean/mesh_boolean_cork.cpp

@@ -0,0 +1,94 @@
+#ifndef IGL_NO_CORK
+#include "mesh_boolean_cork.h"
+#include "to_cork_mesh.h"
+#include "from_cork_mesh.h"
+
+template <
+  typename DerivedVA,
+  typename DerivedFA,
+  typename DerivedVB,
+  typename DerivedFB,
+  typename DerivedVC,
+  typename DerivedFC>
+IGL_INLINE void igl::mesh_boolean_cork(
+  const Eigen::PlainObjectBase<DerivedVA > & VA,
+  const Eigen::PlainObjectBase<DerivedFA > & FA,
+  const Eigen::PlainObjectBase<DerivedVB > & VB,
+  const Eigen::PlainObjectBase<DerivedFB > & FB,
+  const MeshBooleanType & type,
+  Eigen::PlainObjectBase<DerivedVC > & VC,
+  Eigen::PlainObjectBase<DerivedFC > & FC)
+{
+  CorkTriMesh A,B,C;
+  // pointer to output so it's easy to redirect on degenerate cases
+  CorkTriMesh *ret = &C;
+  to_cork_mesh(VA,FA,A);
+  to_cork_mesh(VB,FB,B);
+  switch(type)
+  {
+    case MESH_BOOLEAN_TYPE_UNION:
+      if(A.n_triangles == 0)
+      {
+        ret = &B;
+      }else if(B.n_triangles == 0)
+      {
+        ret = &A;
+      }else
+      {
+        computeUnion(A,B,ret);
+      }
+      break;
+    case MESH_BOOLEAN_TYPE_INTERSECT:
+      if(A.n_triangles == 0 || B.n_triangles == 0)
+      {
+        ret->n_triangles = 0;
+        ret->n_vertices = 0;
+      }else
+      {
+        computeIntersection(A,B,ret);
+      }
+      break;
+    case MESH_BOOLEAN_TYPE_MINUS:
+      if(A.n_triangles == 0)
+      {
+        ret->n_triangles = 0;
+        ret->n_vertices = 0;
+      }else if(B.n_triangles == 0)
+      {
+        ret = &A;
+      }else
+      {
+        computeDifference(A,B,ret);
+      }
+      break;
+    case MESH_BOOLEAN_TYPE_XOR:
+      if(A.n_triangles == 0)
+      {
+        ret = &B;
+      }else if(B.n_triangles == 0)
+      {
+        ret = &A;
+      }else
+      {
+        computeSymmetricDifference(A,B,&C);
+      }
+      break;
+    case MESH_BOOLEAN_TYPE_RESOLVE:
+      resolveIntersections(A,B,&C);
+      break;
+    default:
+      assert(false && "Unknown type");
+      return;
+  }
+  from_cork_mesh(*ret,VC,FC);
+  freeCorkTriMesh(&A);
+  freeCorkTriMesh(&B);
+  freeCorkTriMesh(&C);
+}
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+template void igl::mesh_boolean_cork<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, MeshBooleanType const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::mesh_boolean_cork<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, igl::MeshBooleanType const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
+#endif
+
+#endif

include/igl/boolean/mesh_boolean_cork.h

@@ -0,0 +1,44 @@
+#ifndef MESH_BOOLEAN_CORK_H
+#define MESH_BOOLEAN_CORK_H
+#ifndef IGL_NO_CORK
+#include "MeshBooleanType.h"
+#include <igl/igl_inline.h>
+#include <Eigen/Core>
+#include <cork.h> // for consistent uint
+
+namespace igl
+{
+  // Compute a boolean operation on two input meshes using the cork library.
+  //
+  // Inputs:
+  //   VA  #VA by 3 list of vertex positions of first mesh
+  //   FA  #FA by 3 list of triangle indices into VA
+  //   VB  #VB by 3 list of vertex positions of second mesh
+  //   FB  #FB by 3 list of triangle indices into VB
+  //   type  of boolean operation see MeshBooleanType.h
+  // Outputs:
+  //   VC  #VC by 3 list of vertex positions of output mesh
+  //   FC  #FC by 3 list of triangle indices into VC
+  template <
+    typename DerivedVA,
+    typename DerivedFA,
+    typename DerivedVB,
+    typename DerivedFB,
+    typename DerivedVC,
+    typename DerivedFC>
+  IGL_INLINE void mesh_boolean_cork(
+    const Eigen::PlainObjectBase<DerivedVA > & VA,
+    const Eigen::PlainObjectBase<DerivedFA > & FA,
+    const Eigen::PlainObjectBase<DerivedVB > & VB,
+    const Eigen::PlainObjectBase<DerivedFB > & FB,
+    const MeshBooleanType & type,
+    Eigen::PlainObjectBase<DerivedVC > & VC,
+    Eigen::PlainObjectBase<DerivedFC > & FC);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "mesh_boolean_cork.cpp"
+#endif
+#endif
+
+#endif

include/igl/boolean/to_cork_mesh.cpp

@@ -0,0 +1,38 @@
+#ifndef IGL_NO_CORK
+#include "to_cork_mesh.h"
+template <
+  typename DerivedV,
+  typename DerivedF>
+IGL_INLINE void igl::to_cork_mesh(
+  const Eigen::PlainObjectBase<DerivedV > & V,
+  const Eigen::PlainObjectBase<DerivedF > & F,
+  CorkTriMesh & mesh)
+{
+  using namespace std;
+  assert((F.cols() == 0 || F.cols() == 3) && "Facets should be triangles.");
+  assert((V.cols() == 0 || V.cols() == 3) && "Vertices should be in 3D.");
+  mesh.n_triangles = F.rows();
+  mesh.n_vertices = V.rows();
+  mesh.vertices = new double[mesh.n_vertices*3];
+  mesh.triangles = new uint[mesh.n_triangles*3];
+  for(size_t v = 0;v<mesh.n_vertices;v++)
+  {
+    for(size_t c = 0;c<3;c++)
+    {
+      mesh.vertices[v*3+c] = V(v,c);
+    }
+  }
+  for(size_t f = 0;f<mesh.n_triangles;f++)
+  {
+    for(size_t c = 0;c<3;c++)
+    {
+      mesh.triangles[f*3+c] = F(f,c);
+    }
+  }
+}
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+template void igl::to_cork_mesh<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, CorkTriMesh&);
+template void igl::to_cork_mesh<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, CorkTriMesh&);
+#endif
+#endif

include/igl/boolean/to_cork_mesh.h

@@ -0,0 +1,28 @@
+#ifndef IGL_TO_CORK_MESH_H
+#define IGL_TO_CORK_MESH_H
+#ifndef IGL_NO_CORK
+#include <igl/igl_inline.h>
+#include <cork.h>
+#include <Eigen/Core>
+namespace igl
+{
+  // Convert a (V,F) mesh to a cork's triangle mesh representation.
+  //
+  // Inputs:
+  //   V  #V by 3 list of vertex positions
+  //   F  #F by 3 list of triangle indices into V
+  // Outputs:
+  //   mesh  cork representation of mesh
+  template <
+    typename DerivedV,
+    typename DerivedF>
+  IGL_INLINE void to_cork_mesh(
+    const Eigen::PlainObjectBase<DerivedV > & V,
+    const Eigen::PlainObjectBase<DerivedF > & F,
+    CorkTriMesh & mesh);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "to_cork_mesh.cpp"
+#endif
+#endif
+#endif

include/igl/cgal/SelfIntersectMesh.h

@@ -85,9 +85,11 @@ namespace igl
       // Number of self-intersecting triangle pairs
       typedef typename DerivedF::Index Index;
       Index count;
-      std::vector<std::list<CGAL::Object> > F_objects;
+      typedef std::vector<CGAL::Object> ObjectList;
+      std::vector<ObjectList > F_objects;
       Triangles T;
-      std::list<Index> lIF;
+      typedef std::vector<Index> IndexList;
+      IndexList lIF;
       std::vector<bool> offensive;
       std::vector<Index> offending_index;
       std::vector<Index> offending;
@@ -95,7 +97,7 @@ namespace igl
       typedef std::pair<Index,Index> EMK;
       // Make a short name for the type stored at each edge, the edge map's
       // value
-      typedef std::list<Index> EMV;
+      typedef std::vector<Index> EMV;
       // Make a short name for the edge map
       typedef std::map<EMK,EMV> EdgeMap;
       EdgeMap edge2faces;
@@ -205,11 +207,11 @@ namespace igl
       //   cdt  Contrained delaunay triangulation in projected 2D plane
       inline void projected_delaunay(
           const Triangle_3 & A,
-          const std::list<CGAL::Object> & A_objects_3,
+          const ObjectList & A_objects_3,
           CDT_plus_2 & cdt);
       // Getters:
     public:
-      //const std::list<int>& get_lIF() const{ return lIF;}
+      //const IndexList& get_lIF() const{ return lIF;}
       static inline void box_intersect(
         SelfIntersectMesh * SIM, 
         const SelfIntersectMesh::Box &a, 
@@ -222,6 +224,7 @@ namespace igl
 #include "mesh_to_cgal_triangle_list.h"
 
 #include <igl/REDRUM.h>
+#include <igl/get_seconds.h>
 #include <igl/C_STR.h>
 
 #include <boost/function.hpp>
@@ -322,8 +325,19 @@ inline igl::SelfIntersectMesh<
 {
   using namespace std;
   using namespace Eigen;
+
+  //const auto & tictoc = []()
+  //{
+  //  static double t_start = igl::get_seconds();
+  //  double diff = igl::get_seconds()-t_start;
+  //  t_start += diff;
+  //  return diff;
+  //};
+  //tictoc();
+
   // Compute and process self intersections
   mesh_to_cgal_triangle_list(V,F,T);
+  //cout<<"mesh_to_cgal_triangle_list: "<<tictoc()<<endl;
   // http://www.cgal.org/Manual/latest/doc_html/cgal_manual/Box_intersection_d/Chapter_main.html#Section_63.5 
   // Create the corresponding vector of bounding boxes
   std::vector<Box> boxes;
@@ -338,6 +352,7 @@ inline igl::SelfIntersectMesh<
   // Leapfrog callback
   boost::function<void(const Box &a,const Box &b)> cb
     = boost::bind(&box_intersect, this, _1,_2);
+  //cout<<"boxes and bind: "<<tictoc()<<endl;
   // Run the self intersection algorithm with all defaults
   try{
     CGAL::box_self_intersection_d(boxes.begin(), boxes.end(),cb);
@@ -350,6 +365,7 @@ inline igl::SelfIntersectMesh<
     }
     // Otherwise just fall through
   }
+  //cout<<"box_self_intersection_d: "<<tictoc()<<endl;
 
   // Convert lIF to Eigen matrix
   assert(lIF.size()%2 == 0);
@@ -357,7 +373,7 @@ inline igl::SelfIntersectMesh<
   {
     Index i=0;
     for(
-      typename list<Index>::const_iterator ifit = lIF.begin();
+      typename IndexList::const_iterator ifit = lIF.begin();
       ifit!=lIF.end();
       )
     {
@@ -368,6 +384,7 @@ inline igl::SelfIntersectMesh<
       i++;
     }
   }
+  //cout<<"IF: "<<tictoc()<<endl;
 
   if(params.detect_only)
   {
@@ -380,7 +397,8 @@ inline igl::SelfIntersectMesh<
     typename Eigen::Matrix<typename DerivedFF::Scalar,Dynamic,Dynamic>
     > NF(offending.size());
   // list of new vertices
-  list<Point_3> NV;
+  typedef vector<Point_3> Point_3List;
+  Point_3List NV;
   Index NV_count = 0;
   vector<CDT_plus_2> cdt(offending.size());
   vector<Plane_3> P(offending.size());
@@ -391,7 +409,9 @@ inline igl::SelfIntersectMesh<
   {
     // index in F
     const Index f = offending[o];
-    projected_delaunay(T[f],F_objects[f],cdt[o]);
+    {
+      projected_delaunay(T[f],F_objects[f],cdt[o]);
+    }
     // Q: Is this also delaunay in 3D?
     // A: No, because the projection is affine and delaunay is not affine
     // invariant.
@@ -439,7 +459,7 @@ inline igl::SelfIntersectMesh<
             assert(edge2faces.count(EMK(i,j))==1);
             // loop over neighbors
             for(
-              typename list<Index>::const_iterator nit = edge2faces[EMK(i,j)].begin();
+              typename IndexList::const_iterator nit = edge2faces[EMK(i,j)].begin();
               nit != edge2faces[EMK(i,j)].end() && !found;
               nit++)
             {
@@ -494,6 +514,7 @@ inline igl::SelfIntersectMesh<
       }
     }
   }
+  //cout<<"CDT: "<<tictoc()<<"  "<<t_proj_del<<endl;
   assert(NV_count == (Index)NV.size());
   // Build output
 #ifndef NDEBUG
@@ -535,7 +556,7 @@ inline igl::SelfIntersectMesh<
   {
     Index i = 0;
     for(
-      typename list<Point_3>::const_iterator nvit = NV.begin();
+      typename Point_3List::const_iterator nvit = NV.begin();
       nvit != NV.end();
       nvit++)
     {
@@ -571,7 +592,7 @@ inline igl::SelfIntersectMesh<
   {
     Index v = V.rows();
     for(
-      typename list<Point_3>::const_iterator nvit = NV.begin();
+      typename Point_3List::const_iterator nvit = NV.begin();
       nvit != NV.end();
       nvit++)
     {
@@ -585,6 +606,7 @@ inline igl::SelfIntersectMesh<
       v++;
     }
   }
+  //cout<<"Output + dupes: "<<tictoc()<<endl;
 
   // Q: Does this give the same result as TETGEN?
   // A: For the cow and beast, yes.
@@ -639,7 +661,7 @@ inline void igl::SelfIntersectMesh<
       // Create new list if there is no entry
       if(edge2faces.count(EMK(i,j))==0)
       {
-        edge2faces[EMK(i,j)] = list<Index>();
+        edge2faces[EMK(i,j)] = EMV();
       }
       // append to list
       edge2faces[EMK(i,j)].push_back(f);
@@ -1096,7 +1118,7 @@ inline void igl::SelfIntersectMesh<
   DerivedJ,
   DerivedIM>::projected_delaunay(
   const Triangle_3 & A,
-  const std::list<CGAL::Object> & A_objects_3,
+  const ObjectList & A_objects_3,
   CDT_plus_2 & cdt)
 {
   using namespace std;
@@ -1115,20 +1137,16 @@ inline void igl::SelfIntersectMesh<
     cdt.insert_constraint( corners[(i+1)%3], corners[(i+2)%3]);
   }
   // Insert constraints for intersection objects
-  for(
-    typename list<CGAL::Object>::const_iterator lit = A_objects_3.begin();
-    lit != A_objects_3.end();
-    lit++)
+  for( const auto & obj : A_objects_3)
   {
-    CGAL::Object obj = *lit;
-    if(const Point_3 *ipoint = CGAL::object_cast<Point_3 >(&obj))
-    {
-      // Add point
-      cdt.insert(P.to_2d(*ipoint));
-    } else if(const Segment_3 *iseg = CGAL::object_cast<Segment_3 >(&obj))
+    if(const Segment_3 *iseg = CGAL::object_cast<Segment_3 >(&obj))
     {
       // Add segment constraint
       cdt.insert_constraint(P.to_2d(iseg->vertex(0)),P.to_2d(iseg->vertex(1)));
+    }else if(const Point_3 *ipoint = CGAL::object_cast<Point_3 >(&obj))
+    {
+      // Add point
+      cdt.insert(P.to_2d(*ipoint));
     } else if(const Triangle_3 *itri = CGAL::object_cast<Triangle_3 >(&obj))
     {
       // Add 3 segment constraints

include/igl/cgal/remesh_self_intersections.cpp

@@ -9,6 +9,7 @@
 #include "SelfIntersectMesh.h"
 #include <igl/C_STR.h>
 #include <list>
+#include <iostream>
 
 template <
   typename DerivedV,
@@ -82,4 +83,6 @@ IGL_INLINE void igl::remesh_self_intersections(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 template void igl::remesh_self_intersections<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::RemeshSelfIntersectionsParam const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::remesh_self_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, igl::RemeshSelfIntersectionsParam const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::remesh_self_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -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&, igl::RemeshSelfIntersectionsParam const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/colon.h

@@ -13,6 +13,12 @@ namespace igl
 {
   // Note:
   // This should be potentially replaced with eigen's LinSpaced() function
+  //
+  // If step = 1, it's about 5 times faster to use:
+  //     X = Eigen::VectorXi::LinSpaced(n,0,n-1);
+  // than 
+  //     X = igl::colon<int>(0,n-1);
+  //
 
   // Colon operator like matlab's colon operator. Enumerats values between low
   // and hi with step step.

include/igl/comiso/miq.cpp

@@ -202,12 +202,12 @@ namespace igl {
     const Eigen::PlainObjectBase<DerivedF> &F;
     const Eigen::PlainObjectBase<DerivedF> &TT;
     const Eigen::PlainObjectBase<DerivedF> &TTi;
-    const Eigen::PlainObjectBase<DerivedV> &PD1;
-    const Eigen::PlainObjectBase<DerivedV> &PD2;
+    // const Eigen::PlainObjectBase<DerivedV> &PD1;
+    // const Eigen::PlainObjectBase<DerivedV> &PD2;
 
     const Eigen::Matrix<int, Eigen::Dynamic, 3> &Handle_MMatch;
-    const Eigen::Matrix<int, Eigen::Dynamic, 1> &Handle_Singular; // bool
-    const Eigen::Matrix<int, Eigen::Dynamic, 1> &Handle_SingularDegree; // vertex;
+    // const Eigen::Matrix<int, Eigen::Dynamic, 1> &Handle_Singular; // bool
+    // const Eigen::Matrix<int, Eigen::Dynamic, 1> &Handle_SingularDegree; // vertex;
     const Eigen::Matrix<int, Eigen::Dynamic, 3> &Handle_Seams; // 3 bool
 
 
@@ -232,11 +232,11 @@ namespace igl {
                    const Eigen::PlainObjectBase<DerivedF> &_F,
                    const Eigen::PlainObjectBase<DerivedF> &_TT,
                    const Eigen::PlainObjectBase<DerivedF> &_TTi,
-                   const Eigen::PlainObjectBase<DerivedV> &_PD1,
-                   const Eigen::PlainObjectBase<DerivedV> &_PD2,
+                  //  const Eigen::PlainObjectBase<DerivedV> &_PD1,
+                  //  const Eigen::PlainObjectBase<DerivedV> &_PD2,
                    const Eigen::Matrix<int, Eigen::Dynamic, 3> &_Handle_MMatch,
-                   const Eigen::Matrix<int, Eigen::Dynamic, 1> &_Handle_Singular,
-                   const Eigen::Matrix<int, Eigen::Dynamic, 1> &_Handle_SingularDegree,
+                  //  const Eigen::Matrix<int, Eigen::Dynamic, 1> &_Handle_Singular,
+                  //  const Eigen::Matrix<int, Eigen::Dynamic, 1> &_Handle_SingularDegree,
                    const Eigen::Matrix<int, Eigen::Dynamic, 3> &_Handle_Seams
                    );
 
@@ -308,6 +308,7 @@ namespace igl {
                       bool direct_round=true,
                       int localIter=0,
                       bool _integer_rounding=true,
+                      bool _singularity_rounding=true,
                       std::vector<int> roundVertices = std::vector<int>(),
                       std::vector<std::vector<int> > hardFeatures = std::vector<std::vector<int> >());
 
@@ -543,11 +544,11 @@ namespace igl {
               const Eigen::PlainObjectBase<DerivedF> &F_,
               const Eigen::PlainObjectBase<DerivedV> &PD1_combed,
               const Eigen::PlainObjectBase<DerivedV> &PD2_combed,
-              const Eigen::PlainObjectBase<DerivedV> &BIS1_combed,
-              const Eigen::PlainObjectBase<DerivedV> &BIS2_combed,
+              // const Eigen::PlainObjectBase<DerivedV> &BIS1_combed,
+              // const Eigen::PlainObjectBase<DerivedV> &BIS2_combed,
               const Eigen::Matrix<int, Eigen::Dynamic, 3> &Handle_MMatch,
               const Eigen::Matrix<int, Eigen::Dynamic, 1> &Handle_Singular,
-              const Eigen::Matrix<int, Eigen::Dynamic, 1> &Handle_SingularDegree,
+              // const Eigen::Matrix<int, Eigen::Dynamic, 1> &Handle_SingularDegree,
               const Eigen::Matrix<int, Eigen::Dynamic, 3> &Handle_Seams,
               Eigen::PlainObjectBase<DerivedU> &UV,
               Eigen::PlainObjectBase<DerivedF> &FUV,
@@ -557,6 +558,7 @@ namespace igl {
               int iter = 5,
               int localIter = 5,
               bool DoRound = true,
+              bool SingularityRound=true,
               std::vector<int> roundVertices = std::vector<int>(),
               std::vector<std::vector<int> > hardFeatures = std::vector<std::vector<int> >());
 
@@ -613,11 +615,11 @@ IGL_INLINE igl::VertexIndexing<DerivedV, DerivedF>::VertexIndexing(const Eigen::
                                                         const Eigen::PlainObjectBase<DerivedF> &_F,
                                                         const Eigen::PlainObjectBase<DerivedF> &_TT,
                                                         const Eigen::PlainObjectBase<DerivedF> &_TTi,
-                                                        const Eigen::PlainObjectBase<DerivedV> &_PD1,
-                                                        const Eigen::PlainObjectBase<DerivedV> &_PD2,
+                                                        // const Eigen::PlainObjectBase<DerivedV> &_PD1,
+                                                        // const Eigen::PlainObjectBase<DerivedV> &_PD2,
                                                         const Eigen::Matrix<int, Eigen::Dynamic, 3> &_Handle_MMatch,
-                                                        const Eigen::Matrix<int, Eigen::Dynamic, 1> &_Handle_Singular,
-                                                        const Eigen::Matrix<int, Eigen::Dynamic, 1> &_Handle_SingularDegree,
+                                                        // const Eigen::Matrix<int, Eigen::Dynamic, 1> &_Handle_Singular,
+                                                        // const Eigen::Matrix<int, Eigen::Dynamic, 1> &_Handle_SingularDegree,
                                                         const Eigen::Matrix<int, Eigen::Dynamic, 3> &_Handle_Seams
 
                                                         ):
@@ -625,11 +627,11 @@ V(_V),
 F(_F),
 TT(_TT),
 TTi(_TTi),
-PD1(_PD1),
-PD2(_PD2),
+// PD1(_PD1),
+// PD2(_PD2),
 Handle_MMatch(_Handle_MMatch),
-Handle_Singular(_Handle_Singular),
-Handle_SingularDegree(_Handle_SingularDegree),
+// Handle_Singular(_Handle_Singular),
+// Handle_SingularDegree(_Handle_SingularDegree),
 Handle_Seams(_Handle_Seams)
 {
 
@@ -747,9 +749,9 @@ IGL_INLINE void igl::VertexIndexing<DerivedV, DerivedF>::FindInitialPos(const in
     if (!vertexB)
       //seam=curr_f->IsSeam(next_f);
       seam=IsSeam(curr_f,next_f);
-    if (vertexB)
-      assert(!Handle_Singular(vert));
-    ;
+    // if (vertexB)
+      // assert(!Handle_Singular(vert));
+    // ;
     //assert(!vert->IsSingular());
     possible_split=((on_border)||(seam));
     complete_turn = next_f == f_init;
@@ -963,6 +965,7 @@ IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::SolvePoisson(Eigen::Vect
                                                           bool direct_round,
                                                           int localIter,
                                                           bool _integer_rounding,
+                                                          bool _singularity_rounding,
                                                           std::vector<int> roundVertices,
                                                           std::vector<std::vector<int> > hardFeatures)
 {
@@ -1004,11 +1007,11 @@ IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::SolvePoisson(Eigen::Vect
   if (DEBUGPRINT)
     printf("\n BUILT THE MATRIX \n");
 
-  if (integer_rounding)
-  {
-    AddSingularityRound();
+  if (integer_rounding)    
     AddToRoundVertices(roundVertices);
-  }
+
+  if (_singularity_rounding)
+      AddSingularityRound();
 
   int t1=clock();
   if (DEBUGPRINT) printf("\n time:%d \n",t1-t0);
@@ -1767,6 +1770,7 @@ IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::MixedIntegerSolve(double
   COMISO::ConstrainedSolver solver;
 
   solver.misolver().set_local_iters(localIter);
+
   solver.misolver().set_direct_rounding(direct_round);
 
   std::sort(ids_to_round.begin(),ids_to_round.end());
@@ -1799,12 +1803,13 @@ IGL_INLINE void igl::PoissonSolver<DerivedV, DerivedF>::addSharpEdgeConstraint(i
   int v2 = F(fid,(vid+1)%3);
 
   Eigen::Matrix<typename DerivedV::Scalar, 3, 1> e = V.row(v2) - V.row(v1);
+  e = e.normalized();
 
-  int v1i = GetFirstVertexIndex(v1);
-  int v2i = GetFirstVertexIndex(v2);
+  int v1i = HandleS_Index(fid,vid);//GetFirstVertexIndex(v1);
+  int v2i = HandleS_Index(fid,(vid+1)%3);//GetFirstVertexIndex(v2);
 
-  double d1 = fabs(e.dot(PD1.row(fid)));
-  double d2 = fabs(e.dot(PD2.row(fid)));
+  double d1 = fabs(e.dot(PD1.row(fid).normalized()));
+  double d2 = fabs(e.dot(PD2.row(fid).normalized()));
 
   int offset = 0;
 
@@ -1831,11 +1836,11 @@ IGL_INLINE igl::MIQ_class<DerivedV, DerivedF, DerivedU>::MIQ_class(const Eigen::
                                                         const Eigen::PlainObjectBase<DerivedF> &F_,
                                                         const Eigen::PlainObjectBase<DerivedV> &PD1_combed,
                                                         const Eigen::PlainObjectBase<DerivedV> &PD2_combed,
-                                                        const Eigen::PlainObjectBase<DerivedV> &BIS1_combed,
-                                                        const Eigen::PlainObjectBase<DerivedV> &BIS2_combed,
+                                                        // const Eigen::PlainObjectBase<DerivedV> &BIS1_combed,
+                                                        // const Eigen::PlainObjectBase<DerivedV> &BIS2_combed,
                                                         const Eigen::Matrix<int, Eigen::Dynamic, 3> &Handle_MMatch,
                                                         const Eigen::Matrix<int, Eigen::Dynamic, 1> &Handle_Singular,
-                                                        const Eigen::Matrix<int, Eigen::Dynamic, 1> &Handle_SingularDegree,
+                                                        // const Eigen::Matrix<int, Eigen::Dynamic, 1> &Handle_SingularDegree,
                                                         const Eigen::Matrix<int, Eigen::Dynamic, 3> &Handle_Seams,
                                                         Eigen::PlainObjectBase<DerivedU> &UV,
                                                         Eigen::PlainObjectBase<DerivedF> &FUV,
@@ -1845,6 +1850,7 @@ IGL_INLINE igl::MIQ_class<DerivedV, DerivedF, DerivedU>::MIQ_class(const Eigen::
                                                         int iter,
                                                         int localIter,
                                                         bool DoRound,
+                                                        bool SingularityRound,
                                                         std::vector<int> roundVertices,
                                                         std::vector<std::vector<int> > hardFeatures):
 V(V_),
@@ -1854,7 +1860,7 @@ F(F_)
   igl::triangle_triangle_adjacency(V,F,TT,TTi);
 
   // Prepare indexing for the linear system
-  VertexIndexing<DerivedV, DerivedF> VInd(V, F, TT, TTi, BIS1_combed, BIS2_combed, Handle_MMatch, Handle_Singular, Handle_SingularDegree, Handle_Seams);
+  VertexIndexing<DerivedV, DerivedF> VInd(V, F, TT, TTi, /*BIS1_combed, BIS2_combed,*/ Handle_MMatch, /*Handle_Singular, Handle_SingularDegree,*/ Handle_Seams);
 
   VInd.InitMapping();
   VInd.InitFaceIntegerVal();
@@ -1862,8 +1868,8 @@ F(F_)
 
   Eigen::PlainObjectBase<DerivedV> PD1_combed_for_poisson, PD2_combed_for_poisson;
   // Rotate by 90 degrees CCW
-  PD1_combed_for_poisson.setZero(BIS1_combed.rows(),3);
-  PD2_combed_for_poisson.setZero(BIS2_combed.rows(),3);
+  PD1_combed_for_poisson.setZero(PD1_combed.rows(),3);
+  PD2_combed_for_poisson.setZero(PD2_combed.rows(),3);
   for (unsigned i=0; i<PD1_combed.rows();++i)
   {
     double n1 = PD1_combed.row(i).norm();
@@ -1892,7 +1898,7 @@ F(F_)
                                             PD1_combed_for_poisson,
                                             PD2_combed_for_poisson,
                                             VInd.HandleS_Index,
-                                            VInd.Handle_Singular,
+                                            /*VInd.Handle_Singular*/Handle_Singular,
                                             VInd.Handle_SystemInfo);
   Handle_Stiffness = Eigen::VectorXd::Constant(F.rows(),1);
 
@@ -1901,7 +1907,7 @@ F(F_)
   {
     for (int i=0;i<iter;i++)
     {
-      PSolver.SolvePoisson(Handle_Stiffness, GradientSize,1.f,DirectRound,localIter,DoRound,roundVertices,hardFeatures);
+      PSolver.SolvePoisson(Handle_Stiffness, GradientSize,1.f,DirectRound,localIter,DoRound,SingularityRound,roundVertices,hardFeatures);
       int nflips=NumFlips(PSolver.WUV);
       bool folded = updateStiffeningJacobianDistorsion(GradientSize,PSolver.WUV);
       printf("ITERATION %d FLIPS %d \n",i,nflips);
@@ -1910,7 +1916,7 @@ F(F_)
   }
   else
   {
-    PSolver.SolvePoisson(Handle_Stiffness,GradientSize,1.f,DirectRound,localIter,DoRound,roundVertices,hardFeatures);
+    PSolver.SolvePoisson(Handle_Stiffness,GradientSize,1.f,DirectRound,localIter,DoRound,SingularityRound,roundVertices,hardFeatures);
   }
 
   int nflips=NumFlips(PSolver.WUV);
@@ -2170,11 +2176,11 @@ IGL_INLINE void igl::miq(const Eigen::PlainObjectBase<DerivedV> &V,
                          const Eigen::PlainObjectBase<DerivedF> &F,
                          const Eigen::PlainObjectBase<DerivedV> &PD1_combed,
                          const Eigen::PlainObjectBase<DerivedV> &PD2_combed,
-                         const Eigen::PlainObjectBase<DerivedV> &BIS1_combed,
-                         const Eigen::PlainObjectBase<DerivedV> &BIS2_combed,
+                        //  const Eigen::PlainObjectBase<DerivedV> &BIS1_combed,
+                        //  const Eigen::PlainObjectBase<DerivedV> &BIS2_combed,
                          const Eigen::Matrix<int, Eigen::Dynamic, 3> &Handle_MMatch,
                          const Eigen::Matrix<int, Eigen::Dynamic, 1> &Handle_Singular,
-                         const Eigen::Matrix<int, Eigen::Dynamic, 1> &Handle_SingularDegree,
+                        //  const Eigen::Matrix<int, Eigen::Dynamic, 1> &Handle_SingularDegree,
                          const Eigen::Matrix<int, Eigen::Dynamic, 3> &Handle_Seams,
                          Eigen::PlainObjectBase<DerivedU> &UV,
                          Eigen::PlainObjectBase<DerivedF> &FUV,
@@ -2184,6 +2190,7 @@ IGL_INLINE void igl::miq(const Eigen::PlainObjectBase<DerivedV> &V,
                          int iter,
                          int localIter,
                          bool DoRound,
+                         bool SingularityRound,
                          std::vector<int> roundVertices,
                          std::vector<std::vector<int> > hardFeatures)
 {
@@ -2193,11 +2200,11 @@ IGL_INLINE void igl::miq(const Eigen::PlainObjectBase<DerivedV> &V,
                                                    F,
                                                    PD1_combed,
                                                    PD2_combed,
-                                                   BIS1_combed,
-                                                   BIS2_combed,
+                                                  //  BIS1_combed,
+                                                  //  BIS2_combed,
                                                    Handle_MMatch,
                                                    Handle_Singular,
-                                                   Handle_SingularDegree,
+                                                  //  Handle_SingularDegree,
                                                    Handle_Seams,
                                                    UV,
                                                    FUV,
@@ -2207,6 +2214,7 @@ IGL_INLINE void igl::miq(const Eigen::PlainObjectBase<DerivedV> &V,
                                                    iter,
                                                    localIter,
                                                    DoRound,
+                                                   SingularityRound,
                                                    roundVertices,
                                                    hardFeatures);
 
@@ -2226,6 +2234,7 @@ IGL_INLINE void igl::miq(const Eigen::PlainObjectBase<DerivedV> &V,
                          int iter,
                          int localIter,
                          bool DoRound,
+                         bool SingularityRound,
                          std::vector<int> roundVertices,
                          std::vector<std::vector<int> > hardFeatures)
 {
@@ -2250,7 +2259,7 @@ IGL_INLINE void igl::miq(const Eigen::PlainObjectBase<DerivedV> &V,
   igl::find_cross_field_singularities(V, F, Handle_MMatch, isSingularity, singularityIndex);
 
   Eigen::Matrix<int, Eigen::Dynamic, 3> Handle_Seams;
-  igl::cut_mesh_from_singularities(V, F, Handle_MMatch, isSingularity, singularityIndex, Handle_Seams);
+  igl::cut_mesh_from_singularities(V, F, Handle_MMatch, Handle_Seams);
 
   Eigen::PlainObjectBase<DerivedV> PD1_combed, PD2_combed;
   igl::comb_frame_field(V, F, PD1, PD2, BIS1_combed, BIS2_combed, PD1_combed, PD2_combed);
@@ -2259,11 +2268,11 @@ IGL_INLINE void igl::miq(const Eigen::PlainObjectBase<DerivedV> &V,
            F,
            PD1_combed,
            PD2_combed,
-           BIS1_combed,
-           BIS2_combed,
+          //  BIS1_combed,
+          //  BIS2_combed,
            Handle_MMatch,
            isSingularity,
-           singularityIndex,
+          //  singularityIndex,
            Handle_Seams,
            UV,
            FUV,
@@ -2273,6 +2282,7 @@ IGL_INLINE void igl::miq(const Eigen::PlainObjectBase<DerivedV> &V,
            iter,
            localIter,
            DoRound,
+           SingularityRound,
            roundVertices,
            hardFeatures);
 
@@ -2280,5 +2290,7 @@ IGL_INLINE void igl::miq(const Eigen::PlainObjectBase<DerivedV> &V,
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template void igl::miq<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 2, 0, -1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, double, double, bool, int, int, bool, std::__1::vector<int, std::__1::allocator<int> >, std::__1::vector<std::__1::vector<int, std::__1::allocator<int> >, std::__1::allocator<std::__1::vector<int, std::__1::allocator<int> > > >);
+// template void igl::miq<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 2, 0, -1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, double, double, bool, int, int, bool, std::__1::vector<int, std::__1::allocator<int> >, std::__1::vector<std::__1::vector<int, std::__1::allocator<int> >, std::__1::allocator<std::__1::vector<int, std::__1::allocator<int> > > >);
+// template void igl::miq<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 2, 0, -1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::Matrix<int, -1, 3, 0, -1, 3> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, double, double, bool, int, int, bool, std::__1::vector<int, std::__1::allocator<int> >, std::__1::vector<std::__1::vector<int, std::__1::allocator<int> >, std::__1::allocator<std::__1::vector<int, std::__1::allocator<int> > > >);
+//template void igl::miq<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 2, 0, -1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::Matrix<int, -1, 3, 0, -1, 3> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::Matrix<int, -1, 3, 0, -1, 3> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, double, double, bool, int, int, bool, std::__1::vector<int, std::__1::allocator<int> >, std::__1::vector<std::__1::vector<int, std::__1::allocator<int> >, std::__1::allocator<std::__1::vector<int, std::__1::allocator<int> > > >);
 #endif

include/igl/comiso/miq.h

@@ -54,7 +54,7 @@ namespace igl
                                               bool direct_round = false,
                                               int iter = 5,
                                               int local_iter = 5,
-                                              bool DoRound = true,
+                                              bool DoRound = true,bool SingularityRound=true,
                                               std::vector<int> round_vertices = std::vector<int>(),
                                               std::vector<std::vector<int> > hard_features = std::vector<std::vector<int> >());
 
@@ -72,11 +72,11 @@ namespace igl
                                               const Eigen::PlainObjectBase<DerivedF> &F,
                                               const Eigen::PlainObjectBase<DerivedV> &PD1_combed,
                                               const Eigen::PlainObjectBase<DerivedV> &PD2_combed,
-                                              const Eigen::PlainObjectBase<DerivedV> &BIS1_combed,
-                                              const Eigen::PlainObjectBase<DerivedV> &BIS2_combed,
+                                              // const Eigen::PlainObjectBase<DerivedV> &BIS1_combed,
+                                              // const Eigen::PlainObjectBase<DerivedV> &BIS2_combed,
                                               const Eigen::Matrix<int, Eigen::Dynamic, 3> &MMatch,
                                               const Eigen::Matrix<int, Eigen::Dynamic, 1> &Singular,
-                                              const Eigen::Matrix<int, Eigen::Dynamic, 1> &SingularDegree,
+                                              // const Eigen::Matrix<int, Eigen::Dynamic, 1> &SingularDegree,
                                               const Eigen::Matrix<int, Eigen::Dynamic, 3> &Seams,
                                               Eigen::PlainObjectBase<DerivedU> &UV,
                                               Eigen::PlainObjectBase<DerivedF> &FUV,
@@ -84,7 +84,7 @@ namespace igl
                                               double Stiffness = 5.0,
                                               bool DirectRound = false,
                                               int iter = 5,
-                                              int localIter = 5, bool DoRound = true,
+                                              int localIter = 5, bool DoRound = true,bool SingularityRound=true,
                                               std::vector<int> roundVertices = std::vector<int>(),
                                               std::vector<std::vector<int> > hardFeatures = std::vector<std::vector<int> >());
 };

include/igl/cross_field_missmatch.cpp

@@ -19,7 +19,7 @@
 #include <igl/rotation_matrix_from_directions.h>
 
 namespace igl {
-  template <typename DerivedV, typename DerivedF, typename DerivedO>
+  template <typename DerivedV, typename DerivedF>
   class MissMatchCalculator
   {
   public:
@@ -91,14 +91,14 @@ public:
     igl::triangle_triangle_adjacency(V,F,TT,TTi);
   }
 
-  inline void calculateMissmatch(Eigen::PlainObjectBase<DerivedO> &Handle_MMatch)
+  inline void calculateMissmatch(Eigen::PlainObjectBase<DerivedF> &Handle_MMatch)
   {
     Handle_MMatch.setConstant(F.rows(),3,-1);
-    for (unsigned int i=0;i<F.rows();i++)
+    for (size_t i=0;i<F.rows();i++)
     {
       for (int j=0;j<3;j++)
       {
-        if (i==TT(i,j) || TT(i,j) == -1)
+        if (((int)i)==TT(i,j) || TT(i,j) == -1)
           Handle_MMatch(i,j)=0;
         else
           Handle_MMatch(i,j) = MissMatchByCross(i,TT(i,j));
@@ -108,13 +108,13 @@ public:
 
 };
 }
-template <typename DerivedV, typename DerivedF, typename DerivedO>
+template <typename DerivedV, typename DerivedF>
 IGL_INLINE void igl::cross_field_missmatch(const Eigen::PlainObjectBase<DerivedV> &V,
                                            const Eigen::PlainObjectBase<DerivedF> &F,
                                            const Eigen::PlainObjectBase<DerivedV> &PD1,
                                            const Eigen::PlainObjectBase<DerivedV> &PD2,
                                            const bool isCombed,
-                                           Eigen::PlainObjectBase<DerivedO> &missmatch)
+                                           Eigen::PlainObjectBase<DerivedF> &missmatch)
 {
   Eigen::PlainObjectBase<DerivedV> PD1_combed;
   Eigen::PlainObjectBase<DerivedV> PD2_combed;
@@ -126,11 +126,12 @@ IGL_INLINE void igl::cross_field_missmatch(const Eigen::PlainObjectBase<DerivedV
     PD1_combed = PD1;
     PD2_combed = PD2;
   }
-  igl::MissMatchCalculator<DerivedV, DerivedF, DerivedO> sf(V, F, PD1_combed, PD2_combed);
+  igl::MissMatchCalculator<DerivedV, DerivedF> sf(V, F, PD1_combed, PD2_combed);
   sf.calculateMissmatch(missmatch);
 }
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template void igl::cross_field_missmatch<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
+template void igl::cross_field_missmatch<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
+
 #endif

include/igl/cross_field_missmatch.h

@@ -15,7 +15,7 @@ namespace igl
   // Calculates the missmatch (integer), at each face edge, of a cross field defined on the mesh faces.
   // The integer missmatch is a multiple of pi/2 that transforms the cross on one side of the edge to
   // the cross on the other side. It represents the deviation from a Lie connection across the edge.
-  
+
   // Inputs:
   //   V         #V by 3 eigen Matrix of mesh vertex 3D positions
   //   F         #F by 3 eigen Matrix of face (quad) indices
@@ -27,14 +27,14 @@ namespace igl
   //   Handle_MMatch    #F by 3 eigen Matrix containing the integer missmatch of the cross field
   //                    across all face edges
   //
-  
-  template <typename DerivedV, typename DerivedF, typename DerivedO>
+
+  template <typename DerivedV, typename DerivedF>
   IGL_INLINE void cross_field_missmatch(const Eigen::PlainObjectBase<DerivedV> &V,
                                         const Eigen::PlainObjectBase<DerivedF> &F,
                                         const Eigen::PlainObjectBase<DerivedV> &PD1,
                                         const Eigen::PlainObjectBase<DerivedV> &PD2,
                                         const bool isCombed,
-                                        Eigen::PlainObjectBase<DerivedO> &missmatch);
+                                        Eigen::PlainObjectBase<DerivedF> &missmatch);
 }
 #ifndef IGL_STATIC_LIBRARY
 #include "cross_field_missmatch.cpp"

include/igl/draw_mesh.cpp

@@ -90,7 +90,7 @@ IGL_INLINE void igl::draw_mesh(
   {
     assert(F.maxCoeff() < V.rows());
     assert(V.cols() == 3);
-    assert(rC == rV || rC == rF || rC == rF*3 || C.size() == 0);
+    assert(rC == rV || rC == rF || rC == rF*3 || rC==1 || C.size() == 0);
     assert(C.cols() == 3 || C.size() == 0);
     assert(N.cols() == 3 || N.size() == 0);
     assert(TC.cols() == 2 || TC.size() == 0);

include/igl/exterior_edges.cpp

@@ -42,11 +42,11 @@ IGL_INLINE void igl::exterior_edges(
   all_edges(F,all_E);
   long int n = F.maxCoeff()+1;
   int m = F.minCoeff();
-  const auto & compress = [&n,&m](const int i, const int j)->long int
+  const auto & compress = [n,m](const int i, const int j)->long int
   {
     return n*(i-m)+(j-m);
   };
-  const auto & decompress = [&n,&m](const long int l,int & i, int & j)
+  const auto & decompress = [n,m](const long int l,int & i, int & j)
   {
     i = (l / n) + m;
     j = (l % n) + m;

include/igl/facet_components.cpp

@@ -0,0 +1,71 @@
+#include "facet_components.h"
+#include <igl/triangle_triangle_adjacency.h>
+#include <vector>
+#include <queue>
+template <typename DerivedF, typename DerivedC>
+IGL_INLINE void igl::facet_components(
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedC> & C)
+{
+  using namespace std;
+  using namespace igl;
+  typedef typename DerivedF::Index Index;
+  vector<vector<vector<Index > > > TT;
+  vector<vector<vector<Index > > > TTi;
+  triangle_triangle_adjacency(F,TT,TTi);
+  return facet_components(TT,C);
+}
+
+template <
+  typename TTIndex, 
+  typename DerivedC>
+IGL_INLINE void igl::facet_components(
+  const std::vector<std::vector<std::vector<TTIndex > > > & TT,
+  Eigen::PlainObjectBase<DerivedC> & C)
+{
+  using namespace std;
+  using namespace igl;
+  typedef TTIndex Index;
+  const Index m = TT.size();
+  C.resize(m,1);
+  vector<bool> seen(m,false);
+  Index id = 0;
+  for(Index g = 0;g<m;g++)
+  {
+    if(seen[g])
+    {
+      continue;
+    }
+    queue<Index> Q;
+    Q.push(g);
+    while(!Q.empty())
+    {
+      const Index f = Q.front();
+      Q.pop();
+      if(seen[f])
+      {
+        continue;
+      }
+      seen[f] = true;
+      C(f,0) = id;
+      // Face f's neighbor lists opposite opposite each corner
+      for(const auto & c : TT[f])
+      {
+        // Each neighbor
+        for(const auto & n : c)
+        {
+          if(!seen[n])
+          {
+            Q.push(n);
+          }
+        }
+      }
+    }
+    id++;
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+template void igl::facet_components<long, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(std::vector<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >, std::allocator<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/facet_components.h

@@ -0,0 +1,33 @@
+#ifndef FACET_COMPONENTS_H
+#define FACET_COMPONENTS_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+namespace igl
+{
+  // Compute connected components of facets based on edge-edge adjacency.
+  //
+  // Inputs:
+  //   F  #F by 3 list of triangle indices
+  // Ouputs:
+  //   C  #F list of connected component ids
+  template <typename DerivedF, typename DerivedC>
+  IGL_INLINE void facet_components(
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    Eigen::PlainObjectBase<DerivedC> & C);
+  // Inputs:
+  //   TT  #TT by 3 list of list of adjacency triangles (see
+  //   triangle_triangle_adjacency.h)
+  // Ouputs:
+  //   C  #F list of connected component ids
+  template <
+    typename TTIndex, 
+    typename DerivedC>
+  IGL_INLINE void facet_components(
+    const std::vector<std::vector<std::vector<TTIndex > > > & TT,
+    Eigen::PlainObjectBase<DerivedC> & C);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "facet_components.cpp"
+#endif
+#endif

include/igl/find_cross_field_singularities.cpp

@@ -65,15 +65,22 @@ IGL_INLINE void igl::find_cross_field_singularities(const Eigen::PlainObjectBase
                                                     const Eigen::PlainObjectBase<DerivedV> &PD1,
                                                     const Eigen::PlainObjectBase<DerivedV> &PD2,
                                                     Eigen::PlainObjectBase<DerivedO> &isSingularity,
-                                                    Eigen::PlainObjectBase<DerivedO> &singularityIndex)
+                                                    Eigen::PlainObjectBase<DerivedO> &singularityIndex,
+                                                    bool isCombed)
 {
-
-  Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 3> Handle_MMatch;
-  igl::cross_field_missmatch<DerivedV, DerivedF, DerivedO>(V, F, PD1, PD2, Handle_MMatch);
+  // Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > Handle_MMatch;
+  // Eigen::PlainObjectBase<Eigen::Matrix<int,  Eigen::Dynamic, 3> > Handle_MMatch;
+  // Eigen::Matrix<int, Eigen::Dynamic, 3> Handle_MMatch;
+  Eigen::PlainObjectBase<DerivedF> Handle_MMatch;
+  // Eigen::Matrix<typename DerivedF::Scalar, Eigen::Dynamic, 3> Handle_MMatch;
+  // Eigen::Matrix<typename DerivedF::Scalar, Eigen::Dynamic, Eigen::Dynamic> Handle_MMatch;
+  igl::cross_field_missmatch<DerivedV, DerivedF>(V, F, PD1, PD2, isCombed, Handle_MMatch);
   igl::find_cross_field_singularities(V, F, Handle_MMatch, isSingularity, singularityIndex);
 }
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 template void igl::find_cross_field_singularities<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::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<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::find_cross_field_singularities<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -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, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&,
+Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, bool);
 #endif

include/igl/find_cross_field_singularities.h

@@ -13,8 +13,8 @@
 namespace igl
 {
   // Computes singularities of a cross field, assumed combed
-  
-  
+
+
   // Inputs:
   //   V                #V by 3 eigen Matrix of mesh vertex 3D positions
   //   F                #F by 3 eigen Matrix of face (quad) indices
@@ -47,7 +47,8 @@ namespace igl
                                                  const Eigen::PlainObjectBase<DerivedV> &PD1,
                                                  const Eigen::PlainObjectBase<DerivedV> &PD2,
                                                  Eigen::PlainObjectBase<DerivedO> &isSingularity,
-                                                 Eigen::PlainObjectBase<DerivedO> &singularityIndex);
+                                                 Eigen::PlainObjectBase<DerivedO> &singularityIndex,
+                                                 bool isCombed = false);
 }
 #ifndef IGL_STATIC_LIBRARY
 #include "find_cross_field_singularities.cpp"

include/igl/get_modifiers.cpp

@@ -36,7 +36,7 @@ IGL_INLINE int igl::get_modifiers()
   mod |= (carbon_is_keydown(kVK_Option)?GLUT_ACTIVE_ALT:0);
   mod |= (carbon_is_keydown(kVK_Control)?GLUT_ACTIVE_CTRL:0);
 #else
-#  error "Not supported.
+#  error "Not supported."
 #endif
   return mod;
 }

include/igl/get_seconds.h

@@ -12,6 +12,21 @@
 namespace igl
 {
   // Return the current time in seconds since program start
+  // 
+  // Example:
+  //    const auto & tictoc = []()
+  //    {
+  //      static double t_start = igl::get_seconds();
+  //      double diff = igl::get_seconds()-t_start;
+  //      t_start += diff;
+  //      return diff;
+  //    };
+  //    tictoc();
+  //    ... // part 1
+  //    cout<<"part 1: "<<tictoc()<<endl;
+  //    ... // part 2
+  //    cout<<"part 2: "<<tictoc()<<endl;
+  //    ... // etc
   IGL_INLINE double get_seconds();
 
 }

include/igl/grad.cpp

@@ -105,4 +105,5 @@ IGL_INLINE void igl::grad(const Eigen::PlainObjectBase<DerivedV>&V,
 // Explicit template specialization
 // template void igl::grad<double, int>(Eigen::Matrix<double, -1, -1, 0, -1,-1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&,Eigen::SparseMatrix<double, 0, int>&);
 template void igl::grad<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::SparseMatrix<Eigen::Matrix<double, -1, 3, 0, -1, 3>::Scalar, 0, int>&);
+//template void igl::grad<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::SparseMatrix<Eigen::Matrix<double, -1, 3, 0, -1, 3>::Scalar, 0, int>&);
 #endif

include/igl/matlab/mexStream.cpp

@@ -1,27 +0,0 @@
-// 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 "MexStream.h"
-#include "mex.h"
-
-inline std::streamsize igl::MexStream::xsputn(
-  const char *s, 
-  std::streamsize n) 
-{
-  mexPrintf("%.*s",n,s);
-  mexEvalString("drawnow;"); // to dump string.
-  return n;
-}
-
-inline int igl::MexStream::overflow(int c) 
-{
-    if (c != EOF) {
-      mexPrintf("%.1s",&c);
-      mexEvalString("drawnow;"); // to dump string.
-    }
-    return 1;
-}

include/igl/matlab/mexStream.h

@@ -30,6 +30,24 @@ namespace igl
       inline virtual int overflow(int c = EOF);
   }; 
 }
-// ALWAYS INCLUDE
-#include "MexStream.cpp"
+
+// Implementation 
+#include <mex.h>
+inline std::streamsize igl::MexStream::xsputn(
+  const char *s, 
+  std::streamsize n) 
+{
+  mexPrintf("%.*s",n,s);
+  mexEvalString("drawnow;"); // to dump string.
+  return n;
+}
+
+inline int igl::MexStream::overflow(int c) 
+{
+    if (c != EOF) {
+      mexPrintf("%.1s",&c);
+      mexEvalString("drawnow;"); // to dump string.
+    }
+    return 1;
+}
 #endif

include/igl/matlab/prepare_lhs.cpp

@@ -11,6 +11,17 @@ IGL_INLINE void igl::prepare_lhs_double(
   copy(&V.data()[0],&V.data()[0]+V.size(),Vp);
 }
 
+template <typename DerivedV>
+IGL_INLINE void igl::prepare_lhs_logical(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  mxArray *plhs[])
+{
+  using namespace std;
+  plhs[0] = mxCreateLogicalMatrix(V.rows(),V.cols());
+  mxLogical * Vp = static_cast<mxLogical*>(mxGetData(plhs[0]));
+  copy(&V.data()[0],&V.data()[0]+V.size(),Vp);
+}
+
 template <typename DerivedV>
 IGL_INLINE void igl::prepare_lhs_index(
   const Eigen::PlainObjectBase<DerivedV> & V,
@@ -26,4 +37,5 @@ IGL_INLINE void igl::prepare_lhs_index(
 template void igl::prepare_lhs_index<Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, mxArray_tag**);
 template void igl::prepare_lhs_index<Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, mxArray_tag**);
 template void igl::prepare_lhs_double<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, mxArray_tag**);
+template void igl::prepare_lhs_index<Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, mxArray_tag**);
 #endif

include/igl/matlab/prepare_lhs.h

@@ -15,6 +15,11 @@ namespace igl
   IGL_INLINE void prepare_lhs_double(
     const Eigen::PlainObjectBase<DerivedV> & V,
     mxArray *plhs[]);
+  // Casts to logical
+  template <typename DerivedV>
+  IGL_INLINE void prepare_lhs_logical(
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    mxArray *plhs[]);
   // Writes out a matrix and adds 1
   template <typename DerivedV>
   IGL_INLINE void prepare_lhs_index(

include/igl/n_polyvector.cpp

@@ -12,6 +12,7 @@
 #include <igl/nchoosek.h>
 #include <igl/slice.h>
 #include <igl/polyroots.h>
+#include <igl/igl_inline.h>
 #include <Eigen/Sparse>
 
 #include <iostream>
@@ -39,7 +40,7 @@ namespace igl {
     Eigen::PlainObjectBase<DerivedV> B1, B2, FN;
 
     IGL_INLINE void computek();
-    IGL_INLINE void setFieldFromGeneralCoefficients(const  std::vector<Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1>> &coeffs,
+    IGL_INLINE void setFieldFromGeneralCoefficients(const  std::vector<Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1> > &coeffs,
                                                     std::vector<Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 2> > &pv);
     IGL_INLINE void computeCoefficientLaplacian(int n, Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &D);
     IGL_INLINE void getGeneralCoeffConstraints(const Eigen::VectorXi &isConstrained,
@@ -152,7 +153,7 @@ minQuadWithKnownMini(const Eigen::SparseMatrix<std::complex<typename DerivedV::S
       indu++;
     }
 
-  Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar>> Quu, Quk;
+  Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > Quu, Quk;
 
   igl::slice(Q,unknown, unknown, Quu);
   igl::slice(Q,unknown, known, Quk);
@@ -166,14 +167,14 @@ minQuadWithKnownMini(const Eigen::SparseMatrix<std::complex<typename DerivedV::S
 
   Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > rhs = (Quk*xknown).sparseView()+.5*fu;
 
-  Eigen::SparseLU< Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar>>> solver;
+  Eigen::SparseLU< Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > > solver;
   solver.compute(-Quu);
   if(solver.info()!=Eigen::Success)
   {
     std::cerr<<"Decomposition failed!"<<std::endl;
     return;
   }
-  Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar>>  b  = solver.solve(rhs);
+  Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> >  b  = solver.solve(rhs);
   if(solver.info()!=Eigen::Success)
   {
     std::cerr<<"Solving failed!"<<std::endl;
@@ -207,7 +208,7 @@ IGL_INLINE bool igl::PolyVectorFieldFinder<DerivedV, DerivedF>::
   // ... +
   // (-1)^n c[n-1]
 
-  std::vector<Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1>> coeffs(n,Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1>::Zero(numF, 1));
+  std::vector<Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1> > coeffs(n,Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1>::Zero(numF, 1));
 
   for (int i =0; i<n; ++i)
   {
@@ -241,8 +242,8 @@ IGL_INLINE bool igl::PolyVectorFieldFinder<DerivedV, DerivedF>::
 }
 
 template<typename DerivedV, typename DerivedF>
-IGL_INLINE void igl::PolyVectorFieldFinder<DerivedV, DerivedF>::setFieldFromGeneralCoefficients(const  std::vector<Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1>> &coeffs,
-                                                            std::vector<Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 2>> &pv)
+IGL_INLINE void igl::PolyVectorFieldFinder<DerivedV, DerivedF>::setFieldFromGeneralCoefficients(const  std::vector<Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1> > &coeffs,
+                                                            std::vector<Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 2> > &pv)
 {
   pv.assign(n, Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 2>::Zero(numF, 2));
   for (int i = 0; i <numF; ++i)
@@ -302,7 +303,7 @@ IGL_INLINE void igl::PolyVectorFieldFinder<DerivedV, DerivedF>::setFieldFromGene
 template<typename DerivedV, typename DerivedF>
 IGL_INLINE void igl::PolyVectorFieldFinder<DerivedV, DerivedF>::computeCoefficientLaplacian(int n, Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &D)
 {
-  std::vector<Eigen::Triplet<std::complex<typename DerivedV::Scalar> >> tripletList;
+  std::vector<Eigen::Triplet<std::complex<typename DerivedV::Scalar> > > tripletList;
 
   // For every non-border edge
   for (unsigned eid=0; eid<numE; ++eid)
@@ -343,7 +344,7 @@ IGL_INLINE void igl::PolyVectorFieldFinder<DerivedV, DerivedF>::getGeneralCoeffC
   Ck.resize(numConstrained,1);
   int n = cfW.cols()/3;
 
-  std::vector<std::vector<int>> allCombs;
+  std::vector<std::vector<int> > allCombs;
   igl::nchoosek(0,k+1,n,allCombs);
 
   int ind = 0;

include/igl/nchoosek.cpp

@@ -14,7 +14,7 @@ namespace igl {
   private:
     std::vector<int> combinations;
     void add(const std::vector<int>& v,
-             std::vector<std::vector<int>> &allCombs)
+             std::vector<std::vector<int> > &allCombs)
     {
       allCombs.push_back(v);
     }
@@ -23,7 +23,7 @@ namespace igl {
     void doCombs(int offset,
                  int k,
                  int N,
-                 std::vector<std::vector<int>> &allCombs)
+                 std::vector<std::vector<int> > &allCombs)
     {
       if (k == 0) {
         add(combinations,allCombs);
@@ -44,7 +44,7 @@ namespace igl {
 IGL_INLINE void igl::nchoosek(int offset,
                               int k,
                               int N,
-                              std::vector<std::vector<int>> &allCombs)
+                              std::vector<std::vector<int> > &allCombs)
 {
   CombinationFinder cmbf;
   allCombs.clear();

include/igl/nchoosek.h

@@ -14,19 +14,10 @@
 #include <Eigen/Core>
 
 namespace igl {
-  //todo
-  /// Given 2 vectors centered on origin calculate the rotation matrix from first to the second
-
-  // Inputs:
-  //   v0, v1         the two #3 by 1 vectors
-  //   normalized     boolean, if false, then the vectors are normalized prior to the calculation
-  // Output:
-  //                  3 by 3 rotation matrix that takes v0 to v1
-  //
   IGL_INLINE void nchoosek(int offset,
                            int k,
                            int N,
-                           std::vector<std::vector<int>> &allCombs);
+                           std::vector<std::vector<int> > &allCombs);
 }
 
 

include/igl/outer_facet.cpp

@@ -0,0 +1,71 @@
+#include "outer_facet.h"
+#include "sort.h"
+#include "vertex_triangle_adjacency.h"
+
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedN,
+  typename DerivedI,
+  typename f_type>
+IGL_INLINE void igl::outer_facet(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  const Eigen::PlainObjectBase<DerivedN> & N,
+  const Eigen::PlainObjectBase<DerivedI> & I,
+  f_type & max_f,
+  bool & flip)
+{
+  using namespace std;
+  typedef typename DerivedV::Scalar Scalar;
+  typedef typename DerivedV::Index Index;
+  typedef 
+    typename Eigen::Matrix<Index, DerivedV::RowsAtCompileTime,1> VectorXI;
+  typedef 
+    typename Eigen::Matrix<Scalar, DerivedV::RowsAtCompileTime,1> VectorXS;
+  // "Direct repair of self-intersecting meshes" [Attene 14]
+  const Index mi = I.size();
+  Scalar max_x = -1e26;
+  Index max_v = -1;
+  Scalar max_nx = -1e26;
+  for(Index i = 0;i<mi;i++)
+  {
+    const Index f = I(i);
+    const Scalar nx = N(f,0);
+    if(fabs(nx)>0)
+    {
+      for(Index c = 0;c<3;c++)
+      {
+        const Index v = F(f,c);
+        if(v == max_v)
+        {
+          if(fabs(nx) > max_nx)
+          {
+            // Just update max face and normal
+            max_f = f;
+            max_nx = fabs(nx);
+            flip = nx<0;
+          }
+        }else
+        {
+          const Scalar x = V(v);
+          if(x>max_x)
+          {
+            // update max vertex, face and normal
+            max_v = v;
+            max_x = x;
+            max_f = f;
+            max_nx = fabs(nx);
+            flip = nx<0;
+          }
+        }
+      }
+    }
+  }
+  assert(max_v >=0 && "Very degenerate case, no suitable face found.");
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+template void igl::outer_facet<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, int>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> > const&, int&, bool&);
+#endif

include/igl/outer_facet.h

@@ -0,0 +1,39 @@
+#ifndef IGL_OUTER_FACET_H
+#define IGL_OUTER_FACET_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Compute a single facet which is guaranteed to be part of the "outer hull of
+  // a mesh (V,F). This implementation follows Section 3.6 of "Direct repair of
+  // self-intersecting meshes" [Attene 2014].
+  //
+  // Inputs:
+  //   V  #V by 3 list of vertex positions
+  //   F  #F by 3 list of triangle indices into V
+  //   N  #F by 3 list of face normals
+  //   I  #I list of facets to actually consider
+  // Outputs:
+  //   f  index of facet into V
+  //   flip  whether facet's orientation should be flipped so that
+  //     counter-clockwise normal points outward.
+  //
+  // See also: outer_hull.h
+  template <
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedN,
+    typename DerivedI,
+    typename f_type>
+  IGL_INLINE void outer_facet(
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    const Eigen::PlainObjectBase<DerivedN> & N,
+    const Eigen::PlainObjectBase<DerivedI> & I,
+    f_type & f,
+    bool & flip);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "outer_facet.cpp"
+#endif
+#endif

include/igl/outer_hull.cpp

@@ -0,0 +1,229 @@
+#include "outer_hull.h"
+#include "outer_facet.h"
+#include "facet_components.h"
+
+#include "triangle_triangle_adjacency.h"
+#include "unique_edge_map.h"
+#include "per_face_normals.h"
+#include "all_edges.h"
+#include "colon.h"
+#include "get_seconds.h"
+
+#include <Eigen/Geometry>
+#include <vector>
+#include <map>
+#include <queue>
+#include <iostream>
+
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedG,
+  typename DerivedJ,
+  typename Derivedflip>
+IGL_INLINE void igl::outer_hull(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedG> & G,
+  Eigen::PlainObjectBase<DerivedJ> & J,
+  Eigen::PlainObjectBase<Derivedflip> & flip)
+{
+  using namespace Eigen;
+  using namespace std;
+  using namespace igl;
+  typedef typename DerivedF::Index Index;
+  Matrix<Index,DerivedF::RowsAtCompileTime,1> C;
+  typedef Matrix<typename DerivedF::Scalar,Dynamic,DerivedF::ColsAtCompileTime> MatrixXF;
+  typedef Matrix<typename DerivedG::Scalar,Dynamic,DerivedG::ColsAtCompileTime> MatrixXG;
+  typedef Matrix<typename DerivedJ::Scalar,Dynamic,DerivedJ::ColsAtCompileTime> MatrixXJ;
+  typedef Matrix<typename Derivedflip::Scalar,Dynamic,Derivedflip::ColsAtCompileTime> MatrixXflip;
+  const Index m = F.rows();
+
+  typedef Matrix<typename DerivedF::Scalar,Dynamic,2> MatrixX2I;
+  typedef Matrix<typename DerivedF::Index,Dynamic,1> VectorXI;
+  MatrixX2I E,uE;
+  VectorXI EMAP;
+  vector<vector<typename DerivedF::Index> > uE2E;
+  unique_edge_map(F,E,uE,EMAP,uE2E);
+
+  vector<vector<vector<Index > > > TT,_1;
+  triangle_triangle_adjacency(E,EMAP,uE2E,false,TT,_1);
+  facet_components(TT,C);
+  const Index ncc = C.maxCoeff()+1;
+  G.resize(0,F.cols());
+  J.resize(0,1);
+  flip.setConstant(m,1,false);
+  // precompute face normals
+  typename Eigen::Matrix<
+    typename DerivedV::Scalar,
+    DerivedF::RowsAtCompileTime,
+    3> N;
+  per_face_normals(V,F,N);
+  // H contains list of faces on outer hull;
+  vector<bool> FH(m,false);
+  vector<bool> EH(3*m,false);
+  vector<MatrixXG> vG(ncc);
+  vector<MatrixXJ> vJ(ncc);
+
+  // Total size of G (and J)
+  size_t nG = 0;
+  // This is O( (n+m) * nnc) !
+  for(Index id = 0;id<ncc;id++)
+  {
+    // Determine number of facets in component id
+    Index num_id = 0;
+    for(Index f = 0;f<m;f++)
+    {
+      if(C(f) == id)
+      {
+        num_id++;
+      }
+    }
+    //MatrixXF Fc(num_id,F.cols());
+    MatrixXJ IM(num_id,1);
+    if(C.maxCoeff() == 0)
+    {
+      assert(num_id == F.rows());
+      //Fc = F;
+      IM = MatrixXJ::LinSpaced(F.rows(),0,F.rows()-1);
+    }else
+    {
+      int g = 0;
+      for(Index f = 0;f<m;f++)
+      {
+        if(C(f) == id)
+        {
+          //Fc.row(g) = F.row(f);
+          IM(g) = f;
+          g++;
+        }
+      }
+    }
+    {
+      // starting face that's guaranteed to be on the outer hull and in this
+      // component
+      int f;
+      bool f_flip;
+      outer_facet(V,F,N,IM,f,f_flip);
+      // Q contains list of face edges to continue traversing upong
+      queue<int> Q;
+      Q.push(f+0*m);
+      Q.push(f+1*m);
+      Q.push(f+2*m);
+      flip[f] = f_flip;
+      int FHcount = 0;
+      while(!Q.empty())
+      {
+        // face-edge
+        const int e = Q.front();
+        Q.pop();
+        // face
+        const int f = e%m;
+        // corner
+        const int c = e/m;
+        // Should never see edge again...
+        if(EH[e] == true)
+        {
+          continue;
+        }
+        EH[e] = true;
+        // first time seeing face
+        if(!FH[f])
+        {
+          FH[f] = true;
+          FHcount++;
+        }
+        // find overlapping face-edges
+        const auto & neighbors = uE2E[EMAP(e)];
+        const auto & fN = (flip[f]?-1.:1.)*N.row(f);
+        // source of edge according to f
+        const int fs = flip[f]?F(f,(c+2)%3):F(f,(c+1)%3);
+        // destination of edge according to f
+        const int fd = flip[f]?F(f,(c+1)%3):F(f,(c+2)%3);
+        const auto & eV = (V.row(fd)-V.row(fs)).normalized();
+        // Loop over and find max dihedral angle
+        typename DerivedV::Scalar max_di = -1;
+        int max_ne = -1;
+        typename Eigen::Matrix< typename DerivedV::Scalar, 1, 3> max_nN;
+        for(const auto & ne : neighbors)
+        {
+          const int nf = ne%m;
+          if(nf == f)
+          {
+            continue;
+          }
+          const int nc = ne/m;
+          // are faces consistently oriented
+          //const int ns = F(nf,(nc+1)%3);
+          const int nd = F(nf,(nc+2)%3);
+          const bool cons = (flip[f]?fd:fs) == nd;
+          const auto & nN = (cons? (flip[f]?-1:1.) : (flip[f]?1.:-1.) )*N.row(nf);
+          const auto & ndi = M_PI - atan2( fN.cross(nN).dot(eV), fN.dot(nN));
+          if(ndi>=max_di)
+          {
+            max_ne = ne;
+            max_di = ndi;
+            max_nN = nN;
+          }
+        }
+        if(max_ne>=0)
+        {
+          const int nf = max_ne%m;
+          const int nc = max_ne/m;
+          const int nd = F(nf,(nc+2)%3);
+          const bool cons = (flip[f]?fd:fs) == nd;
+          flip[nf] = (cons ? flip[f] : !flip[f]);
+          const int ne1 = nf+((nc+1)%3)*m;
+          const int ne2 = nf+((nc+2)%3)*m;
+          if(!EH[ne1])
+          {
+            Q.push(ne1);
+          }
+          if(!EH[ne2])
+          {
+            Q.push(ne2);
+          }
+        }
+      }
+      
+      {
+        vG[id].resize(FHcount,3);
+        vJ[id].resize(FHcount,1);
+        nG += FHcount;
+        size_t h = 0;
+        for(int i = 0;i<num_id;i++)
+        {
+          const size_t f = IM(i);
+          //if(f_flip)
+          //{
+          //  flip[f] = !flip[f];
+          //}
+          if(FH[f])
+          {
+            vG[id].row(h) = (flip[f]?F.row(f).reverse().eval():F.row(f));
+            vJ[id](h,0) = f;
+            h++;
+          }
+        }
+        assert(h == FHcount);
+      }
+    }
+  }
+  // collect G and J across components
+  G.resize(nG,3);
+  J.resize(nG,1);
+  {
+    size_t off = 0;
+    for(Index id = 0;id<ncc;id++)
+    {
+      assert(vG[id].rows() == vJ[id].rows());
+      G.block(off,0,vG[id].rows(),vG[id].cols()) = vG[id];
+      J.block(off,0,vJ[id].rows(),vJ[id].cols()) = vJ[id];
+      off += vG[id].rows();
+    }
+  }
+}
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+template void igl::outer_hull<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<bool, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/outer_hull.h

@@ -0,0 +1,40 @@
+#ifndef IGL_OUTER_HULL_H
+#define IGL_OUTER_HULL_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Compute the "outer hull" of a potentially non-manifold mesh (V,F) whose
+  // intersections have been "resolved" (e.g. using `cork` or
+  // `igl::selfintersect`). The outer hull is defined to be all facets
+  // (regardless of orientation) for which there exists some path from infinity
+  // to the face without intersecting any other facets. For solids, this is the
+  // surface of the solid. In general this includes any thin "wings" or "flaps".
+  // This implementation largely follows Section 3.6 of "Direct repair of
+  // self-intersecting meshes" [Attene 2014].
+  //
+  // Inputs:
+  //   V  #V by 3 list of vertex positions
+  //   F  #F by 3 list of triangle indices into V
+  // Outputs:
+  //   G  #G by 3 list of output triangle indices into V
+  //   J  #G list of indices into F
+  //   flip  #F list of whether facet was added to G **and** flipped orientation
+  //     (false for faces not added to G)
+  template <
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedG,
+    typename DerivedJ,
+    typename Derivedflip>
+  IGL_INLINE void outer_hull(
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    Eigen::PlainObjectBase<DerivedG> & G,
+    Eigen::PlainObjectBase<DerivedJ> & J,
+    Eigen::PlainObjectBase<Derivedflip> & flip);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "outer_hull.cpp"
+#endif
+#endif

include/igl/peal_outer_hull_layers.cpp

@@ -0,0 +1,80 @@
+#include <vector>
+#include "peal_outer_hull_layers.h"
+#include "outer_hull.h"
+
+using namespace std;
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename Derivedodd,
+  typename Derivedflip>
+IGL_INLINE void igl::peal_outer_hull_layers(
+  const Eigen::PlainObjectBase<DerivedV > & V,
+  const Eigen::PlainObjectBase<DerivedF > & F,
+  Eigen::PlainObjectBase<Derivedodd > & odd,
+  Eigen::PlainObjectBase<Derivedflip > & flip)
+{
+  using namespace Eigen;
+  using namespace std;
+  typedef typename DerivedF::Index Index;
+  typedef Matrix<typename DerivedF::Scalar,Dynamic,DerivedF::ColsAtCompileTime> MatrixXF;
+  typedef Matrix<Index,Dynamic,1> MatrixXI;
+  typedef Matrix<typename Derivedflip::Scalar,Dynamic,Derivedflip::ColsAtCompileTime> MatrixXflip;
+  const Index m = F.rows();
+
+  // keep track of iteration parity and whether flipped in hull
+  MatrixXF Fr = F;
+  odd.resize(m,1);
+  flip.resize(m,1);
+  // Keep track of index map
+  MatrixXI IM = MatrixXI::LinSpaced(m,0,m-1);
+  // This is O(n * layers)
+  bool odd_iter = true;
+  MatrixXI P(m,1);
+  Index iter = 0;
+  while(Fr.size() > 0)
+  {
+    assert(Fr.rows() == IM.rows());
+    // Compute outer hull of current Fr
+    MatrixXF Fo;
+    MatrixXI Jo;
+    MatrixXflip flipr;
+    outer_hull(V,Fr,Fo,Jo,flipr);
+    assert(Fo.rows() == Jo.rows());
+    // all faces in Fo of Fr
+    vector<bool> in_outer(Fr.rows(),false); 
+    for(Index g = 0;g<Jo.rows();g++)
+    {
+      odd(IM(Jo(g))) = odd_iter;
+      P(IM(Jo(g))) = iter;
+      in_outer[Jo(g)] = true;
+      flip(IM(Jo(g))) = flipr(Jo(g));
+    }
+    // Fr = Fr - Fo
+    // update IM
+    MatrixXF prev_Fr = Fr;
+    MatrixXI prev_IM = IM;
+    Fr.resize(prev_Fr.rows() - Fo.rows(),F.cols());
+    IM.resize(Fr.rows());
+    {
+      Index g = 0;
+      for(Index f = 0;f<prev_Fr.rows();f++)
+      {
+        if(!in_outer[f])
+        {
+          Fr.row(g) = prev_Fr.row(f);
+          IM(g) = prev_IM(f);
+          g++;
+        }
+      }
+    }
+    odd_iter = !odd_iter;
+    iter++;
+  }
+}
+
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+template void igl::peal_outer_hull_layers<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<bool, -1, 1, 0, -1, 1>, Eigen::Matrix<bool, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/peal_outer_hull_layers.h

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

include/igl/per_vertex_normals.cpp

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-//
+// 
 // Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
-//
-// This Source Code Form is subject to the terms of the Mozilla Public License
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "per_vertex_normals.h"
 
@@ -76,7 +76,8 @@ IGL_INLINE void igl::per_vertex_normals(
     // throw normal at each corner
     for(int j = 0; j < 3;j++)
     {
-      // Does this need to be critical?
+      // Q: Does this need to be critical?
+      // A: Yes. Different (i,j)'s could produce the same F(i,j)
 //#pragma omp critical
       N.row(F(i,j)) += W(i,j)*FN.row(i);
     }
@@ -98,7 +99,7 @@ IGL_INLINE void igl::per_vertex_normals(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template void igl::per_vertex_normals<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, igl::PerVertexNormalsWeightingType, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
 template void igl::per_vertex_normals<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
+template void igl::per_vertex_normals<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::per_vertex_normals<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, igl::PerVertexNormalsWeightingType, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
 #endif

include/igl/project.cpp

@@ -90,6 +90,7 @@ IGL_INLINE int igl::project(
   const Eigen::PlainObjectBase<Derivedobj> & obj,
   Eigen::PlainObjectBase<Derivedwin> & win)
 {
+  assert(obj.size() >= 3);
   Eigen::Vector3d dobj(obj(0),obj(1),obj(2));
   Eigen::Vector3d dwin;
   int ret = project(dobj(0),dobj(1),dobj(2),

include/igl/remove_unreferenced.cpp

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

include/igl/slice.cpp

@@ -205,4 +205,5 @@ template void igl::slice<std::__1::complex<double> >(Eigen::SparseMatrix<std::__
 template void igl::slice<Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::Matrix<double, -1, 3, 0, -1, 3> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, int, Eigen::Matrix<double, -1, 3, 0, -1, 3>&);
 template Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > igl::slice<Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&);
 template Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > igl::slice<Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, int);
+template Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> > igl::slice<Eigen::Matrix<double, 1, -1, 1, 1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, int);
 #endif

include/igl/sort.cpp

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-//
+// 
 // Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
-//
-// This Source Code Form is subject to the terms of the Mozilla Public License
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "sort.h"
 
@@ -218,6 +218,8 @@ template void igl::sort<int>(std::vector<int, std::allocator<int> > const&, bool
 template void igl::sort<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 // generated by autoexplicit.sh
 template void igl::sort<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
+template void igl::sort<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::sort<Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::sort<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::sort<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/triangle_triangle_adjacency.cpp

@@ -8,7 +8,8 @@
 #include "triangle_triangle_adjacency.h"
 #include "is_edge_manifold.h"
 #include "all_edges.h"
-#include <map>
+#include "unique_simplices.h"
+#include "unique_edge_map.h"
 #include <algorithm>
 
 template <typename Scalar, typename Index>
@@ -107,60 +108,104 @@ template <
     const Eigen::PlainObjectBase<DerivedF> & F,
     std::vector<std::vector<std::vector<TTIndex> > > & TT,
     std::vector<std::vector<std::vector<TTiIndex> > > & TTi)
+{
+  return triangle_triangle_adjacency(F,true,TT,TTi);
+}
+
+template <
+  typename DerivedF, 
+  typename TTIndex>
+  IGL_INLINE void igl::triangle_triangle_adjacency(
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    std::vector<std::vector<std::vector<TTIndex> > > & TT)
+{
+  std::vector<std::vector<std::vector<TTIndex> > > not_used;
+  return triangle_triangle_adjacency(F,false,TT,not_used);
+}
+
+template <
+  typename DerivedF, 
+  typename TTIndex, 
+  typename TTiIndex>
+  IGL_INLINE void igl::triangle_triangle_adjacency(
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    const bool construct_TTi,
+    std::vector<std::vector<std::vector<TTIndex> > > & TT,
+    std::vector<std::vector<std::vector<TTiIndex> > > & TTi)
 {
   using namespace Eigen;
   using namespace std;
-  using namespace igl;
   assert(F.cols() == 3 && "Faces must be triangles");
-  typedef typename DerivedF::Index Index;
   // number of faces
-  const int m = F.rows();
-  // All occurances of directed edges
-  MatrixXi E;
-  all_edges(F,E);
-  assert(E.rows() == 3*m);
-  // uE2E[i] --> {j,k,...} means unique edge i corresponds to face edges j and
-  // k (where j-edge comes is the j/m edge of face j%m)
-  map<pair<Index,Index>,vector<Index> > uE2E;
-  for(int e = 0;e<E.rows();e++)
-  {
-    Index i = E(e,0);
-    Index j = E(e,1);
-    if(i<j)
-    {
-      uE2E[pair<Index,Index>(i,j)].push_back(e);
-    }else
-    {
-      uE2E[pair<Index,Index>(j,i)].push_back(e);
-    }
-  }
+  typedef typename DerivedF::Index Index;
+  typedef Matrix<typename DerivedF::Scalar,Dynamic,2> MatrixX2I;
+  typedef Matrix<typename DerivedF::Index,Dynamic,1> VectorXI;
+  MatrixX2I E,uE;
+  VectorXI EMAP;
+  vector<vector<Index> > uE2E;
+  unique_edge_map(F,E,uE,EMAP,uE2E);
+  return triangle_triangle_adjacency(E,EMAP,uE2E,construct_TTi,TT,TTi);
+}
+
+template <
+  typename DerivedE, 
+  typename DerivedEMAP,
+  typename uE2EType,
+  typename TTIndex, 
+  typename TTiIndex>
+  IGL_INLINE void igl::triangle_triangle_adjacency(
+    const Eigen::PlainObjectBase<DerivedE> & E,
+    const Eigen::PlainObjectBase<DerivedEMAP> & EMAP,
+    const std::vector<std::vector<uE2EType> > & uE2E,
+    const bool construct_TTi,
+    std::vector<std::vector<std::vector<TTIndex> > > & TT,
+    std::vector<std::vector<std::vector<TTiIndex> > > & TTi)
+{
+  using namespace std;
+  using namespace Eigen;
+  typedef typename DerivedE::Index Index;
+  const size_t m = E.rows()/3;
+  assert((size_t)E.rows() == m*3 && "E should come from list of triangles.");
   // E2E[i] --> {j,k,...} means face edge i corresponds to other faces edges j
   // and k
-  TT.resize (m,vector<vector<Index> >(F.cols()));
-  TTi.resize(m,vector<vector<Index> >(F.cols()));
-  for(int e = 0;e<E.rows();e++)
+  TT.resize (m,vector<vector<TTIndex> >(3));
+  if(construct_TTi)
+  {
+    TTi.resize(m,vector<vector<TTiIndex> >(3));
+  }
+
+  // No race conditions because TT*[f][c]'s are in bijection with e's
+  // Minimum number of iterms per openmp thread
+  //const size_t num_e = E.rows();
+# ifndef IGL_OMP_MIN_VALUE
+#   define IGL_OMP_MIN_VALUE 1000
+# endif
+# pragma omp parallel for if (m>IGL_OMP_MIN_VALUE)
+  // Slightly better memory access than loop over E
+  for(Index f = 0;f<(Index)m;f++)
   {
-    const Index i = E(e,0);
-    const Index j = E(e,1);
-    const Index f = e%m;
-    const Index c = e/m;
-    const vector<Index> & N = 
-      i<j ? uE2E[pair<Index,Index>(i,j)] : uE2E[pair<Index,Index>(j,i)];
-    for(const auto & ne : N)
+    for(Index c = 0;c<3;c++)
     {
-      const Index nf = ne%m;
-      const Index nc = ne/m;
-      TT[f][c].push_back(nf);
-      TTi[f][c].push_back(nc);
+      const Index e = f + m*c;
+      //const Index c = e/m;
+      const vector<Index> & N = uE2E[EMAP(e)];
+      for(const auto & ne : N)
+      {
+        const Index nf = ne%m;
+        TT[f][c].push_back(nf);
+        if(construct_TTi)
+        {
+          const Index nc = ne/m;
+          TTi[f][c].push_back(nc);
+        }
+      }
     }
   }
 }
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template void igl::triangle_triangle_adjacency<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
+template void igl::triangle_triangle_adjacency<Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, long, long, long>(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> > const&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > > const&, bool, std::vector<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >, std::allocator<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > > > >&, std::vector<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >, std::allocator<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > > > >&);
 template void igl::triangle_triangle_adjacency<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
-template void igl::triangle_triangle_adjacency<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::triangle_triangle_adjacency<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
 #endif

include/igl/triangle_triangle_adjacency.h

@@ -76,6 +76,40 @@ namespace igl
       const Eigen::PlainObjectBase<DerivedF> & F,
       std::vector<std::vector<std::vector<TTIndex> > > & TT,
       std::vector<std::vector<std::vector<TTiIndex> > > & TTi);
+  template < typename DerivedF, typename TTIndex>
+    IGL_INLINE void triangle_triangle_adjacency(
+      const Eigen::PlainObjectBase<DerivedF> & F,
+      std::vector<std::vector<std::vector<TTIndex> > > & TT);
+  // Wrapper with bool to choose whether to compute TTi (this prototype should
+  // be "hidden").
+  template <
+    typename DerivedF, 
+    typename TTIndex, 
+    typename TTiIndex>
+    IGL_INLINE void triangle_triangle_adjacency(
+      const Eigen::PlainObjectBase<DerivedF> & F,
+      const bool construct_TTi,
+      std::vector<std::vector<std::vector<TTIndex> > > & TT,
+      std::vector<std::vector<std::vector<TTiIndex> > > & TTi);
+  // Inputs:
+  //   E  #F*3 by 2 list of all of directed edges in order (see `all_edges`)
+  //   EMAP #F*3 list of indices into uE, mapping each directed edge to unique
+  //     undirected edge
+  //   uE2E  #uE list of lists of indices into E of coexisting edges
+  // See also: unique_edge_map, all_edges
+  template <
+    typename DerivedE, 
+    typename DerivedEMAP,
+    typename uE2EType,
+    typename TTIndex, 
+    typename TTiIndex>
+    IGL_INLINE void triangle_triangle_adjacency(
+      const Eigen::PlainObjectBase<DerivedE> & E,
+      const Eigen::PlainObjectBase<DerivedEMAP> & EMAP,
+      const std::vector<std::vector<uE2EType > > & uE2E,
+      const bool construct_TTi,
+      std::vector<std::vector<std::vector<TTIndex> > > & TT,
+      std::vector<std::vector<std::vector<TTiIndex> > > & TTi);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/unique.cpp

@@ -229,4 +229,5 @@ template void igl::unique_rows<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::
 template void igl::unique_rows<Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::unique<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::unique<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::unique_rows<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/unique_edge_map.cpp

@@ -0,0 +1,42 @@
+#include "unique_edge_map.h"
+#include "all_edges.h"
+#include "unique_simplices.h"
+#include <cassert>
+#include <algorithm>
+template <
+  typename DerivedF,
+  typename DerivedE,
+  typename DeriveduE,
+  typename DerivedEMAP,
+  typename uE2EType>
+IGL_INLINE void igl::unique_edge_map(
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedE> & E,
+  Eigen::PlainObjectBase<DeriveduE> & uE,
+  Eigen::PlainObjectBase<DerivedEMAP> & EMAP,
+  std::vector<std::vector<uE2EType> > & uE2E)
+{
+  using namespace Eigen;
+  using namespace std;
+  // All occurances of directed edges
+  all_edges(F,E);
+  const size_t ne = E.rows();
+  // This is 2x faster to create than a map from pairs to lists of edges and 5x
+  // faster to access (actually access is probably assympotically faster O(1)
+  // vs. O(log m)
+  Matrix<typename DerivedEMAP::Scalar,Dynamic,1> IA;
+  unique_simplices(E,uE,IA,EMAP);
+  uE2E.resize(uE.rows());
+  // This does help a little
+  for_each(uE2E.begin(),uE2E.end(),[](vector<uE2EType > & v){v.reserve(2);});
+  assert((size_t)EMAP.size() == ne);
+  for(uE2EType e = 0;e<(uE2EType)ne;e++)
+  {
+    uE2E[EMAP(e)].push_back(e);
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+template void igl::unique_edge_map<Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, long>(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&);
+#endif 

include/igl/unique_edge_map.h

@@ -0,0 +1,37 @@
+#ifndef IGL_EXTERIOR_EDGES_H
+#define IGL_EXTERIOR_EDGES_H
+#include "igl_inline.h"
+#include <Eigen/Dense>
+#include <vector>
+namespace igl
+{
+  // Constuct relationships between facet "half"-(or rather "viewed")-edges E
+  // to unique edges of the mesh seen as a graph.
+  //
+  // Inputs:
+  //   F  #F by 3  list of simplices
+  // Outputs:
+  //   E  #F*3 by 2 list of all of directed edges
+  //   uE  #uE by 2 list of unique undirected edges
+  //   EMAP #F*3 list of indices into uE, mapping each directed edge to unique
+  //     undirected edge
+  //   uE2E  #uE list of lists of indices into E of coexisting edges
+  template <
+    typename DerivedF,
+    typename DerivedE,
+    typename DeriveduE,
+    typename DerivedEMAP,
+    typename uE2EType>
+  IGL_INLINE void unique_edge_map(
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    Eigen::PlainObjectBase<DerivedE> & E,
+    Eigen::PlainObjectBase<DeriveduE> & uE,
+    Eigen::PlainObjectBase<DerivedEMAP> & EMAP,
+    std::vector<std::vector<uE2EType> > & uE2E);
+
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "unique_edge_map.cpp"
+#endif
+
+#endif

include/igl/unique_simplices.cpp

@@ -24,7 +24,7 @@ IGL_INLINE void igl::unique_simplices(
   using namespace igl;
   // Sort each face
   MatrixXi sortF, unusedI;
-  igl::sort(F,2,1,sortF,unusedI);
+  igl::sort(F,2,true,sortF,unusedI);
   // Find unique faces
   MatrixXi C;
   igl::unique_rows(sortF,C,IA,IC);
@@ -50,4 +50,6 @@ IGL_INLINE void igl::unique_simplices(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instanciations
 template void igl::unique_simplices<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::unique_simplices<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
+template void igl::unique_simplices<Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/viewer/OpenGL_state.cpp

@@ -415,16 +415,16 @@ IGL_INLINE void igl::OpenGL_state::init()
   "vec3 Ia = La * Kai;"    // ambient intensity
 
   "vec3 light_position_eye = vec3 (view * vec4 (light_position_world, 1.0));"
-  "vec3 distance_to_light_eye = light_position_eye - position_eye;"
-  "vec3 direction_to_light_eye = normalize (distance_to_light_eye);"
+  "vec3 vector_to_light_eye = light_position_eye - position_eye;"
+  "vec3 direction_to_light_eye = normalize (vector_to_light_eye);"
   "float dot_prod = dot (direction_to_light_eye, normal_eye);"
-  "dot_prod = max (dot_prod, 0.0);"
-  "vec3 Id = Ld * Kdi * dot_prod;"    // Diffuse intensity
+  "float clamped_dot_prod = max (dot_prod, 0.0);"
+  "vec3 Id = Ld * Kdi * clamped_dot_prod;"    // Diffuse intensity
 
   "vec3 reflection_eye = reflect (-direction_to_light_eye, normal_eye);"
   "vec3 surface_to_viewer_eye = normalize (-position_eye);"
   "float dot_prod_specular = dot (reflection_eye, surface_to_viewer_eye);"
-  "dot_prod_specular = max (dot_prod_specular, 0.0);"
+  "dot_prod_specular = float(abs(dot_prod)==dot_prod) * max (dot_prod_specular, 0.0);"
   "float specular_factor = pow (dot_prod_specular, specular_exponent);"
   "vec3 Is = Ls * Ksi * specular_factor;"    // specular intensity
   "vec4 color = vec4(lighting_factor * (Is + Id) + Ia, 1.0) + vec4((1.0-lighting_factor) * Kdi,1.0);"

include/igl/viewer/Viewer.cpp

@@ -911,7 +911,7 @@ namespace igl
     glfwWindowHint(GLFW_OPENGL_PROFILE, GLFW_OPENGL_CORE_PROFILE);
     glfwWindowHint(GLFW_OPENGL_FORWARD_COMPAT, GL_TRUE);
     #endif
-    window = glfwCreateWindow(1280, 800, "IGL Viewer", NULL, NULL);
+    window = glfwCreateWindow(1280, 800, "libigl viewer", NULL, NULL);
     if (!window)
     {
       glfwTerminate();

include/igl/viewer/ViewerData.cpp

@@ -79,7 +79,6 @@ IGL_INLINE void igl::ViewerData::set_mesh(const Eigen::MatrixXd& _V, const Eigen
 
   if (V.rows() == 0 && F.rows() == 0)
   {
-    clear();
     V = V_temp;
     F = _F;
 

include/igl/writeOFF.cpp

@@ -55,4 +55,5 @@ template bool igl::writeOFF<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Mat
 template bool igl::writeOFF<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&);
 template bool igl::writeOFF<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&);
 template bool igl::writeOFF<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&);
+template bool igl::writeOFF<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&);
 #endif

index.html

@@ -0,0 +1,592 @@
+<!DOCTYPE html>
+<html>
+<head>
+	<meta charset="utf-8"/>
+	<title>libigl</title>
+	<meta name="author" content="Alec Jacobson and Daniele Panozzo and others"/>
+	<link type="text/css" rel="stylesheet" href="tutorial/style.css"/>
+<script type='text/javascript' src='http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML'></script>
+<link rel='stylesheet' href='http://yandex.st/highlightjs/7.3/styles/default.min.css'>
+<script src='http://yandex.st/highlightjs/7.3/highlight.min.js'></script>
+<script>hljs.initHighlightingOnLoad();</script>
+</head>
+<body>
+
+<h1 id="libigl-asimplecgeometryprocessinglibrary">libigl - A simple c++ geometry processing library</h1>
+
+<figure>
+<img src="tutorial/images/libigl-logo.jpg" alt="" />
+<figcaption></figcaption></figure>
+
+<p><a href="http://libigl.github.io/libigl/">http://libigl.github.io/libigl/</a>
+<a href="https://github.com/libigl/libigl/">https://github.com/libigl/libigl/</a></p>
+
+<p>Copyright 2014 - Alec Jacobson, Daniele Panozzo, Olga Diamanti, Kenshi
+Takayama, Leo Sacht, Wenzel Jacob, etc.</p>
+
+<p>Libigl is first and foremost a <em>header</em> library. Each header file should
+contain a single function. This function may have multiple overloads and
+prototypes. All functions should use the <code>igl::</code> namespace and should adhere to
+the conventions and styles listed in the <a href="style_guidelines.html">style
+guidelines</a>.</p>
+
+<blockquote>
+<p><strong>New:</strong> As of 1 July 2014, we have released libigl as beta version 1.0.
+There are a number of changes we collected for this release to minimize
+confusion and changes to how you use libigl. See <a href="#version1.0changes">Version 1.0
+Changes</a>.</p>
+</blockquote>
+
+<h2 id="installation">Installation</h2>
+
+<p>Libigl is a <em>header</em> library. You do <strong>not</strong> need to build anything to install.
+Simple add <code>igl/</code> to your include path and include relevant headers. Here&#8217;s a
+small &#8220;Hello, World&#8221; program:</p>
+
+<pre><code class="cpp">#include &lt;igl/cotmatrix.h&gt;
+#include &lt;Eigen/Dense&gt;
+#include &lt;Eigen/Sparse&gt;
+#include &lt;iostream&gt;
+int main()
+{
+  Eigen::MatrixXd V(4,2);
+  V&lt;&lt;0,0,
+     1,0,
+     1,1,
+     0,1;
+  Eigen::MatrixXi F(2,3);
+  F&lt;&lt;0,1,2,
+     0,2,3;
+  Eigen::SparseMatrix&lt;double&gt; L;
+  igl::cotmatrix(V,F,L);
+  std::cout&lt;&lt;&quot;Hello, mesh: &quot;&lt;&lt;std::endl&lt;&lt;L*V&lt;&lt;std::endl;
+  return 0;
+}
+</code></pre>
+
+<p>If you save this in <code>hello.cpp</code>, then on <code>gcc</code> with Eigen installed via
+macports for example you could compile this with:</p>
+
+<pre><code class="bash">gcc -I/opt/local/include/eigen3 -I./igl/ hello.cpp -o hello
+</code></pre>
+
+<p>Running <code>./hello</code> would then produce</p>
+
+<pre><code>Hello, mesh:
+ 0.5  0.5
+-0.5  0.5
+-0.5 -0.5
+ 0.5 -0.5
+</code></pre>
+
+<h2 id="tutorial">Tutorial</h2>
+
+<p>As of version 1.0, libigl includes an introductory
+<a href="http://libigl.github.io/libigl/tutorial/tutorial.html">tutorial</a> that covers
+its basic functionalities.</p>
+
+<h2 id="dependencies">Dependencies</h2>
+
+<ul>
+<li>Eigen3 Last tested with Eigen Version 3.2</li>
+</ul>
+
+<h3 id="optional">Optional</h3>
+
+<ul>
+<li>OpenGL (disable with <code>IGL_NO_OPENGL</code>)
+
+<ul>
+<li>OpenGL &gt;= 4 (enable with <code>IGL_OPENGL_4</code>)</li>
+</ul></li>
+<li>AntTweakBar (disable with <code>IGL_NO_ANTTWEAKBAR</code>) Last tested 1.16 (see
+ <code>libigl/external/AntTweakBar</code>)</li>
+<li>GLEW Windows and Linux</li>
+<li>OpenMP</li>
+<li>libpng libiglpng extra only</li>
+<li>Mosek libiglmosek extra only</li>
+<li>Matlab libiglmatlab extra only</li>
+<li>boost libiglboost, libiglcgal extra only</li>
+<li>SSE/AVX libiglsvd3x3 extra only</li>
+<li>CGAL libiglcgal extra only
+
+<ul>
+<li>boost</li>
+<li>gmp</li>
+<li>mpfr</li>
+</ul></li>
+<li>CoMiSo libcomiso extra only</li>
+</ul>
+
+<h3 id="optionalincludedinexternal">Optional (included in external/)</h3>
+
+<ul>
+<li>TetGen libigltetgen extra only</li>
+<li>Embree libiglembree extra only</li>
+<li>tinyxml2 libiglxml extra only</li>
+<li>glfw libviewer extra only</li>
+<li>LIM liblim extra only</li>
+</ul>
+
+<h2 id="headeronly">Header only</h2>
+
+<p>Libigl is designed to work &#8220;out-of-the-box&#8221; as a headers only library. To
+include libigl in your project. You need only include the libigl/include/
+directory in your include path. To
+compile a hello-word example.cpp:</p>
+
+<pre><code>#include &lt;Eigen/Dense&gt;
+#include &lt;igl/readOBJ.h&gt;
+#include &lt;iostream&gt;
+int main(int argc, char * argv[])
+{
+  if(argc&gt;1)
+  {
+    Eigen::MatrixXd V;
+    Eigen::MatrixXi F;
+    igl::readOBJ(argv[1],V,F);
+    std::cout&lt;&lt;&quot;Hello, mesh with &quot;&lt;&lt;V.rows()&lt;&lt;&quot; vertices!&quot;&lt;&lt;std::endl;
+  }else{
+    std::cout&lt;&lt;&quot;Hello, world!&quot;&lt;&lt;std::endl;
+  }
+  return 0;
+}
+</code></pre>
+
+<p>using gcc (replacing appropriate paths):</p>
+
+<pre><code>g++ -I/usr/local/igl/libigl/include \
+  -I/opt/local/include/eigen3 example.cpp -o example
+</code></pre>
+
+<p>Then run this example with:</p>
+
+<pre><code>./example examples/shared/TinyTorus.obj
+</code></pre>
+
+<h2 id="compilationasastaticlibrary">Compilation as a static library</h2>
+
+<p>Libigl is developed most often on Mac OS X, though has current users in Linux
+and Windows.</p>
+
+<h3 id="linuxmacosxcygwin">Linux/Mac OS X/Cygwin</h3>
+
+<p>Libigl may also be compiled to a static library. This is advantageous when
+building a project with libigl, since when used as an header-only library can
+slow down compile times.</p>
+
+<p>To build the entire libigl library producing lib/libigl.a, issue:</p>
+
+<pre><code>cd build
+make lib
+</code></pre>
+
+<p>You may need to edit Makefile.conf accordingly. Best to give yourself an
+<code>IGL_USERNAME</code> and add a custom install suite for yourself. Then you can enable
+appropriate extras.</p>
+
+<h4 id="extras">Extras</h4>
+
+<p>Once you&#8217;ve set up an <code>IGL_USERNAME</code> and enabled extras within Makefile.conf.
+You can build the extra libraries (into lib/ligiglpng.a, lib/libiglmatlab.a,
+lib/libigltetgen.a, lib/libiglmosek.a, etc.) by issuing:</p>
+
+<pre><code>cd build
+make extras
+</code></pre>
+
+<h4 id="examples">Examples</h4>
+
+<p>You can make a slew of examples by issuing:</p>
+
+<pre><code>cd build
+make examples
+</code></pre>
+
+<h4 id="external">External</h4>
+
+<p>Finally there are a number of external libraries that we include in
+./external/ because they are either difficult to obtain or they have been
+patched for easier use with libigl. Please see the respective readmes in
+those directories.</p>
+
+<h5 id="installinganttweakbar">Installing AntTweakBar</h5>
+
+<p>To build the a static AntTweakBar library on Mac OS X issue:</p>
+
+<pre><code>cd external/AntTweakBar/src
+make -f Makefile.osx.igl
+</code></pre>
+
+<h5 id="installingtetgen">Installing Tetgen</h5>
+
+<p>To build the tetgen library and executable on Mac OS X issue:</p>
+
+<pre><code>cd external/tetgen
+make clean
+rm -f obj/*.o
+make -f Makefile.igl tetgen
+rm -f obj/*.o
+make -f Makefile.igl tetlib
+</code></pre>
+
+<h5 id="installingmedit">Installing medit</h5>
+
+<p>To build the igl version of the medit executable on Mac OS X issue:</p>
+
+<pre><code>cd external/medit
+make -C libmesh
+make -f Makefile.igl medit
+</code></pre>
+
+<h5 id="installingembree2.0">Installing Embree 2.0</h5>
+
+<p>To build the embree library and executables on Mac OS X issue:</p>
+
+<pre><code>cd external/embree
+mkdir build
+cd build
+cmake ..
+# Or using a different compiler
+#cmake .. -DCMAKE_C_COMPILER=/opt/local/bin/gcc -DCMAKE_CXX_COMPILER=/opt/local/bin/g++
+make
+# Could also install embree to your root, but libigl examples don't expect
+# this
+#sudo make install
+</code></pre>
+
+<h5 id="installingtinyxml2">Installing tinyxml2</h5>
+
+<p>To build the a static tinyxml2 library on Mac OS X issue:</p>
+
+<pre><code>cd external/tinyxml2
+cmake .
+make
+</code></pre>
+
+<h5 id="installingyimg">Installing YImg</h5>
+
+<p>To build the a static YImg library on Mac OS X issue:</p>
+
+<pre><code>cd external/yimg
+make
+</code></pre>
+
+<p>You may need to install libpng. Systems with X11 might find this already
+installed at <code>/usr/X11/lib</code>.</p>
+
+<h3 id="windowsexperimental">Windows (Experimental)</h3>
+
+<p>To build a static library (.lib) on windows, open Visual Studio 2010.</p>
+
+<ul>
+<li>New &gt; Project &#8230;</li>
+<li>Visual C++ &gt; Win32</li>
+<li>Win32 Console Application</li>
+<li>Name: libiglVisualStudio</li>
+<li>Uncheck &#8220;Create directory for solution&#8221;</li>
+<li>Then hit OK, and then Next</li>
+<li>Check &#8220;Static Library&#8221;</li>
+<li>Uncheck &#8220;Precompiled headers&#8221;</li>
+<li>Add all include/igl/*.cpp to the sources directory</li>
+<li>Add all include/igl/*.h to the headers directory</li>
+<li>Open Project &gt; libigl Properties&#8230;</li>
+<li>Add the path to eigen3 to the include paths</li>
+<li>Change the target name to libigl</li>
+<li>Build and pray (this should create libigl.lib</li>
+</ul>
+
+<p><a href="http://msdn.microsoft.com/en-us/library/ms235627(v=vs.80).aspx">Source</a></p>
+
+<h2 id="examples">Examples</h2>
+
+<p>To get started, we advise that you take a look at a few examples:</p>
+
+<pre><code>./examples/hello-world/
+
+./examples/meshio/
+
+./examples/basic-topology/
+
+./examples/ReAntTweakBar/
+</code></pre>
+
+<h2 id="extras">Extras</h2>
+
+<p>Libigl compartmentalizes dependences via its organization into a <em>main</em> libigl
+library and &#8220;extras.&#8221;</p>
+
+<h3 id="bbw">bbw</h3>
+
+<p>This library extra contains functions for computing Bounded Biharmonic Weights, can
+be used with and without the <a href="#mosek">mosek</a> extra via the <code>IGL_NO_MOSEK</code>
+macro.</p>
+
+<h3 id="boost">boost</h3>
+
+<p>This library extra utilizes the graph functions in the boost library for find
+connected components and performing breadth-first traversals.</p>
+
+<h3 id="cgal">cgal</h3>
+
+<p>This library extra utilizes CGAL&#8217;s efficient and exact intersection and
+proximity queries.</p>
+
+<h3 id="embree">embree</h3>
+
+<p>This library extra utilizes embree&#8217;s efficient ray tracing queries.</p>
+
+<h3 id="matlab">matlab</h3>
+
+<p>This library extra provides support for reading and writing <code>.mat</code> workspace
+files, interfacing with Matlab at run time and compiling mex functions.</p>
+
+<h3 id="mosek">mosek</h3>
+
+<p>This library extra utilizes mosek&#8217;s efficient interior-point solver for
+quadratic programs.</p>
+
+<h3 id="png">png</h3>
+
+<p>This library extra uses <code>libpng</code> and <code>YImage</code> to read and write <code>.png</code> files.</p>
+
+<h3 id="svd3x3">svd3x3</h3>
+
+<p>This library extra implements &#8220;as-rigid-as-possible&#8221; (ARAP) deformation
+techniques using the fast singular value decomposition routines
+written specifically for 3x3 matrices to use <code>SSE</code> intrinsics. This extra can
+still be compiled without sse support and support should be determined
+automatically at compile time via the <code>__SSE__</code> macro.</p>
+
+<h3 id="tetgen">tetgen</h3>
+
+<p>This library extra provides a simplified wrapper to the tetgen 3d tetrahedral meshing
+library.</p>
+
+<h3 id="viewer">viewer</h3>
+
+<p>This library extra utilizes glfw and glew to open an opengl context and launch
+a simple mesh viewer.</p>
+
+<h3 id="xml">xml</h3>
+
+<p>This library extra utilizes tinyxml2 to read and write serialized classes
+containing Eigen matrices and other standard simple data-structures.</p>
+
+<h2 id="development">Development</h2>
+
+<p>Further documentation for developers is listed in tutorial.html,
+style_guidelines.html</p>
+
+<h2 id="license">License</h2>
+
+<p>See <code>LICENSE.txt</code></p>
+
+<h2 id="zipping">Zipping</h2>
+
+<p>Zip this directory without .git litter and binaries using:</p>
+
+<pre><code>git archive -prefix=libigl/ -o libigl.zip master
+</code></pre>
+
+<h2 id="version1.0changes">Version 1.0 Changes</h2>
+
+<p>Our beta release marks our confidence that this library can be used outside of
+casual experimenting. To maintain order, we have made a few changes which
+current users should read and adapt their code accordingly.</p>
+
+<h3 id="renamedfunctions">Renamed functions</h3>
+
+<p>The following table lists functions which have changed name as of version
+1.0.0:</p>
+
+<table>
+<colgroup>
+<col style="text-align:left;"/>
+<col style="text-align:left;"/>
+</colgroup>
+
+<thead>
+<tr>
+	<th style="text-align:left;">Old</th>
+	<th style="text-align:left;">New</th>
+</tr>
+</thead>
+
+<tbody>
+<tr>
+	<td style="text-align:left;"><code>igl::add_barycenter</code></td>
+	<td style="text-align:left;"><code>igl::false_barycentric_subdivision</code></td>
+</tr>
+<tr>
+	<td style="text-align:left;"><code>igl::areamatrix</code></td>
+	<td style="text-align:left;"><code>igl::vector_area_matrix</code></td>
+</tr>
+<tr>
+	<td style="text-align:left;"><code>igl::barycentric2global</code></td>
+	<td style="text-align:left;"><code>igl::barycentric_to_global</code></td>
+</tr>
+<tr>
+	<td style="text-align:left;"><code>igl::boundary_faces</code></td>
+	<td style="text-align:left;"><code>igl::boundary_facets</code></td>
+</tr>
+<tr>
+	<td style="text-align:left;"><code>igl::boundary_vertices_sorted</code></td>
+	<td style="text-align:left;"><code>igl::boundary_loop</code></td>
+</tr>
+<tr>
+	<td style="text-align:left;"><code>igl::cotangent</code></td>
+	<td style="text-align:left;"><code>igl::cotmatrix_entries</code></td>
+</tr>
+<tr>
+	<td style="text-align:left;"><code>igl::edgetopology</code></td>
+	<td style="text-align:left;"><code>igl::edge_topology</code></td>
+</tr>
+<tr>
+	<td style="text-align:left;"><code>igl::gradMat</code></td>
+	<td style="text-align:left;"><code>igl::grad</code></td>
+</tr>
+<tr>
+	<td style="text-align:left;"><code>igl::is_manifold</code></td>
+	<td style="text-align:left;"><code>igl::is_edge_manifold</code></td>
+</tr>
+<tr>
+	<td style="text-align:left;"><code>igl::mexStream</code></td>
+	<td style="text-align:left;"><code>igl::MexStream</code></td>
+</tr>
+<tr>
+	<td style="text-align:left;"><code>igl::moveFV</code></td>
+	<td style="text-align:left;"><code>igl::average_onto_vertices</code></td>
+</tr>
+<tr>
+	<td style="text-align:left;"><code>igl::moveVF</code></td>
+	<td style="text-align:left;"><code>igl::average_onto_faces</code></td>
+</tr>
+<tr>
+	<td style="text-align:left;"><code>igl::plot_vector</code></td>
+	<td style="text-align:left;"><code>igl::print_vector</code></td>
+</tr>
+<tr>
+	<td style="text-align:left;"><code>igl::pos</code></td>
+	<td style="text-align:left;"><code>igl::HalfEdgeIterator</code></td>
+</tr>
+<tr>
+	<td style="text-align:left;"><code>igl::plane_project</code></td>
+	<td style="text-align:left;"><code>igl::project_isometrically_to_plane</code></td>
+</tr>
+<tr>
+	<td style="text-align:left;"><code>igl::project_points_mesh</code></td>
+	<td style="text-align:left;"><code>igl::line_mesh_intersection</code></td>
+</tr>
+<tr>
+	<td style="text-align:left;"><code>igl::read</code></td>
+	<td style="text-align:left;"><code>igl::read_triangle_mesh</code></td>
+</tr>
+<tr>
+	<td style="text-align:left;"><code>igl::removeDuplicates.cpp</code></td>
+	<td style="text-align:left;"><code>igl::remove_duplicates</code></td>
+</tr>
+<tr>
+	<td style="text-align:left;"><code>igl::removeUnreferenced</code></td>
+	<td style="text-align:left;"><code>igl::remove_unreferenced</code></td>
+</tr>
+<tr>
+	<td style="text-align:left;"><code>igl::tt</code></td>
+	<td style="text-align:left;"><code>igl::triangle_triangle_adjacency</code></td>
+</tr>
+<tr>
+	<td style="text-align:left;"><code>igl::vf</code></td>
+	<td style="text-align:left;"><code>igl::vertex_triangle_adjacency</code></td>
+</tr>
+<tr>
+	<td style="text-align:left;"><code>igl::write</code></td>
+	<td style="text-align:left;"><code>igl::write_triangle_mesh</code></td>
+</tr>
+<tr>
+	<td style="text-align:left;"><code>igl::manifold_patches</code></td>
+	<td style="text-align:left;"><code>igl::orientable_patches</code></td>
+</tr>
+<tr>
+	<td style="text-align:left;"><code>igl::selfintersect</code></td>
+	<td style="text-align:left;"><code>igl::remesh_self_intersections</code></td>
+</tr>
+<tr>
+	<td style="text-align:left;"><code>igl::project_mesh</code></td>
+	<td style="text-align:left;"><code>igl::line_mesh_intersection</code></td>
+</tr>
+<tr>
+	<td style="text-align:left;"><code>igl::triangulate</code></td>
+	<td style="text-align:left;"><code>igl::polygon_mesh_to_triangle_mesh</code></td>
+</tr>
+<tr>
+	<td style="text-align:left;"><code>igl::is_manifold</code></td>
+	<td style="text-align:left;"><code>igl::is_edge_manifold</code></td>
+</tr>
+<tr>
+	<td style="text-align:left;"><code>igl::triangle_wrapper</code></td>
+	<td style="text-align:left;"><code>igl::triangulate</code></td>
+</tr>
+</tbody>
+</table>
+
+<h3 id="miscellaneous">Miscellaneous</h3>
+
+<ul>
+<li>To match interfaces provided by (all) other quadratic optimization
+ libraries, <code>igl::min_quad_with_fixed</code> and <code>igl::active_set</code> now expect as
+ input twice the quadratic coefficients matrix, i.e. the Hessian. For
+ example, <code>igl::min_quad_with_fixed(H,B,...)</code> minimizes <span class="math">\(\frac{1}{2}x^T H
+   x+x^T B\)</span>.</li>
+<li>We have inverted the <code>IGL_HEADER_ONLY</code> macro to <code>IGL_STATIC_LIBRARY</code>. To
+ compile using libigl as a header-only library, simply include headers and
+ libigl in the header search path. To link to libigl, you must define the
+ <code>IGL_STATIC_LIBRARY</code> macro at compile time and link to the <code>libigl*.a</code>
+ libraries.</li>
+<li>Building libigl as a static library is now more organized. There is a
+ <code>build/</code> directory with Makefiles for the main library (<code>Makefile</code>) and each
+ dependency (e.g. <code>Makefile_mosek</code> for <code>libiglmosek.a</code>)</li>
+<li><code>igl::polar_svd</code> now always returns a rotation in <code>R</code>, never a reflection.
+ This mirrors the behavior of <code>igl::polar_svd3x3</code>. Consequently the <code>T</code>
+ part may have negative skews.</li>
+<li>We have organized the static</li>
+<li>The previous <code>igl::grad</code> function, which computed the per-triangle gradient
+ of a per-vertex scalar function has been replaced. Now <code>igl::grad</code> computes
+ the linear operator (previous computed using <code>igl::gradMat</code>). The gradient
+ values can still be recovered by multiplying the operator against the scalar
+ field as a vector and reshaping to have gradients per row.</li>
+<li><code>MASSMATRIX_*</code> has become <code>MASSMATRIX_TYPE_*</code></li>
+<li>The function <code>igl::project_normals</code>, which cast a line for each vertex of
+ mesh <em>A</em> in the normal direction and found the closest intersection along
+ these lines with mesh <em>B</em>, has been removed.</li>
+</ul>
+
+<h2 id="contact">Contact</h2>
+
+<p>Libigl is a group endeavor led by Alec Jacobson and Daniele Panozzo. Please
+contact <a href="&#109;&#97;&#105;&#108;&#x74;&#111;&#x3a;&#x61;&#x6c;&#101;&#x63;&#x6a;&#x61;&#x63;&#111;&#98;&#115;&#111;&#110;&#64;&#103;&#109;&#x61;&#105;&#108;&#x2e;&#x63;&#x6f;&#x6d;">&#97;&#x6c;&#x65;&#99;&#106;&#97;&#x63;&#111;&#x62;&#x73;&#x6f;&#110;&#x40;&#103;&#x6d;&#x61;&#105;&#x6c;&#x2e;&#x63;&#x6f;&#x6d;</a> if you have
+questions or comments. We are happy to get feedback! Enjoy!</p>
+
+<p>If you&#8217;re using libigl in your projects, quickly <a href="&#x6d;&#97;&#x69;&#x6c;&#116;&#x6f;&#x3a;&#97;&#x6c;&#101;&#x63;&#106;&#97;&#x63;&#x6f;&#98;&#115;&#x6f;&#x6e;&#64;&#103;&#x6d;&#97;&#105;&#108;&#x2e;&#x63;&#111;&#x6d;">&#x64;&#x72;&#111;&#112; &#117;&#x73; &#x61;
+&#x6e;&#x6f;&#116;&#101;</a>. Tell us who you are and what you&#8217;re using
+it for. This helps us apply for funding and justify spending time maintaining
+this.</p>
+
+<p>If you find bugs or have problems please use our <a href="https://github.com/libigl/libigl/issues">github issue tracking
+page</a>.</p>
+
+<h2 id="academiccitation">Academic citation</h2>
+
+<p>If you use libigl in your research projects, please cite the papers we
+implement as appropriate. To cite the library in general, you could use this
+BibTeX entry:</p>
+
+<pre><code class="bibtex">@misc{libigl,
+  title = {{libigl}: A simple {C++} geometry processing library},
+  author = {Alec Jacobson and Daniele Panozzo and others},
+  note = {http://libigl.github.io/libigl/},
+  year = {2014},
+}
+</code></pre>
+
+</body>
+</html>

tutorial/403_BoundedBiharmonicWeights/CMakeLists.txt

@@ -3,9 +3,13 @@ project(403_BoundedBiharmonicWeights)
 
 include("../CMakeLists.shared")
 
+find_package(MOSEK REQUIRED)
+
+include_directories( ${MOSEK_INCLUDE_DIR} )
+
 set(SOURCES
 ${PROJECT_SOURCE_DIR}/main.cpp
 )
 
 add_executable(${PROJECT_NAME}_bin ${SOURCES} ${SHARED_SOURCES})
-target_link_libraries(${PROJECT_NAME}_bin ${SHARED_LIBRARIES})
+target_link_libraries(${PROJECT_NAME}_bin ${SHARED_LIBRARIES} ${MOSEK_LIBRARIES})

tutorial/609_Boolean/CMakeLists.txt

@@ -0,0 +1,15 @@
+cmake_minimum_required(VERSION 2.6)
+project(609_Boolean)
+
+find_package(CGAL REQUIRED)
+include(${CGAL_USE_FILE})
+
+# for some reason must come after cgal include
+include("../CMakeLists.shared")
+
+set(SOURCES
+${PROJECT_SOURCE_DIR}/main.cpp
+)
+
+add_executable(${PROJECT_NAME}_bin ${SOURCES} ${SHARED_SOURCES})
+target_link_libraries(${PROJECT_NAME}_bin ${SHARED_LIBRARIES} ${CGAL_LIBRARIES})

tutorial/609_Boolean/main.cpp

@@ -0,0 +1,94 @@
+#include <igl/readOFF.h>
+#define IGL_NO_CORK
+#include <igl/boolean/mesh_boolean.h>
+#include <igl/viewer/Viewer.h>
+
+#include <Eigen/Core>
+#include <iostream>
+
+Eigen::MatrixXd VA,VB,VC;
+Eigen::VectorXi J;
+Eigen::MatrixXi FA,FB,FC;
+igl::MeshBooleanType boolean_type(igl::MESH_BOOLEAN_TYPE_UNION);
+
+const char * MESH_BOOLEAN_TYPE_NAMES[] = 
+{
+  "Union",
+  "Intersect",
+  "Minus",
+  "XOR",
+  "Resolve",
+};
+
+void update(igl::Viewer &viewer)
+{
+  igl::mesh_boolean(VA,FA,VB,FB,boolean_type,VC,FC,J);
+  Eigen::MatrixXd C(FC.rows(),3);
+  for(size_t f = 0;f<C.rows();f++)
+  {
+    if(J(f)<FA.rows())
+    {
+      C.row(f) = Eigen::RowVector3d(1,0,0);
+    }else
+    {
+      C.row(f) = Eigen::RowVector3d(0,1,0);
+    }
+  }
+  viewer.data.clear();
+  viewer.data.set_mesh(VC,FC);
+  viewer.data.set_colors(C);
+}
+
+bool key_down(igl::Viewer &viewer, unsigned char key, int mods)
+{
+  switch(key)
+  {
+    default:
+      return false;
+    case '.':
+      boolean_type = 
+        static_cast<igl::MeshBooleanType>(
+          (boolean_type+1)% igl::NUM_MESH_BOOLEAN_TYPES);
+      break;
+    case ',':
+      boolean_type = 
+        static_cast<igl::MeshBooleanType>(
+          (boolean_type+igl::NUM_MESH_BOOLEAN_TYPES-1)%
+          igl::NUM_MESH_BOOLEAN_TYPES);
+      break;
+    case '[':
+      viewer.core.camera_dnear -= 0.1;
+      return true;
+    case ']':
+      viewer.core.camera_dnear += 0.1;
+      return true;
+  }
+  std::cout<<"A "<<MESH_BOOLEAN_TYPE_NAMES[boolean_type]<<" B."<<std::endl;
+  igl::mesh_boolean(VA,FA,VB,FB,boolean_type,VC,FC);
+  update(viewer);
+  return true;
+}
+
+int main(int argc, char *argv[])
+{
+  using namespace Eigen;
+  using namespace std;
+  igl::readOFF("../shared/cheburashka.off",VA,FA);
+  igl::readOFF("../shared/decimated-knight.off",VB,FB);
+  // Plot the mesh with pseudocolors
+  igl::Viewer viewer;
+
+  // Initialize
+  update(viewer);
+
+  viewer.core.show_lines = true;
+  viewer.callback_key_down = &key_down;
+  viewer.core.camera_dnear = 3.9;
+  cout<<
+    "Press '.' to switch to next boolean operation type."<<endl<<
+    "Press ',' to switch to previous boolean operation type."<<endl<<
+    "Press ']' to push near cutting plane away from camera."<<endl<<
+    "Press '[' to pull near cutting plane closer to camera."<<endl<<
+    "Hint: investigate _inside_ the model to see orientation changes."<<endl;
+  viewer.launch();
+}

tutorial/CMakeLists.shared

@@ -56,6 +56,8 @@ link_directories(
 
 # Disable deprecated opengl code from libigl
 add_definitions(-DIGL_OPENGL_4)
+add_definitions(-DNDEBUG)
+add_definitions(-O3)
 
 
 if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang" OR "${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU")

tutorial/CMakeLists.txt

@@ -12,6 +12,7 @@ find_package(MATLAB QUIET)
 
 ## Check for EMBREE, if not available skip the examples that depends on it
 find_package(EMBREE QUIET)
+find_package(CGAL QUIET)
 
 # Chapter 1
 add_subdirectory("101_FileIO")
@@ -68,3 +69,6 @@ add_subdirectory("606_AmbientOcclusion")
 add_subdirectory("607_Picking")
 endif(EMBREE_FOUND)
 add_subdirectory("608_LIM")
+if(CGAL_FOUND)
+add_subdirectory("609_Boolean")
+endif()

tutorial/cmake/FindCGAL.cmake

@@ -0,0 +1,92 @@
+#
+# The following module is based on FindVTK.cmake
+#
+
+# - Find a CGAL installation or binary tree.
+# The following variables are set if CGAL is found.  If CGAL is not
+# found, CGAL_FOUND is set to false.
+#
+#  CGAL_FOUND         - Set to true when CGAL is found.
+#  CGAL_USE_FILE      - CMake file to use CGAL.
+#
+
+# Construct consitent error messages for use below.
+set(CGAL_DIR_DESCRIPTION "directory containing CGALConfig.cmake. This is either the binary directory where CGAL was configured or PREFIX/lib/CGAL for an installation.")
+set(CGAL_DIR_MESSAGE     "CGAL not found.  Set the CGAL_DIR cmake variable or environment variable to the ${CGAL_DIR_DESCRIPTION}")
+ 
+set(CMAKE_ALLOW_LOOSE_LOOP_CONSTRUCTS true)
+ 
+if ( NOT CGAL_DIR )
+  
+  # Get the system search path as a list.
+  if(UNIX)
+    string(REGEX MATCHALL "[^:]+" CGAL_DIR_SEARCH1 "$ENV{PATH}")
+  else()
+    string(REGEX REPLACE "\\\\" "/" CGAL_DIR_SEARCH1 "$ENV{PATH}")
+  endif()
+  
+  string(REGEX REPLACE "/;" ";" CGAL_DIR_SEARCH2 "${CGAL_DIR_SEARCH1}")
+
+  # Construct a set of paths relative to the system search path.
+  set(CGAL_DIR_SEARCH "")
+  
+  foreach(dir ${CGAL_DIR_SEARCH2})
+  
+    set(CGAL_DIR_SEARCH ${CGAL_DIR_SEARCH} ${dir}/../lib/CGAL )
+      
+  endforeach()
+
+
+  #
+  # Look for an installation or build tree.
+  #
+  find_path(CGAL_DIR CGALConfig.cmake
+
+    # Look for an environment variable CGAL_DIR.
+    $ENV{CGAL_DIR}
+
+    # Look in places relative to the system executable search path.
+    ${CGAL_DIR_SEARCH}
+
+    # Look in standard UNIX install locations.
+    /opt/local/share/CGAL/cmake
+    /usr/local/lib/CGAL
+    /usr/lib/CGAL
+
+    # Read from the CMakeSetup registry entries.  It is likely that
+    # CGAL will have been recently built.
+    [HKEY_CURRENT_USER\\Software\\Kitware\\CMakeSetup\\Settings\\StartPath;WhereBuild1]
+    [HKEY_CURRENT_USER\\Software\\Kitware\\CMakeSetup\\Settings\\StartPath;WhereBuild2]
+    [HKEY_CURRENT_USER\\Software\\Kitware\\CMakeSetup\\Settings\\StartPath;WhereBuild3]
+    [HKEY_CURRENT_USER\\Software\\Kitware\\CMakeSetup\\Settings\\StartPath;WhereBuild4]
+    [HKEY_CURRENT_USER\\Software\\Kitware\\CMakeSetup\\Settings\\StartPath;WhereBuild5]
+    [HKEY_CURRENT_USER\\Software\\Kitware\\CMakeSetup\\Settings\\StartPath;WhereBuild6]
+    [HKEY_CURRENT_USER\\Software\\Kitware\\CMakeSetup\\Settings\\StartPath;WhereBuild7]
+    [HKEY_CURRENT_USER\\Software\\Kitware\\CMakeSetup\\Settings\\StartPath;WhereBuild8]
+    [HKEY_CURRENT_USER\\Software\\Kitware\\CMakeSetup\\Settings\\StartPath;WhereBuild9]
+    [HKEY_CURRENT_USER\\Software\\Kitware\\CMakeSetup\\Settings\\StartPath;WhereBuild10]
+
+    # Help the user find it if we cannot.
+    DOC "The ${CGAL_DIR_DESCRIPTION}"
+  )
+  
+endif()
+
+if ( CGAL_DIR )
+  
+  if ( EXISTS "${CGAL_DIR}/CGALConfig.cmake" )
+    include( "${CGAL_DIR}/CGALConfig.cmake" )
+    set( CGAL_FOUND TRUE )
+  endif()
+
+endif()
+
+if( NOT CGAL_FOUND)
+  if(CGAL_FIND_REQUIRED)
+    MESSAGE(FATAL_ERROR ${CGAL_DIR_MESSAGE})
+  else()
+    if(NOT CGAL_FIND_QUIETLY)
+      MESSAGE(STATUS ${CGAL_DIR_MESSAGE})
+    endif()
+  endif()
+endif()

tutorial/cmake/FindLIBIGL.cmake

@@ -12,6 +12,8 @@ find_package(Mosek QUIET)
 FIND_PATH(LIBIGL_INCLUDE_DIR igl/readOBJ.h
    /usr/include
    /usr/local/include
+   /usr/local/igl/libigl/include
+   $ENV{LIBIGL}/include
    $ENV{LIBIGLROOT}/include
    $ENV{LIBIGL_ROOT}/include
    $ENV{LIBIGL_DIR}/include
@@ -47,17 +49,20 @@ if(LIBIGL_USE_STATIC_LIBRARY)
   FIND_LIBRARY( LIBIGL_LIBRARY NAMES igl PATHS ${LIBIGL_LIB_DIRS})
   if(NOT LIBIGL_LIBRARY)
     set(LIBIGL_FOUND FALSE)
+    message(FATAL_ERROR "could NOT find libigl")
   endif(NOT LIBIGL_LIBRARY)
   set(LIBIGL_LIBRARIES ${LIBIGL_LIBRARIES}  ${LIBIGL_LIBRARY})
 
   FIND_LIBRARY( LIBIGLBBW_LIBRARY NAMES iglbbw PATHS ${LIBIGL_LIB_DIRS})
   if(NOT LIBIGLBBW_LIBRARY)
     set(LIBIGL_FOUND FALSE)
+    message(FATAL_ERROR "could NOT find libiglbbw")
   endif(NOT LIBIGLBBW_LIBRARY)
   set(LIBIGL_LIBRARIES ${LIBIGL_LIBRARIES}  ${LIBIGLBBW_LIBRARY})
   FIND_LIBRARY( LIBIGLMOSEK_LIBRARY NAMES iglmosek PATHS ${LIBIGL_LIB_DIRS})
   if(NOT LIBIGLMOSEK_LIBRARY)
     set(LIBIGL_FOUND FALSE)
+    message(FATAL_ERROR "could NOT find libiglmosek")
   endif(NOT LIBIGLMOSEK_LIBRARY)
   set(LIBIGL_LIBRARIES ${LIBIGL_LIBRARIES}  ${LIBIGLMOSEK_LIBRARY})
   if(MOSEK_FOUND)
@@ -65,11 +70,27 @@ if(LIBIGL_USE_STATIC_LIBRARY)
     set(LIBIGL_LIBRARIES ${LIBIGL_LIBRARIES}  ${MOSEK_LIBRARIES})
   else(MOSEK_FOUND)
     set(LIBIGL_FOUND FALSE)
+    message(FATAL_ERROR "could NOT find mosek")
   endif(MOSEK_FOUND)
 
+  FIND_LIBRARY( LIBIGLCGAL_LIBRARY NAMES iglcgal PATHS ${LIBIGL_LIB_DIRS})
+  if(NOT LIBIGLCGAL_LIBRARY)
+    set(LIBIGL_FOUND FALSE)
+    message(FATAL_ERROR "could NOT find libiglcgal")
+  endif(NOT LIBIGLCGAL_LIBRARY)
+  set(LIBIGL_LIBRARIES ${LIBIGL_LIBRARIES}  ${LIBIGLCGAL_LIBRARY})
+
+  FIND_LIBRARY( LIBIGLBOOLEAN_LIBRARY NAMES iglboolean PATHS ${LIBIGL_LIB_DIRS})
+  if(NOT LIBIGLBOOLEAN_LIBRARY)
+    set(LIBIGL_FOUND FALSE)
+    message(FATAL_ERROR "could NOT find libiglboolean ")
+  endif(NOT LIBIGLBOOLEAN_LIBRARY)
+  set(LIBIGL_LIBRARIES ${LIBIGL_LIBRARIES}  ${LIBIGLBOOLEAN_LIBRARY})
+
   FIND_LIBRARY( LIBIGLEMBREE_LIBRARY NAMES iglembree PATHS ${LIBIGL_LIB_DIRS})
   if(NOT LIBIGLEMBREE_LIBRARY)
     set(LIBIGL_FOUND FALSE)
+    message(FATAL_ERROR "could NOT find libiglembree")
   endif(NOT LIBIGLEMBREE_LIBRARY)
   set(LIBIGL_LIBRARIES ${LIBIGL_LIBRARIES}  ${LIBIGLEMBREE_LIBRARY})
   find_package(Embree REQUIRED)
@@ -77,12 +98,14 @@ if(LIBIGL_USE_STATIC_LIBRARY)
   FIND_LIBRARY( LIBIGLLIM_LIBRARY NAMES igllim PATHS ${LIBIGL_LIB_DIRS})
   if(NOT LIBIGLLIM_LIBRARY)
     set(LIBIGL_FOUND FALSE)
+    message(FATAL_ERROR "could NOT find libigllim")
   endif(NOT LIBIGLLIM_LIBRARY)
   set(LIBIGL_LIBRARIES ${LIBIGL_LIBRARIES}  ${LIBIGLLIM_LIBRARY})
 
   FIND_LIBRARY( LIBIGLMATLAB_LIBRARY NAMES iglmatlab PATHS ${LIBIGL_LIB_DIRS})
   if(NOT LIBIGLMATLAB_LIBRARY)
     set(LIBIGL_FOUND FALSE)
+    message(FATAL_ERROR "could NOT find libiglmatlab")
   endif(NOT LIBIGLMATLAB_LIBRARY)
   set(LIBIGL_LIBRARIES ${LIBIGL_LIBRARIES}  ${LIBIGLMATLAB_LIBRARY})
   find_package(Matlab REQUIRED)
@@ -91,29 +114,34 @@ if(LIBIGL_USE_STATIC_LIBRARY)
     set(LIBIGL_LIBRARIES ${LIBIGL_LIBRARIES}  ${MATLAB_LIBRARIES})
   else(MATLAB_FOUND)
     set(LIBIGL_FOUND FALSE)
+    message(FATAL_ERROR "could NOT find matlab")
   endif(MATLAB_FOUND)
 
   FIND_LIBRARY( LIBIGLSVD3X3_LIBRARY NAMES iglsvd3x3 PATHS ${LIBIGL_LIB_DIRS})
   if(NOT LIBIGLSVD3X3_LIBRARY)
     set(LIBIGL_FOUND FALSE)
+    message(FATAL_ERROR "could NOT find libiglsvd3x3")
   endif(NOT LIBIGLSVD3X3_LIBRARY)
   set(LIBIGL_LIBRARIES ${LIBIGL_LIBRARIES}  ${LIBIGLSVD3X3_LIBRARY})
 
   FIND_LIBRARY( LIBIGLTETGEN_LIBRARY NAMES igltetgen PATHS ${LIBIGL_LIB_DIRS})
   if(NOT LIBIGLTETGEN_LIBRARY)
     set(LIBIGL_FOUND FALSE)
+    message(FATAL_ERROR "could NOT find libigltetgen")
   endif(NOT LIBIGLTETGEN_LIBRARY)
   set(LIBIGL_LIBRARIES ${LIBIGL_LIBRARIES}  ${LIBIGLTETGEN_LIBRARY})
 
   FIND_LIBRARY( LIBIGLTRIANGLE_LIBRARY NAMES igltriangle PATHS ${LIBIGL_LIB_DIRS})
   if(NOT LIBIGLTRIANGLE_LIBRARY)
     set(LIBIGL_FOUND FALSE)
+    message(FATAL_ERROR "could NOT find libigltriangle")
   endif(NOT LIBIGLTRIANGLE_LIBRARY)
   set(LIBIGL_LIBRARIES ${LIBIGL_LIBRARIES}  ${LIBIGLTRIANGLE_LIBRARY})
 
   FIND_LIBRARY( LIBIGLVIEWER_LIBRARY NAMES iglviewer PATHS ${LIBIGL_LIB_DIRS})
   if(NOT LIBIGLVIEWER_LIBRARY)
     set(LIBIGL_FOUND FALSE)
+      message(FATAL_ERROR "could NOT find libiglviewer")
   endif(NOT LIBIGLVIEWER_LIBRARY)
   set(LIBIGL_LIBRARIES ${LIBIGL_LIBRARIES}  ${LIBIGLVIEWER_LIBRARY})
 

tutorial/cmake/FindMATLAB.cmake

@@ -137,7 +137,7 @@ ELSE(WIN32)
     IF((NOT DEFINED MATLAB_ROOT) OR ("${MATLAB_ROOT}" STREQUAL ""))
 
     # Search for a version of Matlab available, starting from the most modern one to older versions
-      FOREACH(MATVER "R2014b" "R2014a" "R2013b" "R2013a" "R2012b" "R2012a" "R2011b" "R2011a" "R2010b" "R2010a" "R2009b" "R2009a" "R2008b")
+      FOREACH(MATVER "R2015b" "R2015a" "R2014b" "R2014a" "R2014a" "R2013b" "R2013a" "R2012b" "R2012a" "R2011b" "R2011a" "R2010b" "R2010a" "R2009b" "R2009a" "R2008b")
         IF((NOT DEFINED MATLAB_ROOT) OR ("${MATLAB_ROOT}" STREQUAL ""))
           IF(EXISTS /Applications/MATLAB_${MATVER}.app)
             SET(MATLAB_ROOT /Applications/MATLAB_${MATVER}.app)
@@ -187,9 +187,14 @@ ELSE(WIN32)
 
 ENDIF(WIN32)
 
+if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU")
+  set(MATLAB_LIBSTDCPP "-lstdc++")
+endif ()
+
 
 # This is common to UNIX and Win32:
 SET(MATLAB_LIBRARIES
+  ${MATLAB_LIBSTDCPP}
   ${MATLAB_MAT_LIBRARY}
   ${MATLAB_MEX_LIBRARY}
   ${MATLAB_MX_LIBRARY}

tutorial/images/cheburashka-knight-boolean.jpg.REMOVED.git-id

@@ -0,0 +1 @@
+687b6cd2769b89898267be8a92eb47e1fe599c02

tutorial/tutorial.html.REMOVED.git-id

@@ -1 +1 @@
-141eb93f575045ab7a950e5674b54ec2fa87b197
+bb9c175cf9dfae345668e73bfd5220e21b28ae0e

tutorial/tutorial.md.REMOVED.git-id

@@ -1 +1 @@
-d0833fffecc624061e2d057e37e4e3e99cc11351
+a7312d6b030f9fcaa361629f047902a139d01d1e