is vertex manifold

Former-commit-id: bf975ae25ac56a44503419c8fe8d397ab733d65d

include/igl/all_edges.cpp

@@ -7,9 +7,10 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "all_edges.h"
 
+template <typename DerivedF, typename DerivedE>
 IGL_INLINE void igl::all_edges(
-  const Eigen::MatrixXi & F,
-  Eigen::MatrixXi & E)
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedE> & E)
 {
   E.resize(F.rows()*F.cols(),F.cols()-1);
   switch(F.cols())

include/igl/all_edges.h

@@ -22,9 +22,10 @@ namespace igl
   // directed edge including repeats (meaning interior edges on a surface will
   // show up once for each direction and non-manifold edges may appear more than
   // once for each direction).
+  template <typename DerivedF, typename DerivedE>
   IGL_INLINE void all_edges(
-    const Eigen::MatrixXi & F,
-    Eigen::MatrixXi & E);
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    Eigen::PlainObjectBase<DerivedE> & E);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/boost/components.h

@@ -23,7 +23,8 @@ namespace igl
   IGL_INLINE void components(
     const Eigen::SparseMatrix<AScalar> & A,
     Eigen::PlainObjectBase<DerivedC> & C);
-  // Ditto but for mesh faces as input
+  // Ditto but for mesh faces as input. This computes connected components of
+  // **vertices** where **edges** establish connectivity.
   //
   // Inputs:
   //   F  n by 3 list of triangle indices

include/igl/is_edge_manifold.cpp

@@ -6,12 +6,56 @@
 // 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 "is_edge_manifold.h"
+#include "all_edges.h"
+#include "unique_simplices.h"
 
 #include <algorithm>
+#include <vector>
+
+template <
+  typename DerivedF, 
+  typename DerivedBF,
+  typename DerivedE,
+  typename DerivedEMAP,
+  typename DerivedBE>
+IGL_INLINE bool igl::is_edge_manifold(
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedBF>& BF,
+  Eigen::PlainObjectBase<DerivedE>& E,
+  Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+  Eigen::PlainObjectBase<DerivedBE>& BE)
+{
+  using namespace Eigen;
+  typedef typename DerivedF::Index Index;
+  typedef Matrix<typename DerivedF::Scalar,Dynamic,1> VectorXF;
+  typedef Matrix<typename DerivedF::Scalar,Dynamic,2> MatrixXF2;
+  MatrixXF2 allE;
+  all_edges(F,allE);
+  // Find unique undirected edges and mapping
+  VectorXF _;
+  unique_simplices(allE,E,_,EMAP);
+  std::vector<typename DerivedF::Index> count(E.rows(),0);
+  for(Index e = 0;e<EMAP.rows();e++)
+  {
+    count[EMAP[e]]++;
+  }
+  const Index m = F.rows();
+  BF.resize(m,3);
+  BE.resize(E.rows(),1);
+  bool all = true;
+  for(Index e = 0;e<EMAP.rows();e++)
+  {
+    const bool manifold = count[EMAP[e]] <= 2;
+    all &= BF(e%m,e/m) = manifold;
+    BE(EMAP[e]) = manifold;
+  }
+  return all;
+}
 
 template <typename DerivedV, typename DerivedF>
-IGL_INLINE bool igl::is_edge_manifold(const Eigen::PlainObjectBase<DerivedV>& /*V*/,
-                                 const Eigen::PlainObjectBase<DerivedF>& F)
+IGL_INLINE bool igl::is_edge_manifold(
+  const Eigen::PlainObjectBase<DerivedV>& /*V*/,
+  const Eigen::PlainObjectBase<DerivedF>& F)
 {
   // List of edges (i,j,f,c) where edge i<j is associated with corner i of face
   // f

include/igl/is_edge_manifold.h

@@ -10,7 +10,6 @@
 #include "igl_inline.h"
 
 #include <Eigen/Core>
-#include <vector>
 
 namespace igl 
 {
@@ -25,6 +24,18 @@ namespace igl
   //  Does not check for non-manifold vertices
   //
   // See also: is_vertex_manifold
+  template <
+    typename DerivedF, 
+    typename DerivedBF,
+    typename DerivedE,
+    typename DerivedEMAP,
+    typename DerivedBE>
+  IGL_INLINE bool is_edge_manifold(
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    Eigen::PlainObjectBase<DerivedBF>& BF,
+    Eigen::PlainObjectBase<DerivedE>& E,
+    Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+    Eigen::PlainObjectBase<DerivedBE>& BE);
   template <typename DerivedV, typename DerivedF>
   IGL_INLINE bool is_edge_manifold(
     const Eigen::PlainObjectBase<DerivedV>& V,

include/igl/is_vertex_manifold.cpp

@@ -0,0 +1,89 @@
+#include "is_vertex_manifold.h"
+#include "triangle_triangle_adjacency.h"
+#include "vertex_triangle_adjacency.h"
+#include "unique.h"
+#include <vector>
+#include <cassert>
+#include <map>
+#include <queue>
+#include <iostream>
+
+template <typename DerivedF,typename DerivedB>
+IGL_INLINE bool igl::is_vertex_manifold(
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedB>& B)
+{
+  using namespace std;
+  using namespace Eigen;
+  assert(F.cols() == 3 && "F must contain triangles");
+  typedef typename DerivedF::Scalar Index;
+  typedef typename DerivedF::Index FIndex;
+  const FIndex m = F.rows();
+  const Index n = F.maxCoeff()+1;
+  vector<vector<vector<FIndex > > > TT;
+  vector<vector<vector<FIndex > > > TTi;
+  triangle_triangle_adjacency(F,TT,TTi);
+
+  vector<vector<FIndex > > V2F,_1;
+  vertex_triangle_adjacency(n,F,V2F,_1);
+
+  const auto & check_vertex = [&](const Index v)->bool
+  {
+    vector<FIndex> uV2Fv;
+    {
+      vector<size_t> _1,_2;
+      unique(V2F[v],uV2Fv,_1,_2);
+    }
+    const FIndex one_ring_size = uV2Fv.size();
+    if(one_ring_size == 0)
+    {
+      return false;
+    }
+    const FIndex g = uV2Fv[0];
+    queue<Index> Q;
+    Q.push(g);
+    map<FIndex,bool> seen;
+    while(!Q.empty())
+    {
+      const FIndex f = Q.front();
+      Q.pop();
+      if(seen.count(f)==1)
+      {
+        continue;
+      }
+      seen[f] = true;
+      // Face f's neighbor lists opposite opposite each corner
+      for(const auto & c : TT[f])
+      {
+        // Each neighbor
+        for(const auto & n : c)
+        {
+          bool contains_v = false;
+          for(Index nc = 0;nc<F.cols();nc++)
+          {
+            if(F(n,nc) == v)
+            {
+              contains_v = true;
+              break;
+            }
+          }
+          if(seen.count(n)==0 && contains_v)
+          {
+            Q.push(n);
+          }
+        }
+      }
+    }
+    return one_ring_size == seen.size();
+  };
+
+  // Unreferenced vertices are considered non-manifold
+  B.setConstant(n,1,false);
+  // Loop over all vertices touched by F
+  bool all = true;
+  for(Index v = 0;v<n;v++)
+  {
+    all &= B(v) = check_vertex(v);
+  }
+  return all;
+}

include/igl/is_vertex_manifold.h

@@ -0,0 +1,39 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 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_IS_VERTEX_MANIFOLD_H
+#define IGL_IS_VERTEX_MANIFOLD_H
+#include "igl_inline.h"
+
+#include <Eigen/Core>
+
+namespace igl 
+{
+  // Check if a mesh is vertex-manifold. This only checks whether the faces
+  // incident on each vertex form exactly one connected component. Vertices
+  // incident on non-manifold edges are not consider non-manifold by this
+  // function (see is_edge_manifold.h). Unreferenced verties are considered
+  // non-manifold (zero components).
+  //
+  // Inputs:
+  //   F  #F by 3 list of triangle indices
+  // Outputs:
+  //   B  #V list whether 
+  // Returns whether mesh is vertex manifold.
+  //
+  // See also: is_edge_manifold
+  template <typename DerivedF,typename DerivedB>
+  IGL_INLINE bool is_vertex_manifold(
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    Eigen::PlainObjectBase<DerivedB>& B);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "is_vertex_manifold.cpp"
+#endif
+
+#endif

include/igl/per_vertex_normals.cpp

@@ -32,13 +32,13 @@ IGL_INLINE void igl::per_vertex_normals(
   return per_vertex_normals(V,F,PER_VERTEX_NORMALS_WEIGHTING_TYPE_DEFAULT,N);
 }
 
-template <typename DerivedV, typename DerivedF>
+template <typename DerivedV, typename DerivedF, typename DerivedFN, typename DerivedN>
 IGL_INLINE void igl::per_vertex_normals(
   const Eigen::PlainObjectBase<DerivedV>& V,
   const Eigen::PlainObjectBase<DerivedF>& F,
   const igl::PerVertexNormalsWeightingType weighting,
-  const Eigen::PlainObjectBase<DerivedV>& FN,
-  Eigen::PlainObjectBase<DerivedV> & N)
+  const Eigen::PlainObjectBase<DerivedFN>& FN,
+  Eigen::PlainObjectBase<DerivedN> & N)
 {
   // Resize for output
   N.setZero(V.rows(),3);
@@ -99,9 +99,4 @@ IGL_INLINE void igl::per_vertex_normals(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 // generated by autoexplicit.sh
-template void igl::per_vertex_normals<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::PerVertexNormalsWeightingType, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
-template void igl::per_vertex_normals<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-template void igl::per_vertex_normals<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&);
-template void igl::per_vertex_normals<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/per_vertex_normals.h

@@ -46,13 +46,13 @@ namespace igl
     Eigen::PlainObjectBase<DerivedV> & N);
   // Inputs:
   //   FN  #F by 3 matrix of face (triangle) normals
-  template <typename DerivedV, typename DerivedF>
+  template <typename DerivedV, typename DerivedF, typename DerivedFN, typename DerivedN>
   IGL_INLINE void per_vertex_normals(
     const Eigen::PlainObjectBase<DerivedV>& V,
     const Eigen::PlainObjectBase<DerivedF>& F,
     const PerVertexNormalsWeightingType weighting,
-    const Eigen::PlainObjectBase<DerivedV>& FN,
-    Eigen::PlainObjectBase<DerivedV> & N);
+    const Eigen::PlainObjectBase<DerivedFN>& FN,
+    Eigen::PlainObjectBase<DerivedN> & N);
   // Without weighting
   template <typename DerivedV, typename DerivedF>
   IGL_INLINE void per_vertex_normals(

include/igl/sort.cpp

@@ -16,12 +16,12 @@
 #include <algorithm>
 #include <iostream>
 
-template <typename DerivedX, typename DerivedIX>
+template <typename DerivedX, typename DerivedY, typename DerivedIX>
 IGL_INLINE void igl::sort(
   const Eigen::PlainObjectBase<DerivedX>& X,
   const int dim,
   const bool ascending,
-  Eigen::PlainObjectBase<DerivedX>& Y,
+  Eigen::PlainObjectBase<DerivedY>& Y,
   Eigen::PlainObjectBase<DerivedIX>& IX)
 {
   // get number of rows (or columns)
@@ -74,12 +74,12 @@ IGL_INLINE void igl::sort(
   }
 }
 
-template <typename DerivedX, typename DerivedIX>
+template <typename DerivedX, typename DerivedY, typename DerivedIX>
 IGL_INLINE void igl::sort_new(
   const Eigen::PlainObjectBase<DerivedX>& X,
   const int dim,
   const bool ascending,
-  Eigen::PlainObjectBase<DerivedX>& Y,
+  Eigen::PlainObjectBase<DerivedY>& Y,
   Eigen::PlainObjectBase<DerivedIX>& IX)
 {
   // get number of rows (or columns)
@@ -137,12 +137,12 @@ IGL_INLINE void igl::sort_new(
   }
 }
 
-template <typename DerivedX, typename DerivedIX>
+template <typename DerivedX, typename DerivedY, typename DerivedIX>
 IGL_INLINE void igl::sort2(
   const Eigen::PlainObjectBase<DerivedX>& X,
   const int dim,
   const bool ascending,
-  Eigen::PlainObjectBase<DerivedX>& Y,
+  Eigen::PlainObjectBase<DerivedY>& Y,
   Eigen::PlainObjectBase<DerivedIX>& IX)
 {
   using namespace Eigen;
@@ -213,24 +213,4 @@ IGL_INLINE void igl::sort(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 // generated by autoexplicit.sh
-template void igl::sort<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
-template void igl::sort<Eigen::Matrix<float, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
-template void igl::sort<Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
-template void igl::sort<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-template void igl::sort<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-template void igl::sort<igl::SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> > >(std::vector<igl::SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> >, std::allocator<igl::SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> > > > const&, bool, std::vector<igl::SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> >, std::allocator<igl::SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> > > >&, std::vector<unsigned long, std::allocator<unsigned long> >&);
-template void igl::sort<int>(std::vector<int, std::allocator<int> > const&, bool, std::vector<int, std::allocator<int> >&, std::vector<unsigned long, std::allocator<unsigned long> >&);
-template void igl::sort<igl::SortableRow<Eigen::Matrix<double, -1, 1, 0, -1, 1> > >(std::vector<igl::SortableRow<Eigen::Matrix<double, -1, 1, 0, -1, 1> >, std::allocator<igl::SortableRow<Eigen::Matrix<double, -1, 1, 0, -1, 1> > > > const&, bool, std::vector<igl::SortableRow<Eigen::Matrix<double, -1, 1, 0, -1, 1> >, std::allocator<igl::SortableRow<Eigen::Matrix<double, -1, 1, 0, -1, 1> > > >&, std::vector<unsigned long, std::allocator<unsigned long> >&);
-template void igl::sort2<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-template void igl::sort2<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-template void igl::sort_new<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-template void igl::sort<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template void igl::sort<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template void igl::sort<Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template void igl::sort<Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 2, 0, -1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&);
-template void igl::sort<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-template void igl::sort<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/sort.h

@@ -30,28 +30,28 @@ namespace igl
   //   Y  m by n matrix whose entries are sorted
   //   IX  m by n matrix of indices so that if dim = 1, then in matlab notation
   //     for j = 1:n, Y(:,j) = X(I(:,j),j); end
-  template <typename DerivedX, typename DerivedIX>
+  template <typename DerivedX, typename DerivedY, typename DerivedIX>
   IGL_INLINE void sort(
     const Eigen::PlainObjectBase<DerivedX>& X,
     const int dim,
     const bool ascending,
-    Eigen::PlainObjectBase<DerivedX>& Y,
+    Eigen::PlainObjectBase<DerivedY>& Y,
     Eigen::PlainObjectBase<DerivedIX>& IX);
-  template <typename DerivedX, typename DerivedIX>
+  template <typename DerivedX, typename DerivedY, typename DerivedIX>
   // Only better if size(X,dim) is small
   IGL_INLINE void sort_new(
     const Eigen::PlainObjectBase<DerivedX>& X,
     const int dim,
     const bool ascending,
-    Eigen::PlainObjectBase<DerivedX>& Y,
+    Eigen::PlainObjectBase<DerivedY>& Y,
     Eigen::PlainObjectBase<DerivedIX>& IX);
   // Special case if size(X,dim) == 2
-  template <typename DerivedX, typename DerivedIX>
+  template <typename DerivedX, typename DerivedY, typename DerivedIX>
   IGL_INLINE void sort2(
     const Eigen::PlainObjectBase<DerivedX>& X,
     const int dim,
     const bool ascending,
-    Eigen::PlainObjectBase<DerivedX>& Y,
+    Eigen::PlainObjectBase<DerivedY>& Y,
     Eigen::PlainObjectBase<DerivedIX>& IX);
 
   // Act like matlab's [Y,I] = SORT(X) for std library vectors

include/igl/triangle_triangle_adjacency.cpp

@@ -6,8 +6,9 @@
 // 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 "triangle_triangle_adjacency.h"
-
-#include <igl/is_edge_manifold.h>
+#include "is_edge_manifold.h"
+#include "all_edges.h"
+#include <map>
 #include <algorithm>
 
 template <typename Scalar, typename Index>
@@ -98,6 +99,63 @@ IGL_INLINE void igl::triangle_triangle_adjacency(const Eigen::PlainObjectBase<Sc
   triangle_triangle_adjacency_extractTTi(F,TTT,TTi);
 }
 
+template <
+  typename DerivedF, 
+  typename TTIndex, 
+  typename TTiIndex>
+  IGL_INLINE void igl::triangle_triangle_adjacency(
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    std::vector<std::vector<std::vector<TTIndex> > > & TT,
+    std::vector<std::vector<std::vector<TTiIndex> > > & TTi)
+{
+  using namespace Eigen;
+  using namespace std;
+  using namespace igl;
+  assert(F.cols() == 3 && "Faces must be triangles");
+  typedef typename DerivedF::Index Index;
+  // number of faces
+  const int m = F.rows();
+  // All occurances of directed edges
+  MatrixXi E;
+  all_edges(F,E);
+  assert(E.rows() == 3*m);
+  // uE2E[i] --> {j,k,...} means unique edge i corresponds to face edges j and
+  // k (where j-edge comes is the j/m edge of face j%m)
+  map<pair<Index,Index>,vector<Index> > uE2E;
+  for(int e = 0;e<E.rows();e++)
+  {
+    Index i = E(e,0);
+    Index j = E(e,1);
+    if(i<j)
+    {
+      uE2E[pair<Index,Index>(i,j)].push_back(e);
+    }else
+    {
+      uE2E[pair<Index,Index>(j,i)].push_back(e);
+    }
+  }
+  // E2E[i] --> {j,k,...} means face edge i corresponds to other faces edges j
+  // and k
+  TT.resize (m,vector<vector<Index> >(F.cols()));
+  TTi.resize(m,vector<vector<Index> >(F.cols()));
+  for(int e = 0;e<E.rows();e++)
+  {
+    const Index i = E(e,0);
+    const Index j = E(e,1);
+    const Index f = e%m;
+    const Index c = e/m;
+    const vector<Index> & N = 
+      i<j ? uE2E[pair<Index,Index>(i,j)] : uE2E[pair<Index,Index>(j,i)];
+    for(const auto & ne : N)
+    {
+      const Index nf = ne%m;
+      const Index nc = ne/m;
+      TT[f][c].push_back(nf);
+      TTi[f][c].push_back(nc);
+    }
+  }
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 template void igl::triangle_triangle_adjacency<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);

include/igl/triangle_triangle_adjacency.h

@@ -29,6 +29,7 @@ namespace igl
   //   TTi  #F by #3 adjacent matrix, the element i,j is the id of edge of the triangle TT(i,j) that is adjacent with triangle i
   // NOTE: the first edge of a triangle is [0,1] the second [1,2] and the third [2,3].
   //       this convention is DIFFERENT from cotmatrix_entries.h
+  // Known bug: this should not need to take V as input.
 
   template <typename Scalar, typename Index>
   IGL_INLINE void triangle_triangle_adjacency(const Eigen::PlainObjectBase<Scalar>& V,
@@ -57,6 +58,24 @@ namespace igl
   IGL_INLINE void triangle_triangle_adjacency_extractTTi(const Eigen::PlainObjectBase<Index>& F,
                                 std::vector<std::vector<int> >& TTT,
                                 Eigen::PlainObjectBase<Index>& TTi);
+  // Adjacency list version, which works with non-manifold meshes
+  //
+  // Inputs:
+  //   F  #F by 3 list of triangle indices
+  // Outputs:
+  //   TT  #F by 3 list of lists so that TT[i][c] --> {j,k,...} means that faces j and
+  //     k etc. are edge-neighbors of face i on face i's edge opposite corner c
+  //   TTj  #F list of lists so that TTj[i][c] --> {j,k,...} means that face
+  //     TT[i][c][0] is an edge-neighbor of face i incident on the edge of face
+  //     TT[i][c][0] opposite corner j, and TT[i][c][1] " corner k, etc.
+  template <
+    typename DerivedF, 
+    typename TTIndex, 
+    typename TTiIndex>
+    IGL_INLINE void triangle_triangle_adjacency(
+      const Eigen::PlainObjectBase<DerivedF> & F,
+      std::vector<std::vector<std::vector<TTIndex> > > & TT,
+      std::vector<std::vector<std::vector<TTiIndex> > > & TTi);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/vertex_triangle_adjacency.cpp

@@ -7,24 +7,23 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "vertex_triangle_adjacency.h"
 
-#include "verbose.h"
-
-template <typename DerivedV, typename DerivedF, typename IndexType>
+template <typename DerivedF, typename VFType, typename VFiType>
 IGL_INLINE void igl::vertex_triangle_adjacency(
-                   const Eigen::PlainObjectBase<DerivedV>& V,
-                   const Eigen::PlainObjectBase<DerivedF>& F,
-                   std::vector<std::vector<IndexType> >& VF,
-                   std::vector<std::vector<IndexType> >& VFi)
+  const typename DerivedF::Scalar n,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  std::vector<std::vector<VFType> >& VF,
+  std::vector<std::vector<VFiType> >& VFi)
 {
   VF.clear();
   VFi.clear();
 
-  VF.resize(V.rows());
-  VFi.resize(V.rows());
+  VF.resize(n);
+  VFi.resize(n);
 
-  for(int fi=0; fi<F.rows(); ++fi)
+  typedef typename DerivedF::Index Index;
+  for(Index fi=0; fi<F.rows(); ++fi)
   {
-    for(int i = 0; i < F.cols(); ++i)
+    for(Index i = 0; i < F.cols(); ++i)
     {
       VF[F(fi,i)].push_back(fi);
       VFi[F(fi,i)].push_back(i);
@@ -32,6 +31,17 @@ IGL_INLINE void igl::vertex_triangle_adjacency(
   }
 }
 
+
+template <typename DerivedV, typename DerivedF, typename IndexType>
+IGL_INLINE void igl::vertex_triangle_adjacency(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  std::vector<std::vector<IndexType> >& VF,
+  std::vector<std::vector<IndexType> >& VFi)
+{
+  return vertex_triangle_adjacency(V.rows(),F,VF,VFi);
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 // generated by autoexplicit.sh

include/igl/vertex_triangle_adjacency.h

@@ -17,7 +17,8 @@ namespace igl
   // vertex_face_adjacency constructs the vertex-face topology of a given mesh (V,F)
   //
   // Inputs:
-  //   V  #V by 3 list of vertex coordinates
+  //   //V  #V by 3 list of vertex coordinates
+  //   n  number of vertices #V (e.g. `F.maxCoeff()+1` or `V.rows()`)
   //   F  #F by dim list of mesh faces (must be triangles)
   // Outputs:
   //   VF  #V list of lists of incident faces (adjacency list)
@@ -27,13 +28,21 @@ namespace igl
   // See also: edges, cotmatrix, diag, vv
   //
   // Known bugs: this should not take V as an input parameter.
-
+  // Known bugs/features: if a facet is combinatorially degenerate then faces
+  // will appear multiple times in VF and correpondingly in VFI (j appears
+  // twice in F.row(i) then i will appear twice in VF[j])
+  template <typename DerivedF, typename VFType, typename VFiType>
+  IGL_INLINE void vertex_triangle_adjacency(
+    const typename DerivedF::Scalar n,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    std::vector<std::vector<VFType> >& VF,
+    std::vector<std::vector<VFiType> >& VFi);
   template <typename DerivedV, typename DerivedF, typename IndexType>
   IGL_INLINE void vertex_triangle_adjacency(
-                     const Eigen::PlainObjectBase<DerivedV>& V,
-                     const Eigen::PlainObjectBase<DerivedF>& F,
-                     std::vector<std::vector<IndexType> >& VF,
-                     std::vector<std::vector<IndexType> >& VFi);
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    std::vector<std::vector<IndexType> >& VF,
+    std::vector<std::vector<IndexType> >& VFi);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/writeDMAT.cpp

@@ -94,6 +94,34 @@ IGL_INLINE bool igl::writeDMAT(
   return true;
 }
 
+template <typename Scalar>
+IGL_INLINE bool igl::writeDMAT(
+  const std::string file_name, 
+  const std::vector<Scalar > W)
+{
+  FILE * fp = fopen(file_name.c_str(),"w");
+  if(fp == NULL)
+  {
+    fprintf(stderr,"IOError: writeDMAT() could not open %s...",file_name.c_str());
+    return false; 
+  }
+  int num_rows = (int)W.size();
+  int num_cols = 0;
+  if(num_rows > 0)
+  {
+    num_cols = 1;
+  }
+  // first line contains number of columns and number of rows
+  fprintf(fp,"%d %d\n",num_cols,num_rows);
+  // loop over rows (down columns) quickly
+  for(int i = 0;i < num_rows;i++)
+  {
+    fprintf(fp,"%0.15lf\n",(double)W[i]);
+  }
+  fclose(fp);
+  return true;
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 // generated by autoexplicit.sh

include/igl/writeDMAT.h

@@ -32,6 +32,10 @@ namespace igl
   IGL_INLINE bool writeDMAT(
     const std::string file_name, 
     const std::vector<std::vector<Scalar> > W);
+  template <typename Scalar>
+  IGL_INLINE bool writeDMAT(
+    const std::string file_name, 
+    const std::vector<Scalar > W);
 }
 
 #ifndef IGL_STATIC_LIBRARY