Merge branch 'master' of github.com:libigl/libigl

Former-commit-id: e99f4b3ab4771c26250399c0ef39a4eec83ebc55

include/igl/bounding_box.cpp

@@ -83,4 +83,5 @@ IGL_INLINE void igl::bounding_box(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 template void igl::bounding_box<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
+template void igl::bounding_box<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -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<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/copyleft/boolean/mesh_boolean.cpp

@@ -344,7 +344,7 @@ IGL_INLINE bool igl::copyleft::boolean::mesh_boolean(
         Eigen::PlainObjectBase<DerivedJ>& J) {
       Eigen::VectorXi I;
       igl::copyleft::cgal::RemeshSelfIntersectionsParam params;
-
+      params.stitch_all = true;
       MatrixXES Vr;
       DerivedFC Fr;
       Eigen::MatrixXi IF;

include/igl/copyleft/cgal/RemeshSelfIntersectionsParam.h

@@ -20,9 +20,14 @@ namespace igl
       {
         bool detect_only;
         bool first_only;
-        RemeshSelfIntersectionsParam():detect_only(false),first_only(false){};
-        RemeshSelfIntersectionsParam(bool _detect_only, bool _first_only):
-          detect_only(_detect_only),first_only(_first_only){};
+        bool stitch_all;
+        RemeshSelfIntersectionsParam(
+          bool _detect_only=false, 
+          bool _first_only=false,
+          bool _stitch_all=false):
+          detect_only(_detect_only),
+          first_only(_first_only),
+          stitch_all(_stitch_all){};
       };
     }
   }

include/igl/copyleft/cgal/SelfIntersectMesh.h

@@ -236,7 +236,6 @@ namespace igl
 
 #include "mesh_to_cgal_triangle_list.h"
 #include "remesh_intersections.h"
-#include "remesh_intersections.h"
 
 #include "../../REDRUM.h"
 #include "../../get_seconds.h"
@@ -427,7 +426,8 @@ inline igl::copyleft::cgal::SelfIntersectMesh<
     return;
   }
 
-  remesh_intersections(V,F,T,offending,edge2faces,true,VV,FF,J,IM);
+  remesh_intersections(
+    V,F,T,offending,edge2faces,params.stitch_all,VV,FF,J,IM);
 
 #ifdef IGL_SELFINTERSECTMESH_DEBUG
   log_time("remesh_intersection");

include/igl/copyleft/cgal/cell_adjacency.cpp

@@ -27,5 +27,5 @@ IGL_INLINE void igl::copyleft::cgal::cell_adjacency(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-template void igl::copyleft::cgal::cell_adjacency<Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, unsigned long, std::__1::vector<std::__1::set<std::__1::tuple<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, bool, unsigned long>, std::__1::less<std::__1::tuple<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, bool, unsigned long> >, std::__1::allocator<std::__1::tuple<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, bool, unsigned long> > >, std::__1::allocator<std::__1::set<std::__1::tuple<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, bool, unsigned long>, std::__1::less<std::__1::tuple<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, bool, unsigned long> >, std::__1::allocator<std::__1::tuple<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, bool, unsigned long> > > > >&);
+template void igl::copyleft::cgal::cell_adjacency<Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, unsigned long, std::vector<std::set<std::tuple<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, bool, unsigned long>, std::less<std::tuple<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, bool, unsigned long> >, std::allocator<std::tuple<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, bool, unsigned long> > >, std::allocator<std::set<std::tuple<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, bool, unsigned long>, std::less<std::tuple<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, bool, unsigned long> >, std::allocator<std::tuple<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, bool, unsigned long> > > > >&);
 #endif

include/igl/copyleft/cgal/outer_hull.cpp

@@ -47,6 +47,7 @@ IGL_INLINE void igl::copyleft::cgal::outer_hull(
   DerivedJ Jr;
   {
     RemeshSelfIntersectionsParam params;
+    params.stitch_all = true;
     Eigen::VectorXi I;
     Eigen::MatrixXi IF;
     remesh_self_intersections(V, F, params, Vr, Fr, IF, Jr, I);

include/igl/copyleft/cgal/piecewise_constant_winding_number.cpp

@@ -22,8 +22,6 @@ IGL_INLINE bool igl::copyleft::cgal::piecewise_constant_winding_number(
   Eigen::Matrix<typename DerivedF::Index,Eigen::Dynamic,1> J;
   Eigen::Matrix<typename DerivedV::Index,Eigen::Dynamic,1> UIM,IM;
   // resolve intersections
-  remesh_self_intersections(V,F,{},VV,FF,IF,J,IM);
-  // _apply_ duplicate vertex mapping IM to FF
-  std::for_each(FF.data(),FF.data()+FF.size(),[&IM](int & a){a=IM(a);});
+  remesh_self_intersections(V,F,{false,false,true},VV,FF,IF,J,IM);
   return igl::piecewise_constant_winding_number(FF);
 }

include/igl/copyleft/cgal/remesh_intersections.cpp

@@ -20,56 +20,59 @@
 //#define REMESH_INTERSECTIONS_TIMING
 
 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::copyleft::cgal::remesh_intersections(
-        const Eigen::PlainObjectBase<DerivedV> & V,
-        const Eigen::PlainObjectBase<DerivedF> & F,
-        const std::vector<CGAL::Triangle_3<Kernel> > & T,
-        const std::map<
-        typename DerivedF::Index,
-        std::vector<
-        std::pair<typename DerivedF::Index, CGAL::Object> > > & offending,
-        const std::map<
-        std::pair<typename DerivedF::Index,typename DerivedF::Index>,
-        std::vector<typename DerivedF::Index> > & edge2faces,
-        Eigen::PlainObjectBase<DerivedVV> & VV,
-        Eigen::PlainObjectBase<DerivedFF> & FF,
-        Eigen::PlainObjectBase<DerivedJ> & J,
-        Eigen::PlainObjectBase<DerivedIM> & IM) {
-    igl::copyleft::cgal::remesh_intersections(
-            V, F, T, offending, edge2faces, false, VV, FF, J, IM);
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  const std::vector<CGAL::Triangle_3<Kernel> > & T,
+  const std::map<
+  typename DerivedF::Index,
+  std::vector<
+  std::pair<typename DerivedF::Index, CGAL::Object> > > & offending,
+  const std::map<
+  std::pair<typename DerivedF::Index,typename DerivedF::Index>,
+  std::vector<typename DerivedF::Index> > & edge2faces,
+  Eigen::PlainObjectBase<DerivedVV> & VV,
+  Eigen::PlainObjectBase<DerivedFF> & FF,
+  Eigen::PlainObjectBase<DerivedJ> & J,
+  Eigen::PlainObjectBase<DerivedIM> & IM) 
+{
+  // by default, no stitching
+  igl::copyleft::cgal::remesh_intersections(
+    V, F, T, offending, edge2faces, false, VV, FF, J, IM);
 }
 
 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::copyleft::cgal::remesh_intersections(
-        const Eigen::PlainObjectBase<DerivedV> & V,
-        const Eigen::PlainObjectBase<DerivedF> & F,
-        const std::vector<CGAL::Triangle_3<Kernel> > & T,
-        const std::map<
-        typename DerivedF::Index,
-        std::vector<
-        std::pair<typename DerivedF::Index, CGAL::Object> > > & offending,
-        const std::map<
-        std::pair<typename DerivedF::Index,typename DerivedF::Index>,
-        std::vector<typename DerivedF::Index> > & /*edge2faces*/,
-        bool stitch_all,
-        Eigen::PlainObjectBase<DerivedVV> & VV,
-        Eigen::PlainObjectBase<DerivedFF> & FF,
-        Eigen::PlainObjectBase<DerivedJ> & J,
-        Eigen::PlainObjectBase<DerivedIM> & IM) {
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  const std::vector<CGAL::Triangle_3<Kernel> > & T,
+  const std::map<
+    typename DerivedF::Index,
+    std::vector<
+      std::pair<typename DerivedF::Index, CGAL::Object> > > & offending,
+  const std::map<
+    std::pair<typename DerivedF::Index,typename DerivedF::Index>,
+    std::vector<typename DerivedF::Index> > & /*edge2faces*/,
+  bool stitch_all,
+  Eigen::PlainObjectBase<DerivedVV> & VV,
+  Eigen::PlainObjectBase<DerivedFF> & FF,
+  Eigen::PlainObjectBase<DerivedJ> & J,
+  Eigen::PlainObjectBase<DerivedIM> & IM)
+{
 
 #ifdef REMESH_INTERSECTIONS_TIMING
     const auto & tictoc = []() -> double
@@ -112,28 +115,35 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
     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;
+    for (const auto itr : offending)
+    {
+      const auto& fid = itr.first;
+      is_offending[fid] = true;
     }
 
+    // Cluster overlaps so that co-planar clusters are resolved only once
     std::unordered_map<Index, std::vector<Index> > intersecting_and_coplanar;
-    for (const auto itr : offending) {
-        const auto& fi = itr.first;
-        const auto P = T[fi].supporting_plane();
-        assert(!P.is_degenerate());
-        for (const auto jtr : itr.second) {
-            const auto& fj = jtr.first;
-            const auto& tj = T[fj];
-            if (P.has_on(tj[0]) && P.has_on(tj[1]) && P.has_on(tj[2])) {
-                auto loc = intersecting_and_coplanar.find(fi);
-                if (loc == intersecting_and_coplanar.end()) {
-                    intersecting_and_coplanar[fi] = {fj};
-                } else {
-                    loc->second.push_back(fj);
-                }
-            }
+    for (const auto itr : offending)
+    {
+      const auto& fi = itr.first;
+      const auto P = T[fi].supporting_plane();
+      assert(!P.is_degenerate());
+      for (const auto jtr : itr.second) 
+      {
+        const auto& fj = jtr.first;
+        const auto& tj = T[fj];
+        if (P.has_on(tj[0]) && P.has_on(tj[1]) && P.has_on(tj[2]))
+        {
+          auto loc = intersecting_and_coplanar.find(fi);
+          if (loc == intersecting_and_coplanar.end())
+          {
+            intersecting_and_coplanar[fi] = {fj};
+          } else
+          {
+            loc->second.push_back(fj);
+          }
         }
+      }
     }
 #ifdef REMESH_INTERSECTIONS_TIMING
     log_time("overlap_analysis");
@@ -144,92 +154,121 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
     std::vector<Point_3> new_vertices;
     EdgeMap edge_vertices;
 
-    /**
-     * Run constraint Delaunay triangulation on the plane.
-     */
-    auto run_delaunay_triangulation = [&offending, &T](const Plane_3& P,
-            const std::vector<Index>& involved_faces,
-            std::vector<Point_3>& vertices,
-            std::vector<std::vector<Index> >& faces) -> void {
-        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& index_obj : itr->second) {
-                //const auto& ofid = index_obj.first;
-                const auto& obj = index_obj.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!");
+    // Run constraint Delaunay triangulation on the plane.
+    // 
+    // Inputs:
+    //   P  plane to triangulate upone
+    //   involved_faces  #F list of indices into triangle of involved faces
+    // Outputs:
+    //   vertices  #V list of vertex positions of output triangulation
+    //   faces   #F list of face indices into vertices of output triangulation
+    //
+    auto run_delaunay_triangulation = [&offending, &T](
+      const Plane_3& P,
+      const std::vector<Index>& involved_faces,
+      std::vector<Point_3>& vertices,
+      std::vector<std::vector<Index> >& faces) -> void 
+    {
+      CDT_plus_2 cdt;
+      // insert each face into a common 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& index_obj : itr->second) 
+        {
+            //const auto& ofid = index_obj.first;
+            const auto& obj = index_obj.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!");
             }
         }
-        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)] });
-        }
+      }
+      // Read off vertices of the cdt, remembering index
+      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++;
+      }
+      // Read off faces and store index triples
+      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 {
-      if (stitch_all) {
+    // Given p on triangle indexed by ori_f, add point to list of vertices return index of p.
+    //
+    // Input:
+    //   p  point to search for
+    //   ori_f  index of triangle p is corner of
+    // Returns global index of vertex (dependent on whether stitch_all flag is
+    // set)
+    //
+    auto find_or_append_point = [&](
+      const Point_3& p, 
+      const size_t ori_f) -> Index 
+    {
+      if (stitch_all) 
+      {
         // No need to check if p shared by multiple triangles because all shared
         // vertices would be merged later on.
         const size_t index = num_base_vertices + new_vertices.size();
         new_vertices.push_back(p);
         return index;
-      } else {
+      } else 
+      {
         // Stitching triangles according to input connectivity.
         // This step is potentially costly.
         const auto& triangle = T[ori_f];
         const auto& f = F.row(ori_f).eval();
 
         // Check if p is one of the triangle corners.
-        for (size_t i=0; i<3; i++) {
+        for (size_t i=0; i<3; i++) 
+        {
           if (p == triangle[i]) return f[i];
         }
 
@@ -272,41 +311,54 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
       }
     };
 
-    /**
-     * Determine the vertex indices for each corner of each output triangle.
-     */
+    // Determine the vertex indices for each corner of each output triangle.
+    // 
+    // Inputs:
+    //   vertices  #V list of vertices of cdt
+    //   faces  #F list of list of face indices into vertices of cdt
+    //   involved_faces  list of involved faces on the plane of cdt
+    // Side effects:
+    //   - add faces to resolved_faces
+    //   - add corresponding original face to source_faces
+    //   - 
     auto post_triangulation_process = [&](
-            const std::vector<Point_3>& vertices,
-            const std::vector<std::vector<Index> >& faces,
-            const std::vector<Index>& involved_faces) -> void {
+      const std::vector<Point_3>& vertices,
+      const std::vector<std::vector<Index> >& faces,
+      const std::vector<Index>& involved_faces) -> void 
+    {
       assert(involved_faces.size() > 0);
-      for (const auto& f : faces) {
+      // for all faces of the cdt
+      for (const auto& f : faces) 
+      {
         const Point_3& v0 = vertices[f[0]];
         const Point_3& v1 = vertices[f[1]];
         const Point_3& v2 = vertices[f[2]];
         Point_3 center(
-            (v0[0] + v1[0] + v2[0]) / 3.0,
-            (v0[1] + v1[1] + v2[1]) / 3.0,
-            (v0[2] + v1[2] + v2[2]) / 3.0);
-        if (involved_faces.size() == 1) {
+          (v0[0] + v1[0] + v2[0]) / 3.0,
+          (v0[1] + v1[1] + v2[1]) / 3.0,
+          (v0[2] + v1[2] + v2[2]) / 3.0);
+        if (involved_faces.size() == 1) 
+        {
           // If only there is only one involved face, all sub-triangles must
           // belong to it and have the correct orientation.
           const auto& ori_f = involved_faces[0];
           std::vector<Index> corners(3);
-          corners[0] = determine_point_index(v0, ori_f);
-          corners[1] = determine_point_index(v1, ori_f);
-          corners[2] = determine_point_index(v2, ori_f);
+          corners[0] = find_or_append_point(v0, ori_f);
+          corners[1] = find_or_append_point(v1, ori_f);
+          corners[2] = find_or_append_point(v2, ori_f);
           resolved_faces.emplace_back(corners);
           source_faces.push_back(ori_f);
-        } else {
-          for (const auto& ori_f : involved_faces) {
+        } else 
+        {
+          for (const auto& ori_f : involved_faces) 
+          {
             const auto& triangle = T[ori_f];
             const Plane_3 P = triangle.supporting_plane();
             if (triangle.has_on(center)) {
               std::vector<Index> corners(3);
-              corners[0] = determine_point_index(v0, ori_f);
-              corners[1] = determine_point_index(v1, ori_f);
-              corners[2] = determine_point_index(v2, ori_f);
+              corners[0] = find_or_append_point(v0, ori_f);
+              corners[1] = find_or_append_point(v1, ori_f);
+              corners[2] = find_or_append_point(v2, ori_f);
               if (CGAL::orientation(
                     P.to_2d(v0), P.to_2d(v1), P.to_2d(v2))
                   == CGAL::RIGHT_TURN) {
@@ -321,48 +373,53 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
     };
 
     // Process un-touched faces.
-    for (size_t i=0; i<num_faces; i++) {
-        if (!is_offending[i] && !T[i].is_degenerate()) {
-            resolved_faces.push_back(
-                    { F(i,0), F(i,1), F(i,2) } );
-            source_faces.push_back(i);
-        }
+    for (size_t i=0; i<num_faces; i++) 
+    {
+      if (!is_offending[i] && !T[i].is_degenerate()) 
+      {
+        resolved_faces.push_back( { F(i,0), F(i,1), F(i,2) } );
+        source_faces.push_back(i);
+      }
     }
 
     // Process self-intersecting faces.
     std::vector<bool> processed(num_faces, false);
     std::vector<std::pair<Plane_3, std::vector<Index> > > cdt_inputs;
-    for (const auto itr : offending) {
-        const auto fid = itr.first;
-        if (processed[fid]) continue;
-        processed[fid] = true;
-
-        const auto loc = intersecting_and_coplanar.find(fid);
-        std::vector<Index> involved_faces;
-        if (loc == intersecting_and_coplanar.end()) {
-            involved_faces.push_back(fid);
-        } else {
-            std::queue<Index> Q;
-            Q.push(fid);
-            while (!Q.empty()) {
-                const auto index = Q.front();
-                involved_faces.push_back(index);
-                Q.pop();
-
-                const auto overlapping_faces =
-                    intersecting_and_coplanar.find(index);
-                assert(overlapping_faces != intersecting_and_coplanar.end());
-
-                for (const auto other_index : overlapping_faces->second) {
-                    if (processed[other_index]) continue;
-                    processed[other_index] = true;
-                    Q.push(other_index);
-                }
-            }
+    for (const auto itr : offending) 
+    {
+      const auto fid = itr.first;
+      if (processed[fid]) continue;
+      processed[fid] = true;
+
+      const auto loc = intersecting_and_coplanar.find(fid);
+      std::vector<Index> involved_faces;
+      if (loc == intersecting_and_coplanar.end()) 
+      {
+        involved_faces.push_back(fid);
+      } else 
+      {
+        std::queue<Index> Q;
+        Q.push(fid);
+        while (!Q.empty()) 
+        {
+          const auto index = Q.front();
+          involved_faces.push_back(index);
+          Q.pop();
+
+          const auto overlapping_faces = intersecting_and_coplanar.find(index);
+          assert(overlapping_faces != intersecting_and_coplanar.end());
+
+          for (const auto other_index : overlapping_faces->second) 
+          {
+            if (processed[other_index]) continue;
+            processed[other_index] = true;
+            Q.push(other_index);
+          }
         }
+      }
 
-        Plane_3 P = T[fid].supporting_plane();
-        cdt_inputs.emplace_back(P, involved_faces);
+      Plane_3 P = T[fid].supporting_plane();
+      cdt_inputs.emplace_back(P, involved_faces);
     }
 #ifdef REMESH_INTERSECTIONS_TIMING
     log_time("preprocess");
@@ -372,8 +429,10 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
     std::vector<std::vector<Point_3> > cdt_vertices(num_cdts);
     std::vector<std::vector<std::vector<Index> > > cdt_faces(num_cdts);
 
-    const auto run_cdt = [&](const size_t first, const size_t last) {
-      for (size_t i=first; i<last; i++) {
+    const auto run_cdt = [&](const size_t first, const size_t last) 
+    {
+      for (size_t i=first; i<last; i++) 
+      {
         auto& vertices = cdt_vertices[i];
         auto& faces = cdt_faces[i];
         const auto& P = cdt_inputs[i].first;
@@ -386,7 +445,8 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
     log_time("cdt");
 #endif
 
-    for (size_t i=0; i<num_cdts; i++) {
+    for (size_t i=0; i<num_cdts; i++) 
+    {
       const auto& vertices = cdt_vertices[i];
       const auto& faces = cdt_faces[i];
       const auto& involved_faces = cdt_inputs[i].second;
@@ -399,23 +459,27 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
     // 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=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));
+
+    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];
+    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);
@@ -427,22 +491,27 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
 
     DerivedVV unique_vv;
     Eigen::VectorXi unique_to_vv, vv_to_unique;
+    // This is not stable... So even if offending is empty V != VV in
+    // general...
     igl::unique_rows(VV, unique_vv, unique_to_vv, vv_to_unique);
-    if (!stitch_all) {
-      // Vertices with the same coordinates would be represented by one vertex.
-      // The IM value of an vertex is the index of the representative vertex.
-      for (Index v=0; v<(Index)VV_size; v++) {
-        IM(v) = unique_to_vv[vv_to_unique[v]];
-      }
-    } else {
-      // Screw IM and representative vertices.  Merge all vertices having the
-      // same coordinates into a single vertex and set IM to identity map.
+    if(stitch_all) 
+    {
+      // Merge all vertices having the same coordinates into a single vertex
+      // and set IM to identity map.
       VV = unique_vv;
       std::transform(FF.data(), FF.data() + FF.rows()*FF.cols(),
           FF.data(), [&vv_to_unique](const typename DerivedFF::Scalar& a)
           { return vv_to_unique[a]; });
       IM.setLinSpaced(unique_vv.rows(), 0, unique_vv.rows()-1);
+    }else 
+    {
+      // Vertices with the same coordinates would be represented by one vertex.
+      // The IM value of a vertex is the index of the representative vertex.
+      for (Index v=0; v<(Index)VV_size; v++) {
+        IM(v) = unique_to_vv[vv_to_unique[v]];
+      }
     }
+
 #ifdef REMESH_INTERSECTIONS_TIMING
     log_time("store_results");
 #endif

include/igl/copyleft/cgal/remesh_intersections.h

@@ -31,13 +31,16 @@ namespace igl
       //     of first finding and list of intersection objects from all
       //     intersections
       //   edge2faces  #edges <= #offending*3 to incident offending faces 
+      //   stitch_all  if true, merge all vertices with thte same coordiante.
       // Outputs:
-      //   VV  #VV by 3 list of vertex positions
+      //   VV  #VV by 3 list of vertex positions, if stitch_all = false then
+      //     first #V vertices will always be V
       //   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.
+      //   IM  / stitch_all = true   #VV list from 0 to #VV-1
+      //       \ stitch_all = false  #VV list of indices into VV of unique vertices.
       //
       template <
         typename DerivedV,
@@ -58,15 +61,12 @@ namespace igl
         const std::map<
           std::pair<typename DerivedF::Index,typename DerivedF::Index>,
           std::vector<typename DerivedF::Index> > & edge2faces,
+        bool stitch_all,
         Eigen::PlainObjectBase<DerivedVV> & VV,
         Eigen::PlainObjectBase<DerivedFF> & FF,
         Eigen::PlainObjectBase<DerivedJ> & J,
         Eigen::PlainObjectBase<DerivedIM> & IM);
-
-      // Same as above except ``stitch_all`` flag:
-      //
-      // Input:
-      //   stitch_all: if true, merge all vertices with thte same coordiante.
+      // Same as above except stitch_all is assumed "false"
       template <
         typename DerivedV,
         typename DerivedF,
@@ -86,7 +86,6 @@ namespace igl
         const std::map<
           std::pair<typename DerivedF::Index,typename DerivedF::Index>,
           std::vector<typename DerivedF::Index> > & edge2faces,
-        bool stitch_all,
         Eigen::PlainObjectBase<DerivedVV> & VV,
         Eigen::PlainObjectBase<DerivedFF> & FF,
         Eigen::PlainObjectBase<DerivedJ> & J,

include/igl/copyleft/tetgen/cdt.cpp

@@ -8,6 +8,7 @@
 #include "cdt.h"
 #include "../../bounding_box.h"
 #include "../../writeOBJ.h"
+#include "tetrahedralize.h"
 
 template <
   typename DerivedV, 
@@ -63,4 +64,5 @@ IGL_INLINE bool igl::copyleft::tetgen::cdt(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
+template bool igl::copyleft::tetgen::cdt<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, 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::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::copyleft::tetgen::CDTParam 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> >&);
 #endif

include/igl/matrix_to_list.cpp

@@ -34,7 +34,7 @@ IGL_INLINE void igl::matrix_to_list(
 {
   using namespace std;
   V.resize(M.size());
-  // loop over rows
+  // loop over cols then rows
   for(int j = 0;j<M.cols();j++)
   {
     for(int i = 0;i<M.rows();i++)

include/igl/setdiff.cpp

@@ -80,4 +80,5 @@ IGL_INLINE void igl::setdiff(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 template void igl::setdiff<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::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::setdiff<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::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/writeDMAT.cpp

@@ -131,4 +131,7 @@ template bool igl::writeDMAT<Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::basic_
 template bool igl::writeDMAT<Eigen::Matrix<double, -1, 1, 0, -1, 1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::Matrix<double, -1, 1, 0, -1, 1> const&, bool);
 template bool igl::writeDMAT<Eigen::Matrix<int, -1, 1, 0, -1, 1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, bool);
 template bool igl::writeDMAT<Eigen::Matrix<double, -1, 3, 0, -1, 3> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::Matrix<double, -1, 3, 0, -1, 3> const&, bool);
+template bool igl::writeDMAT<Eigen::Array<int, -1, 1, 0, -1, 1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::Array<int, -1, 1, 0, -1, 1> const&, bool);
+template bool igl::writeDMAT<Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, bool);
+template bool igl::writeDMAT<Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> > >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> > const&, bool);
 #endif