Merge remote-tracking branch 'upstream/master'

Former-commit-id: 1be72027e8c20e5f7fb8ba58965027f688cd78e2

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

@@ -23,161 +23,56 @@
 //#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)
-    -> 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;
-
-        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) -> void {
-        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 = []() -> double
   {
@@ -191,127 +86,168 @@ IGL_INLINE size_t igl::copyleft::cgal::extract_cells(
       << 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;
-
-    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
-    log_time("extract_single_component_cells");
-#endif
-
-    std::vector<std::vector<std::vector<Index > > > TT,_1;
-    igl::triangle_triangle_adjacency(E, EMAP, uE2E, false, TT, _1);
-#ifdef EXTRACT_CELLS_DEBUG
-    log_time("compute_face_adjacency");
-#endif
-
-    Eigen::VectorXi C, counts;
-    igl::facet_components(TT, C, counts);
+  const size_t num_faces = F.rows();
+  typedef typename DerivedF::Scalar Index;
+  const size_t num_patches = P.maxCoeff()+1;
+
+  // 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");
+
+  // 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");
+
+  // 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();
+  // 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());
+  }
+
+  // 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
-    log_time("form_components");
+  log_time("outer_facet_per_component");
 #endif
 
-    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());
-    }
-
-    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]);
+  // 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));
+        }
+      }
     }
-#ifdef EXTRACT_CELLS_DEBUG
-    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();
+    // 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));
     };
-    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));
-          }
-        }
-        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;
-
-            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 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);
-
-            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);
-                }
-            }
+    
+    // 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);
+      }
+
+      const size_t num_candidate_comps = candidate_comps.size();
+      if (num_candidate_comps == 0) continue;
+
+      // 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]);
+      }
+
+      // 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])
+        {
+          // Awkardly, this doesn't seem to imply that component j is inside
+          // component i. Consider:
+          //
+          //    ________
+          //   /    ___ \       __
+          //  |  2 | i | |  0  / j\
+          //  |     ---  |     \__/
+          //   \________/
+          //
+          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
     log_time("nested_relationship");
 #endif
@@ -387,36 +323,239 @@ IGL_INLINE size_t igl::copyleft::cgal::extract_cells(
 }
 
 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;
-
-    Eigen::MatrixXi E, uE;
-    Eigen::VectorXi EMAP;
-    std::vector<std::vector<size_t> > uE2E;
-    igl::unique_edge_map(F, E, uE, EMAP, uE2E);
-
-    Eigen::VectorXi P;
-    //const size_t num_patches = 
-      igl::extract_manifold_patches(F, EMAP, uE2E, P);
+  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;
+
+    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;
+  };
 
-    DerivedC per_patch_cells;
-    const size_t num_cells =
-        igl::copyleft::cgal::extract_cells(V, F, P, E, uE, uE2E, EMAP, per_patch_cells);
+  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);
+      }
+    }
+  }
+
+  // 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);
+
+  // 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,96 +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.
-        //          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);
+    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.
-        //          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).  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/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/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/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/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/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::__1::vector<std::__1::vector<unsigned long, std::__1::allocator<unsigned long> >, std::__1::allocator<std::__1::vector<unsigned long, std::__1::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/project_to_line_segment.cpp

@@ -28,7 +28,7 @@ IGL_INLINE void igl::project_to_line_segment(
 #pragma omp parallel for if (np>10000)
   for(int p = 0;p<np;p++)
   {
-    const Eigen::PlainObjectBase<DerivedP> Pp = P.row(p);
+    const DerivedP Pp = P.row(p);
     if(t(p)<0)
     {
       sqrD(p) = (Pp-S).squaredNorm();

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