Refactoring and add thread support.

Former-commit-id: 437dc49e8fcbcfe7ee6917d8072b0fb42f17aaf4

include/igl/copyleft/cgal/SelfIntersectMesh.h

@@ -15,6 +15,8 @@
 #include <list>
 #include <map>
 #include <vector>
+#include <thread>
+#include <mutex>
 
 #define IGL_SELFINTERSECTMESH_DEBUG
 #ifndef IGL_FIRST_HIT_EXCEPTION
@@ -29,6 +31,23 @@ namespace igl
   {
     namespace cgal
     {
+      template<typename Kernel>
+      class DeepCoper {
+        public:
+          typename Kernel::FT operator()(const typename Kernel::FT& x) {
+            return typename Kernel::FT(x);
+          }
+      };
+
+      template<>
+      class DeepCoper<CGAL::Exact_predicates_exact_constructions_kernel> {
+        public:
+          typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
+          typename Kernel::FT operator()(const typename Kernel::FT& x) {
+            return Kernel::FT(x.exact() +1) - 1;
+          }
+      };
+
       // Kernel is a CGAL kernel like:
       //     CGAL::Exact_predicates_inexact_constructions_kernel
       // or 
@@ -222,6 +241,8 @@ namespace igl
             SelfIntersectMesh * SIM, 
             const Box &a, 
             const Box &b);
+        private:
+          std::mutex m_offending_lock;
       };
     }
   }
@@ -495,6 +516,7 @@ inline void igl::copyleft::cgal::SelfIntersectMesh<
   const Index fa,
   const Index fb)
 {
+  std::lock_guard<std::mutex> guard(m_offending_lock);
   mark_offensive(fa);
   mark_offensive(fb);
   this->count++;
@@ -538,6 +560,7 @@ inline bool igl::copyleft::cgal::SelfIntersectMesh<
   {
     // Construct intersection
     CGAL::Object result = CGAL::intersection(A,B);
+    std::lock_guard<std::mutex> guard(m_offending_lock);
     offending[fa].push_back({fb, result});
     offending[fb].push_back({fa, result});
   }
@@ -637,6 +660,7 @@ inline bool igl::copyleft::cgal::SelfIntersectMesh<
         A.vertex(va),
         *p));
       count_intersection(fa,fb);
+      std::lock_guard<std::mutex> guard(m_offending_lock);
       offending[fa].push_back({fb, seg});
       offending[fb].push_back({fa, seg});
       return true;
@@ -695,6 +719,8 @@ inline bool igl::copyleft::cgal::SelfIntersectMesh<
 {
   using namespace std;
 
+  {
+    std::lock_guard<std::mutex> guard(m_offending_lock);
   // must be co-planar
   if(
     A.supporting_plane() != B.supporting_plane() &&
@@ -702,6 +728,7 @@ inline bool igl::copyleft::cgal::SelfIntersectMesh<
   {
     return false;
   }
+  }
   // Since A and B are non-degenerate the intersection must be a polygon
   // (triangle). Either
   //   - the vertex of A (B) opposite the shared edge of lies on B (A), or
@@ -784,6 +811,7 @@ inline bool igl::copyleft::cgal::SelfIntersectMesh<
       } else
       {
         // Triangle object
+        std::lock_guard<std::mutex> guard(m_offending_lock);
         offending[fa].push_back({fb, result});
         offending[fb].push_back({fa, result});
         //cerr<<REDRUM("Coplanar at: "<<fa<<" & "<<fb<<" (double shared).")<<endl;
@@ -854,76 +882,99 @@ inline void igl::copyleft::cgal::SelfIntersectMesh<
   DerivedJ,
   DerivedIM>::process_intersecting_boxes()
 {
-  for (const auto& box_pair : candidate_box_pairs) {
-    const auto& a = box_pair.first;
-    const auto& b = box_pair.second;
-    Index fa = a.handle()-T.begin();
-    Index fb = b.handle()-T.begin();
-    const Triangle_3 & A = *a.handle();
-    const Triangle_3 & B = *b.handle();
-
-    // Number of combinatorially shared vertices
-    Index comb_shared_vertices = 0;
-    // Number of geometrically shared vertices (*not* including combinatorially
-    // shared)
-    Index geo_shared_vertices = 0;
-    // Keep track of shared vertex indices
-    std::vector<std::pair<Index,Index> > shared;
-    Index ea,eb;
-    for(ea=0;ea<3;ea++)
-    {
-      for(eb=0;eb<3;eb++)
+  CGAL::Identity_transformation I;
+  std::vector<std::mutex> triangle_locks(T.size());
+  auto process_chunk = [&](size_t first, size_t last) -> void{
+    assert(last >= first);
+
+    for (size_t i=first; i<last; i++) {
+      const auto& box_pair = candidate_box_pairs[i];
+      const auto& a = box_pair.first;
+      const auto& b = box_pair.second;
+      Index fa = a.handle()-T.begin();
+      Index fb = b.handle()-T.begin();
+
+      std::lock_guard<std::mutex> guard_A(triangle_locks[fa]);
+      std::lock_guard<std::mutex> guard_B(triangle_locks[fb]);
+      const Triangle_3& A = *a.handle();
+      const Triangle_3& B = *b.handle();
+
+      // Number of combinatorially shared vertices
+      Index comb_shared_vertices = 0;
+      // Number of geometrically shared vertices (*not* including combinatorially
+      // shared)
+      Index geo_shared_vertices = 0;
+      // Keep track of shared vertex indices
+      std::vector<std::pair<Index,Index> > shared;
+      Index ea,eb;
+      for(ea=0;ea<3;ea++)
       {
-        if(F(fa,ea) == F(fb,eb))
+        for(eb=0;eb<3;eb++)
         {
-          comb_shared_vertices++;
-          shared.emplace_back(ea,eb);
-        }else if(A.vertex(ea) == B.vertex(eb))
-        {
-          geo_shared_vertices++;
-          shared.emplace_back(ea,eb);
+          if(F(fa,ea) == F(fb,eb))
+          {
+            comb_shared_vertices++;
+            shared.emplace_back(ea,eb);
+          }else if(A.vertex(ea) == B.vertex(eb))
+          {
+            geo_shared_vertices++;
+            shared.emplace_back(ea,eb);
+          }
         }
       }
-    }
-    const Index total_shared_vertices = comb_shared_vertices + geo_shared_vertices;
+      const Index total_shared_vertices = comb_shared_vertices + geo_shared_vertices;
 
-    if(comb_shared_vertices== 3)
-    {
-      assert(shared.size() == 3);
-      //// Combinatorially duplicate face, these should be removed by preprocessing
-      //cerr<<REDRUM("Facets "<<fa<<" and "<<fb<<" are combinatorial duplicates")<<endl;
-      continue;
-    }
-    if(total_shared_vertices== 3)
-    {
-      assert(shared.size() == 3);
-      //// Geometrically duplicate face, these should be removed by preprocessing
-      //cerr<<REDRUM("Facets "<<fa<<" and "<<fb<<" are geometrical duplicates")<<endl;
-      continue;
-    }
-    if(total_shared_vertices == 2)
-    {
-      assert(shared.size() == 2);
-      // Q: What about coplanar?
-      //
-      // o    o
-      // |\  /|
-      // | \/ |
-      // | /\ |
-      // |/  \|
-      // o----o
-      double_shared_vertex(A,B,fa,fb,shared);
-      continue;
-    }
-    assert(total_shared_vertices<=1);
-    if(total_shared_vertices==1)
-    {
-      single_shared_vertex(A,B,fa,fb,shared[0].first,shared[0].second);
-    }else
-    {
-      intersect(*a.handle(),*b.handle(),fa,fb);
+      if(comb_shared_vertices== 3)
+      {
+        assert(shared.size() == 3);
+        //// Combinatorially duplicate face, these should be removed by preprocessing
+        //cerr<<REDRUM("Facets "<<fa<<" and "<<fb<<" are combinatorial duplicates")<<endl;
+        continue;
+      }
+      if(total_shared_vertices== 3)
+      {
+        assert(shared.size() == 3);
+        //// Geometrically duplicate face, these should be removed by preprocessing
+        //cerr<<REDRUM("Facets "<<fa<<" and "<<fb<<" are geometrical duplicates")<<endl;
+        continue;
+      }
+      if(total_shared_vertices == 2)
+      {
+        assert(shared.size() == 2);
+        // Q: What about coplanar?
+        //
+        // o    o
+        // |\  /|
+        // | \/ |
+        // | /\ |
+        // |/  \|
+        // o----o
+        double_shared_vertex(A,B,fa,fb,shared);
+        continue;
+      }
+      assert(total_shared_vertices<=1);
+      if(total_shared_vertices==1)
+      {
+        single_shared_vertex(A,B,fa,fb,shared[0].first,shared[0].second);
+      }else
+      {
+        intersect(*a.handle(),*b.handle(),fa,fb);
+      }
     }
+  };
+  const size_t num_threads = std::thread::hardware_concurrency();
+  const size_t num_pairs = candidate_box_pairs.size();
+  const size_t chunk_size = num_pairs / num_threads;
+  std::vector<std::thread> threads;
+  for (size_t i=0; i<num_threads-1; i++) {
+    threads.emplace_back(process_chunk, i*chunk_size, (i+1)*chunk_size);
+  }
+  // Do some work in the master thread.
+  process_chunk((num_threads-1)*chunk_size, num_pairs);
+  for (auto& t : threads) {
+    if (t.joinable()) t.join();
   }
+  //process_chunk(0, candidate_box_pairs.size());
 }
 
 #endif

include/igl/copyleft/cgal/extract_cells.cpp

@@ -19,7 +19,7 @@
 #include <vector>
 #include <queue>
 
-//#define EXTRACT_CELLS_DEBUG
+#define EXTRACT_CELLS_DEBUG
 
 template<
 typename DerivedV,

include/igl/copyleft/cgal/remesh_intersections.cpp

@@ -8,6 +8,7 @@
 //
 #include "remesh_intersections.h"
 #include "assign_scalar.h"
+#include "../../get_seconds.h"
 
 #include <vector>
 #include <map>
@@ -15,6 +16,8 @@
 #include <unordered_map>
 #include <iostream>
 
+#define REMESH_INTERSECTIONS_TIMING
+
 template <
 typename DerivedV,
 typename DerivedF,
@@ -39,6 +42,21 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
         Eigen::PlainObjectBase<DerivedJ> & J,
         Eigen::PlainObjectBase<DerivedIM> & IM) {
 
+#ifdef REMESH_INTERSECTIONS_TIMING
+    const auto & tictoc = []() -> double
+    {
+      static double t_start = igl::get_seconds();
+      double diff = igl::get_seconds()-t_start;
+      t_start += diff;
+      return diff;
+    };
+    const auto log_time = [&](const std::string& label) -> void {
+      std::cout << "remesh_intersections." << label << ": "
+          << tictoc() << std::endl;
+    };
+    tictoc();
+#endif
+
     typedef CGAL::Point_3<Kernel>    Point_3;
     typedef CGAL::Segment_3<Kernel>  Segment_3; 
     typedef CGAL::Triangle_3<Kernel> Triangle_3; 
@@ -60,39 +78,6 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
     };
     typedef std::unordered_map<Edge, std::vector<Index>, EdgeHash > EdgeMap;
 
-    auto normalize_plane_coeff = [](const Plane_3& P) ->
-    std::vector<typename Kernel::FT> {
-        std::vector<typename Kernel::FT> coeffs = {
-            P.a(), P.b(), P.c(), P.d()
-        };
-        const auto max_itr = std::max_element(coeffs.begin(), coeffs.end());
-        const auto min_itr = std::min_element(coeffs.begin(), coeffs.end());
-        typename Kernel::FT max_coeff;
-        if (*max_itr < -1 * *min_itr) {
-            max_coeff = *min_itr;
-        } else {
-            max_coeff = *max_itr;
-        }
-        std::transform(coeffs.begin(), coeffs.end(), coeffs.begin(),
-                [&](const typename Kernel::FT& val)
-                {return val / max_coeff; } );
-        return coeffs;
-    };
-
-    auto plane_comp = [&](const Plane_3& p1, const Plane_3& p2) -> bool{
-        const auto p1_coeffs = normalize_plane_coeff(p1);
-        const auto p2_coeffs = normalize_plane_coeff(p2);
-        if (p1_coeffs[0] != p2_coeffs[0])
-            return p1_coeffs[0] < p2_coeffs[0];
-        if (p1_coeffs[1] != p2_coeffs[1])
-            return p1_coeffs[1] < p2_coeffs[1];
-        if (p1_coeffs[2] != p2_coeffs[2])
-            return p1_coeffs[2] < p2_coeffs[2];
-        if (p1_coeffs[3] != p2_coeffs[3])
-            return p1_coeffs[3] < p2_coeffs[3];
-        return false;
-    };
-
     const size_t num_faces = F.rows();
     const size_t num_base_vertices = V.rows();
     assert(num_faces == T.size());
@@ -121,6 +106,9 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
             }
         }
     }
+#ifdef REMESH_INTERSECTIONS_TIMING
+    log_time("coplanar_analysis");
+#endif
 
     std::vector<std::vector<Index> > resolved_faces;
     std::vector<Index> source_faces;
@@ -130,7 +118,7 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
     /**
      * Run constraint Delaunay triangulation on the plane.
      */
-    auto run_delaunay_triangulation = [&](const Plane_3& P,
+    auto run_delaunay_triangulation = [&offending, &T](const Plane_3& P,
             const std::vector<Index>& involved_faces,
             std::vector<Point_3>& vertices,
             std::vector<std::vector<Index> >& faces) -> void {
@@ -288,9 +276,13 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
             source_faces.push_back(i);
         }
     }
+#ifdef REMESH_INTERSECTIONS_TIMING
+    log_time("copy_untouched_faces");
+#endif
 
     // Process self-intersecting faces.
     std::vector<bool> processed(num_faces, false);
+    std::vector<std::pair<Plane_3, std::vector<Index> > > cdt_inputs;
     for (const auto itr : offending) {
         const auto fid = itr.first;
         if (processed[fid]) continue;
@@ -321,11 +313,40 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
         }
 
         Plane_3 P = T[fid].supporting_plane();
-        std::vector<Point_3> vertices;
-        std::vector<std::vector<Index> > faces;
+        cdt_inputs.emplace_back(P, involved_faces);
+    }
+#ifdef REMESH_INTERSECTIONS_TIMING
+    log_time("prepare_cdt_input");
+#endif
+
+    const size_t num_cdts = cdt_inputs.size();
+    std::vector<std::vector<Point_3> > cdt_vertices(num_cdts);
+    std::vector<std::vector<std::vector<Index> > > cdt_faces(num_cdts);
+
+    const auto run_cdt = [&](const size_t first, const size_t last) {
+      for (size_t i=first; i<last; i++) {
+        auto& vertices = cdt_vertices[i];
+        auto& faces = cdt_faces[i];
+        const auto& P = cdt_inputs[i].first;
+        const auto& involved_faces = cdt_inputs[i].second;
         run_delaunay_triangulation(P, involved_faces, vertices, faces);
+      }
+    };
+    run_cdt(0, num_cdts);
+#ifdef REMESH_INTERSECTIONS_TIMING
+    log_time("cdt");
+#endif
+
+    // Post process
+    for (size_t i=0; i<num_cdts; i++) {
+        const auto& vertices = cdt_vertices[i];
+        const auto& faces = cdt_faces[i];
+        const auto& involved_faces = cdt_inputs[i].second;
         post_triangulation_process(vertices, faces, involved_faces);
     }
+#ifdef REMESH_INTERSECTIONS_TIMING
+    log_time("post_process");
+#endif
 
     // Output resolved mesh.
     const size_t num_out_vertices = new_vertices.size() + num_base_vertices;
@@ -352,12 +373,25 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
     J.resize(num_out_faces);
     std::copy(source_faces.begin(), source_faces.end(), J.data());
 
+    struct PointHash {
+        size_t operator()(const Point_3& p) const {
+            static const std::hash<double> hasher{};
+            double x,y,z;
+            assign_scalar(p.x(), x);
+            assign_scalar(p.y(), y);
+            assign_scalar(p.z(), z);
+            return hasher(x) ^ hasher(y) ^ hasher(z);
+        }
+    };
+
     // Extract unique vertex indices.
-    IM.resize(VV.rows(),1);
-    std::map<Point_3,Index> vv2i;
+    const size_t VV_size = VV.rows();
+    IM.resize(VV_size,1);
+    std::unordered_map<Point_3, Index, PointHash> vv2i;
+    vv2i.reserve(VV_size);
     // 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<VV.rows();v++) {
+    for(Index v = 0;v<VV_size;v++) {
         typename Kernel::FT p0,p1,p2;
         assign_scalar(VV(v,0),p0);
         assign_scalar(VV(v,1),p1);
@@ -369,6 +403,9 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
         assert(vv2i.count(p) == 1);
         IM(v) = vv2i[p];
     }
+#ifdef REMESH_INTERSECTIONS_TIMING
+    log_time("store_results");
+#endif
 }
 
 #ifdef IGL_STATIC_LIBRARY

include/igl/copyleft/cgal/remesh_self_intersections.cpp

@@ -56,7 +56,7 @@ IGL_INLINE void igl::copyleft::cgal::remesh_self_intersections(
         DerivedJ,
         DerivedIM>
       SelfIntersectMeshK;
-    SelfIntersectMeshK SIM = SelfIntersectMeshK(V,F,params,VV,FF,IF,J,IM);
+    SelfIntersectMeshK SIM(V,F,params,VV,FF,IF,J,IM);
 
 //#ifdef __APPLE__
 //    signal(SIGFPE,SIG_DFL);
@@ -76,7 +76,7 @@ IGL_INLINE void igl::copyleft::cgal::remesh_self_intersections(
         DerivedJ,
         DerivedIM>
       SelfIntersectMeshK;
-    SelfIntersectMeshK SIM = SelfIntersectMeshK(V,F,params,VV,FF,IF,J,IM);
+    SelfIntersectMeshK SIM(V,F,params,VV,FF,IF,J,IM);
   }
 }