Fixed a bug (or reimplement) in resolving intersections, which would effect many things.

Former-commit-id: cf6347b1dcb75fd375d3fbb7dc9fa8c3903ed494

include/igl/boolean/mesh_boolean_beta.cpp

@@ -6,6 +6,9 @@
 
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 
+#include "../writePLY.h"
+#include "../writeDMAT.h"
+
 namespace igl {
     namespace boolean {
         namespace mesh_boolean_helper {
@@ -116,19 +119,25 @@ namespace igl {
                         continue;
                     }
                     if (counts[i] == 1) {
+                        bool found = false;
                         for (auto fid : uF2F[i]) {
                             if (fid > 0) {
                                 kept_faces.push_back(abs(fid)-1);
+                                found = true;
                                 break;
                             }
                         }
+                        assert(found);
                     } else if (counts[i] == -1) {
+                        bool found = false;
                         for (auto fid : uF2F[i]) {
                             if (fid < 0) {
                                 kept_faces.push_back(abs(fid)-1);
+                                found = true;
                                 break;
                             }
                         }
+                        assert(found);
                     } else {
                         assert(counts[i] == 0);
                     }
@@ -242,7 +251,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(V, F, labels, W);
+    igl::cgal::propagate_winding_numbers_beta(V, F, labels, W);
     assert(W.rows() == num_faces);
     if (W.cols() == 2) {
         assert(FB.rows() == 0);

include/igl/cgal/SelfIntersectMesh.h

@@ -225,6 +225,7 @@ namespace igl
 
 #include "mesh_to_cgal_triangle_list.h"
 #include "remesh_intersections.h"
+#include "remesh_intersections_beta.h"
 
 #include <igl/REDRUM.h>
 #include <igl/get_seconds.h>
@@ -404,7 +405,8 @@ inline igl::cgal::SelfIntersectMesh<
     return;
   }
 
-  remesh_intersections(V,F,T,offending,edge2faces,VV,FF,J,IM);
+  //remesh_intersections(V,F,T,offending,edge2faces,VV,FF,J,IM);
+  remesh_intersections_beta(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

@@ -25,7 +25,6 @@ IGL_INLINE size_t igl::cgal::extract_cells_single_component(
         const std::vector<std::vector<uE2EType> >& uE2E,
         const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
         Eigen::PlainObjectBase<DerivedC>& cells) {
-
     typedef typename DerivedF::Scalar Index;
     const size_t num_faces = F.rows();
     auto edge_index_to_face_index = [&](size_t index) {

include/igl/cgal/order_facets_around_edge.cpp

@@ -40,7 +40,7 @@ void igl::cgal::order_facets_around_edge(
     size_t s,
     size_t d, 
     const std::vector<int>& adj_faces,
-    Eigen::PlainObjectBase<DerivedI>& order)
+    Eigen::PlainObjectBase<DerivedI>& order, bool debug)
 {
     using namespace igl::cgal::order_facets_around_edges_helper;
 
@@ -196,10 +196,32 @@ void igl::cgal::order_facets_around_edge(
                 throw std::runtime_error("Unknown CGAL state detected.");
         }
     }
+    if (debug) {
+        std::cout << "tie positive: " << std::endl;
+        for (auto& f : tie_positive_oriented) {
+            std::cout << get_face_index(f) << " ";
+        }
+        std::cout << std::endl;
+        std::cout << "positive side: " << std::endl;
+        for (auto& f : positive_side) {
+            std::cout << get_face_index(f) << " ";
+        }
+        std::cout << std::endl;
+        std::cout << "tie negative: " << std::endl;
+        for (auto& f : tie_negative_oriented) {
+            std::cout << get_face_index(f) << " ";
+        }
+        std::cout << std::endl;
+        std::cout << "negative side: " << std::endl;
+        for (auto& f : negative_side) {
+            std::cout << get_face_index(f) << " ";
+        }
+        std::cout << std::endl;
+    }
 
     Eigen::PlainObjectBase<DerivedI> positive_order, negative_order;
-    order_facets_around_edge(V, F, s, d, positive_side, positive_order);
-    order_facets_around_edge(V, F, s, d, negative_side, negative_order);
+    order_facets_around_edge(V, F, s, d, positive_side, positive_order, debug);
+    order_facets_around_edge(V, F, s, d, negative_side, negative_order, debug);
     std::vector<size_t> tie_positive_order = index_sort(tie_positive_oriented);
     std::vector<size_t> tie_negative_order = index_sort(tie_negative_oriented);
 

include/igl/cgal/order_facets_around_edge.h

@@ -43,7 +43,8 @@ namespace igl {
                 const Eigen::PlainObjectBase<DerivedF>& F,
                 size_t s, size_t d, 
                 const std::vector<int>& adj_faces,
-                Eigen::PlainObjectBase<DerivedI>& order);
+                Eigen::PlainObjectBase<DerivedI>& order,
+                bool debug=false);
 
         // This funciton is a wrapper around the one above.  Since the ordering
         // is circular, the pivot point is used to define a starting point.  So

include/igl/cgal/propagate_winding_numbers.cpp

@@ -4,6 +4,7 @@
 #include "../facet_components.h"
 #include "../unique_edge_map.h"
 #include "../writeOBJ.h"
+#include "../writePLY.h"
 #include "order_facets_around_edge.h"
 #include "outer_facet.h"
 #include "closest_facet.h"
@@ -13,6 +14,7 @@
 #include <stdexcept>
 #include <limits>
 #include <vector>
+#include <tuple>
 
 namespace propagate_winding_numbers_helper {
     template<typename DerivedW >
@@ -51,6 +53,45 @@ namespace propagate_winding_numbers_helper {
         }
         return true;
     }
+
+    template<
+        typename DerivedF,
+        typename DeriveduE,
+        typename uE2EType >
+    bool is_orientable(
+            const Eigen::PlainObjectBase<DerivedF>& F,
+            const Eigen::PlainObjectBase<DeriveduE>& uE,
+            const std::vector<std::vector<uE2EType> >& uE2E) {
+        const size_t num_faces = F.rows();
+        const size_t num_edges = uE.rows();
+        auto edge_index_to_face_index = [&](size_t ei) {
+            return ei % num_faces;
+        };
+        auto is_consistent = [&](size_t fid, size_t s, size_t d) {
+            if (F(fid, 0) == s && F(fid, 1) == d) return true;
+            if (F(fid, 1) == s && F(fid, 2) == d) return true;
+            if (F(fid, 2) == s && F(fid, 0) == d) return true;
+
+            if (F(fid, 0) == d && F(fid, 1) == s) return false;
+            if (F(fid, 1) == d && F(fid, 2) == s) return false;
+            if (F(fid, 2) == d && F(fid, 0) == s) return false;
+            throw "Invalid face!!";
+        };
+        for (size_t i=0; i<num_edges; i++) {
+            const size_t s = uE(i,0);
+            const size_t d = uE(i,1);
+            int count=0;
+            for (const auto& ei : uE2E[i]) {
+                const size_t fid = edge_index_to_face_index(ei);
+                if (is_consistent(fid, s, d)) count++;
+                else count--;
+            }
+            if (count != 0) {
+                return false;
+            }
+        }
+        return true;
+    }
 }
 
 template<
@@ -773,3 +814,185 @@ IGL_INLINE void igl::cgal::propagate_winding_numbers(
     }
 }
 
+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;
+
+    Eigen::MatrixXi E, uE;
+    Eigen::VectorXi EMAP;
+    std::vector<std::vector<size_t> > uE2E;
+    igl::unique_edge_map(F, E, uE, EMAP, uE2E);
+    assert(propagate_winding_numbers_helper::is_orientable(F, uE, uE2E));
+
+    Eigen::VectorXi P;
+    const size_t num_patches = igl::extract_manifold_patches(F, EMAP, uE2E, P);
+
+    DerivedW per_patch_cells;
+    const size_t num_cells =
+        igl::cgal::extract_cells(V, F, P, E, uE, uE2E, EMAP, per_patch_cells);
+
+    typedef std::tuple<size_t, bool, size_t> CellConnection;
+    std::vector<std::set<CellConnection> > cell_adjacency(num_cells);
+    for (size_t i=0; i<num_patches; i++) {
+        const int positive_cell = per_patch_cells(i,0);
+        const int negative_cell = per_patch_cells(i,1);
+        cell_adjacency[positive_cell].emplace(negative_cell, false, i);
+        cell_adjacency[negative_cell].emplace(positive_cell, true, i);
+    }
+
+    {
+        auto save_cell = [&](const std::string& filename, size_t cell_id) {
+            std::vector<size_t> faces;
+            for (size_t i=0; i<num_patches; i++) {
+                if ((per_patch_cells.row(i).array() == cell_id).any()) {
+                    for (size_t j=0; j<num_faces; j++) {
+                        if (P[j] == i) {
+                            faces.push_back(j);
+                        }
+                    }
+                }
+            }
+            Eigen::MatrixXi cell_faces(faces.size(), 3);
+            for (size_t i=0; i<faces.size(); i++) {
+                cell_faces.row(i) = F.row(faces[i]);
+            }
+            Eigen::MatrixXd vertices(V.rows(), 3);
+            for (size_t i=0; i<V.rows(); i++) {
+                assign_scalar(V(i,0), vertices(i,0));
+                assign_scalar(V(i,1), vertices(i,1));
+                assign_scalar(V(i,2), vertices(i,2));
+            }
+            writePLY(filename, vertices, cell_faces);
+        };
+
+        // Check for odd cycle.
+        Eigen::VectorXi cell_labels(num_cells);
+        cell_labels.setZero();
+        Eigen::VectorXi parents(num_cells);
+        parents.setConstant(-1);
+        auto trace_parents = [&](size_t idx) {
+            std::list<size_t> path;
+            path.push_back(idx);
+            while (parents[path.back()] != path.back()) {
+                path.push_back(parents[path.back()]);
+            }
+            return path;
+        };
+        for (size_t i=0; i<num_cells; i++) {
+            if (cell_labels[i] == 0) {
+                cell_labels[i] = 1;
+                std::queue<size_t> Q;
+                Q.push(i);
+                parents[i] = i;
+                while (!Q.empty()) {
+                    size_t curr_idx = Q.front();
+                    Q.pop();
+                    int curr_label = cell_labels[curr_idx];
+                    for (const auto& neighbor : cell_adjacency[curr_idx]) {
+                        if (cell_labels[std::get<0>(neighbor)] == 0) {
+                            cell_labels[std::get<0>(neighbor)] = curr_label * -1;
+                            Q.push(std::get<0>(neighbor));
+                            parents[std::get<0>(neighbor)] = curr_idx;
+                        } else {
+                            if (cell_labels[std::get<0>(neighbor)] !=
+                                    curr_label * -1) {
+                                std::cerr << "Odd cell cycle detected!" << std::endl;
+                                auto path = trace_parents(curr_idx);
+                                path.reverse();
+                                auto path2 = trace_parents(std::get<0>(neighbor));
+                                path.insert(path.end(),
+                                        path2.begin(), path2.end());
+                                for (auto cell_id : path) {
+                                    std::cout << cell_id << " ";
+                                    std::stringstream filename;
+                                    filename << "cell_" << cell_id << ".ply";
+                                    save_cell(filename.str(), cell_id);
+                                }
+                                std::cout << std::endl;
+                            }
+                            assert(cell_labels[std::get<0>(neighbor)] == curr_label * -1);
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    size_t outer_facet;
+    bool flipped;
+    Eigen::VectorXi I;
+    I.setLinSpaced(num_faces, 0, num_faces-1);
+    igl::cgal::outer_facet(V, F, I, outer_facet, flipped);
+
+    const size_t outer_patch = P[outer_facet];
+    const size_t infinity_cell = per_patch_cells(outer_patch, flipped?1:0);
+
+    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 per_cell_W(num_cells, num_labels);
+    per_cell_W.setConstant(INVALID);
+    per_cell_W.row(infinity_cell).setZero();
+    std::queue<size_t> Q;
+    Q.push(infinity_cell);
+    while (!Q.empty()) {
+        size_t curr_cell = Q.front();
+        Q.pop();
+        for (const auto& neighbor : cell_adjacency[curr_cell]) {
+            size_t neighbor_cell, patch_idx;
+            bool direction;
+            std::tie(neighbor_cell, direction, patch_idx) = neighbor;
+            if ((per_cell_W.row(neighbor_cell).array() == INVALID).any()) {
+                per_cell_W.row(neighbor_cell) = per_cell_W.row(curr_cell);
+                for (size_t i=0; i<num_labels; i++) {
+                    int inc = (patch_labels[patch_idx] == i) ?
+                        (direction ? -1:1) :0;
+                    per_cell_W(neighbor_cell, i) =
+                        per_cell_W(curr_cell, i) + inc;
+                }
+            } else {
+                for (size_t i=0; i<num_labels; i++) {
+                    if (i == patch_labels[patch_idx]) {
+                        int inc = direction ? -1:1;
+                        assert(per_cell_W(neighbor_cell, i) ==
+                                per_cell_W(curr_cell, i) + inc);
+                    } else {
+                        assert(per_cell_W(neighbor_cell, i) ==
+                                per_cell_W(curr_cell, i));
+                    }
+                }
+            }
+        }
+    }
+
+    W.resize(num_faces, num_labels*2);
+    for (size_t i=0; i<num_faces; i++) {
+        const size_t patch = P[i];
+        const size_t positive_cell = per_patch_cells(patch, 0);
+        const size_t negative_cell = per_patch_cells(patch, 1);
+        for (size_t j=0; j<num_labels; j++) {
+            W(i,j*2  ) = per_cell_W(positive_cell, j);
+            W(i,j*2+1) = per_cell_W(negative_cell, j);
+        }
+    }
+}
+

include/igl/cgal/propagate_winding_numbers.h

@@ -163,6 +163,17 @@ 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_beta.cpp

@@ -0,0 +1,330 @@
+#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

@@ -0,0 +1,62 @@
+#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