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

Former-commit-id: f2d889b069e2204e8171dfcbbb4fe00b6f8938ad

include/igl/point_simplex_squared_distance.cpp

@@ -20,26 +20,29 @@ template <
   typename DerivedV,
   typename DerivedEle,
   typename Derivedsqr_d,
-  typename Derivedc>
+  typename Derivedc,
+  typename Derivedb>
 IGL_INLINE void igl::point_simplex_squared_distance(
   const Eigen::MatrixBase<Derivedp> & p,
   const Eigen::MatrixBase<DerivedV> & V,
   const Eigen::MatrixBase<DerivedEle> & Ele,
   const typename DerivedEle::Index primitive,
   Derivedsqr_d & sqr_d,
-  Eigen::PlainObjectBase<Derivedc> & c)
+  Eigen::PlainObjectBase<Derivedc> & c,
+  Eigen::PlainObjectBase<Derivedb> & bary)
 {
   typedef typename Derivedp::Scalar Scalar;
   typedef typename Eigen::Matrix<Scalar,1,DIM> Vector;
   typedef Vector Point;
+  typedef Eigen::PlainObjectBase<Derivedb> BaryPoint;
 
   const auto & Dot = [](const Point & a, const Point & b)->Scalar
   {
     return a.dot(b);
   };
   // Real-time collision detection, Ericson, Chapter 5
-  const auto & ClosestPtPointTriangle = 
-    [&Dot](Point p, Point a, Point b, Point c)->Point 
+  const auto & ClosestBaryPtPointTriangle = 
+    [&Dot](Point p, Point a, Point b, Point c, BaryPoint& bary_out )->Point 
   {
     // Check if P in vertex region outside A
     Vector ab = b - a;
@@ -47,42 +50,61 @@ IGL_INLINE void igl::point_simplex_squared_distance(
     Vector ap = p - a;
     Scalar d1 = Dot(ab, ap);
     Scalar d2 = Dot(ac, ap);
-    if (d1 <= 0.0 && d2 <= 0.0) return a; // barycentric coordinates (1,0,0)
+    if (d1 <= 0.0 && d2 <= 0.0) {
+      // barycentric coordinates (1,0,0)
+      bary_out << 1, 0, 0;
+      return a;
+    }
     // Check if P in vertex region outside B
     Vector bp = p - b;
     Scalar d3 = Dot(ab, bp);
     Scalar d4 = Dot(ac, bp);
-    if (d3 >= 0.0 && d4 <= d3) return b; // barycentric coordinates (0,1,0)
+    if (d3 >= 0.0 && d4 <= d3) {
+      // barycentric coordinates (0,1,0)
+      bary_out << 0, 1, 0;
+      return b;
+    }
     // Check if P in edge region of AB, if so return projection of P onto AB
     Scalar vc = d1*d4 - d3*d2;
     if( a != b)
     {
       if (vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0) {
         Scalar v = d1 / (d1 - d3);
-        return a + v * ab; // barycentric coordinates (1-v,v,0)
+        // barycentric coordinates (1-v,v,0)
+        bary_out << 1-v, v, 0;
+        return a + v * ab;
       }
     }
     // Check if P in vertex region outside C
     Vector cp = p - c;
     Scalar d5 = Dot(ab, cp);
     Scalar d6 = Dot(ac, cp);
-    if (d6 >= 0.0 && d5 <= d6) return c; // barycentric coordinates (0,0,1)
+    if (d6 >= 0.0 && d5 <= d6) {
+      // barycentric coordinates (0,0,1)
+      bary_out << 0, 0, 1;
+      return c;
+    }
     // Check if P in edge region of AC, if so return projection of P onto AC
     Scalar vb = d5*d2 - d1*d6;
     if (vb <= 0.0 && d2 >= 0.0 && d6 <= 0.0) {
       Scalar w = d2 / (d2 - d6);
-      return a + w * ac; // barycentric coordinates (1-w,0,w)
+      // barycentric coordinates (1-w,0,w)
+      bary_out << 1-w, 0, w;
+      return a + w * ac;
     }
     // Check if P in edge region of BC, if so return projection of P onto BC
     Scalar va = d3*d6 - d5*d4;
     if (va <= 0.0 && (d4 - d3) >= 0.0 && (d5 - d6) >= 0.0) {
       Scalar w = (d4 - d3) / ((d4 - d3) + (d5 - d6));
-      return b + w * (c - b); // barycentric coordinates (0,1-w,w)
+      // barycentric coordinates (0,1-w,w)
+      bary_out << 0, 1-w, w;
+      return b + w * (c - b);
     }
     // P inside face region. Compute Q through its barycentric coordinates (u,v,w)
     Scalar denom = 1.0 / (va + vb + vc);
     Scalar v = vb * denom;
     Scalar w = vc * denom;
+    bary_out << 1.0-v-w, v, w;
     return a + ab * v + ac * w; // = u*a + v*b + w*c, u = va * denom = 1.0-v-w
   };
 
@@ -91,16 +113,38 @@ IGL_INLINE void igl::point_simplex_squared_distance(
   assert(Ele.cols() <= DIM+1);
   assert(Ele.cols() <= 3 && "Only simplices up to triangles are considered");
 
-  c = ClosestPtPointTriangle(
+  bary.resize( DIM, 1 );
+  c = ClosestBaryPtPointTriangle(
     p,
     V.row(Ele(primitive,0)),
     V.row(Ele(primitive,1%Ele.cols())),
-    V.row(Ele(primitive,2%Ele.cols())));
+    V.row(Ele(primitive,2%Ele.cols())),
+    bary);
   sqr_d = (p-c).squaredNorm();
 }
+template <
+  int DIM,
+  typename Derivedp,
+  typename DerivedV,
+  typename DerivedEle,
+  typename Derivedsqr_d,
+  typename Derivedc>
+IGL_INLINE void igl::point_simplex_squared_distance(
+  const Eigen::MatrixBase<Derivedp> & p,
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedEle> & Ele,
+  const typename DerivedEle::Index primitive,
+  Derivedsqr_d & sqr_d,
+  Eigen::PlainObjectBase<Derivedc> & c)
+{
+    Eigen::PlainObjectBase<Derivedc> b;
+    point_simplex_squared_distance<DIM>( p, V, Ele, primitive, sqr_d, c, b );
+}
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instanciation
 template void igl::point_simplex_squared_distance<3, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, double&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&);
 template void igl::point_simplex_squared_distance<2, Eigen::Matrix<double, 1, 2, 1, 1, 2>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, Eigen::Matrix<double, 1, 2, 1, 1, 2> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, double&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> >&);
+template void igl::point_simplex_squared_distance<3, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, double&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&);
+template void igl::point_simplex_squared_distance<2, Eigen::Matrix<double, 1, 2, 1, 1, 2>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, Eigen::Matrix<double, 1, 2, 1, 1, 2> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, double&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> >&);
 #endif

include/igl/point_simplex_squared_distance.h

@@ -36,6 +36,35 @@ namespace igl
     const typename DerivedEle::Index i,
     Derivedsqr_d & sqr_d,
     Eigen::PlainObjectBase<Derivedc> & c);
+  // Determine squared distance from a point to linear simplex.
+  // Also return barycentric coordinate of closest point. 
+  //
+  // Inputs:
+  //   p  d-long query point
+  //   V  #V by d list of vertices
+  //   Ele  #Ele by ss<=d+1 list of simplex indices into V
+  //   i  index into Ele of simplex
+  // Outputs:
+  //   sqr_d  squared distance of Ele(i) to p
+  //   c  closest point on Ele(i) 
+  //   b  barycentric coordinates of closest point on Ele(i) 
+  //
+  template <
+    int DIM,
+    typename Derivedp,
+    typename DerivedV,
+    typename DerivedEle,
+    typename Derivedsqr_d,
+    typename Derivedc,
+    typename Derivedb>
+  IGL_INLINE void point_simplex_squared_distance(
+    const Eigen::MatrixBase<Derivedp> & p,
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedEle> & Ele,
+    const typename DerivedEle::Index i,
+    Derivedsqr_d & sqr_d,
+    Eigen::PlainObjectBase<Derivedc> & c,
+    Eigen::PlainObjectBase<Derivedb> & b);
 }
 #ifndef IGL_STATIC_LIBRARY
 #  include "point_simplex_squared_distance.cpp"

include/igl/remove_duplicates.cpp

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

include/igl/seam_edges.cpp

@@ -11,23 +11,6 @@
 #include <unordered_set>
 #include <cassert>
 
-namespace {
-    // Computes the orientation of `c` relative to the line between `a` and `b`.
-    // Assumes 2D vector input.
-    // Based on: https://www.cs.cmu.edu/~quake/robust.html
-    template< typename scalar_t >
-    inline scalar_t orientation(
-        const Eigen::Matrix< scalar_t, 2, 1 >& a,
-        const Eigen::Matrix< scalar_t, 2, 1 >& b,
-        const Eigen::Matrix< scalar_t, 2, 1 >& c
-        ) {
-        typedef Eigen::Matrix< scalar_t, 2, 1 > Vector2S;
-        const Vector2S row0 = a - c;
-        const Vector2S row1 = b - c;
-        return row0(0)*row1(1) - row1(0)*row0(1);
-    }
-}
-
 // I have verified that this function produces the exact same output as
 // `find_seam_fast.py` for `cow_triangled.obj`.
 template <typename DerivedV, typename DerivedF, typename DerivedT>
@@ -49,6 +32,15 @@ IGL_INLINE void igl::seam_edges(
     // Assume 2D texture coordinates (foldovers tests).
     assert( TC.cols() == 2 );
     typedef Eigen::Matrix< typename DerivedT::Scalar, 2, 1 > Vector2S;
+    // Computes the orientation of `c` relative to the line between `a` and `b`.
+    // Assumes 2D vector input.
+    // Based on: https://www.cs.cmu.edu/~quake/robust.html
+    const auto& Orientation = []( const Vector2S& a, const Vector2S& b, const Vector2S& c ) -> typename DerivedT::Scalar
+    {
+        const Vector2S row0 = a - c;
+        const Vector2S row1 = b - c;
+        return row0(0)*row1(1) - row1(0)*row0(1);
+    };
     
     seams     .setZero( 3*F.rows(), 4 );
     boundaries.setZero( 3*F.rows(), 2 );
@@ -73,7 +65,7 @@ IGL_INLINE void igl::seam_edges(
     typedef std::pair< typename DerivedF::Scalar, typename DerivedF::Scalar > directed_edge;
 	const int numV = V.rows();
 	const int numF = F.rows();
-	auto edge_hasher = [numV]( directed_edge const& e ) { return e.first*numV + e.second; };
+	const auto& edge_hasher = [numV]( directed_edge const& e ) { return e.first*numV + e.second; };
 	// When we pass a hash function object, we also need to specify the number of buckets.
 	// The Euler characteristic says that the number of undirected edges is numV + numF -2*genus.
 	std::unordered_map< directed_edge, std::pair< int, int >, decltype( edge_hasher ) > directed_position_edge2face_position_index( 2*( numV + numF ), edge_hasher );
@@ -84,12 +76,13 @@ IGL_INLINE void igl::seam_edges(
         }
     }
     
-    // First find all undirected position edges (collect both orientations of
+    // First find all undirected position edges (collect a canonical orientation of
     // the directed edges).
     std::unordered_set< directed_edge, decltype( edge_hasher ) > undirected_position_edges( numV + numF, edge_hasher );
     for( const auto& el : directed_position_edge2face_position_index ) {
-        undirected_position_edges.insert( el.first );
-        undirected_position_edges.insert( std::make_pair( el.first.second, el.first.first ) );
+        // The canonical orientation is the one where the smaller of
+        // the two vertex indices is first.
+        undirected_position_edges.insert( std::make_pair( std::min( el.first.first, el.first.second ), std::max( el.first.first, el.first.second ) ) );
     }
     
     // Now we will iterate over all position edges.
@@ -97,6 +90,10 @@ IGL_INLINE void igl::seam_edges(
     // texcoord indices (or one doesn't exist at all in the case of a mesh
     // boundary).
     for( const auto& vp_edge : undirected_position_edges ) {
+        // We should only see canonical edges,
+        // where the first vertex index is smaller.
+        assert( vp_edge.first < vp_edge.second );
+        
         const auto vp_edge_reverse = std::make_pair( vp_edge.second, vp_edge.first );
         // If it and its opposite exist as directed edges, check if their
         // texture coordinate indices match.
@@ -108,12 +105,6 @@ IGL_INLINE void igl::seam_edges(
             // NOTE: They should never be equal.
             assert( forwards != backwards );
             
-            // NOTE: Non-matching seam edges will show up twice, once as
-            //       ( forwards, backwards ) and once as
-            //       ( backwards, forwards ). We don't need to examine both,
-            //       so continue only if forwards < backwards.
-            if( forwards < backwards ) continue;
-
             // If the texcoord indices match (are similarly flipped),
             // this edge is not a seam. It could be a foldover.
             if( std::make_pair( FTC( forwards.first, forwards.second ), FTC( forwards.first, ( forwards.second+1 ) % 3 ) )
@@ -129,8 +120,8 @@ IGL_INLINE void igl::seam_edges(
                 const Vector2S c_backwards = TC.row( FTC( backwards.first, (backwards.second+2) % 3 ) );
                 // If the opposite vertices' texture coordinates fall on the same side
                 // of the edge, we have a UV-space foldover.
-                const auto orientation_forwards = orientation( a, b, c_forwards );
-                const auto orientation_backwards = orientation( a, b, c_backwards );
+                const auto orientation_forwards = Orientation( a, b, c_forwards );
+                const auto orientation_backwards = Orientation( a, b, c_backwards );
                 if( ( orientation_forwards > 0 && orientation_backwards > 0 ) ||
                     ( orientation_forwards < 0 && orientation_backwards < 0 )
                     ) {

include/igl/writeDMAT.cpp

@@ -18,7 +18,7 @@ IGL_INLINE bool igl::writeDMAT(
   const Mat & W,
   const bool ascii)
 {
-  FILE * fp = fopen(file_name.c_str(),"w");
+  FILE * fp = fopen(file_name.c_str(),"wb");
   if(fp == NULL)
   {
     fprintf(stderr,"IOError: writeDMAT() could not open %s...",file_name.c_str());

include/igl/writeOFF.cpp

@@ -43,4 +43,5 @@ template bool igl::writeOFF<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matri
 template bool igl::writeOFF<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&);
 template bool igl::writeOFF<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&);
 template bool igl::writeOFF<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&);
+template bool igl::writeOFF<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&);
 #endif

python/modules/copyleft/py_igl_cgal.cpp

@@ -6,6 +6,8 @@
 #include "../../python_shared.h"
 
 #include <igl/copyleft/cgal/mesh_boolean.h>
+#include <igl/copyleft/cgal/remesh_self_intersections.h>
+#include <igl/copyleft/cgal/RemeshSelfIntersectionsParam.h>
 
 
 void python_export_igl_cgal(py::module &me) {
@@ -14,5 +16,7 @@ void python_export_igl_cgal(py::module &me) {
     "cgal", "Wrappers for libigl functions that use cgal");
 
   #include "../../py_igl/copyleft/cgal/py_mesh_boolean.cpp"
+  #include "../../py_igl/copyleft/cgal/py_remesh_self_intersections.cpp"
+  #include "../../py_igl/copyleft/cgal/py_RemeshSelfIntersectionsParam.cpp"
 
 }

python/py_doc.cpp

@@ -182,6 +182,40 @@ const char *__doc_igl_copyleft_cgal_mesh_boolean = R"igl_Qu8mg5v7(//  MESH_BOOLE
       //
       //  See also: mesh_boolean_cork, intersect_other,
       //  remesh_self_intersections)igl_Qu8mg5v7";
+const char *__doc_igl_copyleft_cgal_remesh_self_intersections = R"igl_Qu8mg5v7(// Given a triangle mesh (V,F) compute a new mesh (VV,FF) which is the same
+      // as (V,F) except that any self-intersecting triangles in (V,F) have been
+      // subdivided (new vertices and face created) so that the self-intersection
+      // contour lies exactly on edges in (VV,FF). New vertices will appear in
+      // original faces or on original edges. New vertices on edges are "merged"
+      // only across original faces sharing that edge. This means that if the input
+      // triangle mesh is a closed manifold the output will be too.
+      //
+      // Inputs:
+      //   V  #V by 3 list of vertex positions
+      //   F  #F by 3 list of triangle indices into V
+      //   params  struct of optional parameters
+      // Outputs:
+      //   VV  #VV by 3 list of vertex positions
+      //   FF  #FF by 3 list of triangle indices into VV
+      //   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.
+      //
+      // Known bugs: If an existing edge in (V,F) lies exactly on another face then
+      // any resulting additional vertices along that edge may not get properly
+      // connected so that the output mesh has the same global topology. This is
+      // because 
+      //
+      // Example:
+      //     // resolve intersections
+      //     igl::copyleft::cgal::remesh_self_intersections(V,F,params,VV,FF,IF,J,IM);
+      //     // _apply_ duplicate vertex mapping IM to FF
+      //     for_each(FF.data(),FF.data()+FF.size(),[&IM](int & a){a=IM(a);});
+      //     // remove any vertices now unreferenced after duplicate mapping.
+      //     igl::remove_unreferenced(VV,FF,SV,SF,UIM);
+      //     // Now (SV,SF) is ready to extract outer hull
+      //     igl::copyleft::cgal::outer_hull(SV,SF,G,J,flip);igl_Qu8mg5v7";
 const char *__doc_igl_copyleft_comiso_miq = R"igl_Qu8mg5v7(// Inputs:
     //   V              #V by 3 list of mesh vertex 3D positions
     //   F              #F by 3 list of faces indices in V
@@ -871,6 +905,25 @@ const char *__doc_igl_read_triangle_mesh = R"igl_Qu8mg5v7(// read mesh from an a
   //   V  eigen double matrix #V by 3
   //   F  eigen int matrix #F by 3
   // Returns true iff success)igl_Qu8mg5v7";
+const char *__doc_igl_remove_duplicate_vertices = R"igl_Qu8mg5v7(  // REMOVE_DUPLICATE_VERTICES Remove duplicate vertices upto a uniqueness
+  // tolerance (epsilon)
+  //
+  // Inputs:
+  //   V  #V by dim list of vertex positions
+  //   epsilon  uniqueness tolerance (significant digit), can probably think of
+  //     this as a tolerance on L1 distance
+  // Outputs:
+  //   SV  #SV by dim new list of vertex positions
+  //   SVI #V by 1 list of indices so SV = V(SVI,:)
+  //   SVJ #SV by 1 list of indices so V = SV(SVJ,:)
+  //
+  // Example:
+  //   % Mesh in (V,F)
+  //   [SV,SVI,SVJ] = remove_duplicate_vertices(V,1e-7);
+  //   % remap faces
+  //   SF = SVJ(F);
+  //
+  //)igl_Qu8mg5v7";
 const char *__doc_igl_rotate_vectors = R"igl_Qu8mg5v7(// Rotate the vectors V by A radiants on the tangent plane spanned by B1 and
   // B2
   //

python/py_doc.h

@@ -13,6 +13,7 @@ extern const char *__doc_igl_comb_cross_field;
 extern const char *__doc_igl_comb_frame_field;
 extern const char *__doc_igl_compute_frame_field_bisectors;
 extern const char *__doc_igl_copyleft_cgal_mesh_boolean;
+extern const char *__doc_igl_copyleft_cgal_remesh_self_intersections;
 extern const char *__doc_igl_copyleft_comiso_miq;
 extern const char *__doc_igl_copyleft_comiso_nrosy;
 extern const char *__doc_igl_copyleft_marching_cubes;
@@ -71,6 +72,7 @@ extern const char *__doc_igl_readMESH;
 extern const char *__doc_igl_readOBJ;
 extern const char *__doc_igl_readOFF;
 extern const char *__doc_igl_read_triangle_mesh;
+extern const char *__doc_igl_remove_duplicate_vertices;
 extern const char *__doc_igl_rotate_vectors;
 extern const char *__doc_igl_setdiff;
 extern const char *__doc_igl_signed_distance;

python/py_igl.cpp

@@ -61,6 +61,7 @@
 #include <igl/readOBJ.h>
 #include <igl/readOFF.h>
 #include <igl/read_triangle_mesh.h>
+#include <igl/remove_duplicate_vertices.h>
 #include <igl/rotate_vectors.h>
 #include <igl/setdiff.h>
 #include <igl/signed_distance.h>
@@ -140,6 +141,7 @@ void python_export_igl(py::module &m)
 #include "py_igl/py_readOBJ.cpp"
 #include "py_igl/py_readOFF.cpp"
 #include "py_igl/py_read_triangle_mesh.cpp"
+#include "py_igl/py_remove_duplicate_vertices.cpp"
 #include "py_igl/py_rotate_vectors.cpp"
 #include "py_igl/py_setdiff.cpp"
 #include "py_igl/py_signed_distance.cpp"

python/py_igl/copyleft/cgal/py_RemeshSelfIntersectionsParam.cpp

@@ -0,0 +1,14 @@
+py::class_<igl::copyleft::cgal::RemeshSelfIntersectionsParam > RemeshSelfIntersectionsParam(m, "RemeshSelfIntersectionsParam");
+
+RemeshSelfIntersectionsParam
+.def("__init__", [](igl::copyleft::cgal::RemeshSelfIntersectionsParam &m)
+{
+    new (&m) igl::copyleft::cgal::RemeshSelfIntersectionsParam();
+    m.detect_only = false;
+    m.first_only = false;
+    m.stitch_all = false;
+})
+.def_readwrite("detect_only", &igl::copyleft::cgal::RemeshSelfIntersectionsParam::detect_only)
+.def_readwrite("first_only", &igl::copyleft::cgal::RemeshSelfIntersectionsParam::first_only)
+.def_readwrite("stitch_all", &igl::copyleft::cgal::RemeshSelfIntersectionsParam::stitch_all)
+;

python/py_igl/copyleft/cgal/py_remesh_self_intersections.cpp

@@ -0,0 +1,22 @@
+m.def("remesh_self_intersections", []
+(
+	const Eigen::MatrixXd& V,	
+	const Eigen::MatrixXi& F,
+	const igl::copyleft::cgal::RemeshSelfIntersectionsParam& params,
+	Eigen::MatrixXd& VV,
+	Eigen::MatrixXi& FF,
+	Eigen::MatrixXi& IF,
+	Eigen::MatrixXi& J,
+	Eigen::MatrixXi& IM
+)
+{
+	assert_is_VectorX("J", J);
+	assert_is_VectorX("IM", IM);
+	Eigen::VectorXi Jt;
+	Eigen::VectorXi IMt;
+	igl::copyleft::cgal::remesh_self_intersections(V, F, params, VV, FF, IF, Jt, IMt);
+	J = Jt;
+	IM = IMt;
+}, __doc_igl_copyleft_cgal_remesh_self_intersections,
+py::arg("V"), py::arg("F"), py::arg("params"), py::arg("VV")
+, py::arg("FF"), py::arg("IF"), py::arg("J"), py::arg("IM")); 

python/py_igl/py_remove_duplicate_vertices.cpp

@@ -0,0 +1,14 @@
+m.def("remove_duplicate_vertices", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  const double epsilon,
+  Eigen::MatrixXd& SV,
+  Eigen::MatrixXi& SVI,
+  Eigen::MatrixXi& SVJ,
+  Eigen::MatrixXi& SF
+)
+{
+  return igl::remove_duplicate_vertices(V, F, epsilon, SV, SVI, SVJ, SF);
+}, __doc_igl_remove_duplicate_vertices,
+py::arg("V"), py::arg("F"), py::arg("epsilon"), py::arg("SV"), py::arg("SVI"), py::arg("SVJ"), py::arg("SF"));

python/python_shared.cpp

@@ -66,6 +66,8 @@ PYBIND11_PLUGIN(pyigl) {
            comb_frame_field
            compute_frame_field_bisectors
            copyleft_cgal_mesh_boolean
+           copyleft_cgal_remesh_self_intersections
+           copyleft_cgal_RemeshSelfIntersectionsParam
            copyleft_comiso_miq
            copyleft_comiso_nrosy
            copyleft_marching_cubes
@@ -117,6 +119,7 @@ PYBIND11_PLUGIN(pyigl) {
            readOBJ
            readOFF
            read_triangle_mesh
+           remove_duplicate_vertices
            rotate_vectors
            setdiff
            signed_distance

shared/cmake/FindMATLAB.cmake

@@ -137,7 +137,7 @@ ELSE(WIN32)
     IF((NOT DEFINED MATLAB_ROOT) OR ("${MATLAB_ROOT}" STREQUAL ""))
 
     # Search for a version of Matlab available, starting from the most modern one to older versions
-      FOREACH(MATVER "R2016a" "R2015b" "R2015a" "R2014b" "R2014a" "R2014a" "R2013b" "R2013a" "R2012b" "R2012a" "R2011b" "R2011a" "R2010b" "R2010a" "R2009b" "R2009a" "R2008b")
+      FOREACH(MATVER "R2016b" "R2016a" "R2015b" "R2015a" "R2014b" "R2014a" "R2014a" "R2013b" "R2013a" "R2012b" "R2012a" "R2011b" "R2011a" "R2010b" "R2010a" "R2009b" "R2009a" "R2008b")
         IF((NOT DEFINED MATLAB_ROOT) OR ("${MATLAB_ROOT}" STREQUAL ""))
           IF(EXISTS /Applications/MATLAB_${MATVER}.app)
             SET(MATLAB_ROOT /Applications/MATLAB_${MATVER}.app)