Speed up CDT for coplanar, overlapping faces.

Former-commit-id: 41da37df694ee51c87f6036799debc44f48541a1

include/igl/copyleft/cgal/SelfIntersectMesh.h

@@ -90,14 +90,14 @@ namespace igl
           // Number of self-intersecting triangle pairs
           typedef typename DerivedF::Index Index;
           Index count;
-          typedef std::vector<CGAL::Object> ObjectList;
+          typedef std::vector<std::pair<Index, CGAL::Object>> ObjectList;
           // Using a vector here makes this **not** output sensitive
           Triangles T;
           typedef std::vector<Index> IndexList;
           IndexList lIF;
           // #F-long list of faces with intersections mapping to the order in
           // which they were first found
-          std::map<Index,std::pair<Index,ObjectList> > offending;
+          std::map<Index,ObjectList> offending;
           // Make a short name for the edge map's key
           typedef std::pair<Index,Index> EMK;
           // Make a short name for the type stored at each edge, the edge map's
@@ -143,7 +143,7 @@ namespace igl
           //   A  triangle in 3D
           //   B  triangle in 3D
           //   fa  index of A in F (and key into offending)
-          //   fb  index of A in F (and key into offending)
+          //   fb  index of B in F (and key into offending)
           // Returns true only if A intersects B
           //
           inline bool intersect(
@@ -455,8 +455,7 @@ inline void igl::copyleft::cgal::SelfIntersectMesh<
   if(offending.count(f) == 0)
   {
     // first time marking, initialize with new id and empty list
-    const Index id = offending.size();
-    offending[f] = {id,{}};
+    offending[f] = {};
     for(Index e = 0; e<3;e++)
     {
       // append face to edge's list
@@ -531,8 +530,8 @@ inline bool igl::copyleft::cgal::SelfIntersectMesh<
   {
     // Construct intersection
     CGAL::Object result = CGAL::intersection(A,B);
-    offending[fa].second.push_back(result);
-    offending[fb].second.push_back(result);
+    offending[fa].push_back({fb, result});
+    offending[fb].push_back({fa, result});
   }
   return true;
 }
@@ -630,8 +629,8 @@ inline bool igl::copyleft::cgal::SelfIntersectMesh<
         A.vertex(va),
         *p));
       count_intersection(fa,fb);
-      offending[fa].second.push_back(seg);
-      offending[fb].second.push_back(seg);
+      offending[fa].push_back({fb, seg});
+      offending[fb].push_back({fa, seg});
       return true;
     }else if(CGAL::object_cast<Segment_3 >(&result))
     {
@@ -777,8 +776,8 @@ inline bool igl::copyleft::cgal::SelfIntersectMesh<
       } else
       {
         // Triangle object
-        offending[fa].second.push_back(result);
-        offending[fb].second.push_back(result);
+        offending[fa].push_back({fb, result});
+        offending[fb].push_back({fa, result});
         //cerr<<REDRUM("Coplanar at: "<<fa<<" & "<<fb<<" (double shared).")<<endl;
         return true;
       }

include/igl/copyleft/cgal/remesh_intersections.cpp

@@ -11,6 +11,7 @@
 
 #include <vector>
 #include <map>
+#include <queue>
 #include <unordered_map>
 #include <iostream>
 
@@ -28,8 +29,8 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
         const std::vector<CGAL::Triangle_3<Kernel> > & T,
         const std::map<
         typename DerivedF::Index,
-        std::pair<typename DerivedF::Index,
-        std::vector<CGAL::Object> > > & offending,
+        std::vector<
+        std::pair<typename DerivedF::Index, CGAL::Object> > > & offending,
         const std::map<
         std::pair<typename DerivedF::Index,typename DerivedF::Index>,
         std::vector<typename DerivedF::Index> > & /*edge2faces*/,
@@ -93,24 +94,45 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
             return p1_coeffs[3] < p2_coeffs[3];
         return false;
     };
-    std::map<Plane_3, std::vector<Index>, decltype(plane_comp)>
-        unique_planes(plane_comp);
+    //std::map<Plane_3, std::vector<Index>, decltype(plane_comp)>
+    //    unique_planes(plane_comp);
+
 
     const size_t num_faces = F.rows();
     const size_t num_base_vertices = V.rows();
     assert(num_faces == T.size());
+
     std::vector<bool> is_offending(num_faces, false);
     for (const auto itr : offending) {
         const auto& fid = itr.first;
         is_offending[fid] = true;
 
-        Plane_3 key = T[fid].supporting_plane();
-        assert(!key.is_degenerate());
-        const auto jtr = unique_planes.find(key);
-        if (jtr == unique_planes.end()) {
-            unique_planes.insert({key, {fid}});
-        } else {
-            jtr->second.push_back(fid);
+        //Plane_3 key = T[fid].supporting_plane();
+        //assert(!key.is_degenerate());
+        //const auto jtr = unique_planes.find(key);
+        //if (jtr == unique_planes.end()) {
+        //    unique_planes.insert({key, {fid}});
+        //} else {
+        //    jtr->second.push_back(fid);
+        //}
+    }
+
+    std::unordered_map<Index, std::vector<Index> > intersecting_and_coplanar;
+    for (const auto itr : offending) {
+        const auto& fi = itr.first;
+        const auto P = T[fi].supporting_plane();
+        assert(!P.is_degenerate());
+        for (const auto jtr : itr.second) {
+            const auto& fj = jtr.first;
+            const auto& tj = T[fj];
+            if (P.has_on(tj[0]) && P.has_on(tj[1]) && P.has_on(tj[2])) {
+                auto loc = intersecting_and_coplanar.find(fi);
+                if (loc == intersecting_and_coplanar.end()) {
+                    intersecting_and_coplanar[fi] = {fj};
+                } else {
+                    loc->second.push_back(fj);
+                }
+            }
         }
     }
 
@@ -135,7 +157,9 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
             cdt.insert_constraint(P.to_2d(triangle[2]), P.to_2d(triangle[0]));
 
             if (itr == offending.end()) continue;
-            for (const auto& obj : itr->second.second) {
+            for (const auto& index_obj : itr->second) {
+                const auto& ofid = index_obj.first;
+                const auto& obj = index_obj.second;
                 if(const Segment_3 *iseg = CGAL::object_cast<Segment_3 >(&obj)) {
                     // Add segment constraint
                     cdt.insert_constraint(
@@ -280,16 +304,53 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
     }
 
     // Process self-intersecting faces.
-    for (const auto itr : unique_planes) {
-        Plane_3 P = itr.first;
-        const auto& involved_faces = itr.second;
+    std::vector<bool> processed(num_faces, false);
+    for (const auto itr : offending) {
+        const auto fid = itr.first;
+        if (processed[fid]) continue;
+        processed[fid] = true;
 
+        const auto loc = intersecting_and_coplanar.find(fid);
+        std::vector<Index> involved_faces;
+        if (loc == intersecting_and_coplanar.end()) {
+            involved_faces.push_back(fid);
+        } else {
+            std::queue<Index> Q;
+            Q.push(fid);
+            while (!Q.empty()) {
+                const auto index = Q.front();
+                involved_faces.push_back(index);
+                Q.pop();
+
+                const auto overlapping_faces =
+                    intersecting_and_coplanar.find(index);
+                assert(overlapping_faces != intersecting_and_coplanar.end());
+
+                for (const auto other_index : overlapping_faces->second) {
+                    if (processed[other_index]) continue;
+                    processed[other_index] = true;
+                    Q.push(other_index);
+                }
+            }
+        }
+
+        Plane_3 P = T[fid].supporting_plane();
         std::vector<Point_3> vertices;
         std::vector<std::vector<Index> > faces;
         run_delaunay_triangulation(P, involved_faces, vertices, faces);
         post_triangulation_process(vertices, faces, involved_faces);
     }
 
+    //for (const auto itr : unique_planes) {
+    //    Plane_3 P = itr.first;
+    //    const auto& involved_faces = itr.second;
+
+    //    std::vector<Point_3> vertices;
+    //    std::vector<std::vector<Index> > faces;
+    //    run_delaunay_triangulation(P, involved_faces, vertices, faces);
+    //    post_triangulation_process(vertices, faces, involved_faces);
+    //}
+
     // Output resolved mesh.
     const size_t num_out_vertices = new_vertices.size() + num_base_vertices;
     VV.resize(num_out_vertices, 3);

include/igl/copyleft/cgal/remesh_intersections.h

@@ -53,8 +53,8 @@ namespace igl
         const std::vector<CGAL::Triangle_3<Kernel> > & T,
         const std::map<
           typename DerivedF::Index,
-          std::pair<typename DerivedF::Index,
-            std::vector<CGAL::Object> > > & offending,
+            std::vector<
+            std::pair<typename DerivedF::Index, CGAL::Object> > > & offending,
         const std::map<
           std::pair<typename DerivedF::Index,typename DerivedF::Index>,
           std::vector<typename DerivedF::Index> > & edge2faces,