Merge pull request #1487 from libigl/connected_components

Connected components

include/igl/connected_components.cpp

@@ -0,0 +1,61 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2020 Alec Jacobson <[email protected]>
+//
+// 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 "connected_components.h"
+#include <queue>
+
+template < typename Atype, typename DerivedC, typename DerivedK>
+IGL_INLINE int igl::connected_components(
+  const Eigen::SparseMatrix<Atype> & A,
+  Eigen::PlainObjectBase<DerivedC> & C,
+  Eigen::PlainObjectBase<DerivedK> & K)
+{
+  typedef typename Eigen::SparseMatrix<Atype>::Index Index;
+  const auto m = A.rows();
+  assert(A.cols() == A.rows() && "A should be square");
+  // 1.1 sec
+  // m  means not yet visited
+  C.setConstant(m,1,m);
+  // Could use amortized dynamic array but didn't see real win.
+  K.setZero(m,1);
+  typename DerivedC::Scalar c = 0;
+  for(Eigen::Index f = 0;f<m;f++)
+  {
+    // already seen
+    if(C(f)<m) continue;
+    // start bfs
+    std::queue<Index> Q;
+    Q.push(f);
+    while(!Q.empty())
+    {
+      const Index g = Q.front();
+      Q.pop();
+      // already seen
+      if(C(g)<m) continue;
+      // see it
+      C(g) = c;
+      K(c)++;
+      for(typename Eigen::SparseMatrix<Atype>::InnerIterator it (A,g); it; ++it)
+      {
+        const Index n = it.index();
+        // already seen
+        if(C(n)<m) continue;
+        Q.push(n);
+      }
+    }
+    c++;
+  }
+  K.conservativeResize(c,1);
+  return c;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template int igl::connected_components<int, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::SparseMatrix<int, 0, int> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/connected_components.h

@@ -0,0 +1,35 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2020 Alec Jacobson <[email protected]>
+//
+// 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_CONNECTED_COMPONENTS_H
+#define IGL_CONNECTED_COMPONENTS_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+namespace igl
+{
+  // Determine the connected components of a graph described by the input
+  // adjacency matrix (similar to MATLAB's graphconncomp).
+  //
+  // Inputs:
+  //    A  #A by #A adjacency matrix (treated as describing an undirected graph)
+  // Outputs:
+  //    C  #A list of component indices into [0,#K-1]
+  //    K  #K list of sizes of each component
+  // Returns number of connected components
+  template < typename Atype, typename DerivedC, typename DerivedK>
+  IGL_INLINE int connected_components(
+    const Eigen::SparseMatrix<Atype> & A,
+    Eigen::PlainObjectBase<DerivedC> & C,
+    Eigen::PlainObjectBase<DerivedK> & K);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "connected_components.cpp"
+#endif
+
+#endif

include/igl/cumsum.cpp

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2013 Alec Jacobson <[email protected]>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "cumsum.h"
 #include <numeric>
@@ -14,43 +14,65 @@ IGL_INLINE void igl::cumsum(
   const Eigen::MatrixBase<DerivedX > & X,
   const int dim,
   Eigen::PlainObjectBase<DerivedY > & Y)
+{
+  return cumsum(X,dim,false,Y);
+}
+
+template <typename DerivedX, typename DerivedY>
+IGL_INLINE void igl::cumsum(
+  const Eigen::MatrixBase<DerivedX > & X,
+  const int dim,
+  const bool zero_prefix,
+  Eigen::PlainObjectBase<DerivedY > & Y)
 {
   using namespace Eigen;
   using namespace std;
-  Y.resizeLike(X);
+  Y.resize(
+    X.rows()+(zero_prefix&&dim==1?1:0),
+    X.cols()+(zero_prefix&&dim==2?1:0));
   // get number of columns (or rows)
-  int num_outer = (dim == 1 ? X.cols() : X.rows() );
+  Eigen::Index num_outer = (dim == 1 ? X.cols() : X.rows() );
   // get number of rows (or columns)
-  int num_inner = (dim == 1 ? X.rows() : X.cols() );
+  Eigen::Index num_inner = (dim == 1 ? X.rows() : X.cols() );
   // This has been optimized so that dim = 1 or 2 is roughly the same cost.
   // (Optimizations assume ColMajor order)
   if(dim == 1)
   {
+    if(zero_prefix)
+    {
+      Y.row(0).setConstant(0);
+    }
 #pragma omp parallel for
-    for(int o = 0;o<num_outer;o++)
+    for(Eigen::Index o = 0;o<num_outer;o++)
     {
       typename DerivedX::Scalar sum = 0;
-      for(int i = 0;i<num_inner;i++)
+      for(Eigen::Index i = 0;i<num_inner;i++)
       {
         sum += X(i,o);
-        Y(i,o) = sum;
+        const Eigen::Index yi = zero_prefix?i+1:i;
+        Y(yi,o) = sum;
       }
     }
   }else
   {
-    for(int i = 0;i<num_inner;i++)
+    if(zero_prefix)
+    {
+      Y.col(0).setConstant(0);
+    }
+    for(Eigen::Index i = 0;i<num_inner;i++)
     {
+      const Eigen::Index yi = zero_prefix?i+1:i;
       // Notice that it is *not* OK to put this above the inner loop
       // Though here it doesn't seem to pay off...
 //#pragma omp parallel for
-      for(int o = 0;o<num_outer;o++)
+      for(Eigen::Index o = 0;o<num_outer;o++)
       {
         if(i == 0)
         {
-          Y(o,i) = X(o,i);
+          Y(o,yi) = X(o,i);
         }else
         {
-          Y(o,i) = Y(o,i-1) + X(o,i);
+          Y(o,yi) = Y(o,yi-1) + X(o,i);
         }
       }
     }
@@ -59,11 +81,8 @@ IGL_INLINE void igl::cumsum(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-// generated by autoexplicit.sh
 template void igl::cumsum<Eigen::Matrix<float, -1, 1, 0, -1, 1>, Eigen::Matrix<float, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> >&);
-// generated by autoexplicit.sh
 template void igl::cumsum<Eigen::Matrix<double, 4, 1, 0, 4, 1>, Eigen::Matrix<double, 4, 1, 0, 4, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, 4, 1, 0, 4, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, 4, 1, 0, 4, 1> >&);
-// generated by autoexplicit.sh
 template void igl::cumsum<Eigen::Matrix<double, 1, 4, 1, 1, 4>, Eigen::Matrix<double, 1, 4, 1, 1, 4> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 4, 1, 1, 4> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 4, 1, 1, 4> >&);
 template void igl::cumsum<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::cumsum<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
@@ -71,6 +90,7 @@ template void igl::cumsum<Eigen::Matrix<double, 3, 1, 0, 3, 1>, Eigen::Matrix<do
 template void igl::cumsum<Eigen::Matrix<unsigned long, 2, 1, 0, 2, 1>, Eigen::Matrix<unsigned long, 2, 1, 0, 2, 1> >(Eigen::MatrixBase<Eigen::Matrix<unsigned long, 2, 1, 0, 2, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<unsigned long, 2, 1, 0, 2, 1> >&);
 template void igl::cumsum<Eigen::Matrix<unsigned long, -1, 1, 0, -1, 1>, Eigen::Matrix<unsigned long, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<unsigned long, -1, 1, 0, -1, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<unsigned long, -1, 1, 0, -1, 1> >&);
 template void igl::cumsum<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::cumsum<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #ifdef WIN32
 template void igl::cumsum<class Eigen::Matrix<unsigned __int64, -1, 1, 0, -1, 1>, class Eigen::Matrix<unsigned __int64, -1, 1, 0, -1, 1>>(class Eigen::MatrixBase<class Eigen::Matrix<unsigned __int64, -1, 1, 0, -1, 1>> const &, int, class Eigen::PlainObjectBase<class Eigen::Matrix<unsigned __int64, -1, 1, 0, -1, 1>> &);
 template void igl::cumsum<class Eigen::Matrix<unsigned __int64, 2, 1, 0, 2, 1>, class Eigen::Matrix<unsigned __int64, 2, 1, 0, 2, 1>>(class Eigen::MatrixBase<class Eigen::Matrix<unsigned __int64, 2, 1, 0, 2, 1>> const &, int, class Eigen::PlainObjectBase<class Eigen::Matrix<unsigned __int64, 2, 1, 0, 2, 1>> &);

include/igl/cumsum.h

@@ -29,6 +29,25 @@ namespace igl
     const Eigen::MatrixBase<DerivedX > & X,
     const int dim,
     Eigen::PlainObjectBase<DerivedY > & Y);
+  // Computes a cumulative sum of the columns of [0;X]
+  //
+  // Inputs:
+  //   X  m by n Matrix to be cumulatively summed.
+  //   dim  dimension to take cumulative sum (1 or 2)
+  //   zero_prefix whe
+  // Output:
+  //   if zero_prefix == false
+  //     Y  m by n Matrix containing cumulative sum
+  //   else
+  //     Y  m+1 by n Matrix containing cumulative sum if dim=1
+  //     or 
+  //     Y  m by n+1 Matrix containing cumulative sum if dim=2
+  template <typename DerivedX, typename DerivedY>
+  IGL_INLINE void cumsum(
+    const Eigen::MatrixBase<DerivedX > & X,
+    const int dim,
+    const bool zero_prefix,
+    Eigen::PlainObjectBase<DerivedY > & Y);
   //template <typename DerivedX, typename DerivedY>
   //IGL_INLINE void cumsum(
   //  const Eigen::MatrixBase<DerivedX > & X,

include/igl/facet_adjacency_matrix.cpp

@@ -0,0 +1,58 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2020 Alec Jacobson <[email protected]>
+//
+// 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 "facet_adjacency_matrix.h"
+#include "unique_edge_map.h"
+
+template <typename DerivedF, typename Atype>
+IGL_INLINE void igl::facet_adjacency_matrix(
+  const Eigen::MatrixBase<DerivedF> & F, Eigen::SparseMatrix<Atype> & A)
+{
+  using namespace Eigen;
+  typedef typename DerivedF::Scalar Index;
+  const auto m = F.rows();
+  Eigen::Matrix<Index,Dynamic,1> EMAP,uEE,uEC;
+  Eigen::Matrix<Index,Dynamic,2> E,uE;
+  igl::unique_edge_map(F,E,uE,EMAP,uEC,uEE);
+  std::vector<Eigen::Triplet<Index> > AIJV;
+  AIJV.reserve(2*uE.rows());
+  const Eigen::Index nu = uE.rows();
+  for(Eigen::Index ue = 0;ue<nu;ue++)
+  {
+    // number of faces incident on this unique edge
+    const Eigen::Index mue = uEC(ue+1)-uEC(ue);
+    // base offset in uEE
+    const Eigen::Index uECue = uEC(ue);
+    assert(uECue < uEE.rows());
+    for(Eigen::Index i = 0;i<mue;i++)
+    {
+      const auto ii = uEE(uECue+i)%m;
+      for(Eigen::Index j = 0;j<mue;j++)
+      {
+        const auto jj = uEE(uECue+j)%m;
+        if(ii != jj){ AIJV.emplace_back(ii,jj,1);}
+      }
+    }
+  }
+  // Facet Adjacency matrix (with (arbitrary) >0
+  A.resize(m,m);
+  A.setFromTriplets(AIJV.begin(),AIJV.end());
+  // Set all non-zerro values to 1
+  for(Eigen::Index g = 0;g<A.outerSize();g++)
+  {
+    for(typename Eigen::SparseMatrix<Atype>::InnerIterator it (A,g); it; ++it)
+    {
+      if(it.value()) it.valueRef() = 1;
+    }
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::facet_adjacency_matrix<Eigen::Matrix<int, -1, -1, 0, -1, -1>, int>(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<int, 0, int>&);
+#endif

include/igl/facet_adjacency_matrix.h

@@ -0,0 +1,33 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2020 Alec Jacobson <[email protected]>
+//
+// 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_FACET_ADJACENCY_MATRIX_H
+#define IGL_FACET_ADJACENCY_MATRIX_H
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+#include "igl_inline.h"
+
+namespace igl
+{
+  // Construct a #F×#F adjacency matrix with A(i,j)>0 indicating that faces i and j
+  // share an edge.
+  //
+  // Inputs:
+  //   F  #F by 3 list of facets
+  // Outputs:
+  //   A  #F by #F adjacency matrix
+  template <typename DerivedF, typename Atype>
+  IGL_INLINE void facet_adjacency_matrix(
+    const Eigen::MatrixBase<DerivedF> & F,
+    Eigen::SparseMatrix<Atype> & A);
+};
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "facet_adjacency_matrix.cpp"
+#endif
+
+#endif

include/igl/facet_components.cpp

@@ -1,26 +1,29 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 //
 // Copyright (C) 2015 Alec Jacobson <[email protected]>
+// Copyright (C) 2020 Alec Jacobson <[email protected]>
 //
 // 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 "facet_components.h"
-#include <igl/triangle_triangle_adjacency.h>
+#include "triangle_triangle_adjacency.h"
+#include "facet_adjacency_matrix.h"
+#include "connected_components.h"
 #include <vector>
 #include <queue>
+
 template <typename DerivedF, typename DerivedC>
 IGL_INLINE void igl::facet_components(
   const Eigen::MatrixBase<DerivedF> & F,
   Eigen::PlainObjectBase<DerivedC> & C)
 {
-  using namespace std;
-  typedef typename DerivedF::Index Index;
-  vector<vector<vector<Index > > > TT;
-  vector<vector<vector<Index > > > TTi;
-  triangle_triangle_adjacency(F,TT,TTi);
-  Eigen::VectorXi counts;
-  return facet_components(TT,C,counts);
+  typedef typename DerivedF::Scalar Index;
+  Eigen::SparseMatrix<Index> A;
+  igl::facet_adjacency_matrix(F,A);
+  Eigen::Matrix<Index,Eigen::Dynamic,1> counts;
+  C = DerivedC::Zero(1,1);
+  connected_components(A,C,counts);
 }
 
 template <

include/igl/matlab/prepare_lhs.cpp

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2015 Alec Jacobson <[email protected]>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "prepare_lhs.h"
 #include <algorithm>
@@ -14,10 +14,10 @@ IGL_INLINE void igl::matlab::prepare_lhs_double(
 {
   using namespace std;
   using namespace Eigen;
-  const int m = V.rows();
-  const int n = V.cols();
+  const auto m = V.rows();
+  const auto n = V.cols();
   plhs[0] = mxCreateDoubleMatrix(m,n, mxREAL);
-  Eigen::Map< Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> > 
+  Eigen::Map< Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> >
     map(mxGetPr(plhs[0]),m,n);
   map = V.template cast<double>();
 }
@@ -29,11 +29,11 @@ IGL_INLINE void igl::matlab::prepare_lhs_logical(
 {
   using namespace std;
   using namespace Eigen;
-  const int m = V.rows();
-  const int n = V.cols();
+  const auto m = V.rows();
+  const auto n = V.cols();
   plhs[0] = mxCreateLogicalMatrix(m,n);
   mxLogical * Vp = static_cast<mxLogical*>(mxGetData(plhs[0]));
-  Eigen::Map< Eigen::Matrix<mxLogical,Eigen::Dynamic,Eigen::Dynamic> > 
+  Eigen::Map< Eigen::Matrix<mxLogical,Eigen::Dynamic,Eigen::Dynamic> >
     map(static_cast<mxLogical*>(mxGetData(plhs[0])),m,n);
   map = V.template cast<mxLogical>();
 }
@@ -54,8 +54,8 @@ IGL_INLINE void igl::matlab::prepare_lhs_double(
   mxArray *plhs[])
 {
   using namespace std;
-  const int m = M.rows();
-  const int n = M.cols();
+  const auto m = M.rows();
+  const auto n = M.cols();
   // THIS WILL NOT WORK FOR ROW-MAJOR
   assert(n==M.outerSize());
   const int nzmax = M.nonZeros();
@@ -84,6 +84,27 @@ IGL_INLINE void igl::matlab::prepare_lhs_double(
 
 }
 
+
+template <typename Vtype>
+IGL_INLINE void igl::prepare_lhs_double(
+  const std::vector<Vtype> & V,
+  mxArray *plhs[])
+{
+  plhs[0] = mxCreateCellMatrix(V.size(), 1);
+  for(int  i=0; i<V.size(); i++)
+  {
+    const auto m = V[i].rows();
+    const auto n = V[i].cols();
+    mxArray * ai = mxCreateDoubleMatrix(m,n, mxREAL);
+    Eigen::Map< Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> >
+      map(mxGetPr(ai),m,n);
+    map = V[i].template cast<double>();
+    mxSetCell(plhs[0],i,ai);
+  }
+}
+
+
+
 #ifdef IGL_STATIC_LIBRARY
 template void igl::matlab::prepare_lhs_double<Eigen::Matrix<double, 3, 3, 0, 3, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 3, 0, 3, 3> > const&, mxArray_tag**);
 template void igl::matlab::prepare_lhs_double<Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, mxArray_tag**);

include/igl/matlab/prepare_lhs.h

@@ -40,6 +40,11 @@ namespace igl
     IGL_INLINE void prepare_lhs_double(
       const Eigen::SparseMatrix<Vtype> & V,
       mxArray *plhs[]);
+    // Vector of matrices -> cell array of matrices
+    template <typename Vtype>
+    IGL_INLINE void prepare_lhs_double(
+      const std::vector<Vtype> & V,
+      mxArray *plhs[]);
   };
 }
 #ifndef IGL_STATIC_LIBRARY

include/igl/slice.h

@@ -22,7 +22,8 @@ namespace igl
   // Output:
   //   Y  #R by #C matrix
   //
-  // See also: slice_mask
+  // See also: slice_mask, and Eigen's unaryExpr
+  //   https://stackoverflow.com/a/49411587/148668
   template <
     typename TX,
     typename TY,

include/igl/unique_edge_map.cpp

@@ -8,8 +8,11 @@
 #include "unique_edge_map.h"
 #include "oriented_facets.h"
 #include "unique_simplices.h"
+#include "cumsum.h"
+#include "accumarray.h"
 #include <cassert>
 #include <algorithm>
+
 template <
   typename DerivedF,
   typename DerivedE,
@@ -22,6 +25,31 @@ IGL_INLINE void igl::unique_edge_map(
   Eigen::PlainObjectBase<DeriveduE> & uE,
   Eigen::PlainObjectBase<DerivedEMAP> & EMAP,
   std::vector<std::vector<uE2EType> > & uE2E)
+{
+  using namespace Eigen;
+  using namespace std;
+  unique_edge_map(F,E,uE,EMAP);
+  uE2E.resize(uE.rows());
+  // This does help a little
+  for_each(uE2E.begin(),uE2E.end(),[](vector<uE2EType > & v){v.reserve(2);});
+  const size_t ne = E.rows();
+  assert((size_t)EMAP.size() == ne);
+  for(uE2EType e = 0;e<(uE2EType)ne;e++)
+  {
+    uE2E[EMAP(e)].push_back(e);
+  }
+}
+
+template <
+  typename DerivedF,
+  typename DerivedE,
+  typename DeriveduE,
+  typename DerivedEMAP>
+IGL_INLINE void igl::unique_edge_map(
+  const Eigen::MatrixBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedE> & E,
+  Eigen::PlainObjectBase<DeriveduE> & uE,
+  Eigen::PlainObjectBase<DerivedEMAP> & EMAP)
 {
   using namespace Eigen;
   using namespace std;
@@ -33,23 +61,54 @@ IGL_INLINE void igl::unique_edge_map(
   // vs. O(log m)
   Matrix<typename DerivedEMAP::Scalar,Dynamic,1> IA;
   unique_simplices(E,uE,IA,EMAP);
-  uE2E.resize(uE.rows());
-  // This does help a little
-  for_each(uE2E.begin(),uE2E.end(),[](vector<uE2EType > & v){v.reserve(2);});
   assert((size_t)EMAP.size() == ne);
-  for(uE2EType e = 0;e<(uE2EType)ne;e++)
+}
+
+template <
+  typename DerivedF,
+  typename DerivedE,
+  typename DeriveduE,
+  typename DerivedEMAP,
+  typename DeriveduEC,
+  typename DeriveduEE>
+IGL_INLINE void igl::unique_edge_map(
+  const Eigen::MatrixBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedE> & E,
+  Eigen::PlainObjectBase<DeriveduE> & uE,
+  Eigen::PlainObjectBase<DerivedEMAP> & EMAP,
+  Eigen::PlainObjectBase<DeriveduEC> & uEC,
+  Eigen::PlainObjectBase<DeriveduEE> & uEE)
+{
+  // Avoid using uE2E
+  igl::unique_edge_map(F,E,uE,EMAP);
+  assert(EMAP.maxCoeff() < uE.rows());
+  // counts of each unique edge
+  typedef Eigen::Matrix<typename DeriveduEC::Scalar,Eigen::Dynamic,1> VectorXI;
+  VectorXI uEK;
+  igl::accumarray(EMAP,1,uEK);
+  assert(uEK.rows() == uE.rows());
+  // base offset in uEE
+  igl::cumsum(uEK,1,true,uEC);
+  assert(uEK.rows()+1 == uEC.rows());
+  // running inner offset in uEE
+  VectorXI uEO = VectorXI::Zero(uE.rows(),1);
+  // flat array of faces incide on each uE
+  uEE.resize(EMAP.rows(),1);
+  for(Eigen::Index e = 0;e<EMAP.rows();e++)
   {
-    uE2E[EMAP(e)].push_back(e);
+    const typename DerivedEMAP::Scalar ue = EMAP(e);
+    const typename DeriveduEC::Scalar i = uEC(ue)+ uEO(ue);
+    uEE(i) = e;
+    uEO(ue)++;
   }
+  assert( (uEK.array()==uEO.array()).all() );
 }
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-// generated by autoexplicit.sh
+template void igl::unique_edge_map<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::unique_edge_map<Eigen::Matrix<int, -1, 3, 1, -1, 3>, 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::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::vector<std::vector<unsigned long, std::allocator<unsigned long> >, std::allocator<std::vector<unsigned long, std::allocator<unsigned long> > > >&);
-// generated by autoexplicit.sh
 template void igl::unique_edge_map<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>, int>(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&);
-// generated by autoexplicit.sh
 template void igl::unique_edge_map<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::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, std::vector<std::vector<unsigned long, std::allocator<unsigned long> >, std::allocator<std::vector<unsigned long, std::allocator<unsigned long> > > >&);
 template void igl::unique_edge_map<Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, long>(Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&);
 template void igl::unique_edge_map<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, long>(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&);

include/igl/unique_edge_map.h

@@ -39,6 +39,37 @@ namespace igl
     Eigen::PlainObjectBase<DeriveduE> & uE,
     Eigen::PlainObjectBase<DerivedEMAP> & EMAP,
     std::vector<std::vector<uE2EType> > & uE2E);
+  template <
+    typename DerivedF,
+    typename DerivedE,
+    typename DeriveduE,
+    typename DerivedEMAP>
+  IGL_INLINE void unique_edge_map(
+    const Eigen::MatrixBase<DerivedF> & F,
+    Eigen::PlainObjectBase<DerivedE> & E,
+    Eigen::PlainObjectBase<DeriveduE> & uE,
+    Eigen::PlainObjectBase<DerivedEMAP> & EMAP);
+  // Outputs:
+  //   uEC  #uEC+1 list of cumulative counts of directed edges sharing each
+  //     unique edge so the uEC(i+1)-uEC(i) is the number of directed edges
+  //     sharing the ith unique edge.
+  //   uEE  #E list of indices into E, so that the consecutive segment of
+  //     indices uEE.segment(uEC(i),uEC(i+1)-uEC(i)) lists all directed edges
+  //     sharing the ith unique edge.
+  template <
+    typename DerivedF,
+    typename DerivedE,
+    typename DeriveduE,
+    typename DerivedEMAP,
+    typename DeriveduEC,
+    typename DeriveduEE>
+  IGL_INLINE void unique_edge_map(
+    const Eigen::MatrixBase<DerivedF> & F,
+    Eigen::PlainObjectBase<DerivedE> & E,
+    Eigen::PlainObjectBase<DeriveduE> & uE,
+    Eigen::PlainObjectBase<DerivedEMAP> & EMAP,
+    Eigen::PlainObjectBase<DeriveduEC> & uEC,
+    Eigen::PlainObjectBase<DeriveduEE> & uEE);
 
 }
 #ifndef IGL_STATIC_LIBRARY

tests/include/igl/facet_components.cpp

@@ -0,0 +1,36 @@
+#include <test_common.h>
+#include <igl/facet_components.h>
+
+TEST_CASE("facet_components: two_triangles", "[igl]")
+{
+  const Eigen::MatrixXi F1 = 
+    (Eigen::MatrixXi(2,3)<<
+     0,1,3,
+     0,3,2).finished();
+  Eigen::VectorXi C1;
+  igl::facet_components(F1,C1);
+  REQUIRE(C1(0) == 0);
+  REQUIRE(C1(0) == C1(1));
+
+  const Eigen::MatrixXi F2 = 
+    (Eigen::MatrixXi(2,3)<<
+     0,1,2,
+     3,4,5).finished();
+  Eigen::VectorXi C2;
+  igl::facet_components(F2,C2);
+  REQUIRE(C2(0) != C2(1));
+  REQUIRE(C2.minCoeff() == 0);
+  REQUIRE(C2.maxCoeff() == 1);
+}
+
+TEST_CASE("facet_components: truck", "[igl]")
+{
+  Eigen::MatrixXi F;
+  {
+    Eigen::MatrixXd V;
+    igl::read_triangle_mesh(test_common::data_path("truck.obj"), V, F);
+  }
+  Eigen::VectorXi C;
+  igl::facet_components(F,C);
+  REQUIRE(C.maxCoeff()+1 == 59);
+}