refactor to extract cdt and assignment functions

Former-commit-id: adeb669961430aae71968356c4f089ae1a596cca

include/igl/cgal/SelfIntersectMesh.h

@@ -95,12 +95,7 @@ namespace igl
         IndexList lIF;
         // #F-long list of faces with intersections mapping to the order in
         // which they were first found
-        struct IDObjectList
-        {
-          Index id;
-          ObjectList objects;
-        };
-        std::map<Index,IDObjectList> offending;
+        std::map<Index,std::pair<Index,ObjectList> > offending;
         // Make a short name for the edge map's key
         typedef std::pair<Index,Index> EMK;
         // Make a short name for the type stored at each edge, the edge map's
@@ -215,22 +210,13 @@ namespace igl
         //   A_objects_3  updated list of intersection objects for A
         // Outputs:
         //   cdt  Contrained delaunay triangulation in projected 2D plane
-        static inline void projected_delaunay(
-            const Triangle_3 & A,
-            const ObjectList & A_objects_3,
-            CDT_plus_2 & cdt);
-        // Getters:
       public:
+        // Getters:
         //const IndexList& get_lIF() const{ return lIF;}
         static inline void box_intersect_static(
           SelfIntersectMesh * SIM, 
           const Box &a, 
           const Box &b);
-        // Annoying wrappers to conver from cgal to double or cgal
-        static inline void to_output_type(const typename Kernel::FT & cgal,double & d);
-        static inline void to_output_type(
-          const typename CGAL::Epeck::FT & cgal,
-          CGAL::Epeck::FT & d);
     };
   }
 }
@@ -238,6 +224,7 @@ namespace igl
 // Implementation
 
 #include "mesh_to_cgal_triangle_list.h"
+#include "remesh_intersections.h"
 
 #include <igl/REDRUM.h>
 #include <igl/get_seconds.h>
@@ -417,262 +404,7 @@ inline igl::cgal::SelfIntersectMesh<
     return;
   }
 
-  const auto & remesh = [](
-    const Eigen::PlainObjectBase<DerivedV> & V,
-    const Eigen::PlainObjectBase<DerivedF> & F,
-    const Triangles & T,
-    const std::map<Index,IDObjectList> & offending,
-    const EdgeMap & edge2faces,
-    Eigen::PlainObjectBase<DerivedVV> & VV,
-    Eigen::PlainObjectBase<DerivedFF> & FF,
-    Eigen::PlainObjectBase<DerivedJ> & J,
-    Eigen::PlainObjectBase<DerivedIM> & IM
-    )
-  {
-    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;
-    // Loop over offending triangles
-    const size_t noff = offending.size();
-#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...
-  //# 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.id;
-      {
-#ifdef IGL_SELFINTERSECTMESH_DEBUG
-        const double t_before = get_seconds();
-#endif
-        projected_delaunay(T[f],okv.second.objects,cdt[o]);
-#ifdef IGL_SELFINTERSECTMESH_DEBUG
-        t_proj_del += (get_seconds()-t_before);
-#endif
-      }
-      // 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);
-              }
-              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.id;
-                // 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;
-                  }
-                }
-              }
-            }
-            if(!found)
-            {
-  //#   pragma omp critical
-              {
-                v2i[vit] = V.rows()+NV_count;
-                NV.push_back(vit_point_3);
-                NV_count++;
-              }
-            }
-          }
-          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++;
-        }
-      }
-    }
-#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++;
-      }
-    }
-    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.id;
-      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
-          to_output_type(p[d], VV(V.rows()+i,d));
-        }
-        i++;
-      }
-    }
-    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++)
-    {
-      const Point_3 p(V(v,0),V(v,1),V(v,2));
-      if(vv2i.count(p)==0)
-      {
-        vv2i[p] = v;
-      }
-      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++;
-      }
-    }
-#ifdef IGL_SELFINTERSECTMESH_DEBUG
-    cout<<"Output + dupes: "<<tictoc()<<endl;
-#endif
-  };
-  remesh(
-    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.
@@ -712,24 +444,13 @@ inline void igl::cgal::SelfIntersectMesh<
   if(offending.count(f) == 0)
   {
     // first time marking, initialize with new id and empty list
-    offending[f] = {(Index)offending.size(),{}};
-    // Add to edge map
+    const Index id = offending.size();
+    offending[f] = {id,{}};
     for(Index e = 0; e<3;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);
-      }
-      // Create new list if there is no entry
-      if(edge2faces.count(EMK(i,j))==0)
-      {
-        edge2faces[EMK(i,j)] = EMV();
-      }
-      // append to list
+      // append face to edge's list
+      Index i = F(f,(e+1)%3) < F(f,(e+2)%3) ? F(f,(e+1)%3) : F(f,(e+2)%3);
+      Index j = F(f,(e+1)%3) < F(f,(e+2)%3) ? F(f,(e+2)%3) : F(f,(e+1)%3);
       edge2faces[EMK(i,j)].push_back(f);
     }
   }
@@ -799,8 +520,8 @@ inline bool igl::cgal::SelfIntersectMesh<
   {
     // Construct intersection
     CGAL::Object result = CGAL::intersection(A,B);
-    offending[fa].objects.push_back(result);
-    offending[fb].objects.push_back(result);
+    offending[fa].second.push_back(result);
+    offending[fb].second.push_back(result);
   }
   return true;
 }
@@ -898,8 +619,8 @@ inline bool igl::cgal::SelfIntersectMesh<
         A.vertex(va),
         *p));
       count_intersection(fa,fb);
-      offending[fa].objects.push_back(seg);
-      offending[fb].objects.push_back(seg);
+      offending[fa].second.push_back(seg);
+      offending[fb].second.push_back(seg);
       return true;
     }else if(CGAL::object_cast<Segment_3 >(&result))
     {
@@ -1045,8 +766,8 @@ inline bool igl::cgal::SelfIntersectMesh<
       } else
       {
         // Triangle object
-        offending[fa].objects.push_back(result);
-        offending[fb].objects.push_back(result);
+        offending[fa].second.push_back(result);
+        offending[fb].second.push_back(result);
         //cerr<<REDRUM("Coplanar at: "<<fa<<" & "<<fb<<" (double shared).")<<endl;
         return true;
       }
@@ -1213,142 +934,5 @@ done:
   return;
 }
 
-// Compute 2D delaunay triangulation of a given 3d triangle and a list of
-// intersection objects (points,segments,triangles). CGAL uses an affine
-// projection rather than an isometric projection, so we're not guaranteed
-// that the 2D delaunay triangulation here will be a delaunay triangulation
-// in 3D.
-//
-// Inputs:
-//   A  triangle in 3D
-//   A_objects_3  updated list of intersection objects for A
-// Outputs:
-//   cdt  Contrained delaunay triangulation in projected 2D plane
-template <
-  typename Kernel,
-  typename DerivedV,
-  typename DerivedF,
-  typename DerivedVV,
-  typename DerivedFF,
-  typename DerivedIF,
-  typename DerivedJ,
-  typename DerivedIM>
-inline void igl::cgal::SelfIntersectMesh<
-  Kernel,
-  DerivedV,
-  DerivedF,
-  DerivedVV,
-  DerivedFF,
-  DerivedIF,
-  DerivedJ,
-  DerivedIM>::projected_delaunay(
-  const Triangle_3 & A,
-  const ObjectList & A_objects_3,
-  CDT_plus_2 & cdt)
-{
-  using namespace std;
-  // http://www.cgal.org/Manual/3.2/doc_html/cgal_manual/Triangulation_2/Chapter_main.html#Section_2D_Triangulations_Constrained_Plus
-  // Plane of triangle A
-  Plane_3 P(A.vertex(0),A.vertex(1),A.vertex(2));
-  // Insert triangle into vertices
-  typename CDT_plus_2::Vertex_handle corners[3];
-  for(Index i = 0;i<3;i++)
-  {
-    corners[i] = cdt.insert(P.to_2d(A.vertex(i)));
-  }
-  // Insert triangle edges as constraints
-  for(Index i = 0;i<3;i++)
-  {
-    cdt.insert_constraint( corners[(i+1)%3], corners[(i+2)%3]);
-  }
-  // Insert constraints for intersection objects
-  for( const auto & obj : A_objects_3)
-  {
-    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
-    {
-      cerr<<REDRUM("What is this object?!")<<endl;
-      assert(false);
-    }
-  }
-}
-
-template <
-  typename Kernel,
-  typename DerivedV,
-  typename DerivedF,
-  typename DerivedVV,
-  typename DerivedFF,
-  typename DerivedIF,
-  typename DerivedJ,
-  typename DerivedIM>
-inline 
-void 
-igl::cgal::SelfIntersectMesh<
-  Kernel,
-  DerivedV,
-  DerivedF,
-  DerivedVV,
-  DerivedFF,
-  DerivedIF,
-  DerivedJ,
-  DerivedIM>::to_output_type(
-    const typename Kernel::FT & cgal,
-    double & d)
-{
-  d = CGAL::to_double(cgal);
-}
-
-template <
-  typename Kernel,
-  typename DerivedV,
-  typename DerivedF,
-  typename DerivedVV,
-  typename DerivedFF,
-  typename DerivedIF,
-  typename DerivedJ,
-  typename DerivedIM>
-inline 
-void 
-igl::cgal::SelfIntersectMesh<
-  Kernel,
-  DerivedV,
-  DerivedF,
-  DerivedVV,
-  DerivedFF,
-  DerivedIF,
-  DerivedJ,
-  DerivedIM>::to_output_type(
-    const typename CGAL::Epeck::FT & cgal,
-    CGAL::Epeck::FT & d)
-{
-  d = cgal;
-}
-
 #endif
 

include/igl/cgal/assign_scalar.cpp

@@ -0,0 +1,15 @@
+#include "assign_scalar.h"
+
+IGL_INLINE void igl::cgal::assign_scalar(
+  const typename CGAL::Epeck::FT & cgal,
+  CGAL::Epeck::FT & d)
+{
+  d = cgal;
+}
+
+IGL_INLINE void igl::cgal::assign_scalar(
+  const typename CGAL::Epeck::FT & cgal,
+  double & d)
+{
+  d = CGAL::to_double(cgal);
+}

include/igl/cgal/assign_scalar.h

@@ -0,0 +1,24 @@
+#ifndef IGL_CGAL_ASSIGN_SCALAR_H
+#define IGL_CGAL_ASSIGN_SCALAR_H
+#include "../igl_inline.h"
+#include "CGAL_includes.hpp"
+namespace igl
+{
+  namespace cgal
+  {
+    // Inputs:
+    //   cgal  cgal scalar
+    // Outputs:
+    //   d  output scalar
+    IGL_INLINE void assign_scalar(
+      const typename CGAL::Epeck::FT & cgal,
+      CGAL::Epeck::FT & d);
+    IGL_INLINE void assign_scalar(
+      const typename CGAL::Epeck::FT & cgal,
+      double & d);
+  }
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "assign_scalar.cpp"
+#endif
+#endif

include/igl/cgal/projected_delaunay.cpp

@@ -0,0 +1,89 @@
+#include "projected_delaunay.h"
+#include "../REDRUM.h"
+#include <iostream>
+#include <cassert>
+
+template <typename Kernel>
+IGL_INLINE void igl::cgal::projected_delaunay(
+  const CGAL::Triangle_3<Kernel> & A,
+  const std::vector<CGAL::Object> & A_objects_3,
+  CGAL::Constrained_triangulation_plus_2<
+    CGAL::Constrained_Delaunay_triangulation_2<
+      Kernel,
+      CGAL::Triangulation_data_structure_2<
+        CGAL::Triangulation_vertex_base_2<Kernel>,
+        CGAL::Constrained_triangulation_face_base_2<Kernel> >,
+      CGAL::Exact_intersections_tag> > & cdt)
+{
+  using namespace std;
+  // 3D Primitives
+  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::Tetrahedron_3<Kernel> Tetrahedron_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;
+
+  // http://www.cgal.org/Manual/3.2/doc_html/cgal_manual/Triangulation_2/Chapter_main.html#Section_2D_Triangulations_Constrained_Plus
+  // Plane of triangle A
+  Plane_3 P(A.vertex(0),A.vertex(1),A.vertex(2));
+  // Insert triangle into vertices
+  typename CDT_plus_2::Vertex_handle corners[3];
+  typedef size_t Index;
+  for(Index i = 0;i<3;i++)
+  {
+    const Point_3 & p3 = A.vertex(i);
+    const Point_2 & p2 = P.to_2d(p3);
+    typename CDT_plus_2::Vertex_handle corner = cdt.insert(p2);
+    corners[i] = corner;
+  }
+  // Insert triangle edges as constraints
+  for(Index i = 0;i<3;i++)
+  {
+    cdt.insert_constraint( corners[(i+1)%3], corners[(i+2)%3]);
+  }
+  // Insert constraints for intersection objects
+  for( const auto & obj : A_objects_3)
+  {
+    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
+    {
+      cerr<<REDRUM("What is this object?!")<<endl;
+      assert(false);
+    }
+  }
+}

include/igl/cgal/projected_delaunay.h

@@ -0,0 +1,36 @@
+#ifndef IGL_CGAL_PROJECTED_DELAUNAY_H
+#define IGL_CGAL_PROJECTED_DELAUNAY_H
+#include "../igl_inline.h"
+#include "CGAL_includes.hpp"
+namespace igl
+{
+  namespace cgal
+  {
+    // Compute 2D delaunay triangulation of a given 3d triangle and a list of
+    // intersection objects (points,segments,triangles). CGAL uses an affine
+    // projection rather than an isometric projection, so we're not guaranteed
+    // that the 2D delaunay triangulation here will be a delaunay triangulation
+    // in 3D.
+    //
+    // Inputs:
+    //   A  triangle in 3D
+    //   A_objects_3  updated list of intersection objects for A
+    // Outputs:
+    //   cdt  Contrained delaunay triangulation in projected 2D plane
+    template <typename Kernel>
+    IGL_INLINE void projected_delaunay(
+      const CGAL::Triangle_3<Kernel> & A,
+      const std::vector<CGAL::Object> & A_objects_3,
+      CGAL::Constrained_triangulation_plus_2<
+        CGAL::Constrained_Delaunay_triangulation_2<
+          Kernel,
+          CGAL::Triangulation_data_structure_2<
+            CGAL::Triangulation_vertex_base_2<Kernel>,
+            CGAL::Constrained_triangulation_face_base_2<Kernel> >,
+          CGAL::Exact_intersections_tag> > & cdt);
+  }
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "projected_delaunay.cpp"
+#endif
+#endif

include/igl/cgal/remesh_intersections.cpp

@@ -0,0 +1,299 @@
+#include "remesh_intersections.h"
+#include "SelfIntersectMesh.h"
+#include "assign_scalar.h"
+#include "projected_delaunay.h"
+#include <iostream>
+#include <cassert>
+
+template <
+  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;
+
+  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;
+  // Loop over offending triangles
+  const size_t noff = offending.size();
+#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...
+//# 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
+    }
+    // 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);
+            }
+            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;
+                }
+              }
+            }
+          }
+          if(!found)
+          {
+//#   pragma omp critical
+            {
+              v2i[vit] = V.rows()+NV_count;
+              NV.push_back(vit_point_3);
+              NV_count++;
+            }
+          }
+        }
+        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++;
+      }
+    }
+  }
+#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++;
+    }
+  }
+  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++;
+    }
+  }
+  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++)
+  {
+    const Point_3 p(V(v,0),V(v,1),V(v,2));
+    if(vv2i.count(p)==0)
+    {
+      vv2i[p] = v;
+    }
+    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++;
+    }
+  }
+#ifdef IGL_SELFINTERSECTMESH_DEBUG
+  cout<<"Output + dupes: "<<tictoc()<<endl;
+#endif
+}

include/igl/cgal/remesh_intersections.h

@@ -0,0 +1,77 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Alec Jacobson <alecjacobson@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 <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
+  {
+    // 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(
+      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.cpp"
+#endif
+  
+#endif
+