Fix split_nonmanifold (#2344)

* failing test

* simply cut along all non-manifold edges

* fix compile

* before rewrite

* after rewrite

---------

Co-authored-by: Alec Jacobson <[email protected]>

include/igl/split_nonmanifold.cpp

@@ -6,10 +6,21 @@
 // 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 "split_nonmanifold.h"
+#include "unique_edge_map.h"
 #include "connected_components.h"
-#include "remove_unreferenced.h"
-#include "find.h"
-#include "ismember_rows.h"
+#include "unique.h"
+#include "sort.h"
+#include "triangle_triangle_adjacency.h"
+#include "is_edge_manifold.h"
+#include <unordered_map>
+#include <cassert>
+#include <type_traits>
+
+#include "is_vertex_manifold.h"
+#include "matlab_format.h"
+#include <iostream>
+#include <unordered_set>
+#include <utility>
 
 template <
   typename DerivedF,
@@ -21,121 +32,406 @@ IGL_INLINE void igl::split_nonmanifold(
   Eigen::PlainObjectBase <DerivedSF> & SF,
   Eigen::PlainObjectBase <DerivedSVI> & SVI)
 {
-  // Number of faces
+  using Scalar = typename DerivedSF::Scalar;
+  // Scalar must allow negative values
+  static_assert(std::is_signed<Scalar>::value,"Scalar must be signed");
+  using MatrixX2I = Eigen::Matrix<Scalar,Eigen::Dynamic,2>;
+  using MatrixX3I = Eigen::Matrix<Scalar,Eigen::Dynamic,3>;
+  using VectorXI = Eigen::Matrix< Scalar,Eigen::Dynamic,1>;
+  MatrixX2I E,uE;
+  VectorXI EMAP,uEC,uEE;
+  igl::unique_edge_map(F,E,uE,EMAP,uEC,uEE);
+
+  // Let's assume the most convenient connectivity data structure and worry
+  // about performance later
+
+  // V[c] = v means that corner c is mapped to vertex v
+  // Start with all corners mapped to singleton vertices[:w
+  Eigen::VectorXi V = Eigen::VectorXi::LinSpaced(F.size(),0,F.size()-1);
+  // Convenience map so that CF(f,i) = V[c] = v where c is the ith corner of
+  // face f.
+  Eigen::Map<Eigen::MatrixXi> CF = Eigen::Map<Eigen::MatrixXi>(V.data(),F.rows(),F.cols());
+ 
+  // C[v][j] = c means that c is the jth corner in the group of corners at i
+  std::vector<std::vector<int> > C(F.size());
+  for(int i = 0;i<F.size();i++) { C[i] = {i}; }
+
   const int m = F.rows();
-  // For moment aact like everything will be split
-  SF.resize(m,3);
+
+  // O(S) where S = |star(v)|
+  const auto star = [&](const int v)->std::vector<int>
+  {
+    std::vector<int> faces(C[v].size());
+    for(int i = 0;i<C[v].size();i++)
+    {
+      faces[i] = C[v][i]%m;
+    }
+    return faces;
+  };
+
+  // O(S) where S = |star(v)|
+  const auto nonmanifold_edge_star = [&](const int v)->std::vector<int>
   {
-    int k =0;
-    for(int j = 0;j<3;j++)
+    std::vector<int> edges;
+    for(int e : C[v])
     {
-      for(int i = 0;i<m;i++)
+      const int f = e%m;
+      for(int j = 1;j<3;j++)
       {
-        SF(i,j) = k++;
+        // next edge
+        const int e1 = (e+j*m)%(3*m);
+        const int u1 = EMAP(e1);
+        if(uEC(u1+1)-uEC(u1) > 2)
+        {
+          edges.push_back(e1);
+        }
       }
     }
-  }
-  // Edges in SF
-  Eigen::MatrixXi E(m*3,2);
-  for(int i = 0;i<m;i++)
+    return edges;
+  };
+
+
+  // O(S) where S = |star(v)|
+  const std::function<void(
+    Eigen::VectorXi &,
+    std::vector<std::vector<int> > &,
+      const int, const int)> merge_vertex = 
+    [&merge_vertex](Eigen::VectorXi & V,
+      std::vector<std::vector<int> > & C,
+        const int u, const int v)
   {
-    E.row(i+0*m) << SF(i,1),SF(i,2);
-    E.row(i+1*m) << SF(i,2),SF(i,0);
-    E.row(i+2*m) << SF(i,0),SF(i,1);
-  }
-  // Reindex E by F
-  Eigen::MatrixXi FE(E.rows(),E.cols());
-  for(int i = 0;i<E.rows();i++)
+    if(u == v) { return; }
+    if(u > v) { merge_vertex(V,C,v,u); return; }
+    assert(u < v);
+    // Consider each corner in v
+    for(const int c : C[v])
+    {
+      V[c] = u;
+    }
+    // Merge C[v] into C[u]
+    C[u].insert(C[u].end(),C[v].begin(),C[v].end());
+    C[v].clear();
+  };
+
+  // O(S) where S is the size of the star of e's first vertex.
+  // This could probably be O(1) with careful bookkeeping
+  const auto is_boundary = [&](const int e)->bool
   {
-    for(int j = 0;j<2;j++)
+    // e----d
+    //  \   |
+    //   \f₁↑
+    //    \ |
+    //      s
+    const int s = (e+1*m)%(3*m);
+    const int d = (e+2*m)%(3*m);
+    const int f = e%m;
+    const int vs = V[s];
+    const int vd = V[d];
+    // Consider every face in the star of s
+    for(const int g : star(vs))
     {
-      const int fi = E(i,j) % m;
-      const int fj = E(i,j) / m;
-      FE(i,j) = F(fi,fj);
+      if(g == f) { continue; }
+      // Consider each edge in g
+      for(int i = 0;i<3;i++)
+      {
+        const int a = (g+(i+1)*m)%(3*m);
+        const int b = (g+(i+2)*m)%(3*m);
+        // Is that edge the same as e?
+        if(V[a] == vd && V[b] == vs) { return false; }
+        if(V[a] == vs && V[b] == vd) { return false; }
+      }
     }
-  }
-  // Flip orientation
-  Eigen::MatrixXi FE_flip = FE.rowwise().reverse();
-  // Find which exist in both directions
-
-
-  Eigen::Array<bool,Eigen::Dynamic,1> I;
-  Eigen::VectorXi J;
-  igl::ismember_rows(FE,FE_flip,I,J);
-  // Just keep those find
-  const auto II = igl::find(I);
-  Eigen::MatrixXi EI = E(II,Eigen::all);
-  Eigen::VectorXi JI = J(II);
-  Eigen::MatrixXi EJI = E(JI,Eigen::all);
-
-  Eigen::MatrixXi EJI_flip = EJI.rowwise().reverse();
-  // Build adjacency matrix
-  std::vector<Eigen::Triplet<bool> > Aijv; 
-  Aijv.reserve(EI.size());
-  for(int i = 0;i<EI.rows();i++)
+    return true;
+  };
+ 
+  // Ω(m) and probably  O(m log m) or worse.
+  // This should take in the candidate merge edge pair, extract the submesh and
+  // just check if that's manifold. Then it would be O(S) where S is the size of
+  // biggest star of the edges' vertices.
+  //
+  // My guess is that is_edge_manifold is O(m) but is_vertex_manifold is
+  // O(max(F))
+  const auto is_manifold = [](Eigen::MatrixXi F)->bool
   {
-    for(int j = 0;j<2;j++)
+    Eigen::Array<bool,Eigen::Dynamic,1> referenced = 
+      Eigen::Array<bool,Eigen::Dynamic,1>::Zero(F.maxCoeff()+1,1);
+    for(int i = 0;i<F.size();i++)
     {
-      Aijv.emplace_back( 
-        EI(i,j), 
-        EJI_flip(i,j), 
-        true);
+      referenced(F(i)) = true;
+    }
+    Eigen::Array<bool,Eigen::Dynamic,1> VM;
+    igl::is_vertex_manifold(F,VM);
+    for(int i = 0;i<VM.size();i++)
+    {
+      if(referenced(i) && !VM(i))
+      {
+        return false;
+      }
+    }
+    return igl::is_edge_manifold(F);
+  };
+
+
+  // Ω(S) where S is the largest star of (vs1,vd2) or (vd1,vs2)
+  // I think that is_vertex/edge_manifold(L) is O(|L| log |L|) so I think that
+  // should make this O(|S| log |S|) with some gross constants because of all
+  // the copying and sorting things into different data structures.
+  //
+  // merging edges (vs1,vd2) and (vd1,vs2) requires merging vertices (vs1→vd1) and
+  // (vd2→vd2).
+  //
+  // Merging vertices (a→b) will change and only change the stars of a and b.
+  // That is, some vertex c ≠ a,b will have the sam star before and after.
+  //
+  // Whether a vertex is singular depends entirely on its star.
+  //
+  // Therefore, the only vertices we need to check for non-manifoldness are
+  // vs=(vs1,vd2) and vd=(vd1,vs2).
+  const auto simulated_merge_is_manifold = 
+    [&](
+        const int vs1, const int vd2,
+        const int vd1, const int vs2)->bool
+  {
+    // all_faces[i] = f means that f is the ith face in the list of stars.
+    std::vector<int> all_faces;
+    for(int v : {vs1,vd2,vd1,vs2})
+    {
+      std::vector<int> star_v = star(v);
+      all_faces.insert(all_faces.end(),star_v.begin(),star_v.end());
+    }
+    // unique_faces[l] = f means that f is the lth unique face in the list of
+    // stars.
+    std::vector<int> unique_faces;
+    std::vector<size_t> _, local;
+    igl::unique(all_faces,unique_faces,_,local);
+    Eigen::MatrixXi L(unique_faces.size(),3);
+    // collect local faces
+    for(int l = 0;l<unique_faces.size();l++)
+    {
+      L.row(l) = CF.row(unique_faces[l]);
+    }
+    {
+      int f = 0;
+      const auto merge_local = [&](const int v1, const int v2)
+      {
+        const int u = std::min(v1,v2);
+        for(const int v : {v1,v2})
+        {
+          for(const int c : C[v])
+          {
+            const int i = c/m;
+            L(local[f++],i) = u;
+          }
+        }
+      };
+      // must match order {vs1,vd2,vd1,vs2} above
+      merge_local(vs1,vd2);
+      merge_local(vd1,vs2);
+    }
+    
+    // remove unreferenced vertices by mapping each index in L to a unique
+    // index between 0 and size(unique(L))
+    std::unordered_map<int,int> M;
+    for(int & i : L.reshaped())
+    {
+      if(M.find(i) == M.end())
+      {
+        M[i] = M.size();
+      }
+      i = M[i];
+    }
+    // Only need to check if the two vertices being merged are manifold
+    Eigen::Array<bool,Eigen::Dynamic,1> VM;
+    const int vs = std::min(vs1,vd2);
+    const int vd = std::min(vd1,vs2);
+    igl::is_vertex_manifold(L,VM);
+    if(!VM(M[vs])) { 
+      return false; 
+    }
+    if(!VM(M[vd])) { 
+      return false; 
+    }
+    // Probably only need to check incident edges in star, but this also
+    // checks link 
+    return igl::is_edge_manifold(L);
+  };
+
+  const auto merge_edge = [&](const int e1, const int e2)
+  {
+    // Ideally we would track whether an edge is a boundary so we can just
+    // assert these. But because of "implied stitches" it's not necessarily just
+    // e1 and e2 which become non-boundary when e1 and e2 are merged.
+    //assert(is_boundary(e1));
+    //assert(is_boundary(e2));
+    if(!is_boundary(e1) || !is_boundary(e2)) { return false; }
+    assert(e1 != e2);
+
+    const bool consistent = E(e1,0) == E(e2,1);
+    // skip if inconsistently oriented
+    if(!consistent) { return false; }
+    // The code below is assuming merging consistently oriented edges
+    assert(E(e1,1) == E(e2,0));
+
+    //
+    // e1--d1  s2--e2
+    //  \   |  |   /
+    //   \f₁↑  ↓f₂/
+    //    \ |  | /
+    //     s1  d2
+    //
+    //
+
+
+    // "Cutting and Stitching: Converting Sets of Polygons to Manifold
+    // Surfaces" [Guéziec et al. 2001]
+    const int s1 = (e1+1*m)%(3*m);
+    const int d1 = (e1+2*m)%(3*m);
+#ifndef NDEBUG
+    {
+      const int f1 = e1 % m;
+      const int i1 = e1 / m;
+      const int s1_test = f1 + ((i1+1)%3)*m;
+      const int d1_test = f1 + ((i1+2)%3)*m;
+      assert(s1 == s1_test);
+      assert(d1 == d1_test);
+    }
+#endif
+    int s2 = (e2+1*m)%(3*m);
+    int d2 = (e2+2*m)%(3*m);
+    const int vs1 = V[s1];
+    const int vd2 = V[d2];
+    const int vd1 = V[d1];
+    const int vs2 = V[s2];
+
+#ifndef IGL_SPLIT_NONMANIFOLD_DEBUG
+    const auto simulated_merge_is_manifold_old = [&]()->bool
+    {
+      Eigen::VectorXi V_copy = V;
+      std::vector<std::vector<int> > C_copy = C;
+      merge_vertex(V_copy,C_copy,vs1,vd2);
+      merge_vertex(V_copy,C_copy,vd1,vs2);
+      Eigen::Map<Eigen::MatrixXi> CF_copy = 
+        Eigen::Map<Eigen::MatrixXi>(V_copy.data(),CF.rows(),CF.cols());
+      if(!is_manifold(CF_copy)) { return false; }
+      return true;
+    };
+    const bool ret_old = simulated_merge_is_manifold_old();
+    const bool ret = simulated_merge_is_manifold(vs1,vd2,vd1,vs2);
+    if(ret != ret_old)
+    {
+      assert(false);
+    }
+#endif
+    // I claim this is completely unnecessary if the unique edge was originally
+    // manifold.
+    //
+    // I also hypothesize that this is unnecessary when conducting depth-first
+    // traversals starting at a successful merge.
+    //
+    // That is, we never need to call this in the current algorithm.
+    const int u = EMAP(e1);
+    assert(EMAP(e1) == EMAP(e2));
+    const int edge_valence = uEC(u+1)-uEC(u);
+    assert(edge_valence >= 2);
+    if(edge_valence>2 && !simulated_merge_is_manifold(vs1,vd2,vd1,vs2))
+    {
+      return false;
     }
-  }
-  // Build A to contain off-diagonals only if both directions are present
-  Eigen::SparseMatrix<bool> A1(m*3,m*3);
-  A1.setFromTriplets(Aijv.begin(),Aijv.end());
-  // For some reason I can't write `A = A1 && A1.transpose();`
-  Eigen::SparseMatrix<bool> A1T = A1.transpose();
-  Eigen::SparseMatrix<bool> A = A1 && A1T;
 
+    // Now we can merge
+    merge_vertex(V,C,vs1,vd2);
+    merge_vertex(V,C,vd1,vs2);
+    return true;
+  };
 
-  Eigen::VectorXi K;
+  // Consider each unique edge in the original mesh
+
+  // number of faces incident on each unique edge
+  Eigen::VectorXi D = uEC.tail(uEC.rows()-1)-uEC.head(uEC.rows()-1);
+  Eigen::VectorXi uI;
   {
-    Eigen::VectorXi _;
-    igl::connected_components(A,K,_);
+    Eigen::VectorXi sD;
+    igl::sort(D,1,true,sD,uI);
   }
 
 
-  // Remap by components
-  for(int j = 0;j<3;j++)
+  const std::function<void(const int)> dfs = [&](const int e)
   {
-    for(int i = 0;i<m;i++)
+    // we just successfully merged e, find all other non-manifold edges sharing
+    // a current vertex with e and try to merge it too.
+    const int s = (e+1*m)%(3*m);
+    const int d = (e+2*m)%(3*m);
+    for(const int c : {s,d})
     {
-      SF(i,j) = K(SF(i,j));
+      const int v = V[c];
+      std::vector<int> nme = nonmanifold_edge_star(v);
+      // My thinking is that this must be size 0 or 2.
+      for(int i = 0;i<nme.size();i++)
+      {
+        const int e1 = nme[i];
+        for(int j = i+1;j<nme.size();j++)
+        {
+          const int e2 = nme[j];
+          if(merge_edge(e1,e2))
+          {
+            dfs(e2);
+          }
+        }
+      }
     }
-  }
-  
-  // Initial mapping
-  Eigen::VectorXi SVI0(m*3);
+  };
+
+  // Every edge starts as a boundary
+  for(auto u : uI)
   {
-    int k =0;
-    for(int j = 0;j<3;j++)
+    // if boundary skip
+    if(uEC(u+1)-uEC(u) == 1) { continue; }
+    for(int j = uEC(u);j<uEC(u+1);j++)
     {
-      for(int i = 0;i<m;i++)
+      const int e1 = uEE(j);
+      for(int k = j+1;k<uEC(u+1);k++)
       {
-        SVI0(k++) = F(i,j);
+        const int e2 = uEE(k);
+        if(merge_edge(e1,e2))
+        { 
+          // for non-manifold edges, launch search from e1 and e2
+          if(uEC(u+1)-uEC(u) > 2)
+          {
+            dfs(e1);
+          }
+          break; 
+        }
       }
     }
   }
-  assert(K.size() == m*3);
-  // Scatter via K
-  // SVI1(K) = SVI(K);
-  Eigen::VectorXi SVI1(m*3);
-  SVI1(K) = SVI0;
 
+
+  
+  // Ideally we'd do this so that all duplicated vertices end up at the end
+  // rather than scrambling the whole mesh.
   {
-    Eigen::VectorXi _,J;
-    igl::remove_unreferenced(SF.maxCoeff()+1,SF,_,J);
-    // Remap by J
-    for(int j = 0;j<3;j++)
+    SVI.resize(F.size());
+    std::vector<bool> marked(F.size());
+    VectorXI J = VectorXI::Constant(F.size(),-1);
+    SF.resize(F.rows(),F.cols());
     {
-      for(int i = 0;i<m;i++)
+      int nv = 0;
+      for(int f = 0;f<m;f++)
       {
-        SF(i,j) = J(SF(i,j));
+        for(int i = 0;i<3;i++)
+        {
+          const int c = CF(f,i);
+          if(J(c) == -1)
+          {
+            J(c) = nv;
+            SVI(nv) = F(f,i);
+            nv++;
+          }
+          SF(f,i) = J(c);
+        }
       }
+      SVI.conservativeResize(nv);
     }
-    SVI = SVI1(J);
   }
 
 }

include/igl/unique_edge_map.h

@@ -75,7 +75,7 @@ namespace igl
   /// 
   /// \code{cpp}
   /// // Using uEC,uEE
-  /// for(int u = 0;u<uE.size();u++)
+  /// for(int u = 0;u<uE.rows();u++)
   /// {
   ///   for(int j = uEC(u);j<uEC(u+1);j++)
   ///   {

tests/include/igl/split_nonmanifold.cpp

@@ -1,11 +1,45 @@
 #include <test_common.h>
 #include <igl/split_nonmanifold.h>
+#include <igl/facet_components.h>
+#include <igl/is_edge_manifold.h>
+#include <igl/is_vertex_manifold.h>
+#include <igl/triangulated_grid.h>
+#include <igl/remove_duplicate_vertices.h>
+
 #include <igl/matlab_format.h>
+#include <igl/writePLY.h>
 #include <iostream>
 
+void check_same_but_manifold(
+    const Eigen::MatrixXi & F,
+    const Eigen::MatrixXi & SF,
+    const Eigen::VectorXi & I)
+{
+  // Check that new mesh has no non-manifold edges
+  REQUIRE(igl::is_edge_manifold(SF));
+  // Check that new mesh has no non-manifold vertices
+  REQUIRE(igl::is_vertex_manifold(SF));
+  // Check the new mesh references original vertex information
+  Eigen::MatrixXi ISF(SF.rows(),SF.cols());
+  for(int i = 0;i<ISF.rows();i++)
+  {
+    for(int j = 0;j<ISF.cols();j++)
+    {
+      ISF(i,j) = I(SF(i,j));
+    }
+  }
+  test_common::assert_eq(F,ISF);
+}
+
 TEST_CASE("split_nonmanifold: edge-fan", "[igl]")
 {
   using namespace igl;
+  Eigen::MatrixXi F(5,3);
+  F<<0,1,3,
+     0,3,2,
+     0,4,3,
+     0,3,5,
+     0,3,6;
   Eigen::MatrixXd V(7,3);
   V << 0,0,0,
        1,0,0,
@@ -14,43 +48,12 @@ TEST_CASE("split_nonmanifold: edge-fan", "[igl]")
        0,0,1,
        0,0,-1,
        1,0,1;
-  Eigen::MatrixXi F(5,3);
-  F<<0,1,3,
-     0,3,2,
-     0,4,3,
-     0,3,5,
-     0,3,6;
+  igl::writePLY("edge-fan.ply",V,F);
   Eigen::MatrixXd SV;
   Eigen::MatrixXi SF;
   Eigen::VectorXi SVI;
-  igl::split_nonmanifold(V,F,SV,SF,SVI);
-  Eigen::MatrixXd SVgt(13,3);
-  SVgt<<
-    0,0,0,
-    0,0,0,
-    0,0,0,
-    0,0,0,
-    1,0,0,
-    0,1,0,
-    0,0,1,
-    0,1,0,
-    0,1,0,
-    -1,0,0,
-    0,1,0,
-    0,0,-1,
-    1,0,1;
-  Eigen::MatrixXi SFgt(5,3);
-  SFgt<<
-    0,4,5,
-    0,5,9,
-    1,6,10,
-    2,7,11,
-    3,8,12;
-  Eigen::VectorXi SVIgt(13);
-  SVIgt<<0,0,0,0,1,3,4,3,3,2,3,5,6;
-  test_common::assert_eq(SV,SVgt);
-  test_common::assert_eq(SF,SFgt);
-  test_common::assert_eq(SVI,SVIgt);
+  igl::split_nonmanifold(F,SF,SVI);
+  check_same_but_manifold(F,SF,SVI);
 }
 
 TEST_CASE("split_nonmanifold: vertex-boundary", "[igl]")
@@ -65,40 +68,21 @@ TEST_CASE("split_nonmanifold: vertex-boundary", "[igl]")
   Eigen::MatrixXi F(2,3);
   F<<0,1,2,
      1,3,4;
-  Eigen::MatrixXd SV;
+  {
+    Eigen::MatrixXd V3(V.rows(),3);
+    V3<<V,Eigen::MatrixXd::Zero(V.rows(),1);
+    igl::writePLY("vertex-boundary.ply",V3,F);
+  }
   Eigen::MatrixXi SF;
   Eigen::VectorXi SVI;
-  igl::split_nonmanifold(V,F,SV,SF,SVI);
-  Eigen::MatrixXd SVgt(6,2);
-  SVgt<<
-    0,0,
-    1,0,
-    1,0,
-    2,0,
-    0,1,
-    1,1;
-  Eigen::MatrixXi SFgt(2,3);
-  SFgt<<
-    0,2,4,
-    1,3,5;
-  Eigen::VectorXi SVIgt(6);
-  SVIgt << 0, 1, 1, 3, 2, 4;
-  test_common::assert_eq(SV,SVgt);
-  test_common::assert_eq(SF,SFgt);
-  test_common::assert_eq(SVI,SVIgt);
+  igl::split_nonmanifold(F,SF,SVI);
+  REQUIRE(SVI.rows() == 6);
+  check_same_but_manifold(F,SF,SVI);
 }
 
 TEST_CASE("split_nonmanifold: edge-disk-flap", "[igl]")
 {
   using namespace igl;
-  Eigen::MatrixXd V(6,3);
-  V<<
-    0,0,0,
-    1,0,0,
-    0,1,0,
-    -1,0,0,
-    0,-1,0,
-    0,0,1;
   Eigen::MatrixXi F(5,3);
   F<<
     0,1,2,
@@ -106,44 +90,24 @@ TEST_CASE("split_nonmanifold: edge-disk-flap", "[igl]")
     0,3,4,
     0,4,1,
     0,5,1;
-  Eigen::MatrixXd SV;
-  Eigen::MatrixXi SF;
-  Eigen::VectorXi SVI;
-  igl::split_nonmanifold(V,F,SV,SF,SVI);
-  Eigen::MatrixXd SVgt(8,3);
-  SVgt<<
-    0,0,0,
+  Eigen::MatrixXd V(6,3);
+  V<<
     0,0,0,
     1,0,0,
     0,1,0,
     -1,0,0,
     0,-1,0,
-    0,0,1,
-    1,0,0;
-  Eigen::MatrixXi SFgt(5,3);
-  SFgt<<
-    0,2,3,
-    0,3,4,
-    0,4,5,
-    0,5,2,
-    1,6,7;
-  Eigen::VectorXi SVIgt(8);
-  SVIgt<<  0,  0,  1,  2,  3,  4,  5,  1;
-  test_common::assert_eq(SV,SVgt);
-  test_common::assert_eq(SF,SFgt);
-  test_common::assert_eq(SVI,SVIgt);
+    0,0,1;
+  igl::writePLY("edge-disk-flap.ply",V,F);
+  Eigen::MatrixXi SF;
+  Eigen::VectorXi SVI;
+  igl::split_nonmanifold(F,SF,SVI);
+  check_same_but_manifold(F,SF,SVI);
 }
 
 TEST_CASE("split_nonmanifold: edge-disk-tent", "[igl]")
 {
   using namespace igl;
-  Eigen::MatrixXd V(5,3);
-  V<<
-    0,0,0,
-    1,0,0,
-    -1,1,0,
-    0,-1,0,
-    0,0,1;
   Eigen::MatrixXi F(5,3);
   F<<
     0,1,2,
@@ -151,46 +115,23 @@ TEST_CASE("split_nonmanifold: edge-disk-tent", "[igl]")
     0,3,1,
     0,4,1,
     1,4,3;
-  Eigen::MatrixXd SV;
-  Eigen::MatrixXi SF;
-  Eigen::VectorXi SVI;
-  igl::split_nonmanifold(V,F,SV,SF,SVI);
-  Eigen::MatrixXd SVgt(8,3);
-  SVgt<<
-    0,0,0,
+  Eigen::MatrixXd V(5,3);
+  V<<
     0,0,0,
     1,0,0,
-    1,0,0,
     -1,1,0,
     0,-1,0,
-    0,0,1,
-    0,-1,0;
-  Eigen::MatrixXi SFgt(5,3);
-  SFgt<<
-  0,3,4,
-  0,4,5,
-  0,5,3,
-  1,6,2,
-  2,6,7;
-  Eigen::VectorXi SVIgt(8);
-  SVIgt<<  0,  0,  1,  1,  2,  3,  4,  3;
-  test_common::assert_eq(SV,SVgt);
-  test_common::assert_eq(SF,SFgt);
-  test_common::assert_eq(SVI,SVIgt);
+    0,0,1;
+  igl::writePLY("edge-disk-tent.ply",V,F);
+  Eigen::MatrixXi SF;
+  Eigen::VectorXi SVI;
+  igl::split_nonmanifold(F,SF,SVI);
+  check_same_but_manifold(F,SF,SVI);
 }
 
 TEST_CASE("split_nonmanifold: vertex-kiss", "[igl]")
 {
   using namespace igl;
-  Eigen::MatrixXd V(7,3);
-  V<<
-    0,0,0,
-    1,0,0,
-    0,1,0,
-    0,0,1,
-    0,0,2,
-    1,0,2,
-    0,1,2;
   Eigen::MatrixXi F(6,3);
   F<<
     0,1,3,
@@ -199,46 +140,25 @@ TEST_CASE("split_nonmanifold: vertex-kiss", "[igl]")
     4,5,3,
     5,6,3,
     6,4,3;
-  Eigen::MatrixXd SV;
-  Eigen::MatrixXi SF;
-  Eigen::VectorXi SVI;
-  igl::split_nonmanifold(V,F,SV,SF,SVI);
-  Eigen::MatrixXd SVgt(8,3);
-  SVgt<<
+  Eigen::MatrixXd V(7,3);
+  V<<
     0,0,0,
     1,0,0,
     0,1,0,
+    0,0,1,
     0,0,2,
     1,0,2,
-    0,1,2,
-    0,0,1,
-    0,0,1;
-  Eigen::MatrixXi SFgt(6,3);
-  SFgt<<
-    0,1,6,
-    1,2,6,
-    2,0,6,
-    3,4,7,
-    4,5,7,
-    5,3,7;
-  Eigen::VectorXi SVIgt(8);
-  SVIgt<<  0,  1,  2,  4,  5,  6,  3,  3;
-  test_common::assert_eq(SV,SVgt);
-  test_common::assert_eq(SF,SFgt);
-  test_common::assert_eq(SVI,SVIgt);
+    0,1,2;
+  igl::writePLY("vertex-kiss.ply",V,F);
+  Eigen::MatrixXi SF;
+  Eigen::VectorXi SVI;
+  igl::split_nonmanifold(F,SF,SVI);
+  check_same_but_manifold(F,SF,SVI);
 }
 
 TEST_CASE("split_nonmanifold: non-orientable", "[igl]")
 {
   using namespace igl;
-  Eigen::MatrixXd V(6,3);
-  V<<
-     6, 0, 0,
-     4, 0, 0,
-    -3, 5, 0,
-    -2, 4, 0,
-    -2,-4, 1,
-    -3,-5,-1;
   Eigen::MatrixXi F(6,3);
   F<<
     0,2,1,
@@ -247,31 +167,106 @@ TEST_CASE("split_nonmanifold: non-orientable", "[igl]")
     4,5,3,
     4,1,5,
     1,0,5;
-  Eigen::MatrixXd SV;
+  Eigen::MatrixXd V(6,3);
+  V<<
+     6, 0, 0,
+     4, 0, 0,
+    -3, 5, 0,
+    -2, 4, 0,
+    -2,-4, 1,
+    -3,-5,-1;
+  igl::writePLY("non-orientable.ply",V,F);
   Eigen::MatrixXi SF;
   Eigen::VectorXi SVI;
-  igl::split_nonmanifold(V,F,SV,SF,SVI);
-  Eigen::MatrixXd SVgt(8,3);
-  SVgt<<
-    6,0,0,
-    -3,5,0,
-    -2,-4,1,
-    4,0,0,
-    -2,4,0,
-    -3,-5,-1,
-    6,0,0,
-    4,0,0;
-  Eigen::MatrixXi SFgt(6,3);
-  SFgt<<
-    0,1,7,
-    1,4,7,
-    1,2,4,
-    2,5,4,
-    2,3,5,
-    3,6,5;
-  Eigen::VectorXi SVIgt(8);
-  SVIgt<<  0,  2,  4,  1,  3,  5,  0,  1;
-  test_common::assert_eq(SV,SVgt);
-  test_common::assert_eq(SF,SFgt);
-  test_common::assert_eq(SVI,SVIgt);
+  igl::split_nonmanifold(F,SF,SVI);
+  check_same_but_manifold(F,SF,SVI);
+}
+
+TEST_CASE("split_nonmanifold: flap", "[igl]")
+{
+  Eigen::MatrixXd V(12,3);
+  V<< 
+    1.5,0,0,
+    0.75,0,0.433,
+    0.75,0,-0.433,
+    0,1.5,0,
+    0,0.75,0.433,
+    0,0.75,-0.433,
+    -1.5,0,0,
+    -0.75,0,-0.433,
+    -0,-1.5,0,
+    -0,-0.75,0.433,
+    1.5,0,1,
+    0,1.5,1;
+  Eigen::MatrixXi F(12,3);
+  const auto check = [&F,&V]()
+  {
+    Eigen::MatrixXi SF;
+    Eigen::VectorXi SVI;
+    igl::split_nonmanifold(F,SF,SVI);
+    igl::writePLY("flap.ply",V,F);
+    check_same_but_manifold(F,SF,SVI);
+    {
+      Eigen::VectorXi C;
+      const int nc = igl::facet_components(SF,C);
+      REQUIRE(nc == 2);
+    }
+  };
+  F<< 
+    0,3,1,
+    3,4,1,
+    2,5,0,
+    5,3,0,
+    3,6,4,
+    5,7,3,
+    7,6,3,
+    8,0,9,
+    0,1,9,
+    2,0,8,
+    0,3,11,
+    0,11,10;
+  check();
+  F<< 
+    0,3,11,
+    0,1,9,
+    2,5,0,
+    7,6,3,
+    0,3,1,
+    3,6,4,
+    5,3,0,
+    5,7,3,
+    8,0,9,
+    3,4,1,
+    2,0,8,
+    0,11,10;
+  check();
+}
+
+TEST_CASE("split_nonmanifold: crosses", "[igl]")
+{
+  for(int p = 1;p<7;p++)
+  {
+    Eigen::MatrixXd V;
+    Eigen::MatrixXi F;
+    // n = 2^p
+    const int n = 1<<p;
+    igl::triangulated_grid(n,3,V,F);
+    V.array() -= 0.5;
+    V.conservativeResize(V.rows(),3);
+    V.col(2).setZero();
+    Eigen::MatrixXd U(V.rows(),3);
+    U << V.col(0),V.col(2),V.col(1);
+    Eigen::MatrixXd VV(V.rows()*2,3);
+    VV << V,U;
+    Eigen::MatrixXi FF(F.rows()*2,3);
+    FF << F,F.array()+V.rows();
+    {
+      Eigen::VectorXi I,J;
+      igl::remove_duplicate_vertices(VV,FF,1e-10,V,I,J,F);
+    }
+    Eigen::MatrixXi SF;
+    Eigen::VectorXi SVI;
+    igl::split_nonmanifold(F,SF,SVI);
+    check_same_but_manifold(F,SF,SVI);
+  }
 }