refactor remeshing

Former-commit-id: ea0a88647b7b9868f39691cee007fc229fbe414d

include/igl/cgal/SelfIntersectMesh.h

@@ -89,10 +89,13 @@ namespace igl
         typedef typename DerivedF::Index Index;
         Index count;
         typedef std::vector<CGAL::Object> ObjectList;
+        // Using a vector here makes this **not** output sensitive
         std::vector<ObjectList > F_objects;
         Triangles T;
         typedef std::vector<Index> IndexList;
         IndexList lIF;
+        // #F-long list of bools revealing whether face is participating in any
+        // intersections
         std::vector<bool> offensive;
         std::vector<Index> offending_index;
         std::vector<Index> offending;
@@ -209,7 +212,7 @@ namespace igl
         //   A_objects_3  updated list of intersection objects for A
         // Outputs:
         //   cdt  Contrained delaunay triangulation in projected 2D plane
-        inline void projected_delaunay(
+        static inline void projected_delaunay(
             const Triangle_3 & A,
             const ObjectList & A_objects_3,
             CDT_plus_2 & cdt);
@@ -414,241 +417,260 @@ inline igl::cgal::SelfIntersectMesh<
     return;
   }
 
-  int NF_count = 0;
-  // list of new faces, we'll fix F later
-  vector<
-    typename Eigen::Matrix<typename DerivedFF::Scalar,Dynamic,Dynamic>
-    > NF(offending.size());
-  // list of new vertices
-  typedef vector<Point_3> Point_3List;
-  Point_3List NV;
-  Index NV_count = 0;
-  vector<CDT_plus_2> cdt(offending.size());
-  vector<Plane_3> P(offending.size());
-  // Use map for *all* faces
-  map<typename CDT_plus_2::Vertex_handle,Index> v2i;
-  // Loop over offending triangles
-  const size_t noff = offending.size();
+  const auto & remesh = [](
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    const std::vector<ObjectList > & F_objects,
+    const Triangles & T,
+    const std::vector<bool> & offensive,
+    const std::vector<Index> & offending_index,
+    const std::vector<Index> & offending,
+    const EdgeMap & edge2faces,
+    Eigen::PlainObjectBase<DerivedVV> & VV,
+    Eigen::PlainObjectBase<DerivedFF> & FF,
+    Eigen::PlainObjectBase<DerivedJ> & J,
+    Eigen::PlainObjectBase<DerivedIM> & IM
+    )
+  {
+    int NF_count = 0;
+    // list of new faces, we'll fix F later
+    vector<
+      typename Eigen::Matrix<typename DerivedFF::Scalar,Dynamic,Dynamic>
+      > NF(offending.size());
+    // list of new vertices
+    typedef vector<Point_3> Point_3List;
+    Point_3List NV;
+    Index NV_count = 0;
+    vector<CDT_plus_2> cdt(offending.size());
+    vector<Plane_3> P(offending.size());
+    // Use map for *all* faces
+    map<typename CDT_plus_2::Vertex_handle,Index> v2i;
+    // Loop over offending triangles
+    const size_t noff = offending.size();
 #ifdef IGL_SELFINTERSECTMESH_DEBUG
-  double t_proj_del = 0;
+    double t_proj_del = 0;
 #endif
-  // Unfortunately it looks like CGAL has trouble allocating memory when
-  // multiple openmp threads are running. Crashes durring CDT...
-//# pragma omp parallel for if (noff>1000)
-  for(Index o = 0;o<(Index)noff;o++)
-  {
-    // index in F
-    const Index f = offending[o];
+    // Unfortunately it looks like CGAL has trouble allocating memory when
+    // multiple openmp threads are running. Crashes durring CDT...
+  //# pragma omp parallel for if (noff>1000)
+    for(Index o = 0;o<(Index)noff;o++)
     {
+      // index in F
+      const Index f = offending[o];
+      {
 #ifdef IGL_SELFINTERSECTMESH_DEBUG
-      const double t_before = get_seconds();
+        const double t_before = get_seconds();
 #endif
-      projected_delaunay(T[f],F_objects[f],cdt[o]);
+        projected_delaunay(T[f],F_objects[f],cdt[o]);
 #ifdef IGL_SELFINTERSECTMESH_DEBUG
-      t_proj_del += (get_seconds()-t_before);
+        t_proj_del += (get_seconds()-t_before);
 #endif
-    }
-    // Q: Is this also delaunay in 3D?
-    // A: No, because the projection is affine and delaunay is not affine
-    // invariant.
-    // Q: Then, can't we first get the 2D delaunay triangulation, then lift it
-    // to 3D and flip any offending edges?
-    // Plane of projection (also used by projected_delaunay)
-    P[o] = Plane_3(T[f].vertex(0),T[f].vertex(1),T[f].vertex(2));
-    // Build index map
-    {
-      Index i=0;
-      for(
-        typename CDT_plus_2::Finite_vertices_iterator vit = cdt[o].finite_vertices_begin();
-        vit != cdt[o].finite_vertices_end();
-        ++vit)
+      }
+      // Q: Is this also delaunay in 3D?
+      // A: No, because the projection is affine and delaunay is not affine
+      // invariant.
+      // Q: Then, can't we first get the 2D delaunay triangulation, then lift it
+      // to 3D and flip any offending edges?
+      // Plane of projection (also used by projected_delaunay)
+      P[o] = Plane_3(T[f].vertex(0),T[f].vertex(1),T[f].vertex(2));
+      // Build index map
       {
-        if(i<3)
+        Index i=0;
+        for(
+          typename CDT_plus_2::Finite_vertices_iterator vit = cdt[o].finite_vertices_begin();
+          vit != cdt[o].finite_vertices_end();
+          ++vit)
         {
-          //cout<<T[f].vertex(i)<<
-          //  (T[f].vertex(i) == P[o].to_3d(vit->point())?" == ":" != ")<<
-          //  P[o].to_3d(vit->point())<<endl;
+          if(i<3)
+          {
+            //cout<<T[f].vertex(i)<<
+            //  (T[f].vertex(i) == P[o].to_3d(vit->point())?" == ":" != ")<<
+            //  P[o].to_3d(vit->point())<<endl;
 #ifndef NDEBUG
-          // I want to be sure that the original corners really show up as the
-          // original corners of the CDT. I.e. I don't trust CGAL to maintain
-          // the order
-          assert(T[f].vertex(i) == P[o].to_3d(vit->point()));
+            // I want to be sure that the original corners really show up as the
+            // original corners of the CDT. I.e. I don't trust CGAL to maintain
+            // the order
+            assert(T[f].vertex(i) == P[o].to_3d(vit->point()));
 #endif
-          // For first three, use original index in F
-//#   pragma omp critical
-          v2i[vit] = F(f,i);
-        }else
-        {
-          const Point_3 vit_point_3 = P[o].to_3d(vit->point());
-          // First look up each edge's neighbors to see if exact point has
-          // already been added (This makes everything a bit quadratic)
-          bool found = false;
-          for(int e = 0; e<3 && !found;e++)
+            // For first three, use original index in F
+  //#   pragma omp critical
+            v2i[vit] = F(f,i);
+          }else
           {
-            // Index of F's eth edge in V
-            Index i = F(f,(e+1)%3);
-            Index j = F(f,(e+2)%3);
-            // Be sure that i<j
-            if(i>j)
+            const Point_3 vit_point_3 = P[o].to_3d(vit->point());
+            // First look up each edge's neighbors to see if exact point has
+            // already been added (This makes everything a bit quadratic)
+            bool found = false;
+            for(int e = 0; e<3 && !found;e++)
             {
-              swap(i,j);
-            }
-            assert(edge2faces.count(EMK(i,j))==1);
-            // loop over neighbors
-            for(
-              typename IndexList::const_iterator nit = edge2faces[EMK(i,j)].begin();
-              nit != edge2faces[EMK(i,j)].end() && !found;
-              nit++)
-            {
-              // don't consider self
-              if(*nit == f)
+              // Index of F's eth edge in V
+              Index i = F(f,(e+1)%3);
+              Index j = F(f,(e+2)%3);
+              // Be sure that i<j
+              if(i>j)
               {
-                continue;
+                swap(i,j);
               }
-              // index of neighbor in offending (to find its cdt)
-              Index no = offending_index[*nit];
-              // Loop over vertices of that neighbor's cdt (might not have been
-              // processed yet, but then it's OK because it'll just be empty)
+              assert(edge2faces.count(EMK(i,j))==1);
+              const EMV & facesij = edge2faces.find(EMK(i,j))->second;
+              // loop over neighbors
               for(
-                  typename CDT_plus_2::Finite_vertices_iterator uit = cdt[no].finite_vertices_begin();
-                  uit != cdt[no].finite_vertices_end() && !found;
-                  ++uit)
+                typename IndexList::const_iterator nit = facesij.begin();
+                nit != facesij.end() && !found;
+                nit++)
               {
-                if(vit_point_3 == P[no].to_3d(uit->point()))
+                // don't consider self
+                if(*nit == f)
+                {
+                  continue;
+                }
+                // index of neighbor in offending (to find its cdt)
+                Index no = offending_index[*nit];
+                // Loop over vertices of that neighbor's cdt (might not have been
+                // processed yet, but then it's OK because it'll just be empty)
+                for(
+                    typename CDT_plus_2::Finite_vertices_iterator uit = cdt[no].finite_vertices_begin();
+                    uit != cdt[no].finite_vertices_end() && !found;
+                    ++uit)
                 {
-                  assert(v2i.count(uit) == 1);
-//#   pragma omp critical
-                  v2i[vit] = v2i[uit];
-                  found = true;
+                  if(vit_point_3 == P[no].to_3d(uit->point()))
+                  {
+                    assert(v2i.count(uit) == 1);
+  //#   pragma omp critical
+                    v2i[vit] = v2i[uit];
+                    found = true;
+                  }
                 }
               }
             }
-          }
-          if(!found)
-          {
-//#   pragma omp critical
+            if(!found)
             {
-              v2i[vit] = V.rows()+NV_count;
-              NV.push_back(vit_point_3);
-              NV_count++;
+  //#   pragma omp critical
+              {
+                v2i[vit] = V.rows()+NV_count;
+                NV.push_back(vit_point_3);
+                NV_count++;
+              }
             }
           }
+          i++;
         }
-        i++;
       }
-    }
-    {
-      Index i = 0;
-      // Resize to fit new number of triangles
-      NF[o].resize(cdt[o].number_of_faces(),3);
-//#   pragma omp atomic
-      NF_count+=NF[o].rows();
-      // Append new faces to NF
-      for(
-        typename CDT_plus_2::Finite_faces_iterator fit = cdt[o].finite_faces_begin();
-        fit != cdt[o].finite_faces_end();
-        ++fit)
       {
-        NF[o](i,0) = v2i[fit->vertex(0)];
-        NF[o](i,1) = v2i[fit->vertex(1)];
-        NF[o](i,2) = v2i[fit->vertex(2)];
-        i++;
+        Index i = 0;
+        // Resize to fit new number of triangles
+        NF[o].resize(cdt[o].number_of_faces(),3);
+  //#   pragma omp atomic
+        NF_count+=NF[o].rows();
+        // Append new faces to NF
+        for(
+          typename CDT_plus_2::Finite_faces_iterator fit = cdt[o].finite_faces_begin();
+          fit != cdt[o].finite_faces_end();
+          ++fit)
+        {
+          NF[o](i,0) = v2i[fit->vertex(0)];
+          NF[o](i,1) = v2i[fit->vertex(1)];
+          NF[o](i,2) = v2i[fit->vertex(2)];
+          i++;
+        }
       }
     }
-  }
 #ifdef IGL_SELFINTERSECTMESH_DEBUG
-  cout<<"CDT: "<<tictoc()<<"  "<<t_proj_del<<endl;
+    cout<<"CDT: "<<tictoc()<<"  "<<t_proj_del<<endl;
 #endif
 
-  assert(NV_count == (Index)NV.size());
-  // Build output
+    assert(NV_count == (Index)NV.size());
+    // Build output
 #ifndef NDEBUG
-  {
-    Index off_count = 0;
-    for(Index f = 0;f<F.rows();f++)
     {
-      off_count+= (offensive[f]?1:0);
+      Index off_count = 0;
+      for(Index f = 0;f<F.rows();f++)
+      {
+        off_count+= (offensive[f]?1:0);
+      }
+      assert(off_count==(Index)offending.size());
     }
-    assert(off_count==(Index)offending.size());
-  }
 #endif
-  // Append faces
-  FF.resize(F.rows()-offending.size()+NF_count,3);
-  J.resize(FF.rows());
-  // First append non-offending original faces
-  // There's an Eigen way to do this in one line but I forget
-  Index off = 0;
-  for(Index f = 0;f<F.rows();f++)
-  {
-    if(!offensive[f])
-    {
-      FF.row(off) = F.row(f);
-      J(off) = f;
-      off++;
-    }
-  }
-  assert(off == (Index)(F.rows()-offending.size()));
-  // Now append replacement faces for offending faces
-  for(Index o = 0;o<(Index)offending.size();o++)
-  {
-    FF.block(off,0,NF[o].rows(),3) = NF[o];
-    J.block(off,0,NF[o].rows(),1).setConstant(offending[o]);
-    off += NF[o].rows();
-  }
-  // Append vertices
-  VV.resize(V.rows()+NV_count,3);
-  VV.block(0,0,V.rows(),3) = V.template cast<typename DerivedVV::Scalar>();
-  {
-    Index i = 0;
-    for(
-      typename Point_3List::const_iterator nvit = NV.begin();
-      nvit != NV.end();
-      nvit++)
+    // Append faces
+    FF.resize(F.rows()-offending.size()+NF_count,3);
+    J.resize(FF.rows());
+    // First append non-offending original faces
+    // There's an Eigen way to do this in one line but I forget
+    Index off = 0;
+    for(Index f = 0;f<F.rows();f++)
     {
-      for(Index d = 0;d<3;d++)
+      if(!offensive[f])
       {
-        const Point_3 & p = *nvit;
-        // Don't convert via double if output type is same as Kernel
-        to_output_type(p[d], VV(V.rows()+i,d));
+        FF.row(off) = F.row(f);
+        J(off) = f;
+        off++;
       }
-      i++;
     }
-  }
-  IM.resize(VV.rows(),1);
-  map<Point_3,Index> vv2i;
-  // Safe to check for duplicates using double for original vertices: if
-  // incoming reps are different then the points are unique.
-  for(Index v = 0;v<V.rows();v++)
-  {
-    const Point_3 p(V(v,0),V(v,1),V(v,2));
-    if(vv2i.count(p)==0)
+    assert(off == (Index)(F.rows()-offending.size()));
+    // Now append replacement faces for offending faces
+    for(Index o = 0;o<(Index)offending.size();o++)
     {
-      vv2i[p] = v;
+      FF.block(off,0,NF[o].rows(),3) = NF[o];
+      J.block(off,0,NF[o].rows(),1).setConstant(offending[o]);
+      off += NF[o].rows();
     }
-    assert(vv2i.count(p) == 1);
-    IM(v) = vv2i[p];
-  }
-  // Must check for duplicates of new vertices using exact.
-  {
-    Index v = V.rows();
-    for(
-      typename Point_3List::const_iterator nvit = NV.begin();
-      nvit != NV.end();
-      nvit++)
+    // Append vertices
+    VV.resize(V.rows()+NV_count,3);
+    VV.block(0,0,V.rows(),3) = V.template cast<typename DerivedVV::Scalar>();
     {
-      const Point_3 & p = *nvit;
+      Index i = 0;
+      for(
+        typename Point_3List::const_iterator nvit = NV.begin();
+        nvit != NV.end();
+        nvit++)
+      {
+        for(Index d = 0;d<3;d++)
+        {
+          const Point_3 & p = *nvit;
+          // Don't convert via double if output type is same as Kernel
+          to_output_type(p[d], VV(V.rows()+i,d));
+        }
+        i++;
+      }
+    }
+    IM.resize(VV.rows(),1);
+    map<Point_3,Index> vv2i;
+    // Safe to check for duplicates using double for original vertices: if
+    // incoming reps are different then the points are unique.
+    for(Index v = 0;v<V.rows();v++)
+    {
+      const Point_3 p(V(v,0),V(v,1),V(v,2));
       if(vv2i.count(p)==0)
       {
         vv2i[p] = v;
       }
       assert(vv2i.count(p) == 1);
       IM(v) = vv2i[p];
-      v++;
     }
-  }
+    // Must check for duplicates of new vertices using exact.
+    {
+      Index v = V.rows();
+      for(
+        typename Point_3List::const_iterator nvit = NV.begin();
+        nvit != NV.end();
+        nvit++)
+      {
+        const Point_3 & p = *nvit;
+        if(vv2i.count(p)==0)
+        {
+          vv2i[p] = v;
+        }
+        assert(vv2i.count(p) == 1);
+        IM(v) = vv2i[p];
+        v++;
+      }
+    }
 #ifdef IGL_SELFINTERSECTMESH_DEBUG
-  cout<<"Output + dupes: "<<tictoc()<<endl;
+    cout<<"Output + dupes: "<<tictoc()<<endl;
 #endif
+  };
+  remesh(
+    V,F,F_objects,T,offensive,offending_index,offending,edge2faces,VV,FF,J,IM);
 
   // Q: Does this give the same result as TETGEN?
   // A: For the cow and beast, yes.
@@ -660,6 +682,7 @@ inline igl::cgal::SelfIntersectMesh<
   // CGAL implementation on the beast takes 98 seconds but tetgen detection
   // takes 14 seconds
 
+
 }