bug fix in boolean, first attempt to handle co-planar facets correctly

Former-commit-id: 7f5ec3cd232c9d0056ad050c64f3753b0ad1c108

include/igl/boolean/mesh_boolean.cpp

@@ -233,16 +233,55 @@ IGL_INLINE void igl::mesh_boolean(
 #ifdef IGL_MESH_BOOLEAN_DEBUG
   cout<<"clean..."<<endl;
 #endif
-  // remove duplicates: cancel out in all cases? assertion in others?
+  // Deal with duplicate faces
   {
-    MatrixXi oldG = G;
-    VectorXi IA,_1;
-    unique_simplices(oldG,G,IA,_1);
-    assert(IA.rows() == G.rows());
-    J.resize(IA.rows(),1);
-    for(size_t j = 0;j<(size_t)J.size();j++)
+    VectorXi IA,IC;
+    MatrixX3I uG;
+    unique_simplices(G,uG,IA,IC);
+    assert(IA.rows() == uG.rows());
+    // faces ontop of each unique face
+    vector<vector<Index> > uG2G(uG.rows());
+    // signed counts
+    VectorXi counts = VectorXi::Zero(uG.rows());
+    // loop over all faces
+    for(Index g = 0;g<gm;g++)
     {
-      J(j) = GJ(IA(j));
+      const int ug = IC(g);
+      assert(ug < uG2G.size());
+      uG2G[ug].push_back(g);
+      // is uG(g,:) just a rotated version of G(g,:) ?
+      const bool consistent = 
+        (G(g,0) == uG(ug,0) && G(g,1) == uG(ug,1) && G(g,2) == uG(ug,2)) ||
+        (G(g,0) == uG(ug,1) && G(g,1) == uG(ug,2) && G(g,2) == uG(ug,0)) ||
+        (G(g,0) == uG(ug,2) && G(g,1) == uG(ug,0) && G(g,2) == uG(ug,1));
+      counts(ug) += consistent ? 1 : -1;
+    }
+    MatrixX3I oldG = G;
+    // Faces of output vG[i] = j means ith face of output should be jth face in
+    // oldG
+    vG.clear();
+    for(size_t ug = 0;ug < uG2G.size();ug++)
+    {
+      // if signed occurrences is zero or ±two then keep none
+      // else if signed occurrences is ±one then keep just one facet
+      if(abs(counts(ug)) == 1)
+      {
+        assert(uG2G.size() > 0);
+        vG.push_back(uG2G[ug][0]);
+      }
+#ifdef IGL_MESH_BOOLEAN_DEBUG
+      else
+      {
+        cout<<"Skipping "<<uG2G.size()<<" facets..."<<endl;
+      }
+#endif
+    }
+    G.resize(vG.size(),3);
+    J.resize(vG.size());
+    for(size_t g = 0;g<vG.size();g++)
+    {
+      G.row(g) = oldG.row(vG[g]);
+      J(g) = GJ(vG[g]);
     }
   }
   // remove unreferenced vertices

include/igl/outer_hull.cpp

@@ -1,15 +1,12 @@
 #include "outer_hull.h"
 #include "outer_facet.h"
-#include "sort.h"
+#include "sortrows.h"
 #include "facet_components.h"
 #include "winding_number.h"
 #include "triangle_triangle_adjacency.h"
 #include "unique_edge_map.h"
 #include "barycenter.h"
 #include "per_face_normals.h"
-#include "all_edges.h"
-#include "colon.h"
-#include "get_seconds.h"
 
 #include <Eigen/Geometry>
 #include <vector>
@@ -64,8 +61,11 @@ IGL_INLINE void igl::outer_hull(
   // uE --> sorted order index 
   // uE --> bool, whether needed to flip face to make "consistent" with unique
   //   edge
+  // Place order of each half-edge in its corresponding sorted list around edge
   VectorXI diIM(3*m);
+  // Whether face's edge used for sorting is consistent with unique edge
   VectorXI dicons(3*m);
+  // dihedral angles of faces around edge with face of edge in dicons
   vector<vector<typename DerivedV::Scalar> > di(uE2E.size());
   // For each list of face-edges incide on a unique edge
   for(size_t ui = 0;ui<(size_t)uE.rows();ui++)
@@ -73,7 +73,7 @@ IGL_INLINE void igl::outer_hull(
     // Base normal vector to orient against
     const auto fe0 = uE2E[ui][0];
     const RowVector3N & eVp = N.row(fe0%m);
-    di[ui].resize(uE2E[ui].size());
+    MatrixXd di_I(uE2E[ui].size(),2);
 
     const typename DerivedF::Scalar d = F(fe0%m,((fe0/m)+2)%3);
     const typename DerivedF::Scalar s = F(fe0%m,((fe0/m)+1)%3);
@@ -94,26 +94,38 @@ IGL_INLINE void igl::outer_hull(
       assert(!cons[fei] || (d == F(f,(c+1)%3)));
       // Angle between n and f
       const RowVector3N & n = N.row(f);
-      di[ui][fei] = M_PI - atan2( eVp.cross(n).dot(eV), eVp.dot(n));
+      di_I(fei,0) = M_PI - atan2( eVp.cross(n).dot(eV), eVp.dot(n));
       if(!cons[fei])
       {
-        di[ui][fei] = di[ui][fei] + M_PI;
-        if(di[ui][fei]>2.*M_PI)
+        di_I(fei,0) = di_I(fei,0) + M_PI;
+        if(di_I(fei,0)>=2.*M_PI)
         {
-          di[ui][fei] = di[ui][fei] - 2.*M_PI;
+          di_I(fei,0) = di_I(fei,0) - 2.*M_PI;
         }
       }
+      di_I(fei,1) = f;
+    }
+    VectorXi IM;
+
+    //igl::sort(di[ui],true,di[ui],IM);
+    // Sort, but break ties using index to ensure that duplicates always show
+    // up in same order.
+    MatrixXd s_di_I;
+    igl::sortrows(di_I,true,s_di_I,IM);
+    di[ui].resize(uE2E[ui].size());
+    for(size_t i = 0;i<di[ui].size();i++)
+    {
+      di[ui][i] = s_di_I(i,0);
     }
-    vector<size_t> IM;
-    igl::sort(di[ui],true,di[ui],IM);
+
     // copy old list
     vector<typename DerivedF::Index> temp = uE2E[ui];
     for(size_t fei = 0;fei<uE2E[ui].size();fei++)
     {
-      uE2E[ui][fei] = temp[IM[fei]];
+      uE2E[ui][fei] = temp[IM(fei)];
       const auto & fe = uE2E[ui][fei];
       diIM(fe) = fei;
-      dicons(fe) = cons[IM[fei]];
+      dicons(fe) = cons[IM(fei)];
     }
 
   }
@@ -206,21 +218,24 @@ IGL_INLINE void igl::outer_hull(
       // first time seeing face
       if(!FH[f])
       {
+//#ifdef IGL_OUTER_HULL_DEBUG
+//      cout<<(f+1)<<endl;
+//#endif
         FH[f] = true;
         FHcount++;
       }
-      // find overlapping face-edges
-      const auto & neighbors = uE2E[EMAP(e)];
-      // normal after possible flipping 
-      const auto & fN = (flip(f)?-1.:1.)*N.row(f);
       // source of edge according to f
       const int fs = flip(f)?F(f,(c+2)%3):F(f,(c+1)%3);
       // destination of edge according to f
       const int fd = flip(f)?F(f,(c+1)%3):F(f,(c+2)%3);
-      // Edge vector according to f's (flipped) orientation.
-      const auto & eV = (V.row(fd)-V.row(fs)).normalized();
       // edge valence
       const size_t val = uE2E[EMAP(e)].size();
+      //// find overlapping face-edges
+      //const auto & neighbors = uE2E[EMAP(e)];
+      //// normal after possible flipping 
+      //const auto & fN = (flip(f)?-1.:1.)*N.row(f);
+      //// Edge vector according to f's (flipped) orientation.
+      ////const auto & eV = (V.row(fd)-V.row(fs)).normalized();
 
 //#warning "EXPERIMENTAL, DO NOT USE"
       //// THIS IS WRONG! The first face is---after sorting---no longer the face
@@ -228,8 +243,34 @@ IGL_INLINE void igl::outer_hull(
       //const auto ui = EMAP(e);
       //const auto fe0 = uE2E[ui][0];
       //const auto es = F(fe0%m,((fe0/m)+1)%3);
+
+      // is edge consistent with edge of face used for sorting
       const int e_cons = (dicons(e) ? 1: -1);
-      const int nfei = (diIM(e) + val + e_cons*(flip(f)?-1:1))%val;
+      int nfei = -1;
+      // Loop once around trying to find suitable next face
+      for(size_t step = 1; step<val+2;step++)
+      {
+        nfei = (diIM(e) + 2*val + e_cons*step*(flip(f)?-1:1))%val;
+#ifdef IGL_OUTER_HULL_DEBUG
+        const int nf = uE2E[EMAP(e)][nfei] % m;
+#endif
+        // Don't consider faces with identical dihedral angles
+        if(di[EMAP(e)][diIM(e)] != di[EMAP(e)][nfei])
+//#warning "THIS IS HACK, FIX ME"
+//        if( abs(di[EMAP(e)][diIM(e)] - di[EMAP(e)][nfei]) < 1e-16 )
+        {
+//#ifdef IGL_OUTER_HULL_DEBUG
+//        cout<<"Next facet: "<<(f+1)<<" --> "<<(nf+1)<<", |"<<
+//          di[EMAP(e)][diIM(e)]<<" - "<<di[EMAP(e)][nfei]<<"| = "<<
+//            abs(di[EMAP(e)][diIM(e)] - di[EMAP(e)][nfei])
+//            <<endl;
+//#endif
+          break;
+        }
+#ifdef IGL_OUTER_HULL_DEBUG
+        cout<<"Skipping co-planar facet: "<<(f+1)<<" --> "<<(nf+1)<<endl;
+#endif
+      }
       const int max_ne_2 = uE2E[EMAP(e)][nfei];
 
       int max_ne = -1;

include/igl/peel_outer_hull_layers.cpp

@@ -6,6 +6,7 @@
 //#define IGL_PEEL_OUTER_HULL_LAYERS_DEBUG
 #ifdef IGL_PEEL_OUTER_HULL_LAYERS_DEBUG
 #include "writePLY.h"
+#include "writeDMAT.h"
 #include "STR.h"
 #endif
 
@@ -59,6 +60,7 @@ IGL_INLINE size_t igl::peel_outer_hull_layers(
 #ifdef IGL_PEEL_OUTER_HULL_LAYERS_DEBUG
   cout<<"calling outer hull..."<<endl;
   writePLY(STR("outer-hull-input-"<<iter<<".ply"),V,Fr);
+  writeDMAT(STR("outer-hull-input-"<<iter<<".dmat"),Nr);
 #endif
     outer_hull(V,Fr,Nr,Fo,Jo,flipr);
 #ifdef IGL_PEEL_OUTER_HULL_LAYERS_DEBUG

include/igl/sort.cpp

@@ -225,4 +225,5 @@ template void igl::sort<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<do
 template void igl::sort_new<Eigen::Matrix<int, 1, 6, 1, 1, 6>, Eigen::Matrix<int, 1, 6, 1, 1, 6>, Eigen::Matrix<int, 1, 6, 1, 1, 6> >(Eigen::PlainObjectBase<Eigen::Matrix<int, 1, 6, 1, 1, 6> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, 1, 6, 1, 1, 6> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, 1, 6, 1, 1, 6> >&);
 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, 4, 0, -1, 4>, Eigen::Matrix<double, -1, 4, 0, -1, 4>, Eigen::Matrix<int, -1, 4, 0, -1, 4> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 4, 0, -1, 4> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 4, 0, -1, 4> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 4, 0, -1, 4> >&);
+template void igl::sort<Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > 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/sortrows.h

@@ -17,20 +17,21 @@ namespace igl
   //
   // Templates:
   //   DerivedX derived scalar type, e.g. MatrixXi or MatrixXd
-  //   DerivedIX derived integer type, e.g. MatrixXi
+  //   DerivedI derived integer type, e.g. MatrixXi
   // Inputs:
   //   X  m by n matrix whose entries are to be sorted
   //   ascending  sort ascending (true, matlab default) or descending (false)
   // Outputs:
-  //   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>
+  //   Y  m by n matrix whose entries are sorted (**should not** be same
+  //     reference as X)
+  //   I  m list of indices so that
+  //     Y = X(I,:);
+  template <typename DerivedX, typename DerivedI>
   IGL_INLINE void sortrows(
     const Eigen::PlainObjectBase<DerivedX>& X,
     const bool ascending,
     Eigen::PlainObjectBase<DerivedX>& Y,
-    Eigen::PlainObjectBase<DerivedIX>& IX);
+    Eigen::PlainObjectBase<DerivedI>& I);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/unique_simplices.cpp

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

include/igl/writeDMAT.cpp

@@ -130,4 +130,5 @@ template bool igl::writeDMAT<double>(std::basic_string<char, std::char_traits<ch
 template bool igl::writeDMAT<Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, bool);
 template bool igl::writeDMAT<Eigen::Matrix<double, -1, 1, 0, -1, 1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::Matrix<double, -1, 1, 0, -1, 1> const&, bool);
 template bool igl::writeDMAT<Eigen::Matrix<int, -1, 1, 0, -1, 1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, bool);
+template bool igl::writeDMAT<Eigen::Matrix<double, -1, 3, 0, -1, 3> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::Matrix<double, -1, 3, 0, -1, 3> const&, bool);
 #endif