mesh booleans with cgal or cork and a lot of template specialization

Former-commit-id: a0086ac516d28416f925f01e2db8e596f35e8ca6

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,14 @@
+#ifndef MESH_BOOLEAN_TYPE_H
+#define MESH_BOOLEAN_TYPE_H
+
+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,36 @@
+#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> >&);
+#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,185 @@
+#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>
+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)
+{
+  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>&)> 
+    empty_fun;
+  return mesh_boolean(VA,FA,VB,FB,type,empty_fun,VC,FC);
+}
+
+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,
+  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>&)> 
+    & resolve_fun,
+  Eigen::PlainObjectBase<DerivedVC > & VC,
+  Eigen::PlainObjectBase<DerivedFC > & FC)
+{
+  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;
+  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)
+  {
+    MatrixX3S SV;
+    MatrixX3I SF;
+    MatrixX2I SIF;
+    VectorXI SJ,SIM,UIM;
+    RemeshSelfIntersectionsParam params;
+    remesh_self_intersections(V,F,params,SV,SF,SIF,SJ,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;
+  if(resolve_fun)
+  {
+    resolve_fun(V,F,CV,CF);
+  }else
+  {
+    libigl_resolve(V,F,CV,CF);
+  }
+
+  if(type == MESH_BOOLEAN_TYPE_RESOLVE)
+  {
+    FC = CF;
+    VC = CV;
+    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(false);
+        }
+        break;
+      default:
+        assert(false && "Unknown type");
+        return;
+    }
+  }
+  const Index gm = vG.size();
+  MatrixX3I G(gm,3);
+  for(Index g = 0;g<gm;g++)
+  {
+    G.row(g) = Gflip[g] ? CF.row(vG[g]).reverse().eval() : CF.row(vG[g]);
+  }
+  // remove duplicates: cancel out in all cases? assertion in others?
+  {
+    MatrixXi oldG = G;
+    unique_simplices(oldG,G);
+  }
+  // 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&, MeshBooleanType const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -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,70 @@
+#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
+  //  
+  //  See also: self_intersect
+  //     
+  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>
+  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> &)> 
+      & resolve_fun,
+    Eigen::PlainObjectBase<DerivedVC > & VC,
+    Eigen::PlainObjectBase<DerivedFC > & FC);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "mesh_boolean.cpp"
+#endif
+
+#endif

include/igl/boolean/mesh_boolean_cork.cpp

@@ -0,0 +1,93 @@
+#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> >&);
+#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,37 @@
+#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&);
+#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/remesh_self_intersections.cpp

@@ -83,4 +83,5 @@ 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> >&);
 #endif

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/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,78 @@
+#include "peal_outer_hull_layers.h"
+#include "outer_hull.h"
+
+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/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/sort.cpp

@@ -218,5 +218,6 @@ 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> >&);
 #endif

include/igl/triangle_triangle_adjacency.cpp

@@ -183,7 +183,7 @@ template <
 # endif
 # pragma omp parallel for if (ne>IGL_OMP_MIN_VALUE)
   // Slightly better memory access than loop over E
-  for(Index f = 0;f<m;f++)
+  for(Index f = 0;f<(Index)m;f++)
   {
     for(Index c = 0;c<3;c++)
     {
@@ -209,4 +209,5 @@ template <
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
+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> > > > > >&);
 #endif

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

@@ -30,8 +30,13 @@ IGL_INLINE void igl::unique_edge_map(
   // This does help a little
   for_each(uE2E.begin(),uE2E.end(),[](vector<uE2EType > & v){v.reserve(2);});
   assert(EMAP.size() == ne);
-  for(uE2EType e = 0;e<ne;e++)
+  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_simplices.cpp

@@ -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