Merge branch 'master' of github.com:libigl/libigl into gh-pages

Former-commit-id: 49fc0b7088280cc68cffc92b07bac35ec053dbe9

README.md

@@ -181,7 +181,7 @@ BibTeX entry:
   title = {{libigl}: A simple {C++} geometry processing library},
   author = {Alec Jacobson and Daniele Panozzo and others},
   note = {http://libigl.github.io/libigl/},
-  year = {2015},
+  year = {2016},
 }
 ```
 
@@ -234,7 +234,7 @@ If you find bugs or have problems please use our [github issue tracking
 page](https://github.com/libigl/libigl/issues).
 
 ## Copyright
-2015 Alec Jacobson, Daniele Panozzo, Christian Schüller, Olga Diamanti, Qingnan
+2016 Alec Jacobson, Daniele Panozzo, Christian Schüller, Olga Diamanti, Qingnan
 Zhou, Nico Pietroni, Stefan Brugger, Kenshi Takayama, Wenzel Jakob, Nikolas De
 Giorgis, Luigi Rocca, Leonardo Sacht, Kevin Walliman, Olga Sorkine-Hornung, and others.
 

include/igl/AABB.h

@@ -657,7 +657,6 @@ igl::AABB<DerivedV,DIM>::squared_distance(
 {
   using namespace Eigen;
   using namespace std;
-  using namespace igl;
   Scalar sqr_d = min_sqr_d;
   //assert(DIM == 3 && "Code has only been tested for DIM == 3");
   assert((Ele.cols() == 3 || Ele.cols() == 2 || Ele.cols() == 1)
@@ -1026,7 +1025,6 @@ inline void igl::AABB<DerivedV,DIM>::leaf_squared_distance(
   RowVectorDIMS & c) const
 {
   using namespace Eigen;
-  using namespace igl;
   using namespace std;
 
   // Simplex size

include/igl/adjacency_matrix.cpp

@@ -11,13 +11,14 @@
 
 #include <vector>
 
-template <typename T>
+template <typename DerivedF, typename T>
 IGL_INLINE void igl::adjacency_matrix(
-  const Eigen::MatrixXi & F, 
+  const Eigen::PlainObjectBase<DerivedF> & F, 
   Eigen::SparseMatrix<T>& A)
 {
   using namespace std;
   using namespace Eigen;
+  typedef typename DerivedF::Scalar Index;
 
   typedef Triplet<T> IJV;
   vector<IJV > ijv;
@@ -29,14 +30,14 @@ IGL_INLINE void igl::adjacency_matrix(
     for(int j = 0;j<F.cols();j++)
     {
       // Get indices of edge: s --> d
-      int s = F(i,j);
-      int d = F(i,(j+1)%F.cols());
+      Index s = F(i,j);
+      Index d = F(i,(j+1)%F.cols());
       ijv.push_back(IJV(s,d,1));
       ijv.push_back(IJV(d,s,1));
     }
   }
 
-  const int n = F.maxCoeff()+1;
+  const Index n = F.maxCoeff()+1;
   A.resize(n,n);
   switch(F.cols())
   {
@@ -65,6 +66,4 @@ IGL_INLINE void igl::adjacency_matrix(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template void igl::adjacency_matrix<int>(Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::SparseMatrix<int, 0, int>&);
-template void igl::adjacency_matrix<double>(Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::SparseMatrix<double, 0, int>&);
 #endif

include/igl/adjacency_matrix.h

@@ -38,9 +38,9 @@ namespace igl
   //   U = A-Adiag;
   //
   // See also: edges, cotmatrix, diag
-  template <typename T>
+  template <typename DerivedF, typename T>
   IGL_INLINE void adjacency_matrix(
-    const Eigen::MatrixXi & F, 
+    const Eigen::PlainObjectBase<DerivedF> & F, 
     Eigen::SparseMatrix<T>& A);
 }
 

include/igl/ambient_occlusion.cpp

@@ -27,7 +27,6 @@ IGL_INLINE void igl::ambient_occlusion(
   Eigen::PlainObjectBase<DerivedS> & S)
 {
   using namespace Eigen;
-  using namespace igl;
   const int n = P.rows();
   // Resize output
   S.resize(n,1);

include/igl/arap.cpp

@@ -173,7 +173,6 @@ IGL_INLINE bool igl::arap_solve(
   ARAPData & data,
   Eigen::PlainObjectBase<DerivedU> & U)
 {
-  using namespace igl;
   using namespace Eigen;
   using namespace std;
   assert(data.b.size() == bc.rows());

include/igl/bbw/bbw.cpp

@@ -59,7 +59,6 @@ IGL_INLINE bool igl::bbw::bbw(
   Eigen::PlainObjectBase<DerivedW> & W
   )
 {
-  using namespace igl;
   using namespace std;
   using namespace Eigen;
 

include/igl/bfs_orient.cpp

@@ -17,7 +17,6 @@ IGL_INLINE void igl::bfs_orient(
   Eigen::PlainObjectBase<DerivedC> & C)
 {
   using namespace Eigen;
-  using namespace igl;
   using namespace std;
   SparseMatrix<int> A;
   orientable_patches(F,C,A);

include/igl/copyleft/boolean/CSGTree.h

@@ -98,7 +98,7 @@ namespace igl
             mesh_boolean(A.V(),A.F(),B.V(),B.F(),type,m_V,m_F,m_J);
             // reindex m_J
             std::for_each(m_J.data(),m_J.data()+m_J.size(),
-              [&](typename VectorJ::Scalar & j)
+              [&](typename VectorJ::Scalar & j) -> void
               {
                 if(j < A.F().rows())
                 {

include/igl/copyleft/boolean/mesh_boolean.cpp

@@ -20,6 +20,8 @@
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 #include <algorithm>
 
+//#define MESH_BOOLEAN_TIMING
+
 template <
   typename DerivedVA,
   typename DerivedFA,
@@ -43,6 +45,21 @@ IGL_INLINE void igl::copyleft::boolean::mesh_boolean(
     Eigen::PlainObjectBase<DerivedFC > & FC,
     Eigen::PlainObjectBase<DerivedJ > & J) {
 
+#ifdef MESH_BOOLEAN_TIMING
+  const auto & tictoc = []() -> double
+  {
+    static double t_start = igl::get_seconds();
+    double diff = igl::get_seconds()-t_start;
+    t_start += diff;
+    return diff;
+  };
+  const auto log_time = [&](const std::string& label) -> void {
+    std::cout << "mesh_boolean." << label << ": "
+      << tictoc() << std::endl;
+  };
+  tictoc();
+#endif
+
   typedef typename DerivedVC::Scalar Scalar;
   //typedef typename DerivedFC::Scalar Index;
   typedef CGAL::Epeck Kernel;
@@ -84,6 +101,9 @@ IGL_INLINE void igl::copyleft::boolean::mesh_boolean(
       //}
       resolve_fun(VV, FF, V, F, CJ);
   }
+#ifdef MESH_BOOLEAN_TIMING
+  log_time("resolve_self_intersection");
+#endif
 
   // Compute winding numbers on each side of each facet.
   const size_t num_faces = F.rows();
@@ -91,7 +111,11 @@ IGL_INLINE void igl::copyleft::boolean::mesh_boolean(
   Eigen::VectorXi labels(num_faces);
   std::transform(CJ.data(), CJ.data()+CJ.size(), labels.data(),
       [&](int i) { return i<FA.rows() ? 0:1; });
-  igl::copyleft::cgal::propagate_winding_numbers(V, F, labels, W);
+  if (num_faces > 0) {
+    igl::copyleft::cgal::propagate_winding_numbers(V, F, labels, W);
+  } else {
+    W.resize(0, 4);
+  }
   assert((size_t)W.rows() == num_faces);
   if (W.cols() == 2) {
     assert(FB.rows() == 0);
@@ -101,6 +125,9 @@ IGL_INLINE void igl::copyleft::boolean::mesh_boolean(
   } else {
     assert(W.cols() == 4);
   }
+#ifdef MESH_BOOLEAN_TIMING
+  log_time("propagate_input_winding_number");
+#endif
 
   // Compute resulting winding number.
   Eigen::MatrixXi Wr(num_faces, 2);
@@ -111,6 +138,9 @@ IGL_INLINE void igl::copyleft::boolean::mesh_boolean(
     Wr(i,0) = wind_num_op(w_out);
     Wr(i,1) = wind_num_op(w_in);
   }
+#ifdef MESH_BOOLEAN_TIMING
+  log_time("compute_output_winding_number");
+#endif
 
   // Extract boundary separating inside from outside.
   auto index_to_signed_index = [&](size_t i, bool ori) -> int{
@@ -141,6 +171,9 @@ IGL_INLINE void igl::copyleft::boolean::mesh_boolean(
     }
     kept_face_indices(i, 0) = CJ[idx];
   }
+#ifdef MESH_BOOLEAN_TIMING
+  log_time("extract_output");
+#endif
 
   // Finally, remove duplicated faces and unreferenced vertices.
   {
@@ -160,6 +193,9 @@ IGL_INLINE void igl::copyleft::boolean::mesh_boolean(
     Eigen::VectorXi newIM;
     igl::remove_unreferenced(Vs,G,VC,FC,newIM);
   }
+#ifdef MESH_BOOLEAN_TIMING
+  log_time("clean_up");
+#endif
 }
 
 template <

include/igl/copyleft/boolean/minkowski_sum.cpp

@@ -1,23 +1,159 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 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 "minkowski_sum.h"
 #include "mesh_boolean.h"
 
 #include "../../slice_mask.h"
 #include "../../unique.h"
+#include "../../get_seconds.h"
+#include "../../edges.h"
+#include "../cgal/assign_scalar.h"
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <cassert>
+#include <vector>
+template <
+  typename DerivedVA,
+  typename DerivedFA,
+  typename DerivedVB,
+  typename DerivedFB,
+  typename DerivedW,
+  typename DerivedG,
+  typename DerivedJ>
+IGL_INLINE void igl::copyleft::boolean::minkowski_sum(
+  const Eigen::PlainObjectBase<DerivedVA> & VA,
+  const Eigen::PlainObjectBase<DerivedFA> & FA,
+  const Eigen::PlainObjectBase<DerivedVB> & VB,
+  const Eigen::PlainObjectBase<DerivedFB> & FB,
+  const bool resolve_overlaps,
+  Eigen::PlainObjectBase<DerivedW> & W,
+  Eigen::PlainObjectBase<DerivedG> & G,
+  Eigen::PlainObjectBase<DerivedJ> & J)
+{
+  using namespace std;
+  using namespace Eigen;
+  assert(FA.cols() == 3 && "FA must contain a closed triangle mesh");
+  assert(FB.cols() <= FA.cols() && 
+    "FB must contain lower diemnsional simplices than FA");
+  const auto tictoc = []()->double
+  {
+    static double t_start;
+    double now = igl::get_seconds();
+    double interval = now-t_start;
+    t_start = now;
+    return interval;
+  };
+  tictoc();
+  Matrix<typename DerivedFB::Scalar,Dynamic,2> EB;
+  edges(FB,EB);
+  Matrix<typename DerivedFA::Scalar,Dynamic,2> EA(0,2);
+  if(FB.cols() == 3)
+  {
+    edges(FA,EA);
+  }
+  // number of copies of A along edges of B
+  const int n_ab = EB.rows();
+  // number of copies of B along edges of A
+  const int n_ba = EA.rows();
+
+  vector<DerivedW> vW(n_ab + n_ba);
+  vector<DerivedG> vG(n_ab + n_ba);
+  vector<DerivedJ> vJ(n_ab + n_ba);
+  vector<int> offsets(n_ab + n_ba + 1);
+  offsets[0] = 0;
+  // sweep A along edges of B
+  for(int e = 0;e<n_ab;e++)
+  {
+    Matrix<typename DerivedJ::Scalar,Dynamic,1> eJ;
+    minkowski_sum(
+      VA,
+      FA,
+      VB.row(EB(e,0)).eval(),
+      VB.row(EB(e,1)).eval(),
+      false,
+      vW[e],
+      vG[e],
+      eJ);
+    assert(vG[e].rows() == eJ.rows());
+    assert(eJ.cols() == 1);
+    vJ[e].resize(vG[e].rows(),2);
+    vJ[e].col(0) = eJ;
+    vJ[e].col(1).setConstant(e);
+    offsets[e+1] = offsets[e] + vW[e].rows();
+  }
+  // sweep B along edges of A
+  for(int e = 0;e<n_ba;e++)
+  {
+    Matrix<typename DerivedJ::Scalar,Dynamic,1> eJ;
+    const int ee = n_ab+e;
+    minkowski_sum(
+      VB,
+      FB,
+      VA.row(EA(e,0)).eval(),
+      VA.row(EA(e,1)).eval(),
+      false,
+      vW[ee],
+      vG[ee],
+      eJ);
+    vJ[ee].resize(vG[ee].rows(),2);
+    vJ[ee].col(0) = eJ.array() + (FA.rows()+1);
+    vJ[ee].col(1).setConstant(ee);
+  }
+  // Combine meshes
+  int n=0,m=0;
+  for_each(vW.begin(),vW.end(),[&n](const DerivedW & w){n+=w.rows();});
+  for_each(vG.begin(),vG.end(),[&m](const DerivedG & g){m+=g.rows();});
+  assert(n == offsets.back());
+
+  W.resize(n,3);
+  G.resize(m,3);
+  J.resize(m,2);
+  {
+    int m_off = 0,n_off = 0;
+    for(int i = 0;i<vG.size();i++)
+    {
+      W.block(n_off,0,vW[i].rows(),3) = vW[i];
+      G.block(m_off,0,vG[i].rows(),3) = vG[i].array()+offsets[i];
+      J.block(m_off,0,vJ[i].rows(),2) = vJ[i];
+      n_off += vW[i].rows();
+      m_off += vG[i].rows();
+    }
+    assert(n == n_off);
+    assert(m == m_off);
+  }
+  if(resolve_overlaps)
+  {
+    Eigen::Matrix<typename DerivedJ::Scalar, Eigen::Dynamic,1> SJ;
+    mesh_boolean(
+      DerivedW(W),
+      DerivedG(G),
+      Matrix<typename DerivedW::Scalar,Dynamic,Dynamic>(),
+      Matrix<typename DerivedG::Scalar,Dynamic,Dynamic>(),
+      MESH_BOOLEAN_TYPE_UNION,
+      W,
+      G,
+      SJ);
+    J = slice(DerivedJ(J),SJ,1);
+  }
+}
 
 template <
   typename DerivedVA,
   typename DerivedFA,
-  typename Deriveds,
-  typename Derivedd,
+  typename sType, int sCols, int sOptions,
+  typename dType, int dCols, int dOptions,
   typename DerivedW,
   typename DerivedG,
   typename DerivedJ>
 IGL_INLINE void igl::copyleft::boolean::minkowski_sum(
   const Eigen::PlainObjectBase<DerivedVA> & VA,
   const Eigen::PlainObjectBase<DerivedFA> & FA,
-  const Eigen::PlainObjectBase<Deriveds> & s,
-  const Eigen::PlainObjectBase<Derivedd> & d,
+  const Eigen::Matrix<sType,1,sCols,sOptions> & s,
+  const Eigen::Matrix<dType,1,dCols,dOptions> & d,
   const bool resolve_overlaps, 
   Eigen::PlainObjectBase<DerivedW> & W,
   Eigen::PlainObjectBase<DerivedG> & G,
@@ -25,6 +161,8 @@ IGL_INLINE void igl::copyleft::boolean::minkowski_sum(
 {
   using namespace Eigen;
   using namespace std;
+  assert(s.cols() == 3 && "s should be a 3d point");
+  assert(d.cols() == 3 && "d should be a 3d point");
   // silly base case
   if(FA.size() == 0)
   {
@@ -41,13 +179,13 @@ IGL_INLINE void igl::copyleft::boolean::minkowski_sum(
   // number of vertices
   const int n = VA.rows();
   // duplicate vertices at s and d, we'll remove unreferernced later
-  DerivedW WW(2*n,dim);
+  W.resize(2*n,dim);
   for(int i = 0;i<n;i++)
   {
     for(int j = 0;j<dim;j++)
     {
-      WW  (i,j) = VA(i,j) + s(j);
-      WW(i+n,j) = VA(i,j) + d(j);
+      W  (i,j) = VA(i,j) + s(j);
+      W(i+n,j) = VA(i,j) + d(j);
     }
   }
   // number of faces
@@ -189,8 +327,8 @@ IGL_INLINE void igl::copyleft::boolean::minkowski_sum(
     assert(nq == k);
   }
 
-  MatrixXI GG(GT.rows()+2*GQ.rows(),3);
-  GG<< 
+  G.resize(GT.rows()+2*GQ.rows(),3);
+  G<< 
     GT,
     GQ.col(0), GQ.col(1), GQ.col(2), 
     GQ.col(0), GQ.col(2), GQ.col(3);
@@ -198,9 +336,9 @@ IGL_INLINE void igl::copyleft::boolean::minkowski_sum(
   J<<JT,DerivedJ::Constant(2*GQ.rows(),1,2*m+1);
   if(resolve_overlaps)
   {
-    DerivedJ SJ;
+    Eigen::Matrix<typename DerivedJ::Scalar, Eigen::Dynamic,1> SJ;
     mesh_boolean(
-      WW,GG,
+      DerivedW(W),DerivedG(G),
       Matrix<typename DerivedVA::Scalar,Dynamic,Dynamic>(),MatrixXI(),
       MESH_BOOLEAN_TYPE_UNION,
       W,G,SJ);
@@ -211,16 +349,16 @@ IGL_INLINE void igl::copyleft::boolean::minkowski_sum(
 template <
   typename DerivedVA,
   typename DerivedFA,
-  typename Deriveds,
-  typename Derivedd,
+  typename sType, int sCols, int sOptions,
+  typename dType, int dCols, int dOptions,
   typename DerivedW,
   typename DerivedG,
   typename DerivedJ>
 IGL_INLINE void igl::copyleft::boolean::minkowski_sum(
   const Eigen::PlainObjectBase<DerivedVA> & VA,
   const Eigen::PlainObjectBase<DerivedFA> & FA,
-  const Eigen::PlainObjectBase<Deriveds> & s,
-  const Eigen::PlainObjectBase<Derivedd> & d,
+  const Eigen::Matrix<sType,1,sCols,sOptions> & s,
+  const Eigen::Matrix<dType,1,dCols,dOptions> & d,
   Eigen::PlainObjectBase<DerivedW> & W,
   Eigen::PlainObjectBase<DerivedG> & G,
   Eigen::PlainObjectBase<DerivedJ> & J)

include/igl/copyleft/boolean/minkowski_sum.h

@@ -1,3 +1,10 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
 #ifndef IGL_COPYLEFT_CGAL_MINKOWSKI_SUM_H
 #define IGL_COPYLEFT_CGAL_MINKOWSKI_SUM_H
 
@@ -11,6 +18,39 @@ namespace igl
     namespace boolean
     {
       // Compute the Minkowski sum of a closed triangle mesh (V,F) and a
+      // set of simplices in 3D.
+      //
+      // Inputs:
+      //   VA  #VA by 3 list of mesh vertices in 3D
+      //   FA  #FA by 3 list of triangle indices into VA
+      //   VB  #VB by 3 list of mesh vertices in 3D
+      //   FB  #FB by ss list of simplex indices into VB, ss<=3
+      //   resolve_overlaps  whether or not to resolve self-union. If false
+      //     then result may contain self-intersections if input mesh is
+      //     non-convex.
+      // Outputs:
+      //   W  #W by 3 list of mesh vertices in 3D
+      //   G  #G by 3 list of triangle indices into W
+      //   J  #G by 2 list of indices into 
+      //   
+      template <
+        typename DerivedVA,
+        typename DerivedFA,
+        typename DerivedVB,
+        typename DerivedFB,
+        typename DerivedW,
+        typename DerivedG,
+        typename DerivedJ>
+      IGL_INLINE void minkowski_sum(
+        const Eigen::PlainObjectBase<DerivedVA> & VA,
+        const Eigen::PlainObjectBase<DerivedFA> & FA,
+        const Eigen::PlainObjectBase<DerivedVB> & VB,
+        const Eigen::PlainObjectBase<DerivedFB> & FB,
+        const bool resolve_overlaps,
+        Eigen::PlainObjectBase<DerivedW> & W,
+        Eigen::PlainObjectBase<DerivedG> & G,
+        Eigen::PlainObjectBase<DerivedJ> & J);
+      // Compute the Minkowski sum of a closed triangle mesh (V,F) and a
       // segment [s,d] in 3D.
       //
       // Inputs:
@@ -29,16 +69,16 @@ namespace igl
       template <
         typename DerivedVA,
         typename DerivedFA,
-        typename Deriveds,
-        typename Derivedd,
+        typename sType, int sCols, int sOptions,
+        typename dType, int dCols, int dOptions,
         typename DerivedW,
         typename DerivedG,
         typename DerivedJ>
       IGL_INLINE void minkowski_sum(
         const Eigen::PlainObjectBase<DerivedVA> & VA,
         const Eigen::PlainObjectBase<DerivedFA> & FA,
-        const Eigen::PlainObjectBase<Deriveds> & s,
-        const Eigen::PlainObjectBase<Derivedd> & d,
+        const Eigen::Matrix<sType,1,sCols,sOptions> & s,
+        const Eigen::Matrix<dType,1,dCols,dOptions> & d,
         const bool resolve_overlaps,
         Eigen::PlainObjectBase<DerivedW> & W,
         Eigen::PlainObjectBase<DerivedG> & G,
@@ -46,20 +86,19 @@ namespace igl
       template <
         typename DerivedVA,
         typename DerivedFA,
-        typename Deriveds,
-        typename Derivedd,
+        typename sType, int sCols, int sOptions,
+        typename dType, int dCols, int dOptions,
         typename DerivedW,
         typename DerivedG,
         typename DerivedJ>
       IGL_INLINE void minkowski_sum(
         const Eigen::PlainObjectBase<DerivedVA> & VA,
         const Eigen::PlainObjectBase<DerivedFA> & FA,
-        const Eigen::PlainObjectBase<Deriveds> & s,
-        const Eigen::PlainObjectBase<Derivedd> & d,
+        const Eigen::Matrix<sType,1,sCols,sOptions> & s,
+        const Eigen::Matrix<dType,1,dCols,dOptions> & d,
         Eigen::PlainObjectBase<DerivedW> & W,
         Eigen::PlainObjectBase<DerivedG> & G,
         Eigen::PlainObjectBase<DerivedJ> & J);
-
     }
   }
 }

include/igl/copyleft/cgal/SelfIntersectMesh.h

@@ -15,6 +15,8 @@
 #include <list>
 #include <map>
 #include <vector>
+#include <thread>
+#include <mutex>
 
 //#define IGL_SELFINTERSECTMESH_DEBUG
 #ifndef IGL_FIRST_HIT_EXCEPTION
@@ -90,14 +92,14 @@ namespace igl
           // Number of self-intersecting triangle pairs
           typedef typename DerivedF::Index Index;
           Index count;
-          typedef std::vector<CGAL::Object> ObjectList;
+          typedef std::vector<std::pair<Index, CGAL::Object>> ObjectList;
           // Using a vector here makes this **not** output sensitive
           Triangles T;
           typedef std::vector<Index> IndexList;
           IndexList lIF;
           // #F-long list of faces with intersections mapping to the order in
           // which they were first found
-          std::map<Index,std::pair<Index,ObjectList> > offending;
+          std::map<Index,ObjectList> offending;
           // Make a short name for the edge map's key
           typedef std::pair<Index,Index> EMK;
           // Make a short name for the type stored at each edge, the edge map's
@@ -107,6 +109,8 @@ namespace igl
           typedef std::map<EMK,EMV> EdgeMap;
           // Maps edges of offending faces to all incident offending faces
           EdgeMap edge2faces;
+          std::vector<std::pair<const Box, const Box> > candidate_box_pairs;
+
         public:
           RemeshSelfIntersectionsParam params;
         public:
@@ -143,7 +147,7 @@ namespace igl
           //   A  triangle in 3D
           //   B  triangle in 3D
           //   fa  index of A in F (and key into offending)
-          //   fb  index of A in F (and key into offending)
+          //   fb  index of B in F (and key into offending)
           // Returns true only if A intersects B
           //
           inline bool intersect(
@@ -200,6 +204,7 @@ namespace igl
           //   a  box containing a triangle
           //   b  box containing a triangle
           inline void box_intersect(const Box& a, const Box& b);
+          inline void process_intersecting_boxes();
         private:
           // Compute 2D delaunay triangulation of a given 3d triangle and a list of
           // intersection objects (points,segments,triangles). CGAL uses an affine
@@ -219,6 +224,8 @@ namespace igl
             SelfIntersectMesh * SIM, 
             const Box &a, 
             const Box &b);
+        private:
+          std::mutex m_offending_lock;
       };
     }
   }
@@ -330,20 +337,24 @@ inline igl::copyleft::cgal::SelfIntersectMesh<
   using namespace Eigen;
 
 #ifdef IGL_SELFINTERSECTMESH_DEBUG
-  const auto & tictoc = []()
+  const auto & tictoc = []() -> double
   {
     static double t_start = igl::get_seconds();
     double diff = igl::get_seconds()-t_start;
     t_start += diff;
     return diff;
   };
+  const auto log_time = [&](const std::string& label) -> void{
+    std::cout << "SelfIntersectMesh." << label << ": "
+      << tictoc() << std::endl;
+  };
   tictoc();
 #endif
 
   // Compute and process self intersections
   mesh_to_cgal_triangle_list(V,F,T);
 #ifdef IGL_SELFINTERSECTMESH_DEBUG
-  cout<<"mesh_to_cgal_triangle_list: "<<tictoc()<<endl;
+  log_time("convert_to_triangle_list");
 #endif
   // 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
@@ -367,7 +378,7 @@ inline igl::copyleft::cgal::SelfIntersectMesh<
       std::placeholders::_1,
       std::placeholders::_2);
 #ifdef IGL_SELFINTERSECTMESH_DEBUG
-  cout<<"boxes and bind: "<<tictoc()<<endl;
+  log_time("box_and_bind");
 #endif
   // Run the self intersection algorithm with all defaults
   try{
@@ -381,8 +392,9 @@ inline igl::copyleft::cgal::SelfIntersectMesh<
     }
     // Otherwise just fall through
   }
+  process_intersecting_boxes();
 #ifdef IGL_SELFINTERSECTMESH_DEBUG
-  cout<<"box_self_intersection_d: "<<tictoc()<<endl;
+  log_time("box_intersection_d");
 #endif
 
   // Convert lIF to Eigen matrix
@@ -403,7 +415,7 @@ inline igl::copyleft::cgal::SelfIntersectMesh<
     }
   }
 #ifdef IGL_SELFINTERSECTMESH_DEBUG
-  cout<<"IF: "<<tictoc()<<endl;
+  log_time("store_intersecting_face_pairs");
 #endif
 
   if(params.detect_only)
@@ -414,7 +426,7 @@ inline igl::copyleft::cgal::SelfIntersectMesh<
   remesh_intersections(V,F,T,offending,edge2faces,VV,FF,J,IM);
 
 #ifdef IGL_SELFINTERSECTMESH_DEBUG
-  cout<<"remesh intersection: "<<tictoc()<<endl;
+  log_time("remesh_intersection");
 #endif
 
   // Q: Does this give the same result as TETGEN?
@@ -455,8 +467,7 @@ inline void igl::copyleft::cgal::SelfIntersectMesh<
   if(offending.count(f) == 0)
   {
     // first time marking, initialize with new id and empty list
-    const Index id = offending.size();
-    offending[f] = {id,{}};
+    offending[f] = {};
     for(Index e = 0; e<3;e++)
     {
       // append face to edge's list
@@ -488,6 +499,7 @@ inline void igl::copyleft::cgal::SelfIntersectMesh<
   const Index fa,
   const Index fb)
 {
+  std::lock_guard<std::mutex> guard(m_offending_lock);
   mark_offensive(fa);
   mark_offensive(fb);
   this->count++;
@@ -531,8 +543,9 @@ inline bool igl::copyleft::cgal::SelfIntersectMesh<
   {
     // Construct intersection
     CGAL::Object result = CGAL::intersection(A,B);
-    offending[fa].second.push_back(result);
-    offending[fb].second.push_back(result);
+    std::lock_guard<std::mutex> guard(m_offending_lock);
+    offending[fa].push_back({fb, result});
+    offending[fb].push_back({fa, result});
   }
   return true;
 }
@@ -630,8 +643,9 @@ inline bool igl::copyleft::cgal::SelfIntersectMesh<
         A.vertex(va),
         *p));
       count_intersection(fa,fb);
-      offending[fa].second.push_back(seg);
-      offending[fb].second.push_back(seg);
+      std::lock_guard<std::mutex> guard(m_offending_lock);
+      offending[fa].push_back({fb, seg});
+      offending[fb].push_back({fa, seg});
       return true;
     }else if(CGAL::object_cast<Segment_3 >(&result))
     {
@@ -688,6 +702,8 @@ inline bool igl::copyleft::cgal::SelfIntersectMesh<
 {
   using namespace std;
 
+  {
+    std::lock_guard<std::mutex> guard(m_offending_lock);
   // must be co-planar
   if(
     A.supporting_plane() != B.supporting_plane() &&
@@ -695,6 +711,7 @@ inline bool igl::copyleft::cgal::SelfIntersectMesh<
   {
     return false;
   }
+  }
   // Since A and B are non-degenerate the intersection must be a polygon
   // (triangle). Either
   //   - the vertex of A (B) opposite the shared edge of lies on B (A), or
@@ -777,8 +794,9 @@ inline bool igl::copyleft::cgal::SelfIntersectMesh<
       } else
       {
         // Triangle object
-        offending[fa].second.push_back(result);
-        offending[fb].second.push_back(result);
+        std::lock_guard<std::mutex> guard(m_offending_lock);
+        offending[fa].push_back({fb, result});
+        offending[fb].push_back({fa, result});
         //cerr<<REDRUM("Coplanar at: "<<fa<<" & "<<fb<<" (double shared).")<<endl;
         return true;
       }
@@ -825,125 +843,129 @@ inline void igl::copyleft::cgal::SelfIntersectMesh<
   const Box& a, 
   const Box& b)
 {
-  using namespace std;
-  // Could we write this as a static function of:
-  //
-  // F.row(fa)
-  // F.row(fb)
-  // A
-  // B
-
-  // index in F and T
-  Index fa = a.handle()-T.begin();
-  Index fb = b.handle()-T.begin();
-  const Triangle_3 & A = *a.handle();
-  const Triangle_3 & B = *b.handle();
-  //// I'm not going to deal with degenerate triangles, though at some point we
-  //// should
-  //assert(!a.handle()->is_degenerate());
-  //assert(!b.handle()->is_degenerate());
-  // Number of combinatorially shared vertices
-  Index comb_shared_vertices = 0;
-  // Number of geometrically shared vertices (*not* including combinatorially
-  // shared)
-  Index geo_shared_vertices = 0;
-  // Keep track of shared vertex indices
-  std::vector<std::pair<Index,Index> > shared;
-  Index ea,eb;
-  for(ea=0;ea<3;ea++)
-  {
-    for(eb=0;eb<3;eb++)
-    {
-      if(F(fa,ea) == F(fb,eb))
+  candidate_box_pairs.push_back({a, b});
+}
+
+template <
+  typename Kernel,
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedVV,
+  typename DerivedFF,
+  typename DerivedIF,
+  typename DerivedJ,
+  typename DerivedIM>
+inline void igl::copyleft::cgal::SelfIntersectMesh<
+  Kernel,
+  DerivedV,
+  DerivedF,
+  DerivedVV,
+  DerivedFF,
+  DerivedIF,
+  DerivedJ,
+  DerivedIM>::process_intersecting_boxes()
+{
+  CGAL::Identity_transformation I;
+  std::vector<std::mutex> triangle_locks(T.size());
+  auto process_chunk = [&](size_t first, size_t last) -> void{
+    assert(last >= first);
+
+    for (size_t i=first; i<last; i++) {
+      const auto& box_pair = candidate_box_pairs[i];
+      const auto& a = box_pair.first;
+      const auto& b = box_pair.second;
+      Index fa = a.handle()-T.begin();
+      Index fb = b.handle()-T.begin();
+
+      std::lock_guard<std::mutex> guard_A(triangle_locks[fa]);
+      std::lock_guard<std::mutex> guard_B(triangle_locks[fb]);
+      const Triangle_3& A = *a.handle();
+      const Triangle_3& B = *b.handle();
+
+      // Number of combinatorially shared vertices
+      Index comb_shared_vertices = 0;
+      // Number of geometrically shared vertices (*not* including combinatorially
+      // shared)
+      Index geo_shared_vertices = 0;
+      // Keep track of shared vertex indices
+      std::vector<std::pair<Index,Index> > shared;
+      Index ea,eb;
+      for(ea=0;ea<3;ea++)
+      {
+        for(eb=0;eb<3;eb++)
+        {
+          if(F(fa,ea) == F(fb,eb))
+          {
+            comb_shared_vertices++;
+            shared.emplace_back(ea,eb);
+          }else if(A.vertex(ea) == B.vertex(eb))
+          {
+            geo_shared_vertices++;
+            shared.emplace_back(ea,eb);
+          }
+        }
+      }
+      const Index total_shared_vertices = comb_shared_vertices + geo_shared_vertices;
+
+      if(comb_shared_vertices== 3)
+      {
+        assert(shared.size() == 3);
+        //// Combinatorially duplicate face, these should be removed by preprocessing
+        //cerr<<REDRUM("Facets "<<fa<<" and "<<fb<<" are combinatorial duplicates")<<endl;
+        continue;
+      }
+      if(total_shared_vertices== 3)
+      {
+        assert(shared.size() == 3);
+        //// Geometrically duplicate face, these should be removed by preprocessing
+        //cerr<<REDRUM("Facets "<<fa<<" and "<<fb<<" are geometrical duplicates")<<endl;
+        continue;
+      }
+      if(total_shared_vertices == 2)
+      {
+        assert(shared.size() == 2);
+        // Q: What about coplanar?
+        //
+        // o    o
+        // |\  /|
+        // | \/ |
+        // | /\ |
+        // |/  \|
+        // o----o
+        double_shared_vertex(A,B,fa,fb,shared);
+        continue;
+      }
+      assert(total_shared_vertices<=1);
+      if(total_shared_vertices==1)
       {
-        comb_shared_vertices++;
-        shared.emplace_back(ea,eb);
-      }else if(A.vertex(ea) == B.vertex(eb))
+        single_shared_vertex(A,B,fa,fb,shared[0].first,shared[0].second);
+      }else
       {
-        geo_shared_vertices++;
-        shared.emplace_back(ea,eb);
+        intersect(*a.handle(),*b.handle(),fa,fb);
       }
     }
+  };
+  size_t num_threads=0;
+  const size_t hardware_limit = std::thread::hardware_concurrency();
+  if (const char* igl_num_threads = std::getenv("LIBIGL_NUM_THREADS")) {
+    num_threads = atoi(igl_num_threads);
   }
-  const Index total_shared_vertices = comb_shared_vertices + geo_shared_vertices;
-  if(comb_shared_vertices== 3)
-  {
-    assert(shared.size() == 3);
-    //// Combinatorially duplicate face, these should be removed by preprocessing
-    //cerr<<REDRUM("Facets "<<fa<<" and "<<fb<<" are combinatorial duplicates")<<endl;
-    goto done;
-  }
-  if(total_shared_vertices== 3)
-  {
-    assert(shared.size() == 3);
-    //// Geometrically duplicate face, these should be removed by preprocessing
-    //cerr<<REDRUM("Facets "<<fa<<" and "<<fb<<" are geometrical duplicates")<<endl;
-    goto done;
+  if (num_threads == 0 || num_threads > hardware_limit) {
+    num_threads = hardware_limit;
   }
-  //// SPECIAL CASES ARE BROKEN FOR COPLANAR TRIANGLES
-  //if(total_shared_vertices > 0)
-  //{
-  //  bool coplanar = 
-  //    CGAL::coplanar(A.vertex(0),A.vertex(1),A.vertex(2),B.vertex(0)) &&
-  //    CGAL::coplanar(A.vertex(0),A.vertex(1),A.vertex(2),B.vertex(1)) &&
-  //    CGAL::coplanar(A.vertex(0),A.vertex(1),A.vertex(2),B.vertex(2));
-  //  if(coplanar)
-  //  {
-  //    cerr<<MAGENTAGIN("Facets "<<fa<<" and "<<fb<<
-  //      " are coplanar and share vertices")<<endl;
-  //    goto full;
-  //  }
-  //}
-
-  if(total_shared_vertices == 2)
-  {
-    assert(shared.size() == 2);
-    // Q: What about coplanar?
-    //
-    // o    o
-    // |\  /|
-    // | \/ |
-    // | /\ |
-    // |/  \|
-    // o----o
-    double_shared_vertex(A,B,fa,fb,shared);
-
-    goto done;
+  assert(num_threads > 0);
+  const size_t num_pairs = candidate_box_pairs.size();
+  const size_t chunk_size = num_pairs / num_threads;
+  std::vector<std::thread> threads;
+  for (size_t i=0; i<num_threads-1; i++) {
+    threads.emplace_back(process_chunk, i*chunk_size, (i+1)*chunk_size);
   }
-  assert(total_shared_vertices<=1);
-  if(total_shared_vertices==1)
-  {
-//#ifndef NDEBUG
-//    CGAL::Object result =
-//#endif
-    single_shared_vertex(A,B,fa,fb,shared[0].first,shared[0].second);
-//#ifndef NDEBUG
-//    if(!CGAL::object_cast<Segment_3 >(&result))
-//    {
-//      const Point_3 * p = CGAL::object_cast<Point_3 >(&result);
-//      assert(p);
-//      for(int ea=0;ea<3;ea++)
-//      {
-//        for(int eb=0;eb<3;eb++)
-//        {
-//          if(F(fa,ea) == F(fb,eb))
-//          {
-//            assert(*p==A.vertex(ea));
-//            assert(*p==B.vertex(eb));
-//          }
-//        }
-//      }
-//    }
-//#endif
-  }else
-  {
-//full:
-    // No geometrically shared vertices, do general intersect
-    intersect(*a.handle(),*b.handle(),fa,fb);
+  // Do some work in the master thread.
+  process_chunk((num_threads-1)*chunk_size, num_pairs);
+  for (auto& t : threads) {
+    if (t.joinable()) t.join();
   }
-done:
-  return;
+  //process_chunk(0, candidate_box_pairs.size());
 }
 
 #endif
-

include/igl/copyleft/cgal/assign_scalar.cpp

@@ -18,7 +18,13 @@ IGL_INLINE void igl::copyleft::cgal::assign_scalar(
   const typename CGAL::Epeck::FT & cgal,
   double & d)
 {
-  d = CGAL::to_double(cgal);
+  const auto interval = CGAL::to_interval(cgal);
+  d = interval.first;
+  do {
+      const double next = nextafter(d, interval.second);
+      if (CGAL::abs(cgal-d) < CGAL::abs(cgal-next)) break;
+      d = next;
+  } while (d < interval.second);
 }
 
 IGL_INLINE void igl::copyleft::cgal::assign_scalar(

include/igl/copyleft/cgal/closest_facet.cpp

@@ -77,7 +77,7 @@ IGL_INLINE void igl::copyleft::cgal::closest_facet(
   Tree tree(triangles.begin(), triangles.end());
   tree.accelerate_distance_queries();
 
-  auto on_the_positive_side = [&](size_t fid, const Point_3& p) {
+  auto on_the_positive_side = [&](size_t fid, const Point_3& p) -> bool{
     const auto& f = F.row(fid).eval();
     Point_3 v0(V(f[0], 0), V(f[0], 1), V(f[0], 2));
     Point_3 v1(V(f[1], 0), V(f[1], 1), V(f[1], 2));
@@ -122,7 +122,7 @@ IGL_INLINE void igl::copyleft::cgal::closest_facet(
 
   enum ElementType { VERTEX, EDGE, FACE };
   auto determine_element_type = [&](const Point_3& p, const size_t fid,
-      size_t& element_index) {
+      size_t& element_index) -> ElementType {
     const auto& tri = triangles[fid];
     const Point_3 p0 = tri[0];
     const Point_3 p1 = tri[1];
@@ -143,7 +143,7 @@ IGL_INLINE void igl::copyleft::cgal::closest_facet(
       size_t query_idx,
       const size_t s, const size_t d,
       size_t preferred_facet,
-      bool& orientation) {
+      bool& orientation) -> size_t {
     Point_3 query_point(
         P(query_idx, 0),
         P(query_idx, 1),
@@ -221,7 +221,7 @@ IGL_INLINE void igl::copyleft::cgal::closest_facet(
 
   auto process_face_case = [&](
       const size_t query_idx, const Point_3& closest_point,
-      const size_t fid, bool& orientation) {
+      const size_t fid, bool& orientation) -> size_t {
     const auto& f = F.row(I(fid, 0));
     return process_edge_case(query_idx, f[0], f[1], I(fid, 0), orientation);
   };
@@ -244,7 +244,7 @@ IGL_INLINE void igl::copyleft::cgal::closest_facet(
     const size_t query_idx, 
     size_t s,
     size_t preferred_facet, 
-    bool& orientation) 
+    bool& orientation) -> size_t
   {
     const Point_3 query_point(
         P(query_idx, 0), P(query_idx, 1), P(query_idx, 2));
@@ -296,7 +296,7 @@ IGL_INLINE void igl::copyleft::cgal::closest_facet(
 
     // A plane is on the exterior if all adj_points lies on or to
     // one side of the plane.
-    auto is_on_exterior = [&](const Plane_3& separator) {
+    auto is_on_exterior = [&](const Plane_3& separator) -> bool{
       size_t positive=0;
       size_t negative=0;
       size_t coplanar=0;

include/igl/copyleft/cgal/complex_to_mesh.cpp

@@ -21,7 +21,6 @@ IGL_INLINE bool igl::copyleft::cgal::complex_to_mesh(
   Eigen::PlainObjectBase<DerivedV> & V, 
   Eigen::PlainObjectBase<DerivedF> & F)
 {
-  using namespace igl;
   using namespace Eigen;
   // CGAL/IO/Complex_2_in_triangulation_3_file_writer.h
   using CGAL::Surface_mesher::number_of_facets_on_surface;

include/igl/copyleft/cgal/extract_cells.cpp

@@ -22,292 +22,228 @@
 //#define EXTRACT_CELLS_DEBUG
 
 template<
-typename DerivedV,
-typename DerivedF,
-typename DerivedP,
-typename DeriveduE,
-typename uE2EType,
-typename DerivedEMAP,
-typename DerivedC >
-IGL_INLINE size_t igl::copyleft::cgal::extract_cells_single_component(
-        const Eigen::PlainObjectBase<DerivedV>& V,
-        const Eigen::PlainObjectBase<DerivedF>& F,
-        const Eigen::PlainObjectBase<DerivedP>& P,
-        const Eigen::PlainObjectBase<DeriveduE>& uE,
-        const std::vector<std::vector<uE2EType> >& uE2E,
-        const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
-        Eigen::PlainObjectBase<DerivedC>& cells) {
-    //typedef typename DerivedF::Scalar Index;
-    const size_t num_faces = F.rows();
-    auto edge_index_to_face_index = [&](size_t index) {
-        return index % num_faces;
-    };
-    auto is_consistent = [&](const size_t fid, const size_t s, const size_t d) {
-        if ((size_t)F(fid, 0) == s && (size_t)F(fid, 1) == d) return false;
-        if ((size_t)F(fid, 1) == s && (size_t)F(fid, 2) == d) return false;
-        if ((size_t)F(fid, 2) == s && (size_t)F(fid, 0) == d) return false;
-
-        if ((size_t)F(fid, 0) == d && (size_t)F(fid, 1) == s) return true;
-        if ((size_t)F(fid, 1) == d && (size_t)F(fid, 2) == s) return true;
-        if ((size_t)F(fid, 2) == d && (size_t)F(fid, 0) == s) return true;
-        throw "Invalid face!";
-        return false;
-    };
-
-    const size_t num_unique_edges = uE.rows();
-    const size_t num_patches = P.maxCoeff() + 1;
-
-    std::vector<std::vector<size_t> > patch_edge_adj(num_patches);
-    std::vector<Eigen::VectorXi> orders(num_unique_edges);
-    std::vector<std::vector<bool> > orientations(num_unique_edges);
-    for (size_t i=0; i<num_unique_edges; i++) {
-        const size_t s = uE(i,0);
-        const size_t d = uE(i,1);
-        const auto adj_faces = uE2E[i];
-        if (adj_faces.size() > 2) {
-            std::vector<int> signed_adj_faces;
-            for (auto ei : adj_faces) {
-                const size_t fid = edge_index_to_face_index(ei);
-                bool cons = is_consistent(fid, s, d);
-                signed_adj_faces.push_back((fid+1)*(cons ? 1:-1));
-            }
-            auto& order = orders[i];
-            igl::copyleft::cgal::order_facets_around_edge(
-                    V, F, s, d, signed_adj_faces, order);
-            auto& orientation = orientations[i];
-            orientation.resize(order.size());
-            std::transform(order.data(), order.data() + order.size(),
-                    orientation.begin(), [&](int index) { return
-                    signed_adj_faces[index] > 0; });
-            std::transform(order.data(), order.data() + order.size(),
-                    order.data(), [&](int index) { return adj_faces[index]; });
-
-            for (auto ei : adj_faces) {
-                const size_t fid = edge_index_to_face_index(ei);
-                patch_edge_adj[P[fid]].push_back(ei);
-            }
-        }
-    }
-
-    const int INVALID = std::numeric_limits<int>::max();
-    cells.resize(num_patches, 2);
-    cells.setConstant(INVALID);
-
-    auto peel_cell_bd = [&](
-            size_t seed_patch_id,
-            short seed_patch_side,
-            size_t cell_idx) {
-        typedef std::pair<size_t, short> PatchSide;
-        std::queue<PatchSide> Q;
-        Q.emplace(seed_patch_id, seed_patch_side);
-        cells(seed_patch_id, seed_patch_side) = cell_idx;
-        while (!Q.empty()) {
-            auto entry = Q.front();
-            Q.pop();
-            const size_t patch_id = entry.first;
-            const short side = entry.second;
-            const auto& adj_edges = patch_edge_adj[patch_id];
-            for (const auto& ei : adj_edges) {
-                const size_t uei = EMAP[ei];
-                const auto& order = orders[uei];
-                const auto& orientation = orientations[uei];
-                const size_t edge_valance = order.size();
-                size_t curr_i = 0;
-                for (curr_i=0; curr_i < edge_valance; curr_i++) {
-                    if ((size_t)order[curr_i] == ei) break;
-                }
-                assert(curr_i < edge_valance);
-
-                const bool cons = orientation[curr_i];
-                size_t next;
-                if (side == 0) {
-                    next = (cons ? (curr_i + 1) :
-                            (curr_i + edge_valance - 1)) % edge_valance;
-                } else {
-                    next = (cons ? (curr_i + edge_valance - 1) :
-                            (curr_i + 1)) % edge_valance;
-                }
-                const size_t next_ei = order[next];
-                const bool next_cons = orientation[next];
-                const size_t next_patch_id = P[next_ei % num_faces];
-                const short next_patch_side = (next_cons != cons) ?
-                    side:abs(side-1);
-                if (cells(next_patch_id, next_patch_side) == INVALID) {
-                    Q.emplace(next_patch_id, next_patch_side);
-                    cells(next_patch_id, next_patch_side) = cell_idx;
-                } else {
-                    assert(
-                      (size_t)cells(next_patch_id, next_patch_side) == 
-                      cell_idx);
-                }
-            }
-        }
-    };
-
-    size_t count=0;
-    for (size_t i=0; i<num_patches; i++) {
-        if (cells(i, 0) == INVALID) {
-            peel_cell_bd(i, 0, count);
-            count++;
-        }
-        if (cells(i, 1) == INVALID) {
-            peel_cell_bd(i, 1, count);
-            count++;
-        }
-    }
-    return count;
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedC >
+IGL_INLINE size_t igl::copyleft::cgal::extract_cells(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedC>& cells)
+{
+  const size_t num_faces = F.rows();
+  // Construct edge adjacency
+  Eigen::MatrixXi E, uE;
+  Eigen::VectorXi EMAP;
+  std::vector<std::vector<size_t> > uE2E;
+  igl::unique_edge_map(F, E, uE, EMAP, uE2E);
+  // Cluster into manifold patches
+  Eigen::VectorXi P;
+  igl::extract_manifold_patches(F, EMAP, uE2E, P);
+  // Extract cells
+  DerivedC per_patch_cells;
+  const size_t num_cells =
+    igl::copyleft::cgal::extract_cells(V,F,P,E,uE,uE2E,EMAP,per_patch_cells);
+  // Distribute per-patch cell information to each face
+  cells.resize(num_faces, 2);
+  for (size_t i=0; i<num_faces; i++) 
+  {
+    cells.row(i) = per_patch_cells.row(P[i]);
+  }
+  return num_cells;
 }
 
+
 template<
-    typename DerivedV,
-    typename DerivedF,
-    typename DerivedP,
-    typename DerivedE,
-    typename DeriveduE,
-    typename uE2EType,
-    typename DerivedEMAP,
-    typename DerivedC >
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedP,
+  typename DerivedE,
+  typename DeriveduE,
+  typename uE2EType,
+  typename DerivedEMAP,
+  typename DerivedC >
 IGL_INLINE size_t igl::copyleft::cgal::extract_cells(
-        const Eigen::PlainObjectBase<DerivedV>& V,
-        const Eigen::PlainObjectBase<DerivedF>& F,
-        const Eigen::PlainObjectBase<DerivedP>& P,
-        const Eigen::PlainObjectBase<DerivedE>& E,
-        const Eigen::PlainObjectBase<DeriveduE>& uE,
-        const std::vector<std::vector<uE2EType> >& uE2E,
-        const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
-        Eigen::PlainObjectBase<DerivedC>& cells) {
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  const Eigen::PlainObjectBase<DerivedP>& P,
+  const Eigen::PlainObjectBase<DerivedE>& E,
+  const Eigen::PlainObjectBase<DeriveduE>& uE,
+  const std::vector<std::vector<uE2EType> >& uE2E,
+  const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+  Eigen::PlainObjectBase<DerivedC>& cells) 
+{
 #ifdef EXTRACT_CELLS_DEBUG
-  const auto & tictoc = []()
+  const auto & tictoc = []() -> double
   {
     static double t_start = igl::get_seconds();
     double diff = igl::get_seconds()-t_start;
     t_start += diff;
     return diff;
   };
+  const auto log_time = [&](const std::string& label) -> void {
+    std::cout << "extract_cells." << label << ": "
+      << tictoc() << std::endl;
+  };
   tictoc();
+#else
+  // no-op
+  const auto log_time = [](const std::string){};
 #endif
-    const size_t num_faces = F.rows();
-    typedef typename DerivedF::Scalar Index;
-    const size_t num_patches = P.maxCoeff()+1;
+  const size_t num_faces = F.rows();
+  typedef typename DerivedF::Scalar Index;
+  const size_t num_patches = P.maxCoeff()+1;
 
-    DerivedC raw_cells;
-    const size_t num_raw_cells =
-        igl::copyleft::cgal::extract_cells_single_component(
-                V, F, P, uE, uE2E, EMAP, raw_cells);
-#ifdef EXTRACT_CELLS_DEBUG
-    std::cout << "Extract single component cells: " << tictoc() << std::endl;
-#endif
+  // Extract all cells...
+  DerivedC raw_cells;
+  const size_t num_raw_cells =
+    extract_cells_single_component(V,F,P,uE,uE2E,EMAP,raw_cells);
+  log_time("extract_single_component_cells");
 
-    std::vector<std::vector<std::vector<Index > > > TT,_1;
-    igl::triangle_triangle_adjacency(E, EMAP, uE2E, false, TT, _1);
-#ifdef EXTRACT_CELLS_DEBUG
-    std::cout << "face adj: " << tictoc() << std::endl;
-#endif
+  // Compute triangle-triangle adjacency data-structure
+  std::vector<std::vector<std::vector<Index > > > TT,_1;
+  igl::triangle_triangle_adjacency(E, EMAP, uE2E, false, TT, _1);
+  log_time("compute_face_adjacency");
 
-    Eigen::VectorXi C, counts;
-    igl::facet_components(TT, C, counts);
-#ifdef EXTRACT_CELLS_DEBUG
-    std::cout << "face comp: " << tictoc() << std::endl;
-#endif
+  // Compute connected components of the mesh
+  Eigen::VectorXi C, counts;
+  igl::facet_components(TT, C, counts);
+  log_time("form_components");
 
-    const size_t num_components = counts.size();
-    std::vector<std::vector<size_t> > components(num_components);
-    for (size_t i=0; i<num_faces; i++) {
-        components[C[i]].push_back(i);
-    }
-    std::vector<Eigen::VectorXi> Is(num_components);
-    for (size_t i=0; i<num_components; i++) {
-        Is[i].resize(components[i].size());
-        std::copy(components[i].begin(), components[i].end(),
-                Is[i].data());
-    }
+  const size_t num_components = counts.size();
+  // components[c] --> list of face indices into F of faces in component c
+  std::vector<std::vector<size_t> > components(num_components);
+  // Loop over all faces
+  for (size_t i=0; i<num_faces; i++) 
+  {
+    components[C[i]].push_back(i);
+  }
+  // Convert vector lists to Eigen lists...
+  std::vector<Eigen::VectorXi> Is(num_components);
+  for (size_t i=0; i<num_components; i++)
+  {
+    Is[i].resize(components[i].size());
+    std::copy(components[i].begin(), components[i].end(),Is[i].data());
+  }
 
-    Eigen::VectorXi outer_facets(num_components);
-    Eigen::VectorXi outer_facet_orientation(num_components);
-    Eigen::VectorXi outer_cells(num_components);
-    for (size_t i=0; i<num_components; i++) {
-        bool flipped;
-        igl::copyleft::cgal::outer_facet(V, F, Is[i], outer_facets[i], flipped);
-        outer_facet_orientation[i] = flipped?1:0;
-        outer_cells[i] = raw_cells(P[outer_facets[i]], outer_facet_orientation[i]);
-    }
+  // Find outer facets, their orientations and cells for each component
+  Eigen::VectorXi outer_facets(num_components);
+  Eigen::VectorXi outer_facet_orientation(num_components);
+  Eigen::VectorXi outer_cells(num_components);
+  for (size_t i=0; i<num_components; i++)
+  {
+    bool flipped;
+    igl::copyleft::cgal::outer_facet(V, F, Is[i], outer_facets[i], flipped);
+    outer_facet_orientation[i] = flipped?1:0;
+    outer_cells[i] = raw_cells(P[outer_facets[i]], outer_facet_orientation[i]);
+  }
 #ifdef EXTRACT_CELLS_DEBUG
-    std::cout << "Per comp outer facet: " << tictoc() << std::endl;
+  log_time("outer_facet_per_component");
 #endif
 
-    auto get_triangle_center = [&](const size_t fid) {
-        return ((V.row(F(fid, 0)) + V.row(F(fid, 1)) + V.row(F(fid, 2)))
-                /3.0).eval();
-    };
-    std::vector<std::vector<size_t> > nested_cells(num_raw_cells);
-    std::vector<std::vector<size_t> > ambient_cells(num_raw_cells);
-    std::vector<std::vector<size_t> > ambient_comps(num_components);
-    if (num_components > 1) {
-        DerivedV bbox_min(num_components, 3);
-        DerivedV bbox_max(num_components, 3);
-        bbox_min.rowwise() = V.colwise().maxCoeff().eval();
-        bbox_max.rowwise() = V.colwise().minCoeff().eval();
-        for (size_t i=0; i<num_faces; i++) {
-          const auto comp_id = C[i];
-          const auto& f = F.row(i);
-          for (size_t j=0; j<3; j++) {
-            bbox_min(comp_id, 0) = std::min(bbox_min(comp_id, 0), V(f[j], 0));
-            bbox_min(comp_id, 1) = std::min(bbox_min(comp_id, 1), V(f[j], 1));
-            bbox_min(comp_id, 2) = std::min(bbox_min(comp_id, 2), V(f[j], 2));
-            bbox_max(comp_id, 0) = std::max(bbox_max(comp_id, 0), V(f[j], 0));
-            bbox_max(comp_id, 1) = std::max(bbox_max(comp_id, 1), V(f[j], 1));
-            bbox_max(comp_id, 2) = std::max(bbox_max(comp_id, 2), V(f[j], 2));
-          }
+  // Compute barycenter of a triangle in mesh (V,F)
+  //
+  // Inputs:
+  //   fid  index into F
+  // Returns row-vector of barycenter coordinates
+  const auto get_triangle_center = [&V,&F](const size_t fid) 
+  {
+    return ((V.row(F(fid,0))+V.row(F(fid,1))+V.row(F(fid,2)))/3.0).eval();
+  };
+  std::vector<std::vector<size_t> > nested_cells(num_raw_cells);
+  std::vector<std::vector<size_t> > ambient_cells(num_raw_cells);
+  std::vector<std::vector<size_t> > ambient_comps(num_components);
+  // Only bother if there's more than one component
+  if(num_components > 1) 
+  {
+    // construct bounding boxes for each component
+    DerivedV bbox_min(num_components, 3);
+    DerivedV bbox_max(num_components, 3);
+    // Why not just initialize to numeric_limits::min, numeric_limits::max?
+    bbox_min.rowwise() = V.colwise().maxCoeff().eval();
+    bbox_max.rowwise() = V.colwise().minCoeff().eval();
+    // Loop over faces
+    for (size_t i=0; i<num_faces; i++)
+    {
+      // component of this face
+      const auto comp_id = C[i];
+      const auto& f = F.row(i);
+      for (size_t j=0; j<3; j++) 
+      {
+        for(size_t d=0;d<3;d++)
+        {
+          bbox_min(comp_id,d) = std::min(bbox_min(comp_id,d), V(f[j],d));
+          bbox_max(comp_id,d) = std::max(bbox_max(comp_id,d), V(f[j],d));
         }
-        auto bbox_intersects = [&](size_t comp_i, size_t comp_j) {
-          return !(
-              bbox_max(comp_i,0) < bbox_min(comp_j,0) ||
-              bbox_max(comp_i,1) < bbox_min(comp_j,1) ||
-              bbox_max(comp_i,2) < bbox_min(comp_j,2) ||
-              bbox_max(comp_j,0) < bbox_min(comp_i,0) ||
-              bbox_max(comp_j,1) < bbox_min(comp_i,1) ||
-              bbox_max(comp_j,2) < bbox_min(comp_i,2));
-        };
-
-        for (size_t i=0; i<num_components; i++) {
-            std::vector<size_t> candidate_comps;
-            candidate_comps.reserve(num_components);
-            for (size_t j=0; j<num_components; j++) {
-                if (i == j) continue;
-                if (bbox_intersects(i,j)) candidate_comps.push_back(j);
-            }
-
-            const size_t num_candidate_comps = candidate_comps.size();
-            if (num_candidate_comps == 0) continue;
+      }
+    }
+    // Return true if box of component ci intersects that of cj
+    const auto bbox_intersects = [&bbox_max,&bbox_min](size_t ci, size_t cj)
+    {
+      return !(
+        bbox_max(ci,0) < bbox_min(cj,0) ||
+        bbox_max(ci,1) < bbox_min(cj,1) ||
+        bbox_max(ci,2) < bbox_min(cj,2) ||
+        bbox_max(cj,0) < bbox_min(ci,0) ||
+        bbox_max(cj,1) < bbox_min(ci,1) ||
+        bbox_max(cj,2) < bbox_min(ci,2));
+    };
+    
+    // Loop over components. This section is O(m²)
+    for (size_t i=0; i<num_components; i++)
+    {
+      // List of components that could overlap with component i
+      std::vector<size_t> candidate_comps;
+      candidate_comps.reserve(num_components);
+      // Loop over components
+      for (size_t j=0; j<num_components; j++) 
+      {
+        if (i == j) continue;
+        if (bbox_intersects(i,j)) candidate_comps.push_back(j);
+      }
 
-            DerivedV queries(num_candidate_comps, 3);
-            for (size_t j=0; j<num_candidate_comps; j++) {
-                const size_t index = candidate_comps[j];
-                queries.row(j) = get_triangle_center(outer_facets[index]);
-            }
+      const size_t num_candidate_comps = candidate_comps.size();
+      if (num_candidate_comps == 0) continue;
 
-            const auto& I = Is[i];
-            Eigen::VectorXi closest_facets, closest_facet_orientations;
-            igl::copyleft::cgal::closest_facet(V, F, I, queries,
-                    uE2E, EMAP, closest_facets, closest_facet_orientations);
+      // Get query points on each candidate component: barycenter of
+      // outer-facet 
+      DerivedV queries(num_candidate_comps, 3);
+      for (size_t j=0; j<num_candidate_comps; j++)
+      {
+        const size_t index = candidate_comps[j];
+        queries.row(j) = get_triangle_center(outer_facets[index]);
+      }
 
-            for (size_t j=0; j<num_candidate_comps; j++) {
-                const size_t index = candidate_comps[j];
-                const size_t closest_patch = P[closest_facets[j]];
-                const size_t closest_patch_side = closest_facet_orientations[j]
-                    ? 0:1;
-                const size_t ambient_cell = raw_cells(closest_patch,
-                        closest_patch_side);
-                if (ambient_cell != (size_t)outer_cells[i]) {
-                    nested_cells[ambient_cell].push_back(outer_cells[index]);
-                    ambient_cells[outer_cells[index]].push_back(ambient_cell);
-                    ambient_comps[index].push_back(i);
-                }
-            }
+      // Gather closest facets in ith component to each query point and their
+      // orientations
+      const auto& I = Is[i];
+      Eigen::VectorXi closest_facets, closest_facet_orientations;
+      closest_facet(V, F, I, queries,
+        uE2E, EMAP, closest_facets, closest_facet_orientations);
+      // Loop over all candidates
+      for (size_t j=0; j<num_candidate_comps; j++)
+      {
+        const size_t index = candidate_comps[j];
+        const size_t closest_patch = P[closest_facets[j]];
+        const size_t closest_patch_side = closest_facet_orientations[j] ? 0:1;
+        // The cell id of the closest patch
+        const size_t ambient_cell =
+          raw_cells(closest_patch,closest_patch_side);
+        if (ambient_cell != (size_t)outer_cells[i])
+        {
+          // ---> component index inside component i, because the cell of the
+          // closest facet on i to component index is **not** the same as the
+          // "outer cell" of component i: component index is **not** outside of
+          // component i (therefore it's inside).
+          nested_cells[ambient_cell].push_back(outer_cells[index]);
+          ambient_cells[outer_cells[index]].push_back(ambient_cell);
+          ambient_comps[index].push_back(i);
         }
+      }
     }
+  }
+
 #ifdef EXTRACT_CELLS_DEBUG
-    std::cout << "Determine nested relaitonship: " << tictoc() << std::endl;
+    log_time("nested_relationship");
 #endif
 
     const size_t INVALID = std::numeric_limits<size_t>::max();
@@ -348,6 +284,10 @@ IGL_INLINE size_t igl::copyleft::cgal::extract_cells(
 
     size_t count = 0;
     std::vector<size_t> mapped_indices(num_raw_cells+1, INVALID);
+    // Always map infinite cell to index 0.
+    mapped_indices[INFINITE_CELL] = count;
+    count++;
+
     for (size_t i=0; i<num_patches; i++) {
         const size_t old_positive_cell_id = raw_cells(i, 0);
         const size_t old_negative_cell_id = raw_cells(i, 1);
@@ -371,42 +311,245 @@ IGL_INLINE size_t igl::copyleft::cgal::extract_cells(
     }
     cells = raw_cells;
 #ifdef EXTRACT_CELLS_DEBUG
-    std::cout << "Finalize and output: " << tictoc() << std::endl;
+    log_time("finalize");
 #endif
     return count;
 }
 
 template<
-typename DerivedV,
-typename DerivedF,
-typename DerivedC >
-IGL_INLINE size_t igl::copyleft::cgal::extract_cells(
-        const Eigen::PlainObjectBase<DerivedV>& V,
-        const Eigen::PlainObjectBase<DerivedF>& F,
-        Eigen::PlainObjectBase<DerivedC>& cells) {
-    const size_t num_faces = F.rows();
-    //typedef typename DerivedF::Scalar Index;
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedP,
+  typename DeriveduE,
+  typename uE2EType,
+  typename DerivedEMAP,
+  typename DerivedC>
+IGL_INLINE size_t igl::copyleft::cgal::extract_cells_single_component(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  const Eigen::PlainObjectBase<DerivedP>& P,
+  const Eigen::PlainObjectBase<DeriveduE>& uE,
+  const std::vector<std::vector<uE2EType> >& uE2E,
+  const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+  Eigen::PlainObjectBase<DerivedC>& cells)
+{
+  const size_t num_faces = F.rows();
+  // Input:
+  //   index  index into #F*3 list of undirect edges
+  // Returns index into face
+  const auto edge_index_to_face_index = [&num_faces](size_t index)
+  {
+    return index % num_faces;
+  };
+  // Determine if a face (containing undirected edge {s,d} is consistently
+  // oriented with directed edge {s,d} (or otherwise it is with {d,s})
+  // 
+  // Inputs:
+  //   fid  face index into F
+  //   s  source index of edge
+  //   d  destination index of edge
+  // Returns true if face F(fid,:) is consistent with {s,d}
+  const auto is_consistent = 
+    [&F](const size_t fid, const size_t s, const size_t d) -> bool
+  {
+    if ((size_t)F(fid, 0) == s && (size_t)F(fid, 1) == d) return false;
+    if ((size_t)F(fid, 1) == s && (size_t)F(fid, 2) == d) return false;
+    if ((size_t)F(fid, 2) == s && (size_t)F(fid, 0) == d) return false;
 
-    Eigen::MatrixXi E, uE;
-    Eigen::VectorXi EMAP;
-    std::vector<std::vector<size_t> > uE2E;
-    igl::unique_edge_map(F, E, uE, EMAP, uE2E);
+    if ((size_t)F(fid, 0) == d && (size_t)F(fid, 1) == s) return true;
+    if ((size_t)F(fid, 1) == d && (size_t)F(fid, 2) == s) return true;
+    if ((size_t)F(fid, 2) == d && (size_t)F(fid, 0) == s) return true;
+    throw "Invalid face!";
+    return false;
+  };
+
+  const size_t num_unique_edges = uE.rows();
+  const size_t num_patches = P.maxCoeff() + 1;
+
+  // patch_edge_adj[p] --> list {e,f,g,...} such that p is a patch id and
+  //   e,f,g etc. are edge indices into 
+  std::vector<std::vector<size_t> > patch_edge_adj(num_patches);
+  // orders[u] --> J  where u is an index into unique edges uE and J is a
+  //   #adjacent-faces list of face-edge indices into F*3 sorted cyclicaly around
+  //   edge u.
+  std::vector<Eigen::VectorXi> orders(num_unique_edges);
+  // orientations[u] ---> list {f1,f2, ...} where u is an index into unique edges uE
+  //   and points to #adj-facets long list of flags whether faces are oriented
+  //   to point their normals clockwise around edge when looking along the
+  //   edge.
+  std::vector<std::vector<bool> > orientations(num_unique_edges);
+  // Loop over unique edges
+  for (size_t i=0; i<num_unique_edges; i++) 
+  {
+    const size_t s = uE(i,0);
+    const size_t d = uE(i,1);
+    const auto adj_faces = uE2E[i];
+    // If non-manifold (more than two incident faces)
+    if (adj_faces.size() > 2) 
+    {
+      // signed_adj_faces[a] --> sid  where a is an index into adj_faces
+      // (list of face edges on {s,d}) and sid is a signed index for resolve
+      // co-planar duplicates consistently (i.e. using simulation of
+      // simplicity).
+      std::vector<int> signed_adj_faces;
+      for (auto ei : adj_faces)
+      {
+        const size_t fid = edge_index_to_face_index(ei);
+        bool cons = is_consistent(fid, s, d);
+        signed_adj_faces.push_back((fid+1)*(cons ? 1:-1));
+      }
+      {
+        // Sort adjacent faces cyclically around {s,d}
+        auto& order = orders[i];
+        // order[f] will reveal the order of face f in signed_adj_faces
+        order_facets_around_edge(V, F, s, d, signed_adj_faces, order);
+        // Determine if each facet is oriented to point its normal clockwise or
+        // not around the {s,d} (normals are not explicitly computed/used)
+        auto& orientation = orientations[i];
+        orientation.resize(order.size());
+        std::transform(
+          order.data(), 
+          order.data() + order.size(),
+          orientation.begin(), 
+          [&signed_adj_faces](const int i){ return signed_adj_faces[i] > 0;});
+        // re-index order from adjacent faces to full face list. Now order
+        // indexes F directly
+        std::transform(
+          order.data(), 
+          order.data() + order.size(),
+          order.data(), 
+          [&adj_faces](const int index){ return adj_faces[index];});
+      }
+      // loop over adjacent faces, remember that patch is adjacent to this edge
+      for(const auto & ei : adj_faces)
+      {
+        const size_t fid = edge_index_to_face_index(ei);
+        patch_edge_adj[P[fid]].push_back(ei);
+      }
+    }
+  }
 
-    Eigen::VectorXi P;
-    //const size_t num_patches = 
-      igl::extract_manifold_patches(F, EMAP, uE2E, P);
+  // Initialize all patch-to-cell indices as "invalid" (unlabeled)
+  const int INVALID = std::numeric_limits<int>::max();
+  cells.resize(num_patches, 2);
+  cells.setConstant(INVALID);
 
-    DerivedC per_patch_cells;
-    const size_t num_cells =
-        igl::copyleft::cgal::extract_cells(V, F, P, E, uE, uE2E, EMAP, per_patch_cells);
+  // Given a "seed" patch id, a cell id, and which side of the patch that cell
+  // lies, identify all other patches bounding this cell (and tell them that
+  // they do)
+  //
+  // Inputs:
+  //   seed_patch_id  index into patches
+  //   cell_idx  index into cells
+  //   seed_patch_side   0 or 1 depending on whether cell_idx is on left or
+  //     right side of seed_patch_id 
+  //   cells  #cells by 2 list of current assignment of cells incident on each
+  //     side of a patch
+  // Outputs:
+  //   cells  udpated (see input)
+  // 
+  const auto & peel_cell_bd = 
+    [&P,&patch_edge_adj,&EMAP,&orders,&orientations,&num_faces](
+    const size_t seed_patch_id, 
+    const short seed_patch_side, 
+    const size_t cell_idx,
+    Eigen::PlainObjectBase<DerivedC>& cells)
+  {
+    typedef std::pair<size_t, short> PatchSide;
+    // Initialize a queue of {p,side} patch id and patch side pairs to BFS over
+    // all patches
+    std::queue<PatchSide> Q;
+    Q.emplace(seed_patch_id, seed_patch_side);
+    // assign cell id of seed patch on appropriate side
+    cells(seed_patch_id, seed_patch_side) = cell_idx;
+    while (!Q.empty())
+    {
+      // Pop patch from Q
+      const auto entry = Q.front();
+      Q.pop();
+      const size_t patch_id = entry.first;
+      const short side = entry.second;
+      // face-edges adjacent to patch
+      const auto& adj_edges = patch_edge_adj[patch_id];
+      // Loop over **ALL EDGES IN THE ENTIRE PATCH** not even just the boundary
+      // edges, all edges... O(n)
+      for (const auto& ei : adj_edges)
+      {
+        // unique edge
+        const size_t uei = EMAP[ei];
+        // ordering of face-edges stored at edge
+        const auto& order = orders[uei];
+        // consistent orientation flags at face-edges stored at edge
+        const auto& orientation = orientations[uei];
+        const size_t edge_valance = order.size();
+        // Search for ei's (i.e. patch_id's) place in the ordering: O(#patches)
+        size_t curr_i = 0;
+        for (curr_i=0; curr_i < edge_valance; curr_i++)
+        {
+          if ((size_t)order[curr_i] == ei) break;
+        }
+        assert(curr_i < edge_valance && "Failed to find edge.");
+        // is the starting face consistent?
+        const bool cons = orientation[curr_i];
+        // Look clockwise or counter-clockwise for the next face, depending on
+        // whether looking to left or right side and whether consistently
+        // oriented or not.
+        size_t next;
+        if (side == 0)
+        {
+          next = (cons ? (curr_i + 1) :
+          (curr_i + edge_valance - 1)) % edge_valance;
+        } else {
+          next = (cons ? (curr_i+edge_valance-1) : (curr_i+1))%edge_valance;
+        }
+        // Determine face-edge index of next face
+        const size_t next_ei = order[next];
+        // Determine whether next is consistently oriented
+        const bool next_cons = orientation[next];
+        // Determine patch of next
+        const size_t next_patch_id = P[next_ei % num_faces];
+        // Determine which side of patch cell_idx is on, based on whether the
+        // consistency of next matches the consistency of this patch and which
+        // side of this patch we're on.
+        const short next_patch_side = (next_cons != cons) ?  side:abs(side-1);
+        // If that side of next's patch is not labeled, then label and add to
+        // queue
+        if (cells(next_patch_id, next_patch_side) == INVALID) 
+        {
+          Q.emplace(next_patch_id, next_patch_side);
+          cells(next_patch_id, next_patch_side) = cell_idx;
+        }else 
+        {
+          assert(
+            (size_t)cells(next_patch_id, next_patch_side) == cell_idx && 
+            "Encountered cell assignment inconsistency");
+        }
+      }
+    }
+  };
 
-    cells.resize(num_faces, 2);
-    for (size_t i=0; i<num_faces; i++) {
-        cells.row(i) = per_patch_cells.row(P[i]);
+  size_t count=0;
+  // Loop over all patches
+  for (size_t i=0; i<num_patches; i++)
+  {
+    // If left side of patch is still unlabeled, create a new cell and find all
+    // patches also on its boundary
+    if (cells(i, 0) == INVALID)
+    {
+      peel_cell_bd(i, 0, count,cells);
+      count++;
     }
-    return num_cells;
+    // Likewise for right side
+    if (cells(i, 1) == INVALID)
+    {
+      peel_cell_bd(i, 1, count,cells);
+      count++;
+    }
+  }
+  return count;
 }
 
+
 #ifdef IGL_STATIC_LIBRARY
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 template unsigned long igl::copyleft::cgal::extract_cells<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -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>, unsigned long, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -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> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<unsigned long, std::allocator<unsigned long> >, std::allocator<std::vector<unsigned long, std::allocator<unsigned long> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);

include/igl/copyleft/cgal/extract_cells.h

@@ -16,94 +16,95 @@
 namespace igl {
   namespace copyleft
   {
-    namespace cgal {
-        // Extract connected 3D space partitioned by mesh (V, F).
-        //
-        // Inputs:
-        //   V  #V by 3 array of vertices.
-        //   F  #F by 3 array of faces.
-        //
-        // Output:
-        //   cells  #F by 2 array of cell indices.  cells(i,0) represents the
-        //          cell index on the positive side of face i, and cells(i,1)
-        //          represents cell index of the negqtive side.
-        template<
-            typename DerivedV,
-            typename DerivedF,
-            typename DerivedC >
-        IGL_INLINE size_t extract_cells(
-                const Eigen::PlainObjectBase<DerivedV>& V,
-                const Eigen::PlainObjectBase<DerivedF>& F,
-                Eigen::PlainObjectBase<DerivedC>& cells);
+    namespace cgal
+    {
+      // Extract connected 3D space partitioned by mesh (V, F).
+      //
+      // Inputs:
+      //   V  #V by 3 array of vertices.
+      //   F  #F by 3 array of faces.
+      //
+      // Output:
+      //   cells  #F by 2 array of cell indices.  cells(i,0) represents the
+      //          cell index on the positive side of face i, and cells(i,1)
+      //          represents cell index of the negqtive side.
+      //          By convension cell with index 0 is the infinite cell.
+      template<
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedC >
+      IGL_INLINE size_t extract_cells(
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        Eigen::PlainObjectBase<DerivedC>& cells);
 
-        // Extract connected 3D space partitioned by mesh (V, F).
-        //
-        // Inputs:
-        //   V  #V by 3 array of vertices.
-        //   F  #F by 3 array of faces.
-        //   P  #F list of patch indices.
-        //   E  #E by 2 array of vertex indices, one edge per row.
-        //  uE    #uE by 2 list of vertex_indices, represents undirected edges.
-        //  uE2E  #uE list of lists that maps uE to E. (a one-to-many map)
-        //  EMAP  #F*3 list of indices into uE.
-        //
-        // Output:
-        //   cells  #P by 2 array of cell indices.  cells(i,0) represents the
-        //          cell index on the positive side of patch i, and cells(i,1)
-        //          represents cell index of the negqtive side.
-        template<
-            typename DerivedV,
-            typename DerivedF,
-            typename DerivedP,
-            typename DerivedE,
-            typename DeriveduE,
-            typename uE2EType,
-            typename DerivedEMAP,
-            typename DerivedC >
-        IGL_INLINE size_t extract_cells(
-                const Eigen::PlainObjectBase<DerivedV>& V,
-                const Eigen::PlainObjectBase<DerivedF>& F,
-                const Eigen::PlainObjectBase<DerivedP>& P,
-                const Eigen::PlainObjectBase<DerivedE>& E,
-                const Eigen::PlainObjectBase<DeriveduE>& uE,
-                const std::vector<std::vector<uE2EType> >& uE2E,
-                const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
-                Eigen::PlainObjectBase<DerivedC>& cells);
-
-        // Extract connected 3D space partitioned by mesh (V, F).  Each
-        // connected component of the mesh partitions the ambient space
-        // separately.
-        //
-        // Inputs:
-        //   V  #V by 3 array of vertices.
-        //   F  #F by 3 array of faces.
-        //   P  #F list of patch indices.
-        //   E  #E by 2 array of vertex indices, one edge per row.
-        //  uE    #uE by 2 list of vertex_indices, represents undirected edges.
-        //  uE2E  #uE list of lists that maps uE to E. (a one-to-many map)
-        //  EMAP  #F*3 list of indices into uE.
-        //
-        // Output:
-        //   cells  #P by 2 array of cell indices.  cells(i,0) represents the
-        //          cell index on the positive side of patch i, and cells(i,1)
-        //          represents cell index of the negqtive side.
-        template<
-            typename DerivedV,
-            typename DerivedF,
-            typename DerivedP,
-            typename DeriveduE,
-            typename uE2EType,
-            typename DerivedEMAP,
-            typename DerivedC >
-        IGL_INLINE size_t extract_cells_single_component(
-                const Eigen::PlainObjectBase<DerivedV>& V,
-                const Eigen::PlainObjectBase<DerivedF>& F,
-                const Eigen::PlainObjectBase<DerivedP>& P,
-                const Eigen::PlainObjectBase<DeriveduE>& uE,
-                const std::vector<std::vector<uE2EType> >& uE2E,
-                const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
-                Eigen::PlainObjectBase<DerivedC>& cells);
+      // Extract connected 3D space partitioned by mesh (V, F).
+      //
+      // Inputs:
+      //   V  #V by 3 array of vertices.
+      //   F  #F by 3 array of faces.
+      //   P  #F list of patch indices.
+      //   E  #E by 2 array of vertex indices, one edge per row.
+      //  uE    #uE by 2 list of vertex_indices, represents undirected edges.
+      //  uE2E  #uE list of lists that maps uE to E. (a one-to-many map)
+      //  EMAP  #F*3 list of indices into uE.
+      //
+      // Output:
+      //   cells  #P by 2 array of cell indices.  cells(i,0) represents the
+      //          cell index on the positive side of patch i, and cells(i,1)
+      //          represents cell index of the negqtive side.
+      //          By convension cell with index 0 is the infinite cell.
+      template<
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedP,
+        typename DerivedE,
+        typename DeriveduE,
+        typename uE2EType,
+        typename DerivedEMAP,
+        typename DerivedC >
+      IGL_INLINE size_t extract_cells(
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        const Eigen::PlainObjectBase<DerivedP>& P,
+        const Eigen::PlainObjectBase<DerivedE>& E,
+        const Eigen::PlainObjectBase<DeriveduE>& uE,
+        const std::vector<std::vector<uE2EType> >& uE2E,
+        const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+        Eigen::PlainObjectBase<DerivedC>& cells);
 
+      // Extract connected 3D space partitioned by mesh (V,F) composed of
+      // **possibly multiple components** (the name of this function is
+      // dubious).
+      //
+      // Inputs:
+      //   V  #V by 3 array of vertices.
+      //   F  #F by 3 array of faces.
+      //   P  #F list of patch indices.
+      //   E  #E by 2 array of vertex indices, one edge per row.
+      //  uE    #uE by 2 list of vertex_indices, represents undirected edges.
+      //  uE2E  #uE list of lists that maps uE to E. (a one-to-many map)
+      //  EMAP  #F*3 list of indices into uE.
+      // Output:
+      //  cells  #P by 2 array of cell indices.  cells(i,0) represents the
+      //    cell index on the positive side of patch i, and cells(i,1)
+      //    represents cell index of the negative side.
+      template<
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedP,
+        typename DeriveduE,
+        typename uE2EType,
+        typename DerivedEMAP,
+        typename DerivedC >
+      IGL_INLINE size_t extract_cells_single_component(
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        const Eigen::PlainObjectBase<DerivedP>& P,
+        const Eigen::PlainObjectBase<DeriveduE>& uE,
+        const std::vector<std::vector<uE2EType> >& uE2E,
+        const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+        Eigen::PlainObjectBase<DerivedC>& cells);
     }
   }
 }

include/igl/copyleft/cgal/intersect_other.cpp

@@ -25,11 +25,12 @@ namespace igl
       static IGL_INLINE void push_result(
         const Eigen::PlainObjectBase<DerivedF> & F,
         const int f,
+        const int f_other,
         const CGAL::Object & result,
         std::map<
           typename DerivedF::Index,
-          std::pair<typename DerivedF::Index,
-            std::vector<CGAL::Object> > > & offending,
+          std::vector<std::pair<typename DerivedF::Index, CGAL::Object> > > &
+          offending,
         std::map<
           std::pair<typename DerivedF::Index,typename DerivedF::Index>,
           std::vector<typename DerivedF::Index> > & edge2faces)
@@ -39,8 +40,7 @@ namespace igl
         if(offending.count(f) == 0)
         {
           // first time marking, initialize with new id and empty list
-          Index id = offending.size();
-          offending[f] = {id,{}};
+          offending[f] = {};
           for(Index e = 0; e<3;e++)
           {
             // append face to edge's list
@@ -49,7 +49,7 @@ namespace igl
             edge2faces[EMK(i,j)].push_back(f);
           }
         }
-        offending[f].second.push_back(result);
+        offending[f].push_back({f_other,result});
       }
       template <
         typename Kernel,
@@ -113,7 +113,7 @@ namespace igl
           CGAL::Box_intersection_d::Box_with_handle_d<double,3,TrianglesIterator> 
           Box;
         typedef 
-          std::map<Index,std::pair<Index,std::vector<CGAL::Object> > > 
+          std::map<Index,std::vector<std::pair<Index,CGAL::Object> > > 
           OffendingMap;
         typedef std::map<std::pair<Index,Index>,std::vector<Index> >  EdgeMap;
         typedef std::pair<Index,Index> EMK;
@@ -125,7 +125,7 @@ namespace igl
         // 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> A_boxes,B_boxes;
-        const auto box_up = [](Triangles & T, std::vector<Box> & boxes)
+        const auto box_up = [](Triangles & T, std::vector<Box> & boxes) -> void
         {
           boxes.reserve(T.size());
           for ( 
@@ -142,7 +142,7 @@ namespace igl
         EdgeMap edge2facesA,edge2facesB;
 
         std::list<int> lIF;
-        const auto cb = [&](const Box &a, const Box &b)
+        const auto cb = [&](const Box &a, const Box &b) -> void
         {
           using namespace std;
           // index in F and T
@@ -163,8 +163,8 @@ namespace igl
             {
               CGAL::Object result = CGAL::intersection(A,B);
 
-              push_result(FA,fa,result,offendingA,edge2facesA);
-              push_result(FB,fb,result,offendingB,edge2facesB);
+              push_result(FA,fa,fb,result,offendingA,edge2facesA);
+              push_result(FB,fb,fa,result,offendingB,edge2facesB);
             }
           }
         };

include/igl/copyleft/cgal/order_facets_around_edge.cpp

@@ -295,7 +295,7 @@ void igl::copyleft::cgal::order_facets_around_edge(
     {
         return abs(signed_idx) -1;
     };
-    auto get_opposite_vertex_index = [&](size_t fid)
+    auto get_opposite_vertex_index = [&](size_t fid) -> typename DerivedF::Scalar
     {
         typedef typename DerivedF::Scalar Index;
         if (F(fid, 0) != (Index)s && F(fid, 0) != (Index)d) return F(fid, 0);

include/igl/copyleft/cgal/outer_element.cpp

@@ -131,7 +131,7 @@ IGL_INLINE void igl::copyleft::cgal::outer_edge(
     Index outer_opp_vid = INVALID;
     bool infinite_slope_detected = false;
     std::vector<Index> incident_faces;
-    auto check_and_update_outer_edge = [&](Index opp_vid, Index fid) {
+    auto check_and_update_outer_edge = [&](Index opp_vid, Index fid) -> void {
         if (opp_vid == outer_opp_vid)
         {
             incident_faces.push_back(fid);

include/igl/copyleft/cgal/points_inside_component.cpp

@@ -169,7 +169,7 @@ namespace igl {
 
                 // A plane is on the exterior if all adj_points lies on or to
                 // one side of the plane.
-                auto is_on_exterior = [&](const Plane_3& separator) {
+                auto is_on_exterior = [&](const Plane_3& separator) -> bool{
                     size_t positive=0;
                     size_t negative=0;
                     size_t coplanar=0;

include/igl/copyleft/cgal/propagate_winding_numbers.cpp

@@ -40,13 +40,17 @@ IGL_INLINE void igl::copyleft::cgal::propagate_winding_numbers(
     const Eigen::PlainObjectBase<DerivedL>& labels,
     Eigen::PlainObjectBase<DerivedW>& W) {
 #ifdef PROPAGATE_WINDING_NUMBER_TIMING
-  const auto & tictoc = []()
+  const auto & tictoc = []() -> double
   {
     static double t_start = igl::get_seconds();
     double diff = igl::get_seconds()-t_start;
     t_start += diff;
     return diff;
   };
+  const auto log_time = [&](const std::string& label) -> void {
+    std::cout << "propagate_winding_num." << label << ": "
+      << tictoc() << std::endl;
+  };
   tictoc();
 #endif
   const size_t num_faces = F.rows();
@@ -63,16 +67,13 @@ IGL_INLINE void igl::copyleft::cgal::propagate_winding_numbers(
 
   Eigen::VectorXi P;
   const size_t num_patches = igl::extract_manifold_patches(F, EMAP, uE2E, P);
-#ifdef PROPAGATE_WINDING_NUMBER_TIMING
-  std::cout << "extract manifold patches: " << tictoc() << std::endl;
-#endif
 
   DerivedW per_patch_cells;
   const size_t num_cells =
     igl::copyleft::cgal::extract_cells(
         V, F, P, E, uE, uE2E, EMAP, per_patch_cells);
 #ifdef PROPAGATE_WINDING_NUMBER_TIMING
-  std::cout << "extract cells: " << tictoc() << std::endl;;
+  log_time("cell_extraction");
 #endif
 
   typedef std::tuple<size_t, bool, size_t> CellConnection;
@@ -84,10 +85,10 @@ IGL_INLINE void igl::copyleft::cgal::propagate_winding_numbers(
     cell_adjacency[negative_cell].emplace(positive_cell, true, i);
   }
 #ifdef PROPAGATE_WINDING_NUMBER_TIMING
-  std::cout << "cell connection: " << tictoc() << std::endl;
+  log_time("cell_connectivity");
 #endif
 
-  auto save_cell = [&](const std::string& filename, size_t cell_id) {
+  auto save_cell = [&](const std::string& filename, size_t cell_id) -> void{
     std::vector<size_t> faces;
     for (size_t i=0; i<num_patches; i++) {
       if ((per_patch_cells.row(i).array() == cell_id).any()) {
@@ -118,7 +119,7 @@ IGL_INLINE void igl::copyleft::cgal::propagate_winding_numbers(
     cell_labels.setZero();
     Eigen::VectorXi parents(num_cells);
     parents.setConstant(-1);
-    auto trace_parents = [&](size_t idx) {
+    auto trace_parents = [&](size_t idx) -> std::list<size_t> {
       std::list<size_t> path;
       path.push_back(idx);
       while ((size_t)parents[path.back()] != path.back()) {
@@ -169,7 +170,7 @@ IGL_INLINE void igl::copyleft::cgal::propagate_winding_numbers(
       }
     }
 #ifdef PROPAGATE_WINDING_NUMBER_TIMING
-    std::cout << "check for odd cycle: " << tictoc() << std::endl;
+    log_time("odd_cycle_check");
 #endif
   }
 #endif
@@ -180,7 +181,7 @@ IGL_INLINE void igl::copyleft::cgal::propagate_winding_numbers(
   I.setLinSpaced(num_faces, 0, num_faces-1);
   igl::copyleft::cgal::outer_facet(V, F, I, outer_facet, flipped);
 #ifdef PROPAGATE_WINDING_NUMBER_TIMING
-  std::cout << "outer facet: " << tictoc() << std::endl;
+  log_time("outer_facet");
 #endif
 
   const size_t outer_patch = P[outer_facet];
@@ -244,7 +245,7 @@ IGL_INLINE void igl::copyleft::cgal::propagate_winding_numbers(
     }
   }
 #ifdef PROPAGATE_WINDING_NUMBER_TIMING
-  std::cout << "propagate winding number: " << tictoc() << std::endl;
+  log_time("propagate_winding_number");
 #endif
 
   W.resize(num_faces, num_labels*2);
@@ -258,7 +259,7 @@ IGL_INLINE void igl::copyleft::cgal::propagate_winding_numbers(
     }
   }
 #ifdef PROPAGATE_WINDING_NUMBER_TIMING
-  std::cout << "save result: " << tictoc() << std::endl;
+  log_time("store_result");
 #endif
 }
 

include/igl/copyleft/cgal/remesh_intersections.h

@@ -53,8 +53,8 @@ namespace igl
         const std::vector<CGAL::Triangle_3<Kernel> > & T,
         const std::map<
           typename DerivedF::Index,
-          std::pair<typename DerivedF::Index,
-            std::vector<CGAL::Object> > > & offending,
+            std::vector<
+            std::pair<typename DerivedF::Index, CGAL::Object> > > & offending,
         const std::map<
           std::pair<typename DerivedF::Index,typename DerivedF::Index>,
           std::vector<typename DerivedF::Index> > & edge2faces,

include/igl/copyleft/cgal/remesh_self_intersections.cpp

@@ -56,7 +56,7 @@ IGL_INLINE void igl::copyleft::cgal::remesh_self_intersections(
         DerivedJ,
         DerivedIM>
       SelfIntersectMeshK;
-    SelfIntersectMeshK SIM = SelfIntersectMeshK(V,F,params,VV,FF,IF,J,IM);
+    SelfIntersectMeshK SIM(V,F,params,VV,FF,IF,J,IM);
 
 //#ifdef __APPLE__
 //    signal(SIGFPE,SIG_DFL);
@@ -76,7 +76,7 @@ IGL_INLINE void igl::copyleft::cgal::remesh_self_intersections(
         DerivedJ,
         DerivedIM>
       SelfIntersectMeshK;
-    SelfIntersectMeshK SIM = SelfIntersectMeshK(V,F,params,VV,FF,IF,J,IM);
+    SelfIntersectMeshK SIM(V,F,params,VV,FF,IF,J,IM);
   }
 }
 

include/igl/copyleft/tetgen/mesh_to_tetgenio.cpp

@@ -64,7 +64,6 @@ IGL_INLINE bool igl::copyleft::tetgen::mesh_to_tetgenio(
   const Eigen::PlainObjectBase<DerivedF>& F,
   tetgenio & in)
 {
-  using namespace igl;
   using namespace std;
   vector<vector<REAL> > vV;
   vector<vector<int> > vF;

include/igl/copyleft/tetgen/mesh_with_skeleton.cpp

@@ -32,7 +32,6 @@ IGL_INLINE bool igl::copyleft::tetgen::mesh_with_skeleton(
   Eigen::MatrixXi & FF)
 {
   using namespace Eigen;
-  using namespace igl;
   using namespace std;
   const string eff_tetgen_flags = 
     (tetgen_flags.length() == 0?DEFAULT_TETGEN_FLAGS:tetgen_flags);

include/igl/copyleft/tetgen/read_into_tetgenio.cpp

@@ -29,7 +29,6 @@ IGL_INLINE bool igl::copyleft::tetgen::read_into_tetgenio(
   const std::string & path,
   tetgenio & in)
 {
-  using namespace igl;
   using namespace std;
   // get file extension
   string dirname,basename,ext,filename;

include/igl/copyleft/tetgen/tetgenio_to_tetmesh.cpp

@@ -100,7 +100,6 @@ IGL_INLINE bool igl::copyleft::tetgen::tetgenio_to_tetmesh(
   Eigen::PlainObjectBase<DerivedT>& T,
   Eigen::PlainObjectBase<DerivedF>& F)
 {
-  using namespace igl;
   using namespace std;
   vector<vector<REAL> > vV;
   vector<vector<int> > vT;

include/igl/copyleft/tetgen/tetrahedralize.cpp

@@ -73,7 +73,6 @@ IGL_INLINE int igl::copyleft::tetgen::tetrahedralize(
   Eigen::PlainObjectBase<DerivedTT>& TT,
   Eigen::PlainObjectBase<DerivedTF>& TF)
 {
-  using namespace igl;
   using namespace std;
   vector<vector<REAL> > vV,vTV;
   vector<vector<int> > vF,vTT,vTF;

include/igl/cut_mesh.h

@@ -5,7 +5,6 @@
 // This Source Code Form is subject to the terms of the Mozilla Public License
 // v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
-
 #ifndef IGL_CUT_MESH
 #define IGL_CUT_MESH
 #include "igl_inline.h"
@@ -13,54 +12,64 @@
 #include <Eigen/Core>
 #include <vector>
 
-namespace igl {
-  // Given a mesh and a list of edges that are to be cut, the function generates a
-  // new disk-topology mesh that has the cuts at its boundary.
+namespace igl
+{
+  // Given a mesh and a list of edges that are to be cut, the function
+  // generates a new disk-topology mesh that has the cuts at its boundary.
+  //
+  // Todo: this combinatorial operation should not depend on the vertex
+  // positions V.
+  //
   // Inputs:
-  //   V                #V by 3 list of the vertex positions
-  //   F                #F by 3 list of the faces (must be triangles)
-  //   VF               #V list of lists of incident faces (adjacency list), e.g.
-  //                    as returned by igl::vertex_triangle_adjacency
-  //   VFi               #V list of lists of index of incidence within incident faces listed
-  //                    in VF, e.g. as returned by igl::vertex_triangle_adjacency
-  //   TT               #F by 3 triangle to triangle adjacent matrix (e.g. computed
-  //                    via igl:triangle_triangle_adjacency)
-  //   TTi              #F by 3 adjacent matrix, the element i,j is the id of edge of the
-  //                    triangle TT(i,j) that is adjacent with triangle i (e.g. computed
-  //                    via igl:triangle_triangle_adjacency)
-  //   V_border         #V by 1 list of booleans, indicating if the corresponging vertex is
-  //                    at the mesh boundary, e.g. as returned by igl::is_border_vertex
-  //   cuts             #F by 3 list of boolean flags, indicating the edges that need to be cut
+  //   V  #V by 3 list of the vertex positions
+  //   F  #F by 3 list of the faces (must be triangles)
+  //   VF  #V list of lists of incident faces (adjacency list), e.g.  as
+  //     returned by igl::vertex_triangle_adjacency
+  //   VFi  #V list of lists of index of incidence within incident faces listed
+  //     in VF, e.g. as returned by igl::vertex_triangle_adjacency
+  //   TT  #F by 3 triangle to triangle adjacent matrix (e.g. computed via
+  //     igl:triangle_triangle_adjacency)
+  //   TTi  #F by 3 adjacent matrix, the element i,j is the id of edge of the
+  //     triangle TT(i,j) that is adjacent with triangle i (e.g. computed via
+  //     igl:triangle_triangle_adjacency)
+  //   V_border  #V by 1 list of booleans, indicating if the corresponging
+  //     vertex is at the mesh boundary, e.g. as returned by
+  //     igl::is_border_vertex
+  //   cuts  #F by 3 list of boolean flags, indicating the edges that need to
+  //     be cut (has 1 at the face edges that are to be cut, 0 otherwise)
   // Outputs:
-  //                    (has 1 at the face edges that are to be cut, 0 otherwise)
-  //   Vcut             #V by 3 list of the vertex positions of the cut mesh. This matrix will be
-  //                    similar to the original vertices except some rows will be duplicated.
-  //   Fcut             #F by 3 list of the faces of the cut mesh(must be triangles). This matrix
-  //                    will be similar to the original face matrix except some indices will be
-  //                    redirected to point to the newly duplicated vertices.
+  //   Vcut  #V by 3 list of the vertex positions of the cut mesh. This matrix
+  //     will be similar to the original vertices except some rows will be
+  //     duplicated.
+  //   Fcut  #F by 3 list of the faces of the cut mesh(must be triangles). This
+  //     matrix will be similar to the original face matrix except some indices
+  //     will be redirected to point to the newly duplicated vertices.
   //
-  template <typename DerivedV, typename DerivedF, typename VFType, typename DerivedTT, typename DerivedC>
-  IGL_INLINE void cut_mesh(const Eigen::PlainObjectBase<DerivedV> &V,
-                           const Eigen::PlainObjectBase<DerivedF> &F,
-                           const std::vector<std::vector<VFType> >& VF,
-                           const std::vector<std::vector<VFType> >& VFi,
-                           const Eigen::PlainObjectBase<DerivedTT>& TT,
-                           const Eigen::PlainObjectBase<DerivedTT>& TTi,
-                           const std::vector<bool> &V_border,
-                           const Eigen::PlainObjectBase<DerivedC> &cuts,
-                           Eigen::PlainObjectBase<DerivedV> &Vcut,
-                           Eigen::PlainObjectBase<DerivedF> &Fcut);
-  
-  
+  template <
+    typename DerivedV, 
+    typename DerivedF, 
+    typename VFType, 
+    typename DerivedTT, 
+    typename DerivedC>
+  IGL_INLINE void cut_mesh(
+    const Eigen::PlainObjectBase<DerivedV> &V,
+    const Eigen::PlainObjectBase<DerivedF> &F,
+    const std::vector<std::vector<VFType> >& VF,
+    const std::vector<std::vector<VFType> >& VFi,
+    const Eigen::PlainObjectBase<DerivedTT>& TT,
+    const Eigen::PlainObjectBase<DerivedTT>& TTi,
+    const std::vector<bool> &V_border,
+    const Eigen::PlainObjectBase<DerivedC> &cuts,
+    Eigen::PlainObjectBase<DerivedV> &Vcut,
+    Eigen::PlainObjectBase<DerivedF> &Fcut);
   //Wrapper of the above with only vertices and faces as mesh input
   template <typename DerivedV, typename DerivedF, typename DerivedC>
   IGL_INLINE void cut_mesh(
-                           const Eigen::PlainObjectBase<DerivedV> &V,
-                           const Eigen::PlainObjectBase<DerivedF> &F,
-                           const Eigen::PlainObjectBase<DerivedC> &cuts,
-                           Eigen::PlainObjectBase<DerivedV> &Vcut,
-                           Eigen::PlainObjectBase<DerivedF> &Fcut);
-  
+    const Eigen::PlainObjectBase<DerivedV> &V,
+    const Eigen::PlainObjectBase<DerivedF> &F,
+    const Eigen::PlainObjectBase<DerivedC> &cuts,
+    Eigen::PlainObjectBase<DerivedV> &Vcut,
+    Eigen::PlainObjectBase<DerivedF> &Fcut);
 };
 
 
@@ -69,4 +78,4 @@ namespace igl {
 #endif
 
 
-#endif /* defined(IGL_CUT_MESH) */
+#endif

include/igl/cut_mesh_from_singularities.h

@@ -12,21 +12,28 @@
 #include <Eigen/Core>
 namespace igl
 {
-  // Given a mesh (V,F) and the integer mismatch of a cross field per edge (MMatch),
-  // finds the cut_graph connecting the singularities (seams) and the degree of the singularities
-  // singularity_index
+  // Given a mesh (V,F) and the integer mismatch of a cross field per edge
+  // (MMatch), finds the cut_graph connecting the singularities (seams) and the
+  // degree of the singularities singularity_index
   //
   // Input:
-  // V:                 #V by 3 list of mesh vertex positions
-  // F:                 #F by 3 list of faces
-  // MMatch:            #F by 3 list of per corner integer mismatch
-  // seams:             #F by 3 list of per corner booleans that denotes if an edge is a seam or not
+  //   V  #V by 3 list of mesh vertex positions
+  //   F  #F by 3 list of faces
+  //   MMatch  #F by 3 list of per corner integer mismatch
+  // Outputs:
+  //   seams  #F by 3 list of per corner booleans that denotes if an edge is a
+  //     seam or not
   //
-  template <typename DerivedV, typename DerivedF, typename DerivedM, typename DerivedO>
-  IGL_INLINE void cut_mesh_from_singularities(const Eigen::PlainObjectBase<DerivedV> &V,
-                                                   const Eigen::PlainObjectBase<DerivedF> &F,
-                                                   const Eigen::PlainObjectBase<DerivedM> &MMatch,
-                                                   Eigen::PlainObjectBase<DerivedO> &seams);
+  template <
+    typename DerivedV, 
+    typename DerivedF, 
+    typename DerivedM, 
+    typename DerivedO> 
+  IGL_INLINE void cut_mesh_from_singularities(
+    const Eigen::PlainObjectBase<DerivedV> &V, 
+    const Eigen::PlainObjectBase<DerivedF> &F, 
+    const Eigen::PlainObjectBase<DerivedM> &MMatch,
+    Eigen::PlainObjectBase<DerivedO> &seams);
 }
 #ifndef IGL_STATIC_LIBRARY
 #include "cut_mesh_from_singularities.cpp"

include/igl/diag.cpp

@@ -95,7 +95,6 @@ IGL_INLINE void igl::diag(
   for(int i = 0;i<V.size();i++)
   {
     dyn_X.coeffRef(i,i) += V[i];
-    i++;
   }
   X = Eigen::SparseMatrix<T>(dyn_X);
 }

include/igl/edge_lengths.cpp

@@ -96,4 +96,7 @@ template void igl::edge_lengths<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::M
 template void igl::edge_lengths<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, 3, 1, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&);
 template void igl::edge_lengths<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 6, 0, -1, 6> >(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, 6, 0, -1, 6> >&);
 template void igl::edge_lengths<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 4, 0, -1, 4>, 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, 4, 0, -1, 4> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+template void igl::edge_lengths<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+template void igl::edge_lengths<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> >&);
+template void igl::edge_lengths<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<float, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> >&);
 #endif

include/igl/edge_lengths.h

@@ -32,7 +32,7 @@ namespace igl
   //     for edges, column of lengths
   //     for triangles, columns correspond to edges [1,2],[2,0],[0,1]
   //     for tets, columns correspond to edges
-  //     [1,2],[2,0],[0,1],[3,0],[3,1],[3,2]
+  //     [3 0],[3 1],[3 2],[1 2],[2 0],[0 1]
   //
   template <typename DerivedV, typename DerivedF, typename DerivedL>
   IGL_INLINE void edge_lengths(

include/igl/edges.cpp

@@ -6,13 +6,16 @@
 // 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 "edges.h"
-
 #include "adjacency_matrix.h"
 
-IGL_INLINE void igl::edges( const Eigen::MatrixXi& F, Eigen::MatrixXi& E)
+template <typename DerivedF, typename DerivedE>
+IGL_INLINE void igl::edges(
+  const Eigen::PlainObjectBase<DerivedF> & F, 
+  Eigen::PlainObjectBase<DerivedE> & E)
 {
   // build adjacency matrix
-  Eigen::SparseMatrix<int> A;
+  typedef typename DerivedF::Scalar Index;
+  Eigen::SparseMatrix<Index> A;
   igl::adjacency_matrix(F,A);
   // Number of non zeros should be twice number of edges
   assert(A.nonZeros()%2 == 0);
@@ -23,7 +26,7 @@ IGL_INLINE void igl::edges( const Eigen::MatrixXi& F, Eigen::MatrixXi& E)
   for(int k=0; k<A.outerSize(); ++k)
   {
     // Iterate over inside
-    for(Eigen::SparseMatrix<int>::InnerIterator it (A,k); it; ++it)
+    for(typename Eigen::SparseMatrix<Index>::InnerIterator it (A,k); it; ++it)
     {
       // only add edge in one direction
       if(it.row()<it.col())

include/igl/edges.h

@@ -24,7 +24,10 @@ namespace igl
   //   E #E by 2 list of edges in no particular order
   //
   // See also: adjacency_matrix
-  IGL_INLINE void edges( const Eigen::MatrixXi& F, Eigen::MatrixXi& E);
+  template <typename DerivedF, typename DerivedE>
+  IGL_INLINE void edges(
+    const Eigen::PlainObjectBase<DerivedF> & F, 
+    Eigen::PlainObjectBase<DerivedE> & E);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/embree/ambient_occlusion.cpp

@@ -46,7 +46,6 @@ IGL_INLINE void igl::embree::ambient_occlusion(
   const int num_samples,
   Eigen::PlainObjectBase<DerivedS> & S)
 {
-  using namespace igl;
   using namespace Eigen;
   EmbreeIntersector ei;
   ei.init(V.template cast<float>(),F.template cast<int>());

include/igl/embree/bone_visible.cpp

@@ -47,7 +47,6 @@ IGL_INLINE void igl::embree::bone_visible(
   const Eigen::PlainObjectBase<DerivedSD> & d,
   Eigen::PlainObjectBase<Derivedflag>  & flag)
 {
-  using namespace igl;
   using namespace std;
   using namespace Eigen;
   flag.resize(V.rows());

include/igl/embree/unproject_in_mesh.cpp

@@ -21,7 +21,6 @@ IGL_INLINE int igl::embree::unproject_in_mesh(
   Eigen::PlainObjectBase<Derivedobj> & obj,
   std::vector<igl::Hit > & hits)
 {
-  using namespace igl;
   using namespace std;
   using namespace Eigen;
   const auto & shoot_ray = [&ei](

include/igl/extract_manifold_patches.cpp

@@ -20,7 +20,7 @@ IGL_INLINE size_t igl::extract_manifold_patches(
     auto face_and_corner_index_to_edge_index = [&](size_t fi, size_t ci) {
         return ci*num_faces + fi;
     };
-    auto is_manifold_edge = [&](size_t fi, size_t ci) {
+    auto is_manifold_edge = [&](size_t fi, size_t ci) -> bool {
         const size_t ei = face_and_corner_index_to_edge_index(fi, ci);
         return uE2E[EMAP(ei, 0)].size() == 2;
     };

include/igl/face_areas.h

@@ -15,22 +15,28 @@ namespace igl
 {
   // Constructs a list of face areas of faces opposite each index in a tet list
   //
-  // Templates:
-  //   DerivedV derived from vertex positions matrix type: i.e. MatrixXd
-  //   DerivedT derived from face indices matrix type: i.e. MatrixXi
-  //   DerivedL derived from edge lengths matrix type: i.e. MatrixXd
   // Inputs:
-  //   V  eigen matrix #V by 3
-  //   T  #T by 3 list of mesh faces (must be triangles)
+  //   V  #V by 3 list of mesh vertex positions
+  //   T  #T by 3 list of tet mesh indices into V
   // Outputs:
-  //   E #E by 2 list of edges in no particular order
+  //   A   #T by 4 list of face areas corresponding to faces opposite vertices
+  //     0,1,2,3
   //
-  // See also: adjacency_matrix
   template <typename DerivedV, typename DerivedT, typename DerivedA>
   IGL_INLINE void face_areas(
     const Eigen::PlainObjectBase<DerivedV>& V,
     const Eigen::PlainObjectBase<DerivedT>& T,
     Eigen::PlainObjectBase<DerivedA>& A);
+  // Compute tet-mesh face areas from edge lengths.
+  //
+  // Inputs:
+  //   L  #T by 6 list of tet-mesh edge lengths corresponding to edges
+  //     [3 0],[3 1],[3 2],[1 2],[2 0],[0 1]
+  // Outputs:
+  //   A   #T by 4 list of face areas corresponding to faces opposite vertices 
+  //     0,1,2,3: i.e. made of edges [123],[024],[015],[345]
+  //    
+  //
   template <typename DerivedL, typename DerivedA>
   IGL_INLINE void face_areas(
     const Eigen::PlainObjectBase<DerivedL>& L,

include/igl/in_element.cpp

@@ -41,7 +41,6 @@ IGL_INLINE void igl::in_element(
 {
   using namespace std;
   using namespace Eigen;
-  using namespace igl;
   const int Qr = Q.rows();
   std::vector<Triplet<Scalar> > IJV;
   IJV.reserve(Qr);

include/igl/lim/lim.cpp

@@ -8,9 +8,8 @@
 #include "lim.h"
 #include <LIMSolverInterface.h>
 
-using namespace igl::lim;
 
-IGL_INLINE State lim(
+IGL_INLINE igl::lim::State igl::lim::lim(
   Eigen::Matrix<double,Eigen::Dynamic,3>& vertices,
   const Eigen::Matrix<double,Eigen::Dynamic,3>& initialVertices,
   const Eigen::Matrix<int,Eigen::Dynamic,Eigen::Dynamic>& elements,
@@ -34,7 +33,7 @@ IGL_INLINE State lim(
     );
 }
 
-IGL_INLINE State igl::lim::lim(
+IGL_INLINE igl::lim::State igl::lim::lim(
   Eigen::Matrix<double,Eigen::Dynamic,3>& vertices,
   const Eigen::Matrix<double,Eigen::Dynamic,3>& initialVertices,
   const Eigen::Matrix<int,Eigen::Dynamic,Eigen::Dynamic>& elements,
@@ -68,7 +67,7 @@ IGL_INLINE State igl::lim::lim(
     );
 }
 
-IGL_INLINE State igl::lim::lim(
+IGL_INLINE igl::lim::State igl::lim::lim(
   Eigen::Matrix<double,Eigen::Dynamic,3>& vertices,
   const Eigen::Matrix<double,Eigen::Dynamic,3>& initialVertices,
   const Eigen::Matrix<int,Eigen::Dynamic,Eigen::Dynamic>& elements,
@@ -96,7 +95,7 @@ IGL_INLINE State igl::lim::lim(
     );
 }
 
-IGL_INLINE State igl::lim::lim(
+IGL_INLINE igl::lim::State igl::lim::lim(
   Eigen::Matrix<double,Eigen::Dynamic,3>& vertices,
   const Eigen::Matrix<double,Eigen::Dynamic,3>& initialVertices,
   const Eigen::Matrix<int,Eigen::Dynamic,Eigen::Dynamic>& elements,

include/igl/list_to_matrix.cpp

@@ -134,5 +134,5 @@ template bool igl::list_to_matrix<bool, Eigen::Array<bool, -1, 1, 0, -1, 1> >(st
 template bool igl::list_to_matrix<bool, Eigen::Array<bool, -1, -1, 0, -1, -1> >(std::vector<std::vector<bool, std::allocator<bool> >, std::allocator<std::vector<bool, std::allocator<bool> > > > const&, Eigen::PlainObjectBase<Eigen::Array<bool, -1, -1, 0, -1, -1> >&);
 template bool igl::list_to_matrix<bool, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::vector<bool, std::allocator<bool> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template bool igl::list_to_matrix<bool, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::vector<std::vector<bool, std::allocator<bool> >, std::allocator<std::vector<bool, std::allocator<bool> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-template bool igl::list_to_matrix<double, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(std::__1::vector<double, std::__1::allocator<double> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template bool igl::list_to_matrix<double, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(std::vector<double, std::allocator<double> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/mosek/mosek_quadprog.cpp

@@ -281,7 +281,6 @@ IGL_INLINE bool igl::mosek::mosek_quadprog(
 {
   using namespace Eigen;
   using namespace std;
-  using namespace igl;
 
   // Q should be square
   assert(Q.rows() == Q.cols());

include/igl/opengl2/MouseController.h

@@ -293,7 +293,6 @@ inline bool igl::opengl2::MouseController::up(const int x, const int y)
 
 inline void igl::opengl2::MouseController::draw() const
 {
-  using namespace igl;
   if(any_selection())
   {
     m_widget.draw();
@@ -342,7 +341,6 @@ inline void igl::opengl2::MouseController::set_selection_from_last_drag(
 {
   using namespace Eigen;
   using namespace std;
-  using namespace igl;
   m_rotations_at_selection = m_rotations;
   assert(BE.rows() == P.rows());
   m_selection = VectorXb::Zero(BE.rows());
@@ -378,7 +376,6 @@ inline void igl::opengl2::MouseController::set_selection(
     const Eigen::VectorXi & P,
     const Eigen::VectorXi & RP)
 {
-  using namespace igl;
   using namespace Eigen;
   using namespace std;
   vector<Quaterniond,aligned_allocator<Quaterniond> > & vQ = 

include/igl/opengl2/lens_flare.cpp

@@ -113,7 +113,6 @@ IGL_INLINE void igl::opengl2::lens_flare_draw(
   glBlendFunc(GL_ONE, GL_ONE);
 
   using namespace Eigen;
-  using namespace igl;
   using namespace std;
 
   //// view_dir  direction from eye to position is is looking at

include/igl/piecewise_constant_winding_number.cpp

@@ -68,3 +68,7 @@ IGL_INLINE bool igl::piecewise_constant_winding_number(
   unique_edge_map(F, E, uE, EMAP, uE2E);
   return piecewise_constant_winding_number(F,uE,uE2E);
 }
+
+#ifdef IGL_STATIC_LIBRARY
+template bool igl::piecewise_constant_winding_number<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, unsigned long>(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<unsigned long, std::allocator<unsigned long> >, std::allocator<std::vector<unsigned long, std::allocator<unsigned long> > > > const&);
+#endif

include/igl/png/texture_from_file.cpp

@@ -19,7 +19,6 @@
 
 IGL_INLINE bool igl::png::texture_from_file(const std::string filename, GLuint & id)
 {
-  using namespace igl;
   using namespace igl::opengl;
   using namespace std;
   // dirname, basename, extension and filename

include/igl/ray_mesh_intersect.cpp

@@ -18,7 +18,6 @@ IGL_INLINE bool igl::ray_mesh_intersect(
   std::vector<igl::Hit> & hits)
 {
   using namespace Eigen;
-  using namespace igl;
   using namespace std;
   // Should be but can't be const 
   Vector3d s_d = s.template cast<double>();

include/igl/remove_duplicate_vertices.cpp

@@ -55,14 +55,6 @@ IGL_INLINE void igl::remove_duplicate_vertices(
   using namespace Eigen;
   using namespace std;
   remove_duplicate_vertices(V,epsilon,SV,SVI,SVJ);
-  // remap faces
-// #ifndef _WIN32
-//   SF = F.unaryExpr(bind1st(mem_fun( 
-//     static_cast<VectorXi::Scalar&(VectorXi::*)(VectorXi::Index)>
-//       (&VectorXi::operator())), &SVJ)).eval();
-// #else
-  // Why doesn't the above compile on windows?
-  // Daniele: it also does not compile with CLANG
   SF.resize(F.rows(),F.cols());
   for(int f = 0;f<F.rows();f++)
   {
@@ -71,7 +63,6 @@ IGL_INLINE void igl::remove_duplicate_vertices(
       SF(f,c) = SVJ(F(f,c));
     }
   }
-// #endif
 }
 
 #ifdef IGL_STATIC_LIBRARY

include/igl/signed_distance.cpp

@@ -303,7 +303,6 @@ IGL_INLINE void igl::signed_distance_winding_number(
 {
   using namespace Eigen;
   using namespace std;
-  using namespace igl;
   sqrd = tree.squared_distance(V,F,q,i,c);
   const double w = hier.winding_number(q.transpose());
   s = 1.-2.*w;
@@ -321,7 +320,6 @@ IGL_INLINE void igl::signed_distance_winding_number(
 {
   using namespace Eigen;
   using namespace std;
-  using namespace igl;
   sqrd = tree.squared_distance(V,F,q,i,c);
   double w;
   winding_number_2(V.data(), V.rows(), F.data(), F.rows(), q.data(), 1, &w);

include/igl/slice.cpp

@@ -267,4 +267,7 @@ template void igl::slice<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix
 template void igl::slice<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<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::slice<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> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, Eigen::Matrix<int, -1, -1, 0, -1, -1>&);
 template void igl::slice<Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >, Eigen::Matrix<long, -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> > const&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::slice<Eigen::PlainObjectBase<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> > >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+template void igl::slice<Eigen::PlainObjectBase<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> > >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::slice<Eigen::PlainObjectBase<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> > >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/snap_to_fixed_up.cpp

@@ -7,19 +7,26 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "snap_to_fixed_up.h"
 
+template <typename Qtype>
 IGL_INLINE void igl::snap_to_fixed_up(
-  const Eigen::Quaterniond & q,
-  Eigen::Quaterniond & s)
+  const Eigen::Quaternion<Qtype> & q,
+  Eigen::Quaternion<Qtype> & s)
 {
   using namespace Eigen;
-  const Vector3d up = q.matrix() * Vector3d(0,1,0);
-  Vector3d proj_up(0,up(1),up(2));
+  typedef Eigen::Matrix<Qtype,3,1> Vector3Q;
+  const Vector3Q up = q.matrix() * Vector3Q(0,1,0);
+  Vector3Q proj_up(0,up(1),up(2));
   if(proj_up.norm() == 0)
   {
-    proj_up = Vector3d(0,1,0);
+    proj_up = Vector3Q(0,1,0);
   }
   proj_up.normalize();
-  Quaterniond dq;
-  dq = Quaterniond::FromTwoVectors(up,proj_up);
+  Quaternion<Qtype> dq;
+  dq = Quaternion<Qtype>::FromTwoVectors(up,proj_up);
   s = dq * q;
 }
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instanciations
+template void igl::snap_to_fixed_up<float>(Eigen::Quaternion<float, 0> const&, Eigen::Quaternion<float, 0>&);
+#endif

include/igl/snap_to_fixed_up.h

@@ -24,9 +24,10 @@ namespace igl
   //   s the resulting rotation (as quaternion)
   //
   // See also: two_axis_valuator_fixed_up
+  template <typename Qtype>
   IGL_INLINE void snap_to_fixed_up(
-    const Eigen::Quaterniond & q,
-    Eigen::Quaterniond & s);
+    const Eigen::Quaternion<Qtype> & q,
+    Eigen::Quaternion<Qtype> & s);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/unproject_in_mesh.cpp

@@ -24,7 +24,6 @@ template < typename Derivedobj>
       Eigen::PlainObjectBase<Derivedobj> & obj,
       std::vector<igl::Hit > & hits)
 {
-  using namespace igl;
   using namespace std;
   using namespace Eigen;
   Vector3f s,dir;
@@ -65,7 +64,6 @@ template < typename DerivedV, typename DerivedF, typename Derivedobj>
       Eigen::PlainObjectBase<Derivedobj> & obj,
       std::vector<igl::Hit > & hits)
 {
-  using namespace igl;
   using namespace std;
   using namespace Eigen;
   const auto & shoot_ray = [&V,&F](

include/igl/unproject_ray.cpp

@@ -23,7 +23,6 @@ IGL_INLINE void igl::unproject_ray(
   Eigen::PlainObjectBase<Deriveds> & s,
   Eigen::PlainObjectBase<Deriveddir> & dir)
 {
-  using namespace igl;
   using namespace std;
   using namespace Eigen;
   // Source and direction on screen

include/igl/unzip_corners.cpp

@@ -0,0 +1,41 @@
+#include "unzip_corners.h"
+
+#include "unique.h"
+#include "slice.h"
+
+template < typename DerivedA, typename DerivedU, typename DerivedG, typename DerivedJ >
+IGL_INLINE void igl::unzip_corners(
+  const std::vector<std::reference_wrapper<DerivedA> > & A,
+  Eigen::PlainObjectBase<DerivedU> & U,
+  Eigen::PlainObjectBase<DerivedG> & G,
+  Eigen::PlainObjectBase<DerivedJ> & J)
+{
+  if(A.size() == 0)
+  {
+    U.resize(0,0);
+    G.resize(0,3);
+    J.resize(0,0);
+    return;
+  }
+  const size_t num_a = A.size();
+  const typename DerivedA::Index m = A[0].get().rows();
+  DerivedU C(m*3,num_a);
+  for(int a = 0;a<num_a;a++)
+  {
+    assert(A[a].get().rows() == m && "All attributes should be same size");
+    assert(A[a].get().cols() == 3 && "Must be triangles");
+    C.block(0*m,a,m,1) = A[a].get().col(0);
+    C.block(1*m,a,m,1) = A[a].get().col(1);
+    C.block(2*m,a,m,1) = A[a].get().col(2);
+  }
+  DerivedJ I;
+  igl::unique_rows(C,U,I,J);
+  G.resize(m,3);
+  for(int f = 0;f<m;f++)
+  {
+    for(int c = 0;c<3;c++)
+    {
+      G(f,c) = J(f+c*m);
+    }
+  }
+}

include/igl/unzip_corners.h

@@ -0,0 +1,54 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_UNZIP_CORNERS_H
+#define IGL_UNZIP_CORNERS_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+#include <functional>
+
+namespace igl
+{
+  // UNZIP_CORNERS Given a triangle mesh where corners of each triangle index
+  // different matrices of attributes (e.g. read from an OBJ file), unzip the
+  // corners into unique efficiently: attributes become properly vertex valued
+  // (usually creating greater than #V but less than #F*3 vertices).
+  //
+  // To pass a list of attributes this function takes an std::vector of
+  // std::reference_wrapper of an Eigen::... type. This allows you to use list
+  // initializers **without** incurring a copy, but means you'll need to
+  // provide the derived type of A as an explicit template parameter:
+  //
+  //     unzip_corners<Eigen::MatrixXi>({F,FTC,FN},U,G,J);
+  //
+  // Inputs:
+  //   A  #A list of #F by 3 attribute indices, typically {F,FTC,FN}
+  // Outputs:
+  //   U  #U by #A list of indices into each attribute for each unique mesh
+  //     vertex: U(v,a) is the attribute index of vertex v in attribute a.
+  //   G  #F by 3 list of triangle indices into U
+  // Example:
+  //   [V,F,TC,FTC] = readOBJ('~/Downloads/kiwis/kiwi.obj');
+  //   [U,G] = unzip_corners(cat(3,F,FTC));
+  //   % display mesh
+  //   tsurf(G,V(U(:,1),:));
+  //   % display texture coordinates
+  //   tsurf(G,TC(U(:,2),:));
+  //
+  template < typename DerivedA, typename DerivedU, typename DerivedG, typename DerivedJ>
+  IGL_INLINE void unzip_corners(
+    const std::vector<std::reference_wrapper<DerivedA> > & A,
+    Eigen::PlainObjectBase<DerivedU> & U,
+    Eigen::PlainObjectBase<DerivedG> & G,
+    Eigen::PlainObjectBase<DerivedJ> & J);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "unzip_corners.cpp"
+#endif
+#endif

include/igl/viewer/Viewer.cpp

@@ -458,6 +458,69 @@ namespace viewer
     for (unsigned int i = 0; i<plugins.size(); ++i)
       if (plugins[i]->key_pressed(unicode_key, modifiers))
         return true;
+    const std::string usage(R"(igl::viewer::Viewer:
+A,a  Toggle animation (tight draw loop)
+I,i  Toggle invert normals
+L,l  Toggle wireframe
+O,o  Toggle orthographic/perspective projection
+T,t  Toggle filled faces
+Z    Snap to canonical view
+[,]  Toggle between rotation control types (e.g. trackball, two-axis valuator
+     with fixed up)
+)");
+    switch(unicode_key)
+    {
+      case 'A':
+      case 'a':
+      {
+        core.is_animating = !core.is_animating;
+        return true;
+      }
+      case 'I':
+      case 'i':
+      {
+        data.dirty |= ViewerData::DIRTY_NORMAL;
+        core.invert_normals = !core.invert_normals;
+        return true;
+      }
+      case 'L':
+      case 'l':
+      {
+        core.show_lines = !core.show_lines;
+        return true;
+      }
+      case 'O':
+      case 'o':
+      {
+        core.orthographic = !core.orthographic;
+        return true;
+      }
+      case 'T':
+      case 't':
+      {
+        core.show_faces = !core.show_faces;
+        return true;
+      }
+      case 'Z':
+      {
+        snap_to_canonical_quaternion();
+        return true;
+      }
+      case '[':
+      case ']':
+      {
+        if(core.rotation_type == ViewerCore::ROTATION_TYPE_TRACKBALL)
+        {
+          core.set_rotation_type(
+            ViewerCore::ROTATION_TYPE_TWO_AXIS_VALUATOR_FIXED_UP);
+        }else
+        {
+          core.set_rotation_type(ViewerCore::ROTATION_TYPE_TRACKBALL);
+        }
+        return true;
+      }
+      default: break;//do nothing
+    }
     return false;
   }
 

include/igl/viewer/Viewer.h

@@ -127,6 +127,8 @@ namespace viewer
     IGL_INLINE void open_dialog_save_mesh();
 
     // C++-style functions
+    //
+    // Returns **true** if action should be cancelled.
     std::function<bool(Viewer& viewer)> callback_init;
     std::function<bool(Viewer& viewer)> callback_pre_draw;
     std::function<bool(Viewer& viewer)> callback_post_draw;

include/igl/viewer/ViewerCore.cpp

@@ -8,6 +8,7 @@
 
 #include "ViewerCore.h"
 #include <igl/quat_to_mat.h>
+#include <igl/snap_to_fixed_up.h>
 #include <igl/look_at.h>
 #include <igl/frustum.h>
 #include <igl/ortho.h>
@@ -457,6 +458,20 @@ IGL_INLINE void igl::viewer::ViewerCore::draw_buffer(ViewerData& data,
   glDeleteFramebuffers(1, &frameBuffer);
 }
 
+IGL_INLINE void igl::viewer::ViewerCore::set_rotation_type(
+  const igl::viewer::ViewerCore::RotationType & value)
+{
+  using namespace Eigen;
+  using namespace std;
+  const RotationType old_rotation_type = rotation_type;
+  rotation_type = value;
+  if(rotation_type == ROTATION_TYPE_TWO_AXIS_VALUATOR_FIXED_UP &&
+    old_rotation_type != ROTATION_TYPE_TWO_AXIS_VALUATOR_FIXED_UP)
+  {
+    snap_to_fixed_up(Quaternionf(trackball_angle),trackball_angle);
+  }
+}
+
 
 IGL_INLINE igl::viewer::ViewerCore::ViewerCore()
 {
@@ -473,7 +488,7 @@ IGL_INLINE igl::viewer::ViewerCore::ViewerCore()
 
   // Default trackball
   trackball_angle = Eigen::Quaternionf::Identity();
-  rotation_type = ViewerCore::ROTATION_TYPE_TRACKBALL;
+  set_rotation_type(ViewerCore::ROTATION_TYPE_TWO_AXIS_VALUATOR_FIXED_UP);
 
   // Defalut model viewing parameters
   model_zoom = 1.0f;

include/igl/viewer/ViewerCore.h

@@ -74,13 +74,23 @@ public:
 
   // Draw everything
   IGL_INLINE void draw(ViewerData& data, OpenGL_state& opengl, bool update_matrices = true);
-  IGL_INLINE void draw_buffer(ViewerData& data,
-                              OpenGL_state& opengl,
-                              bool update_matrices,
-                              Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic>& R,
-                              Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic>& G,
-                              Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic>& B,
-                              Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic>& A);
+  IGL_INLINE void draw_buffer(
+    ViewerData& data,
+    OpenGL_state& opengl,
+    bool update_matrices,
+    Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic>& R,
+    Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic>& G,
+    Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic>& B,
+    Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic>& A);
+
+  // Trackball angle (quaternion)
+  enum RotationType
+  {
+    ROTATION_TYPE_TRACKBALL = 0,
+    ROTATION_TYPE_TWO_AXIS_VALUATOR_FIXED_UP = 1,
+    NUM_ROTATION_TYPES = 2
+  };
+  IGL_INLINE void set_rotation_type(const RotationType & value);
 
   // ------------------- Properties
 
@@ -100,13 +110,8 @@ public:
   Eigen::Vector3f light_position;
   float lighting_factor;
 
-  // Trackball angle (quaternion)
-  enum RotationType
-  {
-    ROTATION_TYPE_TRACKBALL = 0,
-    ROTATION_TYPE_TWO_AXIS_VALUATOR_FIXED_UP = 1,
-    NUM_ROTATION_TYPES = 2
-  } rotation_type;
+  RotationType rotation_type;
+
   Eigen::Quaternionf trackball_angle;
 
   // Model viewing parameters

include/igl/viewer/ViewerData.h

@@ -60,12 +60,28 @@ public:
   // Inputs:
   //   C  #V|#F|1 by 3 list of colors
   IGL_INLINE void set_colors(const Eigen::MatrixXd &C);
+  // Set per-vertex UV coordinates
+  //
+  // Inputs:
+  //   UV  #V by 2 list of UV coordinates (indexed by F)
   IGL_INLINE void set_uv(const Eigen::MatrixXd& UV);
+  // Set per-corner UV coordinates
+  //
+  // Inputs:
+  //   UV_V  #UV by 2 list of UV coordinates
+  //   UV_F  #F by 3 list of UV indices into UV_V
   IGL_INLINE void set_uv(const Eigen::MatrixXd& UV_V, const Eigen::MatrixXi& UV_F);
+  // Set the texture associated with the mesh.
+  //
+  // Inputs:
+  //   R  width by height image matrix of red channel
+  //   G  width by height image matrix of green channel
+  //   B  width by height image matrix of blue channel
+  //
   IGL_INLINE void set_texture(
-                    const Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic>& R,
-                    const Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic>& G,
-                    const Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic>& B);
+    const Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic>& R,
+    const Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic>& G,
+    const Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic>& B);
 
   // Sets points given a list of point vertices. In constrast to `set_points`
   // this will (purposefully) clober existing points.