Merge pull request #193 from qnzhou/master

Bug fix

Former-commit-id: fa19da132540a279cf21b215045ff7564b3046a3

include/igl/copyleft/boolean/mesh_boolean.cpp

@@ -16,11 +16,13 @@
 #include "../../unique_simplices.h"
 #include "../../slice.h"
 #include "../../resolve_duplicated_faces.h"
+#include "../../get_seconds.h"
 
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 #include <algorithm>
 
 //#define MESH_BOOLEAN_TIMING
+//#define DOUBLE_CHECK_EXACT_OUTPUT
 
 template <
   typename DerivedVA,
@@ -182,6 +184,30 @@ IGL_INLINE void igl::copyleft::boolean::mesh_boolean(
     igl::resolve_duplicated_faces(kept_faces, G, JJ);
     igl::slice(kept_face_indices, JJ, 1, J);
 
+#ifdef DOUBLE_CHECK_EXACT_OUTPUT
+    {
+      // Sanity check on exact output.
+      igl::copyleft::cgal::RemeshSelfIntersectionsParam params;
+      params.detect_only = true;
+      params.first_only = true;
+      MatrixXES dummy_VV;
+      DerivedFC dummy_FF, dummy_IF;
+      Eigen::VectorXi dummy_J, dummy_IM;
+      igl::copyleft::cgal::SelfIntersectMesh<
+        Kernel,
+        MatrixXES, DerivedFC,
+        MatrixXES, DerivedFC,
+        DerivedFC,
+        Eigen::VectorXi,
+        Eigen::VectorXi
+      > checker(V, G, params,
+          dummy_VV, dummy_FF, dummy_IF, dummy_J, dummy_IM);
+      if (checker.count != 0) {
+        throw "Self-intersection not fully resolved.";
+      }
+    }
+#endif
+
     MatrixX3S Vs(V.rows(), V.cols());
     for (size_t i=0; i<(size_t)V.rows(); i++)
     {

include/igl/copyleft/cgal/SelfIntersectMesh.h

@@ -10,6 +10,7 @@
 
 #include "CGAL_includes.hpp"
 #include "RemeshSelfIntersectionsParam.h"
+#include "../../unique.h"
 
 #include <Eigen/Dense>
 #include <list>
@@ -423,7 +424,7 @@ inline igl::copyleft::cgal::SelfIntersectMesh<
     return;
   }
 
-  remesh_intersections(V,F,T,offending,edge2faces,VV,FF,J,IM);
+  remesh_intersections(V,F,T,offending,edge2faces,true,VV,FF,J,IM);
 
 #ifdef IGL_SELFINTERSECTMESH_DEBUG
   log_time("remesh_intersection");
@@ -702,8 +703,6 @@ 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() &&
@@ -711,7 +710,6 @@ 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
@@ -865,22 +863,46 @@ inline void igl::copyleft::cgal::SelfIntersectMesh<
   DerivedJ,
   DerivedIM>::process_intersecting_boxes()
 {
-  CGAL::Identity_transformation I;
   std::vector<std::mutex> triangle_locks(T.size());
-  auto process_chunk = [&](size_t first, size_t last) -> void{
+  std::vector<std::mutex> vertex_locks(V.rows());
+  std::mutex index_lock;
+  auto process_chunk = [&](const size_t first, const 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();
+      Index fa=T.size(), fb=T.size();
+      {
+        // Before knowing which triangles are involved, we need to lock
+        // everything to prevent race condition in updating reference counters.
+        std::lock_guard<std::mutex> guard(index_lock);
+        const auto& box_pair = candidate_box_pairs[i];
+        const auto& a = box_pair.first;
+        const auto& b = box_pair.second;
+        fa = a.handle()-T.begin();
+        fb = b.handle()-T.begin();
+      }
+      assert(fa < T.size());
+      assert(fb < T.size());
 
+      // Lock triangles
       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();
+
+      // Lock vertices
+      std::list<std::lock_guard<std::mutex> > guard_vertices;
+      {
+        std::vector<typename DerivedF::Scalar> unique_vertices;
+        std::vector<size_t> tmp1, tmp2;
+        igl::unique({F(fa,0), F(fa,1), F(fa,2), F(fb,0), F(fb,1), F(fb,2)},
+            unique_vertices, tmp1, tmp2);
+        std::for_each(unique_vertices.begin(), unique_vertices.end(),
+            [&](const typename DerivedF::Scalar& vi) {
+            guard_vertices.emplace_back(vertex_locks[vi]);
+            });
+      }
+
+      const Triangle_3& A = T[fa];
+      const Triangle_3& B = T[fb];
 
       // Number of combinatorially shared vertices
       Index comb_shared_vertices = 0;
@@ -941,7 +963,7 @@ inline void igl::copyleft::cgal::SelfIntersectMesh<
         single_shared_vertex(A,B,fa,fb,shared[0].first,shared[0].second);
       }else
       {
-        intersect(*a.handle(),*b.handle(),fa,fb);
+        intersect(A,B,fa,fb);
       }
     }
   };

include/igl/copyleft/cgal/extract_cells.cpp

@@ -16,6 +16,7 @@
 #include "order_facets_around_edge.h"
 #include "outer_facet.h"
 
+#include <iostream>
 #include <vector>
 #include <queue>
 
@@ -366,186 +367,107 @@ IGL_INLINE size_t igl::copyleft::cgal::extract_cells_single_component(
   const size_t num_unique_edges = uE.rows();
   const size_t num_patches = P.maxCoeff() + 1;
 
-  // patch_edge_adj[p] --> list {e,f,g,...} such that p is a patch id and
-  //   e,f,g etc. are edge indices into 
-  std::vector<std::vector<size_t> > patch_edge_adj(num_patches);
-  // orders[u] --> J  where u is an index into unique edges uE and J is a
-  //   #adjacent-faces list of face-edge indices into F*3 sorted cyclicaly around
-  //   edge u.
-  std::vector<Eigen::VectorXi> orders(num_unique_edges);
-  // orientations[u] ---> list {f1,f2, ...} where u is an index into unique edges uE
-  //   and points to #adj-facets long list of flags whether faces are oriented
-  //   to point their normals clockwise around edge when looking along the
-  //   edge.
-  std::vector<std::vector<bool> > orientations(num_unique_edges);
-  // Loop over unique edges
-  for (size_t i=0; i<num_unique_edges; i++) 
-  {
+  // Build patch-patch adjacency list.
+  std::vector<std::map<size_t, size_t> > patch_adj(num_patches);
+  for (size_t i=0; i<num_unique_edges; i++) {
     const size_t s = uE(i,0);
     const size_t d = uE(i,1);
     const auto adj_faces = uE2E[i];
-    // If non-manifold (more than two incident faces)
-    if (adj_faces.size() > 2) 
-    {
-      // signed_adj_faces[a] --> sid  where a is an index into adj_faces
-      // (list of face edges on {s,d}) and sid is a signed index for resolve
-      // co-planar duplicates consistently (i.e. using simulation of
-      // simplicity).
+    const size_t num_adj_faces = adj_faces.size();
+    if (num_adj_faces > 2) {
+      for (size_t j=0; j<num_adj_faces; j++) {
+        const size_t patch_j = P[edge_index_to_face_index(adj_faces[j])];
+        for (size_t k=j+1; k<num_adj_faces; k++) {
+          const size_t patch_k = P[edge_index_to_face_index(adj_faces[k])];
+          if (patch_adj[patch_j].find(patch_k) == patch_adj[patch_j].end()) {
+            patch_adj[patch_j].insert({patch_k, i});
+          }
+          if (patch_adj[patch_k].find(patch_j) == patch_adj[patch_k].end()) {
+            patch_adj[patch_k].insert({patch_j, i});
+          }
+        }
+      }
+    }
+  }
+
+
+  const int INVALID = std::numeric_limits<int>::max();
+  std::vector<size_t> cell_labels(num_patches * 2);
+  for (size_t i=0; i<num_patches; i++) cell_labels[i] = i;
+  std::vector<std::set<size_t> > equivalent_cells(num_patches*2);
+  std::vector<bool> processed(num_unique_edges, false);
+
+  size_t label_count=0;
+  for (size_t i=0; i<num_patches; i++) {
+    for (const auto& entry : patch_adj[i]) {
+      const size_t neighbor_patch = entry.first;
+      const size_t uei = entry.second;
+      if (processed[uei]) continue;
+      processed[uei] = true;
+
+      const auto& adj_faces = uE2E[uei];
+      const size_t num_adj_faces = adj_faces.size();
+      assert(num_adj_faces > 2);
+
+      const size_t s = uE(uei,0);
+      const size_t d = uE(uei,1);
+
       std::vector<int> signed_adj_faces;
-      for (auto ei : adj_faces)
+      for (auto ej : adj_faces)
       {
-        const size_t fid = edge_index_to_face_index(ei);
+        const size_t fid = edge_index_to_face_index(ej);
         bool cons = is_consistent(fid, s, d);
         signed_adj_faces.push_back((fid+1)*(cons ? 1:-1));
       }
       {
         // Sort adjacent faces cyclically around {s,d}
-        auto& order = orders[i];
+        Eigen::VectorXi order;
         // order[f] will reveal the order of face f in signed_adj_faces
         order_facets_around_edge(V, F, s, d, signed_adj_faces, order);
-        // Determine if each facet is oriented to point its normal clockwise or
-        // not around the {s,d} (normals are not explicitly computed/used)
-        auto& orientation = orientations[i];
-        orientation.resize(order.size());
-        std::transform(
-          order.data(), 
-          order.data() + order.size(),
-          orientation.begin(), 
-          [&signed_adj_faces](const int i){ return signed_adj_faces[i] > 0;});
-        // re-index order from adjacent faces to full face list. Now order
-        // indexes F directly
-        std::transform(
-          order.data(), 
-          order.data() + order.size(),
-          order.data(), 
-          [&adj_faces](const int index){ return adj_faces[index];});
-      }
-      // loop over adjacent faces, remember that patch is adjacent to this edge
-      for(const auto & ei : adj_faces)
-      {
-        const size_t fid = edge_index_to_face_index(ei);
-        patch_edge_adj[P[fid]].push_back(ei);
+        for (size_t j=0; j<num_adj_faces; j++) {
+          const size_t curr_idx = j;
+          const size_t next_idx = (j+1)%num_adj_faces;
+          const size_t curr_patch_idx =
+            P[edge_index_to_face_index(adj_faces[order[curr_idx]])];
+          const size_t next_patch_idx =
+            P[edge_index_to_face_index(adj_faces[order[next_idx]])];
+          const bool curr_cons = signed_adj_faces[order[curr_idx]] > 0;
+          const bool next_cons = signed_adj_faces[order[next_idx]] > 0;
+          const size_t curr_cell_idx = curr_patch_idx*2 + (curr_cons?0:1);
+          const size_t next_cell_idx = next_patch_idx*2 + (next_cons?1:0);
+          equivalent_cells[curr_cell_idx].insert(next_cell_idx);
+          equivalent_cells[next_cell_idx].insert(curr_cell_idx);
+        }
       }
     }
   }
 
-  // Initialize all patch-to-cell indices as "invalid" (unlabeled)
-  const int INVALID = std::numeric_limits<int>::max();
+  size_t count=0;
   cells.resize(num_patches, 2);
   cells.setConstant(INVALID);
-
-  // Given a "seed" patch id, a cell id, and which side of the patch that cell
-  // lies, identify all other patches bounding this cell (and tell them that
-  // they do)
-  //
-  // Inputs:
-  //   seed_patch_id  index into patches
-  //   cell_idx  index into cells
-  //   seed_patch_side   0 or 1 depending on whether cell_idx is on left or
-  //     right side of seed_patch_id 
-  //   cells  #cells by 2 list of current assignment of cells incident on each
-  //     side of a patch
-  // Outputs:
-  //   cells  udpated (see input)
-  // 
-  const auto & peel_cell_bd = 
-    [&P,&patch_edge_adj,&EMAP,&orders,&orientations,&num_faces](
-    const size_t seed_patch_id, 
-    const short seed_patch_side, 
-    const size_t cell_idx,
-    Eigen::PlainObjectBase<DerivedC>& cells)
-  {
-    typedef std::pair<size_t, short> PatchSide;
-    // Initialize a queue of {p,side} patch id and patch side pairs to BFS over
-    // all patches
-    std::queue<PatchSide> Q;
-    Q.emplace(seed_patch_id, seed_patch_side);
-    // assign cell id of seed patch on appropriate side
-    cells(seed_patch_id, seed_patch_side) = cell_idx;
-    while (!Q.empty())
-    {
-      // Pop patch from Q
-      const auto entry = Q.front();
+  const auto extract_equivalent_cells = [&](size_t i) {
+    if (cells(i/2, i%2) != INVALID) return;
+    std::queue<size_t> Q;
+    Q.push(i);
+    cells(i/2, i%2) = count;
+    while (!Q.empty()) {
+      const size_t index = Q.front();
       Q.pop();
-      const size_t patch_id = entry.first;
-      const short side = entry.second;
-      // face-edges adjacent to patch
-      const auto& adj_edges = patch_edge_adj[patch_id];
-      // Loop over **ALL EDGES IN THE ENTIRE PATCH** not even just the boundary
-      // edges, all edges... O(n)
-      for (const auto& ei : adj_edges)
-      {
-        // unique edge
-        const size_t uei = EMAP[ei];
-        // ordering of face-edges stored at edge
-        const auto& order = orders[uei];
-        // consistent orientation flags at face-edges stored at edge
-        const auto& orientation = orientations[uei];
-        const size_t edge_valance = order.size();
-        // Search for ei's (i.e. patch_id's) place in the ordering: O(#patches)
-        size_t curr_i = 0;
-        for (curr_i=0; curr_i < edge_valance; curr_i++)
-        {
-          if ((size_t)order[curr_i] == ei) break;
-        }
-        assert(curr_i < edge_valance && "Failed to find edge.");
-        // is the starting face consistent?
-        const bool cons = orientation[curr_i];
-        // Look clockwise or counter-clockwise for the next face, depending on
-        // whether looking to left or right side and whether consistently
-        // oriented or not.
-        size_t next;
-        if (side == 0)
-        {
-          next = (cons ? (curr_i + 1) :
-          (curr_i + edge_valance - 1)) % edge_valance;
-        } else {
-          next = (cons ? (curr_i+edge_valance-1) : (curr_i+1))%edge_valance;
-        }
-        // Determine face-edge index of next face
-        const size_t next_ei = order[next];
-        // Determine whether next is consistently oriented
-        const bool next_cons = orientation[next];
-        // Determine patch of next
-        const size_t next_patch_id = P[next_ei % num_faces];
-        // Determine which side of patch cell_idx is on, based on whether the
-        // consistency of next matches the consistency of this patch and which
-        // side of this patch we're on.
-        const short next_patch_side = (next_cons != cons) ?  side:abs(side-1);
-        // If that side of next's patch is not labeled, then label and add to
-        // queue
-        if (cells(next_patch_id, next_patch_side) == INVALID) 
-        {
-          Q.emplace(next_patch_id, next_patch_side);
-          cells(next_patch_id, next_patch_side) = cell_idx;
-        }else 
-        {
-          assert(
-            (size_t)cells(next_patch_id, next_patch_side) == cell_idx && 
-            "Encountered cell assignment inconsistency");
+      for (const auto j : equivalent_cells[index]) {
+        if (cells(j/2, j%2) == INVALID) {
+          cells(j/2, j%2) = count;
+          Q.push(j);
         }
       }
     }
+    count++;
   };
-
-  size_t count=0;
-  // Loop over all patches
-  for (size_t i=0; i<num_patches; i++)
-  {
-    // If left side of patch is still unlabeled, create a new cell and find all
-    // patches also on its boundary
-    if (cells(i, 0) == INVALID)
-    {
-      peel_cell_bd(i, 0, count,cells);
-      count++;
-    }
-    // Likewise for right side
-    if (cells(i, 1) == INVALID)
-    {
-      peel_cell_bd(i, 1, count,cells);
-      count++;
-    }
+  for (size_t i=0; i<num_patches; i++) {
+    extract_equivalent_cells(i*2);
+    extract_equivalent_cells(i*2+1);
   }
+
+  assert((cells.array() != INVALID).all());
   return count;
 }
 

include/igl/copyleft/cgal/remesh_intersections.cpp

@@ -19,6 +19,33 @@
 
 //#define REMESH_INTERSECTIONS_TIMING
 
+template <
+typename DerivedV,
+typename DerivedF,
+typename Kernel,
+typename DerivedVV,
+typename DerivedFF,
+typename DerivedJ,
+typename DerivedIM>
+IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::PlainObjectBase<DerivedF> & F,
+        const std::vector<CGAL::Triangle_3<Kernel> > & T,
+        const std::map<
+        typename DerivedF::Index,
+        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,
+        Eigen::PlainObjectBase<DerivedVV> & VV,
+        Eigen::PlainObjectBase<DerivedFF> & FF,
+        Eigen::PlainObjectBase<DerivedJ> & J,
+        Eigen::PlainObjectBase<DerivedIM> & IM) {
+    igl::copyleft::cgal::remesh_intersections(
+            V, F, T, offending, edge2faces, false, VV, FF, J, IM);
+}
+
 template <
 typename DerivedV,
 typename DerivedF,
@@ -38,6 +65,7 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
         const std::map<
         std::pair<typename DerivedF::Index,typename DerivedF::Index>,
         std::vector<typename DerivedF::Index> > & /*edge2faces*/,
+        bool stitch_all,
         Eigen::PlainObjectBase<DerivedVV> & VV,
         Eigen::PlainObjectBase<DerivedFF> & FF,
         Eigen::PlainObjectBase<DerivedJ> & J,
@@ -187,51 +215,61 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
      * Given p on triangle indexed by ori_f, determine the index of p.
      */
     auto determine_point_index = [&](
-            const Point_3& p, const size_t ori_f) -> Index {
+        const Point_3& p, const size_t ori_f) -> Index {
+      if (stitch_all) {
+        // No need to check if p shared by multiple triangles because all shared
+        // vertices would be merged later on.
+        const size_t index = num_base_vertices + new_vertices.size();
+        new_vertices.push_back(p);
+        return index;
+      } else {
+        // Stitching triangles according to input connectivity.
+        // This step is potentially costly.
         const auto& triangle = T[ori_f];
         const auto& f = F.row(ori_f).eval();
 
         // Check if p is one of the triangle corners.
         for (size_t i=0; i<3; i++) {
-            if (p == triangle[i]) return f[i];
+          if (p == triangle[i]) return f[i];
         }
 
         // Check if p is on one of the edges.
         for (size_t i=0; i<3; i++) {
-            const Point_3 curr_corner = triangle[i];
-            const Point_3 next_corner = triangle[(i+1)%3];
-            const Segment_3 edge(curr_corner, next_corner);
-            if (edge.has_on(p)) {
-                const Index curr = f[i];
-                const Index next = f[(i+1)%3];
-                Edge key;
-                key.first = curr<next?curr:next;
-                key.second = curr<next?next:curr;
-                auto itr = edge_vertices.find(key);
-                if (itr == edge_vertices.end()) {
-                    const Index index =
-                        num_base_vertices + new_vertices.size();
-                    edge_vertices.insert({key, {index}});
-                    new_vertices.push_back(p);
-                    return index;
-                } else {
-                    for (const auto vid : itr->second) {
-                        if (p == new_vertices[vid - num_base_vertices]) {
-                            return vid;
-                        }
-                    }
-                    const size_t index = num_base_vertices + new_vertices.size();
-                    new_vertices.push_back(p);
-                    itr->second.push_back(index);
-                    return index;
+          const Point_3 curr_corner = triangle[i];
+          const Point_3 next_corner = triangle[(i+1)%3];
+          const Segment_3 edge(curr_corner, next_corner);
+          if (edge.has_on(p)) {
+            const Index curr = f[i];
+            const Index next = f[(i+1)%3];
+            Edge key;
+            key.first = curr<next?curr:next;
+            key.second = curr<next?next:curr;
+            auto itr = edge_vertices.find(key);
+            if (itr == edge_vertices.end()) {
+              const Index index =
+                num_base_vertices + new_vertices.size();
+              edge_vertices.insert({key, {index}});
+              new_vertices.push_back(p);
+              return index;
+            } else {
+              for (const auto vid : itr->second) {
+                if (p == new_vertices[vid - num_base_vertices]) {
+                  return vid;
                 }
+              }
+              const size_t index = num_base_vertices + new_vertices.size();
+              new_vertices.push_back(p);
+              itr->second.push_back(index);
+              return index;
             }
+          }
         }
 
         // p must be in the middle of the triangle.
         const size_t index = num_base_vertices + new_vertices.size();
         new_vertices.push_back(p);
         return index;
+      }
     };
 
     /**
@@ -241,32 +279,45 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
             const std::vector<Point_3>& vertices,
             const std::vector<std::vector<Index> >& faces,
             const std::vector<Index>& involved_faces) -> void {
-        for (const auto& f : faces) {
-            const Point_3& v0 = vertices[f[0]];
-            const Point_3& v1 = vertices[f[1]];
-            const Point_3& v2 = vertices[f[2]];
-            Point_3 center(
-                    (v0[0] + v1[0] + v2[0]) / 3.0,
-                    (v0[1] + v1[1] + v2[1]) / 3.0,
-                    (v0[2] + v1[2] + v2[2]) / 3.0);
-            for (const auto& ori_f : involved_faces) {
-                const auto& triangle = T[ori_f];
-                const Plane_3 P = triangle.supporting_plane();
-                if (triangle.has_on(center)) {
-                    std::vector<Index> corners(3);
-                    corners[0] = determine_point_index(v0, ori_f);
-                    corners[1] = determine_point_index(v1, ori_f);
-                    corners[2] = determine_point_index(v2, ori_f);
-                    if (CGAL::orientation(
-                                P.to_2d(v0), P.to_2d(v1), P.to_2d(v2))
-                            == CGAL::RIGHT_TURN) {
-                        std::swap(corners[0], corners[1]);
-                    }
-                    resolved_faces.emplace_back(corners);
-                    source_faces.push_back(ori_f);
-                }
+      assert(involved_faces.size() > 0);
+      for (const auto& f : faces) {
+        const Point_3& v0 = vertices[f[0]];
+        const Point_3& v1 = vertices[f[1]];
+        const Point_3& v2 = vertices[f[2]];
+        Point_3 center(
+            (v0[0] + v1[0] + v2[0]) / 3.0,
+            (v0[1] + v1[1] + v2[1]) / 3.0,
+            (v0[2] + v1[2] + v2[2]) / 3.0);
+        if (involved_faces.size() == 1) {
+          // If only there is only one involved face, all sub-triangles must
+          // belong to it and have the correct orientation.
+          const auto& ori_f = involved_faces[0];
+          std::vector<Index> corners(3);
+          corners[0] = determine_point_index(v0, ori_f);
+          corners[1] = determine_point_index(v1, ori_f);
+          corners[2] = determine_point_index(v2, ori_f);
+          resolved_faces.emplace_back(corners);
+          source_faces.push_back(ori_f);
+        } else {
+          for (const auto& ori_f : involved_faces) {
+            const auto& triangle = T[ori_f];
+            const Plane_3 P = triangle.supporting_plane();
+            if (triangle.has_on(center)) {
+              std::vector<Index> corners(3);
+              corners[0] = determine_point_index(v0, ori_f);
+              corners[1] = determine_point_index(v1, ori_f);
+              corners[2] = determine_point_index(v2, ori_f);
+              if (CGAL::orientation(
+                    P.to_2d(v0), P.to_2d(v1), P.to_2d(v2))
+                  == CGAL::RIGHT_TURN) {
+                std::swap(corners[0], corners[1]);
+              }
+              resolved_faces.emplace_back(corners);
+              source_faces.push_back(ori_f);
             }
+          }
         }
+      }
     };
 
     // Process un-touched faces.
@@ -335,12 +386,11 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
     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);
+      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("stitching");
@@ -371,28 +421,27 @@ 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);
-        }
-    };
-
-    // TODO: use igl::unique()
     // Extract unique vertex indices.
-
     const size_t VV_size = VV.rows();
     IM.resize(VV_size,1);
 
     DerivedVV unique_vv;
     Eigen::VectorXi unique_to_vv, vv_to_unique;
     igl::unique_rows(VV, unique_vv, unique_to_vv, vv_to_unique);
-    for (Index v=0; v<VV_size; v++) {
-      IM(v) = unique_to_vv[vv_to_unique[v]];
+    if (!stitch_all) {
+      // Vertices with the same coordinates would be represented by one vertex.
+      // The IM value of an vertex is the index of the representative vertex.
+      for (Index v=0; v<VV_size; v++) {
+        IM(v) = unique_to_vv[vv_to_unique[v]];
+      }
+    } else {
+      // Screw IM and representative vertices.  Merge all vertices having the
+      // same coordinates into a single vertex and set IM to identity map.
+      VV = unique_vv;
+      std::transform(FF.data(), FF.data() + FF.rows()*FF.cols(),
+          FF.data(), [&vv_to_unique](const typename DerivedFF::Scalar& a)
+          { return vv_to_unique[a]; });
+      IM.setLinSpaced(unique_vv.rows(), 0, unique_vv.rows()-1);
     }
 #ifdef REMESH_INTERSECTIONS_TIMING
     log_time("store_results");

include/igl/copyleft/cgal/remesh_intersections.h

@@ -62,6 +62,35 @@ namespace igl
         Eigen::PlainObjectBase<DerivedFF> & FF,
         Eigen::PlainObjectBase<DerivedJ> & J,
         Eigen::PlainObjectBase<DerivedIM> & IM);
+
+      // Same as above except ``stitch_all`` flag:
+      //
+      // Input:
+      //   stitch_all: if true, merge all vertices with thte same coordiante.
+      template <
+        typename DerivedV,
+        typename DerivedF,
+        typename Kernel,
+        typename DerivedVV,
+        typename DerivedFF,
+        typename DerivedJ,
+        typename DerivedIM>
+      IGL_INLINE void remesh_intersections(
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::PlainObjectBase<DerivedF> & F,
+        const std::vector<CGAL::Triangle_3<Kernel> > & T,
+        const std::map<
+          typename DerivedF::Index,
+            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,
+        bool stitch_all,
+        Eigen::PlainObjectBase<DerivedVV> & VV,
+        Eigen::PlainObjectBase<DerivedFF> & FF,
+        Eigen::PlainObjectBase<DerivedJ> & J,
+        Eigen::PlainObjectBase<DerivedIM> & IM);
     }
   }
 }