Merge branch 'develop'

Former-commit-id: 48b2190caedea3cff7c0a97224b34d48a505c585

include/igl/boolean/mesh_boolean.cpp

@@ -211,7 +211,7 @@ IGL_INLINE void igl::boolean::per_face_winding_number_binary_operation(
     Eigen::VectorXi labels(num_faces);
     std::transform(CJ.data(), CJ.data()+CJ.size(), labels.data(),
             [&](int i) { return i<FA.rows() ? 0:1; });
-    igl::cgal::propagate_winding_numbers_beta(V, F, labels, W);
+    igl::cgal::propagate_winding_numbers(V, F, labels, W);
     assert(W.rows() == num_faces);
     if (W.cols() == 2) {
         assert(FB.rows() == 0);

include/igl/cgal/SelfIntersectMesh.h

@@ -225,7 +225,7 @@ namespace igl
 
 #include "mesh_to_cgal_triangle_list.h"
 #include "remesh_intersections.h"
-#include "remesh_intersections_beta.h"
+#include "remesh_intersections.h"
 
 #include <igl/REDRUM.h>
 #include <igl/get_seconds.h>
@@ -405,8 +405,7 @@ inline igl::cgal::SelfIntersectMesh<
     return;
   }
 
-  //remesh_intersections(V,F,T,offending,edge2faces,VV,FF,J,IM);
-  remesh_intersections_beta(V,F,T,offending,edge2faces,VV,FF,J,IM);
+  remesh_intersections(V,F,T,offending,edge2faces,VV,FF,J,IM);
 
   // Q: Does this give the same result as TETGEN?
   // A: For the cow and beast, yes.

include/igl/cgal/extract_cells.cpp

@@ -215,31 +215,33 @@ IGL_INLINE size_t igl::cgal::extract_cells(
     std::vector<std::vector<size_t> > nested_cells(num_raw_cells);
     std::vector<std::vector<size_t> > ambient_cells(num_raw_cells);
     std::vector<std::vector<size_t> > ambient_comps(num_components);
-    for (size_t i=0; i<num_components; i++) {
-        DerivedV queries(num_components-1, 3);
-        for (size_t j=0; j<num_components; j++) {
-            if (i == j) continue;
-            size_t index = j<i ? j:j-1;
-            queries.row(index) = get_triangle_center(outer_facets[j]);
-        }
+    if (num_components > 1) {
+        for (size_t i=0; i<num_components; i++) {
+            DerivedV queries(num_components-1, 3);
+            for (size_t j=0; j<num_components; j++) {
+                if (i == j) continue;
+                size_t index = j<i ? j:j-1;
+                queries.row(index) = get_triangle_center(outer_facets[j]);
+            }
 
-        const auto& I = Is[i];
-        Eigen::VectorXi closest_facets, closest_facet_orientations;
-        igl::cgal::closest_facet(V, F, I, queries,
-                closest_facets, closest_facet_orientations);
+            const auto& I = Is[i];
+            Eigen::VectorXi closest_facets, closest_facet_orientations;
+            igl::cgal::closest_facet(V, F, I, queries,
+                    closest_facets, closest_facet_orientations);
 
-        for (size_t j=0; j<num_components; j++) {
-            if (i == j) continue;
-            size_t index = j<i ? j:j-1;
-            const size_t closest_patch = P[closest_facets[index]];
-            const size_t closest_patch_side = closest_facet_orientations[index]
-                ? 0:1;
-            const size_t ambient_cell = raw_cells(closest_patch,
-                    closest_patch_side);
-            if (ambient_cell != outer_cells[i]) {
-                nested_cells[ambient_cell].push_back(outer_cells[j]);
-                ambient_cells[outer_cells[j]].push_back(ambient_cell);
-                ambient_comps[j].push_back(i);
+            for (size_t j=0; j<num_components; j++) {
+                if (i == j) continue;
+                size_t index = j<i ? j:j-1;
+                const size_t closest_patch = P[closest_facets[index]];
+                const size_t closest_patch_side = closest_facet_orientations[index]
+                    ? 0:1;
+                const size_t ambient_cell = raw_cells(closest_patch,
+                        closest_patch_side);
+                if (ambient_cell != outer_cells[i]) {
+                    nested_cells[ambient_cell].push_back(outer_cells[j]);
+                    ambient_cells[outer_cells[j]].push_back(ambient_cell);
+                    ambient_comps[j].push_back(i);
+                }
             }
         }
     }

include/igl/cgal/extract_cells.h

@@ -15,23 +15,40 @@
 
 namespace igl {
     namespace cgal {
+        // Extract connected 3D space partitioned by mesh (V, F).
+        //
+        // Inputs:
+        //   V  #V by 3 array of vertices.
+        //   F  #F by 3 array of faces.
+        //
+        // Output:
+        //   cells  #F by 2 array of cell indices.  cells(i,0) represents the
+        //          cell index on the positive side of face i, and cells(i,1)
+        //          represents cell index of the negqtive side.
         template<
             typename DerivedV,
             typename DerivedF,
-            typename DerivedP,
-            typename DeriveduE,
-            typename uE2EType,
-            typename DerivedEMAP,
             typename DerivedC >
-        IGL_INLINE size_t extract_cells_single_component(
+        IGL_INLINE size_t extract_cells(
                 const Eigen::PlainObjectBase<DerivedV>& V,
                 const Eigen::PlainObjectBase<DerivedF>& F,
-                const Eigen::PlainObjectBase<DerivedP>& P,
-                const Eigen::PlainObjectBase<DeriveduE>& uE,
-                const std::vector<std::vector<uE2EType> >& uE2E,
-                const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
                 Eigen::PlainObjectBase<DerivedC>& cells);
 
+        // Extract connected 3D space partitioned by mesh (V, F).
+        //
+        // Inputs:
+        //   V  #V by 3 array of vertices.
+        //   F  #F by 3 array of faces.
+        //   P  #F list of patch indices.
+        //   E  #E by 2 array of vertex indices, one edge per row.
+        //  uE    #uE by 2 list of vertex_indices, represents undirected edges.
+        //  uE2E  #uE list of lists that maps uE to E. (a one-to-many map)
+        //  EMAP  #F*3 list of indices into uE.
+        //
+        // Output:
+        //   cells  #P by 2 array of cell indices.  cells(i,0) represents the
+        //          cell index on the positive side of patch i, and cells(i,1)
+        //          represents cell index of the negqtive side.
         template<
             typename DerivedV,
             typename DerivedF,
@@ -51,14 +68,40 @@ namespace igl {
                 const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
                 Eigen::PlainObjectBase<DerivedC>& cells);
 
+        // Extract connected 3D space partitioned by mesh (V, F).  Each
+        // connected component of the mesh partitions the ambient space
+        // separately.
+        //
+        // Inputs:
+        //   V  #V by 3 array of vertices.
+        //   F  #F by 3 array of faces.
+        //   P  #F list of patch indices.
+        //   E  #E by 2 array of vertex indices, one edge per row.
+        //  uE    #uE by 2 list of vertex_indices, represents undirected edges.
+        //  uE2E  #uE list of lists that maps uE to E. (a one-to-many map)
+        //  EMAP  #F*3 list of indices into uE.
+        //
+        // Output:
+        //   cells  #P by 2 array of cell indices.  cells(i,0) represents the
+        //          cell index on the positive side of patch i, and cells(i,1)
+        //          represents cell index of the negqtive side.
         template<
             typename DerivedV,
             typename DerivedF,
+            typename DerivedP,
+            typename DeriveduE,
+            typename uE2EType,
+            typename DerivedEMAP,
             typename DerivedC >
-        IGL_INLINE size_t extract_cells(
+        IGL_INLINE size_t extract_cells_single_component(
                 const Eigen::PlainObjectBase<DerivedV>& V,
                 const Eigen::PlainObjectBase<DerivedF>& F,
+                const Eigen::PlainObjectBase<DerivedP>& P,
+                const Eigen::PlainObjectBase<DeriveduE>& uE,
+                const std::vector<std::vector<uE2EType> >& uE2E,
+                const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
                 Eigen::PlainObjectBase<DerivedC>& cells);
+
     }
 }
 

include/igl/cgal/propagate_winding_numbers.cpp

@@ -26,43 +26,6 @@
 #include <queue>
 
 namespace propagate_winding_numbers_helper {
-    template<typename DerivedW >
-    bool winding_number_assignment_is_consistent(
-            const std::vector<Eigen::VectorXi>& orders,
-            const std::vector<std::vector<bool> >& orientations,
-            const Eigen::PlainObjectBase<DerivedW>& per_patch_winding_number) {
-        const size_t num_edge_curves = orders.size();
-        const size_t num_labels = per_patch_winding_number.cols() / 2;
-        for (size_t i=0; i<num_edge_curves; i++) {
-            const auto& order = orders[i];
-            const auto& orientation = orientations[i];
-            assert(order.size() == orientation.size());
-            const size_t order_size = order.size();
-            for (size_t j=0; j<order_size; j++) {
-                const size_t curr = j;
-                const size_t next = (j+1) % order_size;
-
-                for (size_t k=0; k<num_labels; k++) {
-                    // Retrieve the winding numbers of the current partition from
-                    // the current patch and the next patch.  If the patches forms
-                    // the boundary of a 3D volume, the winding number assignments
-                    // should be consistent.
-                    int curr_winding_number = per_patch_winding_number(
-                            order[curr], k*2+(orientation[curr]? 0:1));
-                    int next_winding_number = per_patch_winding_number(
-                            order[next], k*2+(orientation[next]? 1:0));
-                    if (curr_winding_number != next_winding_number) {
-                        std::cout << "edge: " << i << std::endl;
-                        std::cout << curr_winding_number << " != " <<
-                            next_winding_number << std::endl;
-                        return false;
-                    }
-                }
-            }
-        }
-        return true;
-    }
-
     template<
         typename DerivedF,
         typename DeriveduE,
@@ -103,583 +66,6 @@ namespace propagate_winding_numbers_helper {
     }
 }
 
-template<
-    typename DerivedV,
-    typename DerivedF,
-    typename DeriveduE,
-    typename uE2EType,
-    typename DerivedC,
-    typename DerivedP,
-    typename DerivedW >
-IGL_INLINE bool igl::cgal::propagate_winding_numbers_single_component_patch_wise(
-        const Eigen::PlainObjectBase<DerivedV>& V,
-        const Eigen::PlainObjectBase<DerivedF>& F,
-        const Eigen::PlainObjectBase<DeriveduE>& uE,
-        const std::vector<std::vector<uE2EType> >& uE2E,
-        const Eigen::PlainObjectBase<DerivedC>& labels,
-        const Eigen::PlainObjectBase<DerivedP>& P,
-        const std::vector<std::vector<size_t> >& intersection_curves,
-        Eigen::PlainObjectBase<DerivedW>& patch_W) {
-    const size_t num_faces = F.rows();
-    const size_t num_patches = P.maxCoeff() + 1;
-    assert(labels.size() == num_patches);
-    // Utility functions.
-    auto edge_index_to_face_index = [&](size_t ei) { return ei % num_faces; };
-    auto edge_index_to_corner_index = [&](size_t ei) { return ei / num_faces; };
-    auto face_and_corner_index_to_edge_index = [&](size_t fi, size_t ci) {
-        return ci*num_faces + fi;
-    };
-    auto get_edge_end_points = [&](size_t ei, size_t& s, size_t& d) {
-        const size_t fi = edge_index_to_face_index(ei);
-        const size_t ci = edge_index_to_corner_index(ei);
-        s = F(fi, (ci+1)%3);
-        d = F(fi, (ci+2)%3);
-    };
-    auto is_positively_orientated =
-        [&](size_t fi, size_t s, size_t d) -> bool{
-            const Eigen::Vector3i f = F.row(fi);
-            if (f[0] == d && f[1] == s) return true;
-            if (f[1] == d && f[2] == s) return true;
-            if (f[2] == d && f[0] == s) return true;
-            if (f[0] == s && f[1] == d) return false;
-            if (f[1] == s && f[2] == d) return false;
-            if (f[2] == s && f[0] == d) return false;
-            throw std::runtime_error("Edge does not belong to face.");
-            return false;
-        };
-
-    auto compute_signed_index = [&](size_t fi, size_t s, size_t d) -> int{
-        return int(fi+1) * (is_positively_orientated(fi, s, d) ? 1:-1);
-    };
-    auto compute_unsigned_index = [&](int signed_idx) -> size_t {
-        return abs(signed_idx) - 1;
-    };
-
-    // Order patches around each intersection curve.
-    const size_t num_edge_curves = intersection_curves.size();
-    std::vector<Eigen::VectorXi> orders(num_edge_curves);
-    std::vector<std::vector<bool> > orientations(num_edge_curves);
-    std::vector<std::vector<size_t> > patch_curve_adjacency(num_patches);
-    for (size_t i=0; i<num_edge_curves; i++) {
-        const auto& curve = intersection_curves[i];
-        const size_t uei = curve[0];
-        size_t s = uE(uei, 0);
-        size_t d = uE(uei, 1);
-        std::vector<int> adj_faces;
-        for (size_t ei : uE2E[uei]) {
-            const size_t fi = edge_index_to_face_index(ei);
-            const size_t signed_fi =
-                compute_signed_index(fi, s, d);
-            adj_faces.push_back(signed_fi);
-            patch_curve_adjacency[P[fi]].push_back(i);
-        }
-
-        auto& order = orders[i];
-        igl::cgal::order_facets_around_edge(
-                V, F, s, d, adj_faces, order);
-        assert(order.minCoeff() == 0);
-        assert(order.maxCoeff() == adj_faces.size() - 1);
-
-        auto& orientation = orientations[i];
-        orientation.resize(order.size());
-        std::transform(order.data(), order.data()+order.size(),
-                orientation.begin(), 
-                [&](size_t index) { return adj_faces[index] > 0; });
-        std::transform(order.data(), order.data()+order.size(),
-                order.data(), [&](size_t index) {
-                return P[compute_unsigned_index(adj_faces[index])];
-                } );
-    }
-
-    // Propagate winding number from infinity.
-    // Assuming infinity has winding number 0.
-    const size_t num_labels = labels.maxCoeff() + 1;
-    const int INVALID = std::numeric_limits<int>::max();
-    patch_W.resize(num_patches, 2*num_labels);
-    patch_W.setConstant(INVALID);
-
-    size_t outer_facet_idx;
-    bool outer_facet_is_flipped;
-    Eigen::VectorXi face_indices(num_faces);
-    face_indices.setLinSpaced(num_faces, 0, num_faces-1);
-    igl::cgal::outer_facet(V, F, face_indices,
-            outer_facet_idx, outer_facet_is_flipped);
-    size_t outer_patch_idx = P[outer_facet_idx];
-    size_t outer_patch_label = labels[outer_patch_idx];
-    patch_W.row(outer_patch_idx).setZero();
-    if (outer_facet_is_flipped) {
-        patch_W(outer_patch_idx, outer_patch_label*2) = -1;
-    } else {
-        patch_W(outer_patch_idx, outer_patch_label*2+1) = 1;
-    }
-
-    auto winding_num_assigned = [&](size_t patch_idx) -> bool{
-        return (patch_W.row(patch_idx).array() != INVALID).all();
-    };
-
-    std::queue<size_t> Q;
-    Q.push(outer_patch_idx);
-    while (!Q.empty()) {
-        const size_t curr_patch_idx = Q.front();
-        const size_t curr_patch_label = labels[curr_patch_idx];
-        Q.pop();
-
-        const auto& adj_curves = patch_curve_adjacency[curr_patch_idx];
-        for (size_t curve_idx : adj_curves) {
-            const auto& order = orders[curve_idx];
-            const auto& orientation = orientations[curve_idx];
-            const size_t num_adj_patches = order.size();
-            assert(num_adj_patches == orientation.size());
-
-            size_t curr_i = std::numeric_limits<size_t>::max();
-            for (size_t i=0; i<num_adj_patches; i++) {
-                if (order[i] == curr_patch_idx) {
-                    curr_i = i;
-                    break;
-                }
-            }
-            assert(curr_i < num_adj_patches);
-            const bool curr_ori = orientation[curr_i];
-
-            //for (size_t i=0; i<num_adj_patches; i++) {
-            //    const size_t next_i = (curr_i + 1) % num_adj_patches;
-            //    const size_t num_patch_idx = order[next_i];
-            //    const bool next_ori = orientation[next_i];
-            //    const size_t next_patch_label = labels[next_patch_idx];
-            //    // TODO
-            //}
-
-            const size_t next_i = curr_ori ? (curr_i+1) % num_adj_patches
-                : (curr_i+num_adj_patches-1)%num_adj_patches;
-            const size_t prev_i = curr_ori ?
-                (curr_i+num_adj_patches-1)%num_adj_patches
-                : (curr_i+1) % num_adj_patches;
-            const size_t next_patch_idx = order[next_i];
-            const size_t prev_patch_idx = order[prev_i];
-
-            if (!winding_num_assigned(next_patch_idx)) {
-                const bool next_ori = orientation[next_i];
-                const bool next_cons = next_ori != curr_ori;
-                const size_t next_patch_label = labels[next_patch_idx];
-                for (size_t i=0; i<num_labels; i++) {
-                    const int shared_winding_number =
-                        patch_W(curr_patch_idx, i*2);
-
-                    if (i == next_patch_label) {
-                        // Truth table:
-                        // curr_ori  next_ori  wind_# inc
-                        // True      True      -1
-                        // True      False      1
-                        // False     True       1
-                        // False     False     -1
-
-                        patch_W(next_patch_idx, i*2+(next_cons ?0:1)) =
-                            shared_winding_number;
-                        patch_W(next_patch_idx, i*2+(next_cons ?1:0)) = 
-                            shared_winding_number + (next_cons ? 1:-1);
-                    } else {
-                        patch_W(next_patch_idx, i*2  ) = shared_winding_number;
-                        patch_W(next_patch_idx, i*2+1) = shared_winding_number;
-                    }
-                }
-                Q.push(next_patch_idx);
-            }
-            if (!winding_num_assigned(prev_patch_idx)) {
-                const bool prev_ori = orientation[prev_i];
-                const bool prev_cons = prev_ori != curr_ori;
-                const size_t prev_patch_label = labels[prev_patch_idx];
-
-                for (size_t i=0; i<num_labels; i++) {
-                    const int shared_winding_number =
-                        patch_W(curr_patch_idx, i*2+1);
-
-                    if (i == prev_patch_label) {
-                        // Truth table:
-                        // curr_ori  next_ori  wind_# inc
-                        // True      True       1
-                        // True      False     -1
-                        // False     True      -1
-                        // False     False      1
-
-                        patch_W(prev_patch_idx, i*2+(prev_cons ?1:0)) =
-                            shared_winding_number;
-                        patch_W(prev_patch_idx, i*2+(prev_cons ?0:1)) =
-                            shared_winding_number + (prev_cons ? -1:1);
-                    } else {
-                        patch_W(prev_patch_idx, i*2  ) = shared_winding_number; 
-                        patch_W(prev_patch_idx, i*2+1) = shared_winding_number; 
-                    }
-                }
-                Q.push(prev_patch_idx);
-            }
-        }
-    }
-    assert((patch_W.array() != INVALID).all());
-
-    using namespace propagate_winding_numbers_helper;
-    return winding_number_assignment_is_consistent(orders, orientations, patch_W);
-}
-
-#if 0
-template<
-    typename DerivedV,
-    typename DerivedF,
-    typename DeriveduE,
-    typename uE2EType,
-    typename DerivedEMAP,
-    typename DerivedC,
-    typename DerivedP,
-    typename DerivedW >
-IGL_INLINE bool igl::cgal::propagate_winding_numbers_single_component_patch_wise(
-        const Eigen::PlainObjectBase<DerivedV>& V,
-        const Eigen::PlainObjectBase<DerivedF>& F,
-        const Eigen::PlainObjectBase<DeriveduE>& uE,
-        const std::vector<std::vector<uE2EType> >& uE2E,
-        const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
-        const Eigen::PlainObjectBase<DerivedC>& labels,
-        const Eigen::PlainObjectBase<DerivedP>& P,
-        Eigen::PlainObjectBase<DerivedW>& patch_W) {
-    const size_t num_faces = F.rows();
-    const size_t num_patches = P.maxCoeff() + 1;
-    assert(labels.size() == num_patches);
-    // Utility functions.
-    auto edge_index_to_face_index = [&](size_t ei) { return ei % num_faces; };
-    auto edge_index_to_corner_index = [&](size_t ei) { return ei / num_faces; };
-    auto face_and_corner_index_to_edge_index = [&](size_t fi, size_t ci) {
-        return ci*num_faces + fi;
-    };
-    auto get_edge_end_points = [&](size_t ei, size_t& s, size_t& d) {
-        const size_t fi = edge_index_to_face_index(ei);
-        const size_t ci = edge_index_to_corner_index(ei);
-        s = F(fi, (ci+1)%3);
-        d = F(fi, (ci+2)%3);
-    };
-    auto is_positively_orientated =
-        [&](size_t fi, size_t s, size_t d) -> bool{
-            const Eigen::Vector3i f = F.row(fi);
-            if (f[0] == d && f[1] == s) return true;
-            if (f[1] == d && f[2] == s) return true;
-            if (f[2] == d && f[0] == s) return true;
-            if (f[0] == s && f[1] == d) return false;
-            if (f[1] == s && f[2] == d) return false;
-            if (f[2] == s && f[0] == d) return false;
-            throw std::runtime_error("Edge does not belong to face.");
-            return false;
-        };
-
-    auto compute_signed_index = [&](size_t fi, size_t s, size_t d) -> int{
-        return int(fi+1) * (is_positively_orientated(fi, s, d) ? 1:-1);
-    };
-    auto compute_unsigned_index = [&](int signed_idx) -> size_t {
-        return abs(signed_idx) - 1;
-    };
-
-    // Order patches around each non-manifold edge.
-    const size_t num_edges = uE.rows();
-    std::vector<Eigen::VectorXi> orders(num_edges);
-    std::vector<std::vector<bool> > orientations(num_edges);
-    std::vector<std::vector<size_t> > patch_edge_adjacency(num_patches);
-    for (size_t uei=0; uei<num_edges; uei++) {
-        const auto& edges = uE2E[uei];
-        if (edges.size() <= 2) continue;
-        size_t s = uE(uei, 0);
-        size_t d = uE(uei, 1);
-        std::vector<int> adj_faces;
-        for (size_t ei : uE2E[uei]) {
-            const size_t fi = edge_index_to_face_index(ei);
-            const size_t signed_fi =
-                compute_signed_index(fi, s, d);
-            adj_faces.push_back(signed_fi);
-            patch_edge_adjacency[P[fi]].push_back(ei);
-        }
-
-        auto& order = orders[uei];
-        igl::cgal::order_facets_around_edge(
-                V, F, s, d, adj_faces, order);
-        assert(order.minCoeff() == 0);
-        assert(order.maxCoeff() == adj_faces.size() - 1);
-
-        auto& orientation = orientations[uei];
-        orientation.resize(order.size());
-        std::transform(order.data(), order.data()+order.size(),
-                orientation.begin(), 
-                [&](size_t index) { return adj_faces[index] > 0; });
-        std::transform(order.data(), order.data()+order.size(),
-                order.data(), [&](size_t index) {
-                return compute_unsigned_index(adj_faces[index]);
-                } );
-    }
-
-    // Propagate winding number from infinity.
-    // Assuming infinity has winding number 0.
-    const size_t num_labels = labels.maxCoeff() + 1;
-    const int INVALID = std::numeric_limits<int>::max();
-    patch_W.resize(num_patches, 2*num_labels);
-    patch_W.setConstant(INVALID);
-
-    size_t outer_facet_idx;
-    bool outer_facet_is_flipped;
-    Eigen::VectorXi face_indices(num_faces);
-    face_indices.setLinSpaced(num_faces, 0, num_faces-1);
-    igl::cgal::outer_facet(V, F, face_indices,
-            outer_facet_idx, outer_facet_is_flipped);
-    size_t outer_patch_idx = P[outer_facet_idx];
-    size_t outer_patch_label = labels[outer_patch_idx];
-    patch_W.row(outer_patch_idx).setZero();
-    if (outer_facet_is_flipped) {
-        patch_W(outer_patch_idx, outer_patch_label*2) = -1;
-    } else {
-        patch_W(outer_patch_idx, outer_patch_label*2+1) = 1;
-    }
-
-    auto winding_num_assigned = [&](size_t patch_idx) -> bool{
-        return (patch_W.row(patch_idx).array() != INVALID).all();
-    };
-
-    std::queue<size_t> Q;
-    Q.push(outer_patch_idx);
-    while (!Q.empty()) {
-        const size_t curr_patch_idx = Q.front();
-        const size_t curr_patch_label = labels[curr_patch_idx];
-        Q.pop();
-
-        const auto& adj_edges = patch_edge_adjacency[curr_patch_idx];
-        for (size_t ei: adj_edges) {
-            const size_t uei = EMAP[ei];
-            const size_t fid = edge_index_to_face_index(ei);
-            const auto& order = orders[uei];
-            const auto& orientation = orientations[uei];
-            const size_t num_adj_patches = order.size();
-            assert(num_adj_patches == orientation.size());
-
-            size_t curr_i = std::numeric_limits<size_t>::max();
-            for (size_t i=0; i<num_adj_patches; i++) {
-                std::cout << (orientation[i]?"+":"-") << order[i] << " ";
-            }
-            std::cout << std::endl;
-            for (size_t i=0; i<num_adj_patches; i++) {
-                if (order[i] == fid) {
-                    curr_i = i;
-                    break;
-                }
-            }
-            assert(curr_i < num_adj_patches);
-            const bool curr_ori = orientation[curr_i];
-
-            const size_t next_i = curr_ori ? (curr_i+1) % num_adj_patches
-                : (curr_i+num_adj_patches-1)%num_adj_patches;
-            const size_t prev_i = curr_ori ?
-                (curr_i+num_adj_patches-1)%num_adj_patches
-                : (curr_i+1) % num_adj_patches;
-            const size_t next_patch_idx = P[order[next_i]];
-            const size_t prev_patch_idx = P[order[prev_i]];
-
-            if (!winding_num_assigned(next_patch_idx)) {
-                const bool next_ori = orientation[next_i];
-                const bool next_cons = next_ori != curr_ori;
-                const size_t next_patch_label = labels[next_patch_idx];
-                for (size_t i=0; i<num_labels; i++) {
-                    const int shared_winding_number =
-                        patch_W(curr_patch_idx, i*2);
-
-                    if (i == next_patch_label) {
-                        // Truth table:
-                        // curr_ori  next_ori  wind_# inc
-                        // True      True      -1
-                        // True      False      1
-                        // False     True       1
-                        // False     False     -1
-
-                        patch_W(next_patch_idx, i*2+(next_cons ?0:1)) =
-                            shared_winding_number;
-                        patch_W(next_patch_idx, i*2+(next_cons ?1:0)) = 
-                            shared_winding_number + (next_cons ? 1:-1);
-                    } else {
-                        patch_W(next_patch_idx, i*2  ) = shared_winding_number;
-                        patch_W(next_patch_idx, i*2+1) = shared_winding_number;
-                    }
-                }
-                Q.push(next_patch_idx);
-            } else {
-                const bool next_ori = orientation[next_i];
-                const bool next_cons = next_ori != curr_ori;
-                const size_t next_patch_label = labels[next_patch_idx];
-                for (size_t i=0; i<num_labels; i++) {
-                    const int shared_winding_number =
-                        patch_W(curr_patch_idx, i*2);
-
-                    if (i == next_patch_label) {
-                        // Truth table:
-                        // curr_ori  next_ori  wind_# inc
-                        // True      True      -1
-                        // True      False      1
-                        // False     True       1
-                        // False     False     -1
-
-                        assert(patch_W(next_patch_idx, i*2+(next_cons ?0:1)) ==
-                            shared_winding_number);
-                        assert(patch_W(next_patch_idx, i*2+(next_cons ?1:0)) == 
-                            shared_winding_number + (next_cons ? 1:-1));
-                    } else {
-                        assert(patch_W(next_patch_idx, i*2) ==
-                                shared_winding_number);
-                        assert(patch_W(next_patch_idx, i*2+1) ==
-                                shared_winding_number);
-                    }
-                }
-            }
-            if (!winding_num_assigned(prev_patch_idx)) {
-                const bool prev_ori = orientation[prev_i];
-                const bool prev_cons = prev_ori != curr_ori;
-                const size_t prev_patch_label = labels[prev_patch_idx];
-
-                for (size_t i=0; i<num_labels; i++) {
-                    const int shared_winding_number =
-                        patch_W(curr_patch_idx, i*2+1);
-
-                    if (i == prev_patch_label) {
-                        // Truth table:
-                        // curr_ori  next_ori  wind_# inc
-                        // True      True       1
-                        // True      False     -1
-                        // False     True      -1
-                        // False     False      1
-
-                        patch_W(prev_patch_idx, i*2+(prev_cons ?1:0)) =
-                            shared_winding_number;
-                        patch_W(prev_patch_idx, i*2+(prev_cons ?0:1)) =
-                            shared_winding_number + (prev_cons ? -1:1);
-                    } else {
-                        patch_W(prev_patch_idx, i*2  ) = shared_winding_number; 
-                        patch_W(prev_patch_idx, i*2+1) = shared_winding_number; 
-                    }
-                }
-                Q.push(prev_patch_idx);
-            } else {
-                const bool prev_ori = orientation[prev_i];
-                const bool prev_cons = prev_ori != curr_ori;
-                const size_t prev_patch_label = labels[prev_patch_idx];
-
-                for (size_t i=0; i<num_labels; i++) {
-                    const int shared_winding_number =
-                        patch_W(curr_patch_idx, i*2+1);
-
-                    if (i == prev_patch_label) {
-                        // Truth table:
-                        // curr_ori  next_ori  wind_# inc
-                        // True      True       1
-                        // True      False     -1
-                        // False     True      -1
-                        // False     False      1
-
-                        assert(patch_W(prev_patch_idx, i*2+(prev_cons ?1:0)) ==
-                            shared_winding_number);
-                        assert(patch_W(prev_patch_idx, i*2+(prev_cons ?0:1)) ==
-                            shared_winding_number + (prev_cons ? -1:1));
-                    } else {
-                        assert(patch_W(prev_patch_idx, i*2  ) ==
-                                shared_winding_number); 
-                        assert(patch_W(prev_patch_idx, i*2+1) ==
-                                shared_winding_number); 
-                    }
-                }
-            }
-        }
-    }
-    assert((patch_W.array() != INVALID).all());
-
-    using namespace propagate_winding_numbers_helper;
-    if (!winding_number_assignment_is_consistent(orders, orientations, patch_W)) {
-        Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic,
-            DerivedV::IsRowMajor> Vd(V.rows(), V.cols());
-        for (size_t j=0; j<V.rows(); j++) {
-            for (size_t k=0; k<V.cols(); k++) {
-                igl::cgal::assign_scalar(V(j,k), Vd(j,k));
-            }
-        }
-        //for (size_t j=0; j<num_faces; j++) {
-        //    std::cout << patch_W(P[j], 1) << std::endl;
-        //}
-        igl::writeOBJ("debug_wn.obj", Vd, F);
-        assert(false);
-    }
-    return winding_number_assignment_is_consistent(orders, orientations, patch_W);
-}
-#endif
-
-template<
-typename DerivedV,
-typename DerivedF,
-typename DerivedC,
-typename DerivedW>
-IGL_INLINE bool igl::cgal::propagate_winding_numbers_single_component(
-        const Eigen::PlainObjectBase<DerivedV>& V,
-        const Eigen::PlainObjectBase<DerivedF>& F,
-        const Eigen::PlainObjectBase<DerivedC>& labels,
-        Eigen::PlainObjectBase<DerivedW>& W) {
-    typedef typename DerivedF::Scalar Index;
-    const size_t num_faces = F.rows();
-
-    // Extract unique edges.
-    std::vector<std::vector<size_t> > uE2E;
-    Eigen::Matrix<Index, Eigen::Dynamic, Eigen::Dynamic> E, uE;
-    Eigen::Matrix<Index, Eigen::Dynamic, 1> EMAP;
-    igl::unique_edge_map(F, E, uE, EMAP, uE2E);
-
-    // Extract manifold patches and intersection curves.
-    Eigen::VectorXi P;
-    std::vector<std::vector<size_t> > intersection_curves;
-    size_t num_patches =
-        igl::extract_manifold_patches(F, EMAP, uE2E, P);
-    igl::extract_non_manifold_edge_curves(F, EMAP, uE2E,
-            intersection_curves);
-    assert(P.size() == num_faces);
-    assert(P.maxCoeff() + 1 == num_patches);
-
-    Eigen::VectorXi patch_labels(num_patches);
-    const int INVALID = std::numeric_limits<int>::max();
-    patch_labels.setConstant(INVALID);
-    for (size_t i=0; i<num_faces; i++) {
-        if (patch_labels[P[i]] == INVALID) {
-            patch_labels[P[i]] = labels[i];
-        } else {
-            assert(patch_labels[P[i]] == labels[i]);
-        }
-    }
-    assert((patch_labels.array() != INVALID).all());
-    const size_t num_labels = patch_labels.maxCoeff()+1;
-
-    Eigen::MatrixXi winding_numbers;
-    bool is_consistent =
-        igl::cgal::propagate_winding_numbers_single_component_patch_wise(
-            V, F, uE, uE2E, patch_labels, P, intersection_curves, winding_numbers);
-    assert(winding_numbers.rows() == num_patches);
-    assert(winding_numbers.cols() == 2 * num_labels);
-
-    W.resize(num_faces, 2*num_labels);
-    W.setConstant(INVALID);
-    for (size_t i=0; i<num_faces; i++) {
-        W.row(i) = winding_numbers.row(P[i]);
-    }
-    assert((W.array() != INVALID).any());
-
-    return is_consistent;
-}
-
-template<
-typename DerivedV,
-typename DerivedF,
-typename DerivedW>
-IGL_INLINE bool igl::cgal::propagate_winding_numbers_single_component(
-        const Eigen::PlainObjectBase<DerivedV>& V,
-        const Eigen::PlainObjectBase<DerivedF>& F,
-        Eigen::PlainObjectBase<DerivedW>& W) {
-    const size_t num_faces = F.rows();
-    Eigen::VectorXi labels(num_faces);
-    labels.setZero();
-    return igl::cgal::propagate_winding_numbers_single_component(V, F, labels, W);
-}
-
 template<
 typename DerivedV,
 typename DerivedF,
@@ -690,149 +76,6 @@ IGL_INLINE void igl::cgal::propagate_winding_numbers(
         const Eigen::PlainObjectBase<DerivedF>& F,
         const Eigen::PlainObjectBase<DerivedL>& labels,
         Eigen::PlainObjectBase<DerivedW>& W) {
-    typedef typename DerivedF::Scalar Index;
-    const size_t num_faces = F.rows();
-
-    // Extract unique edges.
-    std::vector<std::vector<size_t> > uE2E;
-    Eigen::Matrix<Index, Eigen::Dynamic, Eigen::Dynamic> E, uE;
-    Eigen::Matrix<Index, Eigen::Dynamic, 1> EMAP;
-    igl::unique_edge_map(F, E, uE, EMAP, uE2E);
-
-    // Check to make sure there is no boundaries and no non-manifold edges with
-    // odd number of adjacent faces.
-    for (const auto& adj_faces : uE2E) {
-        if (adj_faces.size() % 2 == 1) {
-            std::stringstream err_msg;
-            err_msg << "Input mesh contains odd number of faces "
-                << "sharing a single edge" << std::endl;
-            err_msg << "Indicating the input mesh does not represent a valid volume"
-                << ", and winding number cannot be propagated." << std::endl;
-            throw std::runtime_error(err_msg.str());
-        }
-    }
-
-    // Gather connected components.
-    std::vector<std::vector<std::vector<Index > > > TT,_1;
-    triangle_triangle_adjacency(E,EMAP,uE2E,false,TT,_1);
-    Eigen::VectorXi counts;
-    Eigen::VectorXi C;
-    igl::facet_components(TT,C,counts);
-
-    const size_t num_components = counts.size();
-    std::vector<std::vector<size_t> > components(num_components);
-    for (size_t i=0; i<num_faces; i++) {
-        components[C[i]].push_back(i);
-    }
-    std::vector<Eigen::MatrixXi> comp_faces(num_components);
-    std::vector<Eigen::VectorXi> comp_labels(num_components);
-    for (size_t i=0; i<num_components; i++) {
-        const auto& comp = components[i];
-        auto& faces = comp_faces[i];
-        auto& c_labels = comp_labels[i];
-        const size_t comp_size = comp.size();
-        faces.resize(comp_size, 3);
-        c_labels.resize(comp_size);
-        for (size_t j=0; j<comp_size; j++) {
-            faces.row(j) = F.row(comp[j]);
-            c_labels[j] = labels[comp[j]];
-        }
-    }
-
-    // Compute winding number for each component.
-    const size_t num_labels = labels.maxCoeff()+1;
-    W.resize(num_faces, 2*num_labels);
-    W.setZero();
-    std::vector<Eigen::MatrixXi> Ws(num_components);
-    for (size_t i=0; i<num_components; i++) {
-        bool is_consistent =
-            propagate_winding_numbers_single_component(V, comp_faces[i], comp_labels[i], Ws[i]);
-        const size_t num_labels_in_i = comp_labels[i].maxCoeff()+1;
-        const size_t num_faces_in_comp = comp_faces[i].rows();
-        assert(Ws[i].cols() == num_labels_in_i*2);
-        assert(Ws[i].rows() == num_faces_in_comp);
-        const auto& comp = components[i];
-        for (size_t j=0; j<num_faces_in_comp; j++) {
-            const size_t fid = comp[j];
-            W.block(fid, 0, 1, num_labels_in_i*2) = Ws[i].row(j);
-        }
-
-        if (!is_consistent) {
-            Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic,
-                DerivedV::IsRowMajor> Vd(V.rows(), V.cols());
-            for (size_t j=0; j<V.rows(); j++) {
-                for (size_t k=0; k<V.cols(); k++) {
-                    igl::cgal::assign_scalar(V(j,k), Vd(j,k));
-                }
-            }
-            igl::writeOBJ("debug_wn.obj", Vd, comp_faces[i]);
-            std::cout << Ws[i].col(1) << std::endl;
-            std::stringstream err_msg;
-            err_msg << "Component " << i
-                << " has inconsistant winding number assignment." << std::endl;
-            //throw std::runtime_error(err_msg.str());
-            std::cout << err_msg.str() << std::endl;
-        }
-    }
-
-    auto sample_component = [&](size_t cid) {
-        const auto& f = comp_faces[cid].row(0).eval();
-        return ((V.row(f[0]) + V.row(f[1]) + V.row(f[2])) / 3.0).eval();
-    };
-
-    std::vector<Eigen::MatrixXi> nested_Ws = Ws;
-    Eigen::MatrixXi ambient_correction(num_components, 2*num_labels);
-    ambient_correction.setZero();
-    if (num_components > 1) {
-        for (size_t i=0; i<num_components; i++) {
-            DerivedV samples(num_components-1, 3);
-            Eigen::VectorXi is_inside;
-            auto index_without_i = [&](size_t index) {
-                return index < i ? index:index-1;
-            };
-            for (size_t j=0; j<num_components; j++) {
-                if (i == j) continue;
-                samples.row(index_without_i(j)) = sample_component(j);
-            }
-            Eigen::VectorXi fids;
-            Eigen::Matrix<bool, Eigen::Dynamic, 1> orientation;
-            igl::cgal::closest_facet(V, comp_faces[i], samples,
-                    fids, orientation);
-
-            const auto& comp = components[i];
-            for (size_t j=0; j<num_components; j++) {
-                if (i == j) continue;
-                const size_t index = index_without_i(j);
-                const size_t fid = fids(index, 0);
-                const bool ori = orientation(index, 0);
-                for (size_t k=0; k<num_labels; k++) {
-                    const int correction = W(comp[fid], k*2+(ori?0:1));
-                    ambient_correction(j, k*2  ) += correction;
-                    ambient_correction(j, k*2+1) += correction;
-                }
-            }
-        }
-    }
-
-    for (size_t i=0; i<num_components; i++) {
-        const auto& comp = components[i];
-        const auto& correction = ambient_correction.row(i).eval();
-        for (const auto& fid : comp) {
-            W.row(fid) += correction;
-        }
-    }
-}
-
-template<
-typename DerivedV,
-typename DerivedF,
-typename DerivedL,
-typename DerivedW>
-IGL_INLINE void igl::cgal::propagate_winding_numbers_beta(
-        const Eigen::PlainObjectBase<DerivedV>& V,
-        const Eigen::PlainObjectBase<DerivedF>& F,
-        const Eigen::PlainObjectBase<DerivedL>& labels,
-        Eigen::PlainObjectBase<DerivedW>& W) {
     const size_t num_faces = F.rows();
     typedef typename DerivedF::Scalar Index;
 

include/igl/cgal/propagate_winding_numbers.h

@@ -21,129 +21,6 @@
 
 namespace igl {
     namespace cgal {
-        // Compute winding number on each side of each patch.  All patches must
-        // be connected.  The per-patch label partitions the patches into
-        // different components, and winding number is computed for each
-        // component.  For example, ``patch_W(i, j*2)`` for ``i`` in [0, #P)
-        // and ``j`` in [0, k) is the winding number on the positive side of
-        // patch ``i`` with respect to component ``j``.  Similarly,
-        // ``patch_W(i, j*2+1)`` is the winding number on the negative side of
-        // patch ``i`` with respect to component ``j``.
-        //
-        // Patch and intersection curve can be generated using:
-        // ``igl::extract_manifold_patches(...)`` and
-        // ``igl::extract_non_manifold_edge_curves(...)``
-        //
-        //
-        // Inputs:
-        //   V  #V by 3 list of vertices.
-        //   F  #F by 3 list of trianlges.
-        //   uE #uE by 2 list of undirected edges.
-        //   uE2E  #uE list of lists mapping each undirected edge to directed
-        //         edges.
-        //   labels  #P list of labels, ranging from 0 to k-1.
-        //           These labels partition patches into k components, and
-        //           winding number is computed for each component.
-        //   P  #F list of patch indices.
-        //   intersection_curves  A list of non-manifold curves that separates
-        //                        the mesh into patches.
-        //
-        // Outputs:
-        //   patch_W  #P by k*2 list of winding numbers on each side of a
-        //            patch for each component.
-        //
-        // Returns:
-        //   True iff integer winding numbers can be consistently assigned.
-        template<
-            typename DerivedV,
-            typename DerivedF,
-            typename DeriveduE,
-            typename uE2EType,
-            typename DerivedC,
-            typename DerivedP,
-            typename DerivedW >
-        IGL_INLINE bool propagate_winding_numbers_single_component_patch_wise(
-                const Eigen::PlainObjectBase<DerivedV>& V,
-                const Eigen::PlainObjectBase<DerivedF>& F,
-                const Eigen::PlainObjectBase<DeriveduE>& uE,
-                const std::vector<std::vector<uE2EType> >& uE2E,
-                const Eigen::PlainObjectBase<DerivedC>& labels,
-                const Eigen::PlainObjectBase<DerivedP>& P,
-                const std::vector<std::vector<size_t> >& intersection_curves,
-                Eigen::PlainObjectBase<DerivedW>& patch_W);
-
-#if 0
-        template<
-            typename DerivedV,
-            typename DerivedF,
-            typename DeriveduE,
-            typename uE2EType,
-            typename DerivedEMAP,
-            typename DerivedC,
-            typename DerivedP,
-            typename DerivedW >
-        IGL_INLINE bool propagate_winding_numbers_single_component_patch_wise(
-                const Eigen::PlainObjectBase<DerivedV>& V,
-                const Eigen::PlainObjectBase<DerivedF>& F,
-                const Eigen::PlainObjectBase<DeriveduE>& uE,
-                const std::vector<std::vector<uE2EType> >& uE2E,
-                const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
-                const Eigen::PlainObjectBase<DerivedC>& labels,
-                const Eigen::PlainObjectBase<DerivedP>& P,
-                Eigen::PlainObjectBase<DerivedW>& patch_W);
-#endif
-
-        // Compute winding number on each side of the face.  The input mesh must
-        // be a single edge-connected component.
-        //
-        // Inputs:
-        //   V  #V by 3 list of vertex positions.
-        //   F  #F by 3 list of triangle indices into V.
-        //   labels  #F list of facet labels ranging from 0 to k-1.
-        //
-        // Output:
-        //   W  #F by 2*k list of winding numbers.
-        //      ``W(i, 2*j)`` is the winding number on the positive side of
-        //      facet ``i`` with respect to facets labelled ``j``.
-        //      ``W(i, 2*j+1)`` is the winding number on the negative side of
-        //      facet ``i`` with respect to facets labelled ``j``.
-        //
-        // Returns:
-        //   True iff integer winding number can be consistently assigned.
-        template<
-            typename DerivedV,
-            typename DerivedF,
-            typename DerivedC,
-            typename DerivedW>
-        IGL_INLINE bool propagate_winding_numbers_single_component(
-                const Eigen::PlainObjectBase<DerivedV>& V,
-                const Eigen::PlainObjectBase<DerivedF>& F,
-                const Eigen::PlainObjectBase<DerivedC>& labels,
-                Eigen::PlainObjectBase<DerivedW>& W);
-
-        // Compute the winding number on each side of the face.
-        // This method assumes the input mesh (V, F) forms a single connected
-        // component.
-        //
-        // Inputs:
-        //   V  #V by 3 list of vertex positions.
-        //   F  #F by 3 list of triangle indices into V.
-        //
-        // Output:
-        //   W  #F by 2 list of winding numbers.  W(i,0) is the winding number
-        //   on the positive side of facet i, and W(i, 1) is the winding
-        //   number on the negative side of facet I[i].
-        //
-        // Returns:
-        //   True iff integer winding number can be consistently assigned.
-        template<
-            typename DerivedV,
-            typename DerivedF,
-            typename DerivedW>
-        IGL_INLINE bool propagate_winding_numbers_single_component(
-                const Eigen::PlainObjectBase<DerivedV>& V,
-                const Eigen::PlainObjectBase<DerivedF>& F,
-                Eigen::PlainObjectBase<DerivedW>& W);
 
         // Compute winding number on each side of the face.  The input mesh
         // could contain multiple connected components.  The input mesh must
@@ -171,17 +48,6 @@ namespace igl {
                 const Eigen::PlainObjectBase<DerivedF>& F,
                 const Eigen::PlainObjectBase<DerivedL>& labels,
                 Eigen::PlainObjectBase<DerivedW>& W);
-
-        template<
-            typename DerivedV,
-            typename DerivedF,
-            typename DerivedL,
-            typename DerivedW>
-        IGL_INLINE void propagate_winding_numbers_beta(
-                const Eigen::PlainObjectBase<DerivedV>& V,
-                const Eigen::PlainObjectBase<DerivedF>& F,
-                const Eigen::PlainObjectBase<DerivedL>& labels,
-                Eigen::PlainObjectBase<DerivedW>& W);
     }
 }
 

include/igl/cgal/remesh_intersections.cpp

@@ -1,328 +1,338 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 // 
-// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
+// Copyright (C) 2015 Qingnan Zhou <qnzhou@gmail.com>
 // 
 // This Source Code Form is subject to the terms of the Mozilla Public License 
 // v. 2.0. If a copy of the MPL was not distributed with this file, You can 
 // obtain one at http://mozilla.org/MPL/2.0/.
+//
 #include "remesh_intersections.h"
-#include "SelfIntersectMesh.h"
-#include "assign_scalar.h"
-#include "projected_delaunay.h"
+
+#include <vector>
+#include <map>
+#include <unordered_map>
 #include <iostream>
-#include <cassert>
 
 template <
-  typename DerivedV,
-  typename DerivedF,
-  typename Kernel,
-  typename DerivedVV,
-  typename DerivedFF,
-  typename DerivedJ,
-  typename DerivedIM>
+typename DerivedV,
+typename DerivedF,
+typename Kernel,
+typename DerivedVV,
+typename DerivedFF,
+typename DerivedJ,
+typename DerivedIM>
 IGL_INLINE void igl::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::pair<typename DerivedF::Index,
-      std::vector<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)
-{
-  using namespace std;
-  using namespace Eigen;
-  typedef typename DerivedF::Index          Index;
-  typedef CGAL::Point_3<Kernel>    Point_3;
-  //typedef CGAL::Segment_3<Kernel>  Segment_3; 
-  //typedef CGAL::Triangle_3<Kernel> Triangle_3; 
-  typedef CGAL::Plane_3<Kernel>    Plane_3;
-  //typedef CGAL::Point_2<Kernel>    Point_2;
-  //typedef CGAL::Segment_2<Kernel>  Segment_2; 
-  //typedef CGAL::Triangle_2<Kernel> Triangle_2; 
-  typedef CGAL::Triangulation_vertex_base_2<Kernel>  TVB_2;
-  typedef CGAL::Constrained_triangulation_face_base_2<Kernel> CTFB_2;
-  typedef CGAL::Triangulation_data_structure_2<TVB_2,CTFB_2> TDS_2;
-  typedef CGAL::Exact_intersections_tag Itag;
-  typedef CGAL::Constrained_Delaunay_triangulation_2<Kernel,TDS_2,Itag> 
-    CDT_2;
-  typedef CGAL::Constrained_triangulation_plus_2<CDT_2> CDT_plus_2;
-  typedef std::pair<Index,Index> EMK;
-  typedef std::vector<Index> EMV;
-  //typedef std::map<EMK,EMV> EdgeMap;
-  typedef std::pair<Index,Index> EMK;
-  //typedef std::vector<CGAL::Object> ObjectList;
-  typedef std::vector<Index> IndexList;
+        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::pair<typename DerivedF::Index,
+        std::vector<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) {
+
+    typedef CGAL::Point_3<Kernel>    Point_3;
+    typedef CGAL::Segment_3<Kernel>  Segment_3; 
+    typedef CGAL::Triangle_3<Kernel> Triangle_3; 
+    typedef CGAL::Plane_3<Kernel>    Plane_3;
+    typedef CGAL::Point_2<Kernel>    Point_2;
+    typedef CGAL::Segment_2<Kernel>  Segment_2; 
+    typedef CGAL::Triangle_2<Kernel> Triangle_2; 
+    typedef CGAL::Triangulation_vertex_base_2<Kernel>  TVB_2;
+    typedef CGAL::Constrained_triangulation_face_base_2<Kernel> CTFB_2;
+    typedef CGAL::Triangulation_data_structure_2<TVB_2,CTFB_2> TDS_2;
+    typedef CGAL::Exact_intersections_tag Itag;
+    typedef CGAL::Constrained_Delaunay_triangulation_2<Kernel,TDS_2,Itag> 
+        CDT_2;
+    typedef CGAL::Constrained_triangulation_plus_2<CDT_2> CDT_plus_2;
+
+    typedef typename DerivedF::Index Index;
+    typedef std::pair<Index, Index> Edge;
+    struct EdgeHash {
+        size_t operator()(const Edge& e) const {
+            return (e.first * 805306457) ^ (e.second * 201326611);
+        }
+    };
+    typedef std::unordered_map<Edge, std::vector<Index>, EdgeHash > EdgeMap;
+
+    auto normalize_plane_coeff = [](const Plane_3& P) {
+        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) {
+        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;
+    };
+    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;
 
-  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;
-#ifdef IGL_SELFINTERSECTMESH_DEBUG
-  double t_proj_del = 0;
-#endif
-  // Unfortunately it looks like CGAL has trouble allocating memory when
-  // multiple openmp threads are running. Crashes durring CDT...
-  //// Loop over offending triangles
-  //const size_t noff = offending.size();
-//# pragma omp parallel for if (noff>1000)
-  for(const auto & okv : offending)
-  {
-    // index in F
-    const Index f = okv.first;
-    const Index o = okv.second.first;
-    {
-#ifdef IGL_SELFINTERSECTMESH_DEBUG
-      const double t_before = get_seconds();
-#endif
-      CDT_plus_2 cdt_o;
-      projected_delaunay(T[f],okv.second.second,cdt_o);
-      cdt[o] = cdt_o;
-#ifdef IGL_SELFINTERSECTMESH_DEBUG
-      t_proj_del += (get_seconds()-t_before);
-#endif
+        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);
+        }
     }
-    // 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)
-      {
-        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()));
-#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++)
-          {
-            // 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)
-            {
-              swap(i,j);
+
+    const size_t INVALID = std::numeric_limits<size_t>::max();
+    std::vector<std::vector<Index> > resolved_faces;
+    std::vector<Index> source_faces;
+    std::vector<Point_3> new_vertices;
+    EdgeMap edge_vertices;
+
+    /**
+     * Run constraint Delaunay triangulation on the plane.
+     */
+    auto run_delaunay_triangulation = [&](const Plane_3& P,
+            const std::vector<Index>& involved_faces,
+            std::vector<Point_3>& vertices,
+            std::vector<std::vector<Index> >& faces) {
+        CDT_plus_2 cdt;
+        for (const auto& fid : involved_faces) {
+            const auto itr = offending.find(fid);
+            const auto& triangle = T[fid];
+            cdt.insert_constraint(P.to_2d(triangle[0]), P.to_2d(triangle[1]));
+            cdt.insert_constraint(P.to_2d(triangle[1]), P.to_2d(triangle[2]));
+            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) {
+                if(const Segment_3 *iseg = CGAL::object_cast<Segment_3 >(&obj)) {
+                    // Add segment constraint
+                    cdt.insert_constraint(
+                            P.to_2d(iseg->vertex(0)),P.to_2d(iseg->vertex(1)));
+                }else if(const Point_3 *ipoint = CGAL::object_cast<Point_3 >(&obj)) {
+                    // Add point
+                    cdt.insert(P.to_2d(*ipoint));
+                } else if(const Triangle_3 *itri = CGAL::object_cast<Triangle_3 >(&obj)) {
+                    // Add 3 segment constraints
+                    cdt.insert_constraint(
+                            P.to_2d(itri->vertex(0)),P.to_2d(itri->vertex(1)));
+                    cdt.insert_constraint(
+                            P.to_2d(itri->vertex(1)),P.to_2d(itri->vertex(2)));
+                    cdt.insert_constraint(
+                            P.to_2d(itri->vertex(2)),P.to_2d(itri->vertex(0)));
+                } else if(const std::vector<Point_3 > *polyp = 
+                        CGAL::object_cast< std::vector<Point_3 > >(&obj)) {
+                    //cerr<<REDRUM("Poly...")<<endl;
+                    const std::vector<Point_3 > & poly = *polyp;
+                    const Index m = poly.size();
+                    assert(m>=2);
+                    for(Index p = 0;p<m;p++)
+                    {
+                        const Index np = (p+1)%m;
+                        cdt.insert_constraint(P.to_2d(poly[p]),P.to_2d(poly[np]));
+                    }
+                }else {
+                    throw std::runtime_error("Unknown intersection object!");
+                }
             }
-            assert(edge2faces.count(EMK(i,j))==1);
-            const EMV & facesij = edge2faces.find(EMK(i,j))->second;
-            // loop over neighbors
-            for(
-              typename IndexList::const_iterator nit = facesij.begin();
-              nit != facesij.end() && !found;
-              nit++)
-            {
-              // don't consider self
-              if(*nit == f)
-              {
-                continue;
-              }
-              // index of neighbor in offending (to find its cdt)
-              assert(offending.count(*nit) == 1);
-              Index no = offending.find(*nit)->second.first;
-              // 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)
-              {
-                if(vit_point_3 == P[no].to_3d(uit->point()))
-                {
-                  assert(v2i.count(uit) == 1);
-//#   pragma omp critical
-                  v2i[vit] = v2i[uit];
-                  found = true;
+        }
+        const size_t num_vertices = cdt.number_of_vertices();
+        const size_t num_faces = cdt.number_of_faces();
+        std::map<typename CDT_plus_2::Vertex_handle,Index> v2i;
+        size_t count=0;
+        for (auto itr = cdt.finite_vertices_begin();
+                itr != cdt.finite_vertices_end(); itr++) {
+            vertices.push_back(P.to_3d(itr->point()));
+            v2i[itr] = count;
+            count++;
+        }
+        for (auto itr = cdt.finite_faces_begin();
+                itr != cdt.finite_faces_end(); itr++) {
+            faces.push_back( {
+                    v2i[itr->vertex(0)],
+                    v2i[itr->vertex(1)],
+                    v2i[itr->vertex(2)] });
+        }
+    };
+
+    /**
+     * 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 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];
+        }
+
+        // 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;
                 }
-              }
             }
-          }
-          if(!found)
-          {
-//#   pragma omp critical
-            {
-              v2i[vit] = V.rows()+NV_count;
-              NV.push_back(vit_point_3);
-              NV_count++;
+        }
+
+        // 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;
+    };
+
+    /**
+     * Determine the vertex indices for each corner of each output triangle.
+     */
+    auto post_triangulation_process = [&](
+            const std::vector<Point_3>& vertices,
+            const std::vector<std::vector<Index> >& faces,
+            const std::vector<Index>& involved_faces) {
+        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);
+                }
             }
-          }
         }
-        i++;
-      }
+    };
+
+    // Process un-touched faces.
+    for (size_t i=0; i<num_faces; i++) {
+        if (!is_offending[i]) {
+            resolved_faces.push_back(
+                    { F(i,0), F(i,1), F(i,2) } );
+            source_faces.push_back(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++;
-      }
+
+    // Process self-intersecting 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);
     }
-  }
-#ifdef IGL_SELFINTERSECTMESH_DEBUG
-  cout<<"CDT: "<<tictoc()<<"  "<<t_proj_del<<endl;
-#endif
 
-  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);
-  //  }
-  //  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(!offending.count(f))
-    {
-      FF.row(off) = F.row(f);
-      J(off) = f;
-      off++;
+    // Output resolved mesh.
+    const size_t num_out_vertices = new_vertices.size() + num_base_vertices;
+    VV.resize(num_out_vertices, 3);
+    for (size_t i=0; i<num_base_vertices; i++) {
+        assign_scalar(V(i,0), VV(i,0));
+        assign_scalar(V(i,1), VV(i,1));
+        assign_scalar(V(i,2), VV(i,2));
     }
-  }
-  assert(off == (Index)(F.rows()-offending.size()));
-  // Now append replacement faces for offending faces
-  for(const auto & okv : offending)
-  {
-    // index in F
-    const Index f = okv.first;
-    const Index o = okv.second.first;
-    FF.block(off,0,NF[o].rows(),3) = NF[o];
-    J.block(off,0,NF[o].rows(),1).setConstant(f);
-    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++)
-    {
-      for(Index d = 0;d<3;d++)
-      {
-        const Point_3 & p = *nvit;
-        // Don't convert via double if output type is same as Kernel
-        assign_scalar(p[d], VV(V.rows()+i,d));
-      }
-      i++;
+    for (size_t i=num_base_vertices; i<num_out_vertices; i++) {
+        assign_scalar(new_vertices[i-num_base_vertices][0], VV(i,0));
+        assign_scalar(new_vertices[i-num_base_vertices][1], VV(i,1));
+        assign_scalar(new_vertices[i-num_base_vertices][2], VV(i,2));
     }
-  }
-  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++)
-  {
-    typename Kernel::FT p0,p1,p2;
-    assign_scalar(V(v,0),p0);
-    assign_scalar(V(v,1),p1);
-    assign_scalar(V(v,2),p2);
-    const Point_3 p(p0,p1,p2);
-    if(vv2i.count(p)==0)
-    {
-      vv2i[p] = v;
+
+    const size_t num_out_faces = resolved_faces.size();
+    FF.resize(num_out_faces, 3);
+    for (size_t i=0; i<num_out_faces; i++) {
+        FF(i,0) = resolved_faces[i][0];
+        FF(i,1) = resolved_faces[i][1];
+        FF(i,2) = resolved_faces[i][2];
     }
-    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++)
-    {
-      const Point_3 & p = *nvit;
-      if(vv2i.count(p)==0)
-      {
-        vv2i[p] = v;
-      }
-      assert(vv2i.count(p) == 1);
-      IM(v) = vv2i[p];
-      v++;
+
+    J.resize(num_out_faces);
+    std::copy(source_faces.begin(), source_faces.end(), J.data());
+
+    // Extract unique vertex indices.
+    IM.resize(VV.rows(),1);
+    std::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<VV.rows();v++) {
+        typename Kernel::FT p0,p1,p2;
+        assign_scalar(VV(v,0),p0);
+        assign_scalar(VV(v,1),p1);
+        assign_scalar(VV(v,2),p2);
+        const Point_3 p(p0,p1,p2);
+        if(vv2i.count(p)==0) {
+            vv2i[p] = v;
+        }
+        assert(vv2i.count(p) == 1);
+        IM(v) = vv2i[p];
     }
-  }
-#ifdef IGL_SELFINTERSECTMESH_DEBUG
-  cout<<"Output + dupes: "<<tictoc()<<endl;
-#endif
 }
 
-#ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
-template void igl::cgal::remesh_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epeck, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -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<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epeck>, std::allocator<CGAL::Triangle_3<CGAL::Epeck> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > > > > > const&, std::map<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::less<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::allocator<std::pair<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> const, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, 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> >&);
-template void igl::cgal::remesh_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epeck, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epeck>, std::allocator<CGAL::Triangle_3<CGAL::Epeck> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > > > > > const&, std::map<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::less<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::allocator<std::pair<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> const, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template void igl::cgal::remesh_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epeck, Eigen::Matrix<double, -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<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epeck>, std::allocator<CGAL::Triangle_3<CGAL::Epeck> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > > > > > const&, std::map<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::less<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::allocator<std::pair<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> const, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, 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> >&);
-template void igl::cgal::remesh_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epeck, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epeck>, std::allocator<CGAL::Triangle_3<CGAL::Epeck> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > > > > > const&, std::map<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::less<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::allocator<std::pair<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> const, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template void igl::cgal::remesh_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epick, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -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<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > > > > > const&, std::map<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::less<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::allocator<std::pair<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> const, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, 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> >&);
-template void igl::cgal::remesh_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epick, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > > > > > const&, std::map<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::less<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::allocator<std::pair<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> const, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template void igl::cgal::remesh_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epick, Eigen::Matrix<double, -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<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > > > > > const&, std::map<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::less<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::allocator<std::pair<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> const, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, 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> >&);
-template void igl::cgal::remesh_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epick, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > > > > > const&, std::map<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::less<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::allocator<std::pair<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> const, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template void igl::cgal::remesh_intersections<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, CGAL::Epeck, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<CGAL::Triangle_3<CGAL::Epeck>, std::allocator<CGAL::Triangle_3<CGAL::Epeck> > > const&, std::map<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index const, std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > > > > > const&, std::map<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index>, std::vector<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, std::allocator<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index> >, std::less<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index> >, std::allocator<std::pair<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index> const, std::vector<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, std::allocator<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index> > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template void igl::cgal::remesh_intersections<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, CGAL::Epick, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > > const&, std::map<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index const, std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > > > > > const&, std::map<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index>, std::vector<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, std::allocator<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index> >, std::less<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index> >, std::allocator<std::pair<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index> const, std::vector<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, std::allocator<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index> > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template void igl::cgal::remesh_intersections<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epeck, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epeck>, std::allocator<CGAL::Triangle_3<CGAL::Epeck> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > > > > > const&, std::map<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::less<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::allocator<std::pair<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> const, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template void igl::cgal::remesh_intersections<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epick, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > > > > > const&, std::map<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::less<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::allocator<std::pair<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> const, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template void igl::cgal::remesh_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epeck, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epeck>, std::allocator<CGAL::Triangle_3<CGAL::Epeck> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > > > > > const&, std::map<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::less<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::allocator<std::pair<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> const, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
-template void igl::cgal::remesh_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epick, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<CGAL::Object, std::allocator<CGAL::Object> > > > > > const&, std::map<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::less<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> >, std::allocator<std::pair<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> const, std::vector<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::allocator<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index> > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
-#endif

include/igl/cgal/remesh_intersections.h

@@ -1,24 +1,18 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 // 
-// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
+// Copyright (C) 2015 Qingnan Zhou <qnzhou@gmail.com>
 // 
 // This Source Code Form is subject to the terms of the Mozilla Public License 
 // v. 2.0. If a copy of the MPL was not distributed with this file, You can 
 // obtain one at http://mozilla.org/MPL/2.0/.
+//
 #ifndef IGL_CGAL_REMESH_INTERSECTIONS_H
 #define IGL_CGAL_REMESH_INTERSECTIONS_H
-#include <igl/igl_inline.h>
 
+#include <igl/igl_inline.h>
 #include <Eigen/Dense>
 #include "CGAL_includes.hpp"
 
-#ifdef MEX
-#  include <mex.h>
-#  include <cassert>
-#  undef assert
-#  define assert( isOK ) ( (isOK) ? (void)0 : (void) mexErrMsgTxt(C_STR(__FILE__<<":"<<__LINE__<<": failed assertion `"<<#isOK<<"'"<<std::endl) ) )
-#endif
-  
 namespace igl
 {
   namespace cgal
@@ -72,6 +66,5 @@ namespace igl
 #ifndef IGL_STATIC_LIBRARY
 #  include "remesh_intersections.cpp"
 #endif
-  
-#endif
 
+#endif

include/igl/cgal/remesh_intersections_beta.cpp

@@ -1,338 +0,0 @@
-// This file is part of libigl, a simple c++ geometry processing library.
-// 
-// Copyright (C) 2015 Qingnan Zhou <qnzhou@gmail.com>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
-// obtain one at http://mozilla.org/MPL/2.0/.
-//
-#include "remesh_intersections_beta.h"
-
-#include <vector>
-#include <map>
-#include <unordered_map>
-#include <iostream>
-
-template <
-typename DerivedV,
-typename DerivedF,
-typename Kernel,
-typename DerivedVV,
-typename DerivedFF,
-typename DerivedJ,
-typename DerivedIM>
-IGL_INLINE void igl::cgal::remesh_intersections_beta(
-        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::pair<typename DerivedF::Index,
-        std::vector<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) {
-
-    typedef CGAL::Point_3<Kernel>    Point_3;
-    typedef CGAL::Segment_3<Kernel>  Segment_3; 
-    typedef CGAL::Triangle_3<Kernel> Triangle_3; 
-    typedef CGAL::Plane_3<Kernel>    Plane_3;
-    typedef CGAL::Point_2<Kernel>    Point_2;
-    typedef CGAL::Segment_2<Kernel>  Segment_2; 
-    typedef CGAL::Triangle_2<Kernel> Triangle_2; 
-    typedef CGAL::Triangulation_vertex_base_2<Kernel>  TVB_2;
-    typedef CGAL::Constrained_triangulation_face_base_2<Kernel> CTFB_2;
-    typedef CGAL::Triangulation_data_structure_2<TVB_2,CTFB_2> TDS_2;
-    typedef CGAL::Exact_intersections_tag Itag;
-    typedef CGAL::Constrained_Delaunay_triangulation_2<Kernel,TDS_2,Itag> 
-        CDT_2;
-    typedef CGAL::Constrained_triangulation_plus_2<CDT_2> CDT_plus_2;
-
-    typedef typename DerivedF::Index Index;
-    typedef std::pair<Index, Index> Edge;
-    struct EdgeHash {
-        size_t operator()(const Edge& e) const {
-            return (e.first * 805306457) ^ (e.second * 201326611);
-        }
-    };
-    typedef std::unordered_map<Edge, std::vector<Index>, EdgeHash > EdgeMap;
-
-    auto normalize_plane_coeff = [](const Plane_3& P) {
-        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) {
-        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;
-    };
-    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);
-        }
-    }
-
-    const size_t INVALID = std::numeric_limits<size_t>::max();
-    std::vector<std::vector<Index> > resolved_faces;
-    std::vector<Index> source_faces;
-    std::vector<Point_3> new_vertices;
-    EdgeMap edge_vertices;
-
-    /**
-     * Run constraint Delaunay triangulation on the plane.
-     */
-    auto run_delaunay_triangulation = [&](const Plane_3& P,
-            const std::vector<Index>& involved_faces,
-            std::vector<Point_3>& vertices,
-            std::vector<std::vector<Index> >& faces) {
-        CDT_plus_2 cdt;
-        for (const auto& fid : involved_faces) {
-            const auto itr = offending.find(fid);
-            const auto& triangle = T[fid];
-            cdt.insert_constraint(P.to_2d(triangle[0]), P.to_2d(triangle[1]));
-            cdt.insert_constraint(P.to_2d(triangle[1]), P.to_2d(triangle[2]));
-            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) {
-                if(const Segment_3 *iseg = CGAL::object_cast<Segment_3 >(&obj)) {
-                    // Add segment constraint
-                    cdt.insert_constraint(
-                            P.to_2d(iseg->vertex(0)),P.to_2d(iseg->vertex(1)));
-                }else if(const Point_3 *ipoint = CGAL::object_cast<Point_3 >(&obj)) {
-                    // Add point
-                    cdt.insert(P.to_2d(*ipoint));
-                } else if(const Triangle_3 *itri = CGAL::object_cast<Triangle_3 >(&obj)) {
-                    // Add 3 segment constraints
-                    cdt.insert_constraint(
-                            P.to_2d(itri->vertex(0)),P.to_2d(itri->vertex(1)));
-                    cdt.insert_constraint(
-                            P.to_2d(itri->vertex(1)),P.to_2d(itri->vertex(2)));
-                    cdt.insert_constraint(
-                            P.to_2d(itri->vertex(2)),P.to_2d(itri->vertex(0)));
-                } else if(const std::vector<Point_3 > *polyp = 
-                        CGAL::object_cast< std::vector<Point_3 > >(&obj)) {
-                    //cerr<<REDRUM("Poly...")<<endl;
-                    const std::vector<Point_3 > & poly = *polyp;
-                    const Index m = poly.size();
-                    assert(m>=2);
-                    for(Index p = 0;p<m;p++)
-                    {
-                        const Index np = (p+1)%m;
-                        cdt.insert_constraint(P.to_2d(poly[p]),P.to_2d(poly[np]));
-                    }
-                }else {
-                    throw std::runtime_error("Unknown intersection object!");
-                }
-            }
-        }
-        const size_t num_vertices = cdt.number_of_vertices();
-        const size_t num_faces = cdt.number_of_faces();
-        std::map<typename CDT_plus_2::Vertex_handle,Index> v2i;
-        size_t count=0;
-        for (auto itr = cdt.finite_vertices_begin();
-                itr != cdt.finite_vertices_end(); itr++) {
-            vertices.push_back(P.to_3d(itr->point()));
-            v2i[itr] = count;
-            count++;
-        }
-        for (auto itr = cdt.finite_faces_begin();
-                itr != cdt.finite_faces_end(); itr++) {
-            faces.push_back( {
-                    v2i[itr->vertex(0)],
-                    v2i[itr->vertex(1)],
-                    v2i[itr->vertex(2)] });
-        }
-    };
-
-    /**
-     * 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 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];
-        }
-
-        // 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;
-                }
-            }
-        }
-
-        // 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;
-    };
-
-    /**
-     * Determine the vertex indices for each corner of each output triangle.
-     */
-    auto post_triangulation_process = [&](
-            const std::vector<Point_3>& vertices,
-            const std::vector<std::vector<Index> >& faces,
-            const std::vector<Index>& involved_faces) {
-        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);
-                }
-            }
-        }
-    };
-
-    // Process un-touched faces.
-    for (size_t i=0; i<num_faces; i++) {
-        if (!is_offending[i]) {
-            resolved_faces.push_back(
-                    { F(i,0), F(i,1), F(i,2) } );
-            source_faces.push_back(i);
-        }
-    }
-
-    // Process self-intersecting 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);
-    for (size_t i=0; i<num_base_vertices; i++) {
-        assign_scalar(V(i,0), VV(i,0));
-        assign_scalar(V(i,1), VV(i,1));
-        assign_scalar(V(i,2), VV(i,2));
-    }
-    for (size_t i=num_base_vertices; i<num_out_vertices; i++) {
-        assign_scalar(new_vertices[i-num_base_vertices][0], VV(i,0));
-        assign_scalar(new_vertices[i-num_base_vertices][1], VV(i,1));
-        assign_scalar(new_vertices[i-num_base_vertices][2], VV(i,2));
-    }
-
-    const size_t num_out_faces = resolved_faces.size();
-    FF.resize(num_out_faces, 3);
-    for (size_t i=0; i<num_out_faces; i++) {
-        FF(i,0) = resolved_faces[i][0];
-        FF(i,1) = resolved_faces[i][1];
-        FF(i,2) = resolved_faces[i][2];
-    }
-
-    J.resize(num_out_faces);
-    std::copy(source_faces.begin(), source_faces.end(), J.data());
-
-    // Extract unique vertex indices.
-    IM.resize(VV.rows(),1);
-    std::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<VV.rows();v++) {
-        typename Kernel::FT p0,p1,p2;
-        assign_scalar(VV(v,0),p0);
-        assign_scalar(VV(v,1),p1);
-        assign_scalar(VV(v,2),p2);
-        const Point_3 p(p0,p1,p2);
-        if(vv2i.count(p)==0) {
-            vv2i[p] = v;
-        }
-        assert(vv2i.count(p) == 1);
-        IM(v) = vv2i[p];
-    }
-}
-

include/igl/cgal/remesh_intersections_beta.h

@@ -1,70 +0,0 @@
-// This file is part of libigl, a simple c++ geometry processing library.
-// 
-// Copyright (C) 2015 Qingnan Zhou <qnzhou@gmail.com>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
-// obtain one at http://mozilla.org/MPL/2.0/.
-//
-#ifndef IGL_CGAL_REMESH_INTERSECTIONS_BETA_H
-#define IGL_CGAL_REMESH_INTERSECTIONS_BETA_H
-
-#include <igl/igl_inline.h>
-#include <Eigen/Dense>
-#include "CGAL_includes.hpp"
-
-namespace igl
-{
-  namespace cgal
-  {
-    // Remesh faces according to results of intersection detection and
-    // construction (e.g. from `igl::cgal::intersect_other` or
-    // `igl::cgal::SelfIntersectMesh`)
-    //
-    // Inputs:
-    //   V  #V by 3 list of vertex positions
-    //   F  #F by 3 list of triangle indices into V
-    //   T  #F list of cgal triangles
-    //   offending #offending map taking face indices into F to pairs of order
-    //     of first finding and list of intersection objects from all
-    //     intersections
-    //   edge2faces  #edges <= #offending*3 to incident offending faces 
-    // Outputs:
-    //   VV  #VV by 3 list of vertex positions
-    //   FF  #FF by 3 list of triangle indices into V
-    //   IF  #intersecting face pairs by 2  list of intersecting face pairs,
-    //     indexing F
-    //   J  #FF list of indices into F denoting birth triangle
-    //   IM  #VV list of indices into VV of unique vertices.
-    //
-    template <
-      typename DerivedV,
-      typename DerivedF,
-      typename Kernel,
-      typename DerivedVV,
-      typename DerivedFF,
-      typename DerivedJ,
-      typename DerivedIM>
-    IGL_INLINE void remesh_intersections_beta(
-      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::pair<typename DerivedF::Index,
-          std::vector<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);
-  }
-}
-
-#ifndef IGL_STATIC_LIBRARY
-#  include "remesh_intersections_beta.cpp"
-#endif
-
-#endif