Merge branch 'dev' of https://github.com/libigl/libigl into dev

Former-commit-id: f6be209bad67412327be32e04653f2db29a448b5

include/igl/adjacency_list.cpp

@@ -165,4 +165,5 @@ template void igl::adjacency_list<Eigen::Matrix<int, -1, 2, 0, -1, 2>, int>(Eige
 // generated by autoexplicit.sh
 template void igl::adjacency_list<Eigen::Matrix<int, -1, -1, 0, -1, -1>, int>(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&, bool);
 template void igl::adjacency_list<Eigen::Matrix<int, -1, 3, 0, -1, 3>, int>(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&, bool);
+template void igl::adjacency_list<class Eigen::Matrix<int, -1, -1, 0, -1, -1>, unsigned int>(class Eigen::PlainObjectBase<class Eigen::Matrix<int, -1, -1, 0, -1, -1> > const &, class std::vector<class std::vector<unsigned int, class std::allocator<unsigned int> >, class std::allocator<class std::vector<unsigned int, class std::allocator<unsigned int> > > > &, bool);
 #endif

include/igl/centroid.cpp

@@ -64,4 +64,5 @@ IGL_INLINE void igl::centroid(
 // generated by autoexplicit.sh
 template void igl::centroid<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<float, 1, 3, 1, 1, 3> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, 1, 3, 1, 1, 3> >&);
 template void igl::centroid<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, 3, 1, 0, 3, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 1, 0, 3, 1> >&);
+template void igl::centroid<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&);
 #endif

include/igl/copyleft/marching_cubes.cpp

@@ -33,44 +33,45 @@ extern const int edgeTable[256];
 extern const int triTable[256][2][17];
 extern const int polyTable[8][16];
 
-struct EdgeKey
-{
-  EdgeKey(unsigned i0, unsigned i1) : i0_(i0), i1_(i1) {}
 
-  bool operator==(const EdgeKey& _rhs) const
+template <typename DerivedValues, typename DerivedPoints, typename DerivedVertices, typename DerivedIndices, typename DerivedFaces>
+class MarchingCubes
+{
+  struct EdgeKey
   {
-    return i0_ == _rhs.i0_ && i1_ == _rhs.i1_;
-  }
+    EdgeKey(unsigned i0, unsigned i1) : i0_(i0), i1_(i1) {}
 
-  unsigned i0_, i1_;
-};
+    bool operator==(const EdgeKey& _rhs) const
+    {
+      return i0_ == _rhs.i0_ && i1_ == _rhs.i1_;
+    }
 
-struct EdgeHash
-{
+    unsigned i0_, i1_;
+  };
+
+  struct EdgeHash
+  {
     std::size_t operator()(const EdgeKey& key) const {
-        std::size_t seed = 0;
-        seed ^= key.i0_ + 0x9e3779b9 + (seed<<6) + (seed>>2); // Copied from boost::hash_combine
-        seed ^= key.i1_ + 0x9e3779b9 + (seed<<6) + (seed>>2);
-        return std::hash<std::size_t>()(seed);
+      std::size_t seed = 0;
+      seed ^= key.i0_ + 0x9e3779b9 + (seed<<6) + (seed>>2); // Copied from boost::hash_combine
+      seed ^= key.i1_ + 0x9e3779b9 + (seed<<6) + (seed>>2);
+      return std::hash<std::size_t>()(seed);
     }
-};
-
+  };
 
-template <typename Derivedvalues, typename Derivedpoints,typename Derivedvertices, typename DerivedF>
-class MarchingCubes
-{
   typedef std::unordered_map<EdgeKey, unsigned, EdgeHash> MyMap;
   typedef typename MyMap::const_iterator                  MyMapIterator;
 
 public:
-  MarchingCubes(
-                const Eigen::PlainObjectBase<Derivedvalues> &values,
-                const Eigen::PlainObjectBase<Derivedpoints> &points,
+  // Dense index grid version
+  MarchingCubes(const Eigen::MatrixBase<DerivedValues> &values,
+                const Eigen::MatrixBase<DerivedPoints> &points,
                 const unsigned x_res,
                 const unsigned y_res,
                 const unsigned z_res,
-                Eigen::PlainObjectBase<Derivedvertices> &vertices,
-                Eigen::PlainObjectBase<DerivedF> &faces)
+                const double isovalue,
+                Eigen::PlainObjectBase<DerivedVertices> &vertices,
+                Eigen::PlainObjectBase<DerivedFaces> &faces)
   {
     assert(values.cols() == 1);
     assert(points.cols() == 3);
@@ -101,7 +102,7 @@ public:
     {
 
       unsigned         corner[8];
-      typename DerivedF::Scalar samples[12];
+      typename DerivedFaces::Scalar samples[12];
       unsigned char    cubetype(0);
       unsigned int     i;
 
@@ -126,7 +127,7 @@ public:
 
       // determine cube type
       for (i=0; i<8; ++i)
-        if (values(corner[i]) > 0.0)
+        if (values(corner[i]) > isovalue)
           cubetype |= (1<<i);
 
 
@@ -182,77 +183,220 @@ public:
     vertices.conservativeResize(num_vertices, Eigen::NoChange);
     faces.conservativeResize(num_faces, Eigen::NoChange);
 
-  };
+  }
 
-  static typename DerivedF::Scalar  add_vertex(const Eigen::PlainObjectBase<Derivedvalues> &values,
-                                               const Eigen::PlainObjectBase<Derivedpoints> &points,
-                                               unsigned int i0,
-                                               unsigned int i1,
-                                               Eigen::PlainObjectBase<Derivedvertices> &vertices,
-                                               int &num_vertices,
-                                               MyMap &edge2vertex)
+  // Sparse index grid version
+  MarchingCubes(const Eigen::MatrixBase<DerivedValues> &values,
+                const Eigen::MatrixBase<DerivedPoints> &points,
+                const Eigen::MatrixBase<DerivedIndices> &cubes,
+                const double isovalue,
+                Eigen::PlainObjectBase<DerivedVertices> &vertices,
+                Eigen::PlainObjectBase<DerivedFaces> &faces)
+  {
+    assert(values.cols() == 1);
+    assert(points.cols() == 3);
+    assert(cubes.cols() == 8);
+
+    if(cubes.rows() == 0)
+    {
+      return;
+    }
+
+    faces.resize(10000,3);
+    int num_faces = 0;
+
+    vertices.resize(10000,3);
+    int num_vertices = 0;
+
+
+    unsigned n_cubes  = cubes.rows();
+
+    for (unsigned cube_it =0 ; cube_it < n_cubes; ++cube_it)
+    {
+      typedef Eigen::Matrix<typename DerivedIndices::Scalar, 1, 8> CubeIndexVector;
+      typedef typename DerivedFaces::Scalar SampleScalar;
+
+      CubeIndexVector  cube = cubes.row(cube_it);
+      SampleScalar     samples[12];
+      unsigned char    cubetype(0);
+
+      // determine cube type
+      for (int i=0; i<8; ++i)
+      {
+        if (values[cube[i]] > isovalue)
+        {
+          cubetype |= (1<<i);
+        }
+      }
+
+
+      // trivial reject ?
+      if (cubetype == 0 || cubetype == 255)
+      {
+        continue;
+      }
+
+      // compute samples on cube's edges
+      if (edgeTable[cubetype]&1)
+        samples[0]  = add_vertex(values, points, cube[0], cube[1], vertices, num_vertices, edge2vertex);
+      if (edgeTable[cubetype]&2)
+        samples[1]  = add_vertex(values, points, cube[1], cube[2], vertices, num_vertices, edge2vertex);
+      if (edgeTable[cubetype]&4)
+        samples[2]  = add_vertex(values, points, cube[3], cube[2], vertices, num_vertices, edge2vertex);
+      if (edgeTable[cubetype]&8)
+        samples[3]  = add_vertex(values, points, cube[0], cube[3], vertices, num_vertices, edge2vertex);
+      if (edgeTable[cubetype]&16)
+        samples[4]  = add_vertex(values, points, cube[4], cube[5], vertices, num_vertices, edge2vertex);
+      if (edgeTable[cubetype]&32)
+        samples[5]  = add_vertex(values, points, cube[5], cube[6], vertices, num_vertices, edge2vertex);
+      if (edgeTable[cubetype]&64)
+        samples[6]  = add_vertex(values, points, cube[7], cube[6], vertices, num_vertices, edge2vertex);
+      if (edgeTable[cubetype]&128)
+        samples[7]  = add_vertex(values, points, cube[4], cube[7], vertices, num_vertices, edge2vertex);
+      if (edgeTable[cubetype]&256)
+        samples[8]  = add_vertex(values, points, cube[0], cube[4], vertices, num_vertices, edge2vertex);
+      if (edgeTable[cubetype]&512)
+        samples[9]  = add_vertex(values, points, cube[1], cube[5], vertices, num_vertices, edge2vertex);
+      if (edgeTable[cubetype]&1024)
+        samples[10] = add_vertex(values, points, cube[2], cube[6], vertices, num_vertices, edge2vertex);
+      if (edgeTable[cubetype]&2048)
+        samples[11] = add_vertex(values, points, cube[3], cube[7], vertices, num_vertices, edge2vertex);
+
+      // connect samples by triangles
+      for (int i=0; triTable[cubetype][0][i] != -1; i+=3 )
+      {
+        num_faces++;
+        if (num_faces > faces.rows())
+        {
+          faces.conservativeResize(faces.rows()+10000, Eigen::NoChange);
+        }
+        faces.row(num_faces-1) <<
+        samples[triTable[cubetype][0][i  ]],
+        samples[triTable[cubetype][0][i+1]],
+        samples[triTable[cubetype][0][i+2]];
+
+      }
+
+    }
+
+    vertices.conservativeResize(num_vertices, Eigen::NoChange);
+    faces.conservativeResize(num_faces, Eigen::NoChange);
+
+  }
+
+  static typename DerivedFaces::Scalar add_vertex(const Eigen::MatrixBase<DerivedValues> &values,
+                                                  const Eigen::MatrixBase<DerivedPoints> &points,
+                                                  unsigned int i0,
+                                                  unsigned int i1,
+                                                  Eigen::PlainObjectBase<DerivedVertices> &vertices,
+                                                  int &num_vertices,
+                                                  MyMap &edge2vertex)
   {
     // find vertex if it has been computed already
     MyMapIterator it = edge2vertex.find(EdgeKey(i0, i1));
     if (it != edge2vertex.end())
+    {
       return it->second;
-    ;
+    }
 
     // generate new vertex
-    const Eigen::Matrix<typename Derivedpoints::Scalar, 1, 3> & p0 = points.row(i0);
-    const Eigen::Matrix<typename Derivedpoints::Scalar, 1, 3> & p1 = points.row(i1);
-
-    typename Derivedvalues::Scalar s0 = fabs(values(i0));
-    typename Derivedvalues::Scalar s1 = fabs(values(i1));
-    typename Derivedvalues::Scalar t  = s0 / (s0+s1);
+    const Eigen::Matrix<typename DerivedPoints::Scalar, 1, 3> & p0 = points.row(i0);
+    const Eigen::Matrix<typename DerivedPoints::Scalar, 1, 3> & p1 = points.row(i1);
 
+    typename DerivedValues::Scalar s0 = fabs(values[i0]);
+    typename DerivedValues::Scalar s1 = fabs(values[i1]);
+    typename DerivedValues::Scalar t  = s0 / (s0+s1);
 
     num_vertices++;
     if (num_vertices > vertices.rows())
+    {
       vertices.conservativeResize(vertices.rows()+10000, Eigen::NoChange);
+    }
 
     // Linear interpolation based on linearly interpolating values
-    vertices.row(num_vertices-1)  = ((1.0f-t)*p0 + t*p1).template cast<typename Derivedvertices::Scalar>();
+    vertices.row(num_vertices-1)  = ((1.0f-t)*p0 + t*p1).template cast<typename DerivedVertices::Scalar>();
     edge2vertex[EdgeKey(i0, i1)] = num_vertices-1;
 
     return num_vertices-1;
   }
-  ;
 
   // maps an edge to the sample vertex generated on it
   MyMap  edge2vertex;
 };
 
 
-template <typename Derivedvalues, typename Derivedpoints, typename Derivedvertices, typename DerivedF>
+template <typename DerivedValues, typename DerivedPoints, typename DerivedVertices, typename DerivedFaces>
+IGL_INLINE void igl::copyleft::marching_cubes(
+  const Eigen::MatrixBase<DerivedValues> &values,
+  const Eigen::MatrixBase<DerivedPoints> &points,
+  const unsigned x_res,
+  const unsigned y_res,
+  const unsigned z_res,
+  const double isovalue,
+  Eigen::PlainObjectBase<DerivedVertices> &vertices,
+  Eigen::PlainObjectBase<DerivedFaces> &faces)
+{
+  typedef Eigen::MatrixXi Shim; /* DerivedIndices shim type is unused in this instantiation*/
+  MarchingCubes<DerivedValues, DerivedPoints, DerivedVertices, Shim, DerivedFaces>
+          mc(values, points, x_res, y_res, z_res, isovalue, vertices, faces);
+}
+
+template <typename DerivedValues, typename DerivedPoints, typename DerivedVertices, typename DerivedFaces>
 IGL_INLINE void igl::copyleft::marching_cubes(
-  const Eigen::PlainObjectBase<Derivedvalues> &values,
-  const Eigen::PlainObjectBase<Derivedpoints> &points,
+  const Eigen::MatrixBase<DerivedValues> &values,
+  const Eigen::MatrixBase<DerivedPoints> &points,
   const unsigned x_res,
   const unsigned y_res,
   const unsigned z_res,
-  Eigen::PlainObjectBase<Derivedvertices> &vertices,
-  Eigen::PlainObjectBase<DerivedF> &faces)
+  Eigen::PlainObjectBase<DerivedVertices> &vertices,
+  Eigen::PlainObjectBase<DerivedFaces> &faces)
 {
-  MarchingCubes<Derivedvalues, Derivedpoints, Derivedvertices, DerivedF> mc(values,
-                                       points,
-                                       x_res,
-                                       y_res,
-                                       z_res,
-                                       vertices,
-                                       faces);
+  typedef Eigen::MatrixXi Shim; /* DerivedIndices shim type is unused in this instantiation*/
+  MarchingCubes<DerivedValues, DerivedPoints, DerivedVertices, Shim, DerivedFaces>
+          mc(values, points, x_res, y_res, z_res, 0.0 /*isovalue*/, vertices, faces);
 }
+
+template <typename DerivedValues, typename DerivedPoints, typename DerivedVertices, typename DerivedIndices, typename DerivedFaces>
+IGL_INLINE void igl::copyleft::marching_cubes(
+  const Eigen::MatrixBase<DerivedValues>& values,
+  const Eigen::MatrixBase<DerivedPoints>& points,
+  const Eigen::MatrixBase<DerivedIndices>& indices,
+  const double isovalue,
+  Eigen::PlainObjectBase<DerivedVertices>& vertices,
+  Eigen::PlainObjectBase<DerivedFaces> &faces)
+{
+  MarchingCubes<DerivedValues, DerivedPoints, DerivedVertices, DerivedIndices, DerivedFaces> mc(values, points, indices, isovalue, vertices, faces);
+}
+
+template <typename DerivedValues, typename DerivedPoints, typename DerivedVertices, typename DerivedIndices, typename DerivedFaces>
+IGL_INLINE void igl::copyleft::marching_cubes(
+  const Eigen::MatrixBase<DerivedValues> &values,
+  const Eigen::MatrixBase<DerivedPoints> &points,
+  const Eigen::MatrixBase<DerivedIndices> & indices,
+  Eigen::PlainObjectBase<DerivedVertices> &vertices,
+  Eigen::PlainObjectBase<DerivedFaces> &faces)
+{
+  MarchingCubes<DerivedValues, DerivedPoints, DerivedVertices, DerivedIndices, DerivedFaces> mc(values, points, indices, 0.0 /*isovalue*/, vertices, faces);
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+
 // generated by autoexplicit.sh
-template void igl::copyleft::marching_cubes<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, unsigned int, unsigned int, unsigned int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> >&);
-// generated by autoexplicit.sh
-template void igl::copyleft::marching_cubes<Eigen::Matrix<float, -1, 1, 0, -1, 1>, Eigen::Matrix<float, -1, 3, 0, -1, 3>, Eigen::Matrix<float, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> > const&, unsigned int, unsigned int, unsigned int, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
+template void igl::copyleft::marching_cubes<Eigen::Matrix<float, -1, 1, 0, -1, 1>, Eigen::Matrix<float, -1, 3, 0, -1, 3>, Eigen::Matrix<float, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> > const&, unsigned int, unsigned int, unsigned int, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
 // generated by autoexplicit.sh
-template void igl::copyleft::marching_cubes<Eigen::Matrix<float, -1, 1, 0, -1, 1>, Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, unsigned int, unsigned int, unsigned int, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> >&);
+template void igl::copyleft::marching_cubes<Eigen::Matrix<float, -1, 1, 0, -1, 1>, Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, unsigned int, unsigned int, unsigned int, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> >&);
 // generated by autoexplicit.sh
-template void igl::copyleft::marching_cubes<Eigen::Matrix<float, -1, 1, 0, -1, 1>, Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, unsigned int, unsigned int, unsigned int, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> >&);
+template void igl::copyleft::marching_cubes<Eigen::Matrix<float, -1, 1, 0, -1, 1>, Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, unsigned int, unsigned int, unsigned int, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> >&);
 // generated by autoexplicit.sh
-template void igl::copyleft::marching_cubes<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, unsigned int, unsigned int, unsigned int, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> >&);
-template void igl::copyleft::marching_cubes< Eigen::Matrix<double, -1, 1, 0, -1, 1>, 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> > const&, unsigned int, unsigned int, unsigned int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::copyleft::marching_cubes<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, unsigned int, unsigned int, unsigned int, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> >&);
+template void igl::copyleft::marching_cubes<Eigen::Matrix<double, -1, 1, 0, -1, 1>, 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::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, unsigned int, unsigned int, unsigned int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::copyleft::marching_cubes<Eigen::Matrix<double, -1, 1, 0, -1, 1>, 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::Matrix<int, -1, -1, 0, -1, -1> >(const Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, const Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::copyleft::marching_cubes<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(const Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, const Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, const Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::copyleft::marching_cubes<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(const Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, const Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::copyleft::marching_cubes<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(const Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, const Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, const Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::copyleft::marching_cubes<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 8, 0, -1, 8>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(const Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, const Eigen::MatrixBase<Eigen::Matrix<int, -1, 8, 0, -1, 8> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::copyleft::marching_cubes<Eigen::Matrix<double, -1, 1, 0, -1, 1>,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::Matrix<int, -1, 3, 0, -1, 3> >(const Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&,const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&,const Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&,Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&,Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
+template void igl::copyleft::marching_cubes<Eigen::Matrix<double, -1, 1, 0, -1, 1>,Eigen::Matrix<double, -1, -1, 0, -1, -1>,Eigen::Matrix<double, -1, -1, 0, -1, -1>,Eigen::Matrix<int, -1, 8, 0, -1, 8>,Eigen::Matrix<int, -1, 3, 0, -1, 3> >(const Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&,const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&,const Eigen::MatrixBase<Eigen::Matrix<int, -1, 8, 0, -1, 8> >&,Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&,Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
+template void igl::copyleft::marching_cubes<Eigen::Matrix<double, -1, 1, 0, -1, 1>,Eigen::Matrix<double, -1, 3, 0, -1, 3>,Eigen::Matrix<double, -1, 3, 0, -1, 3>,Eigen::Matrix<int, -1, 8, 0, -1, 8>,Eigen::Matrix<int, -1, 3, 0, -1, 3> >(const Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&,const Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&,const Eigen::MatrixBase<Eigen::Matrix<int, -1, 8, 0, -1, 8> >&,Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&,Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
 #endif

include/igl/copyleft/marching_cubes.h

@@ -14,10 +14,10 @@ namespace igl
 {
   namespace copyleft
   {
-    // marching_cubes( values, points, x_res, y_res, z_res, vertices, faces )
+    // marching_cubes( values, points, x_res, y_res, z_res, isovalue, vertices, faces )
     //
-    // performs marching cubes reconstruction on the grid defined by values, and
-    // points, and generates vertices and faces
+    // performs marching cubes reconstruction on a grid defined by values, and
+    // points, and generates a mesh defined by vertices and faces
     //
     // Input:
     //  values  #number_of_grid_points x 1 array -- the scalar values of an
@@ -34,24 +34,75 @@ namespace igl
     //  xres  resolutions of the grid in x dimension
     //  yres  resolutions of the grid in y dimension
     //  zres  resolutions of the grid in z dimension
+    //  isovalue  the isovalue of the surface to reconstruct
     // Output:
     //   vertices  #V by 3 list of mesh vertex positions
     //   faces  #F by 3 list of mesh triangle indices
     //
-    template <
-      typename Derivedvalues, 
-      typename Derivedpoints, 
-      typename Derivedvertices, 
-      typename DerivedF>
-      IGL_INLINE void marching_cubes(
-        const Eigen::PlainObjectBase<Derivedvalues> &values,
-        const Eigen::PlainObjectBase<Derivedpoints> &points,
+    template <typename DerivedValues, typename DerivedPoints, typename DerivedVertices, typename DerivedFaces>
+    IGL_INLINE void marching_cubes(
+        const Eigen::MatrixBase<DerivedValues> &values,
+        const Eigen::MatrixBase<DerivedPoints> &points,
         const unsigned x_res,
         const unsigned y_res,
         const unsigned z_res,
-        Eigen::PlainObjectBase<Derivedvertices> &vertices,
-        Eigen::PlainObjectBase<DerivedF> &faces);
+        const double isovalue,
+        Eigen::PlainObjectBase<DerivedVertices> &vertices,
+        Eigen::PlainObjectBase<DerivedFaces> &faces);
+
+    // Overload of the above function where the isovalue defaults to 0.0
+    template <typename DerivedValues, typename DerivedPoints, typename DerivedVertices, typename DerivedFaces>
+    IGL_INLINE void marching_cubes(
+      const Eigen::MatrixBase<DerivedValues> &values,
+      const Eigen::MatrixBase<DerivedPoints> &points,
+      const unsigned x_res,
+      const unsigned y_res,
+      const unsigned z_res,
+      Eigen::PlainObjectBase<DerivedVertices> &vertices,
+      Eigen::PlainObjectBase<DerivedFaces> &faces);
+
+
+
+    // marching_cubes( values, points, indices, vertices, faces )
+    //
+    // Perform marching cubes reconstruction on the grid cells defined by indices.
+    // The indices parameter is an nx8 dense array of index values into the points and values arrays.
+    // Each row of indices represents a cube for which to generate vertices and faces over.
+    //
+    // Input:
+    //  values  #number_of_grid_points x 1 array -- the scalar values of an
+    //    implicit function defined on the grid points (<0 in the inside of the
+    //    surface, 0 on the border, >0 outside)
+    //  points  #number_of_grid_points x 3 array -- 3-D positions of the grid
+    //    points, ordered in x,y,z order:
+    //  indices  #cubes x 8 array -- one row for each cube where each value is
+    //    the index of a vertex in points and a scalar in values.
+    //    i.e. points[indices[i, j]] = the position of the j'th vertex of the i'th cube
+    // Output:
+    //   vertices  #V by 3 list of mesh vertex positions
+    //   faces  #F by 3 list of mesh triangle indices
+    // Note: The winding direction of the cube indices will affect the output winding of the faces
+    //
+    template <typename DerivedValues, typename DerivedPoints, typename DerivedVertices, typename DerivedIndices, typename DerivedFaces>
+    IGL_INLINE void marching_cubes(
+      const Eigen::MatrixBase<DerivedValues> &values,
+      const Eigen::MatrixBase<DerivedPoints> &points,
+      const Eigen::MatrixBase<DerivedIndices> &indices,
+      const double isovalue,
+      Eigen::PlainObjectBase<DerivedVertices> &vertices,
+      Eigen::PlainObjectBase<DerivedFaces> &faces);
+
+    // Overload of the above function where the isovalue defaults to 0.0
+    template <typename DerivedValues, typename DerivedPoints, typename DerivedVertices, typename DerivedIndices, typename DerivedFaces>
+    IGL_INLINE void marching_cubes(
+      const Eigen::MatrixBase<DerivedValues> &values,
+      const Eigen::MatrixBase<DerivedPoints> &points,
+      const Eigen::MatrixBase<DerivedIndices> &indices,
+      Eigen::PlainObjectBase<DerivedVertices> &vertices,
+      Eigen::PlainObjectBase<DerivedFaces> &faces);
+
   }
+
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/cut_to_disk.cpp

@@ -0,0 +1,330 @@
+#include "cut_to_disk.h"
+
+#include <map>
+#include <set>
+#include <deque>
+#include <algorithm>
+
+namespace igl {
+  template <typename DerivedF, typename Index>
+  void cut_to_disk(
+    const Eigen::MatrixBase<DerivedF> &F,
+    std::vector<std::vector<Index> > &cuts) 
+  {
+    cuts.clear();
+
+    Index nfaces = F.rows();
+
+    if (nfaces == 0)
+        return;
+
+    std::map<std::pair<Index, Index>, std::vector<Index> > edges;
+    // build edges
+    
+    for (Index i = 0; i < nfaces; i++)
+    {
+        for (int j = 0; j < 3; j++)
+        {
+            Index v0 = F(i, j);
+            Index v1 = F(i, (j + 1) % 3);
+            std::pair<Index, Index> e;
+            e.first = std::min(v0, v1);
+            e.second = std::max(v0, v1);
+            edges[e].push_back(i);
+        }
+    }
+
+    int nedges = edges.size();
+    Eigen::Matrix<Index, -1, -1> edgeVerts(nedges,2);
+    Eigen::Matrix<Index, -1, -1> edgeFaces(nedges,2);
+    Eigen::Matrix<Index, -1, -1> faceEdges(nfaces, 3);
+    std::set<Index> boundaryEdges;
+    std::map<std::pair<Index, Index>, Index> edgeidx;
+    Index idx = 0;
+    for (auto it : edges)
+    {
+        edgeidx[it.first] = idx;
+        edgeVerts(idx, 0) = it.first.first;
+        edgeVerts(idx, 1) = it.first.second;
+        edgeFaces(idx, 0) = it.second[0];
+        if (it.second.size() > 1)
+        {
+            edgeFaces(idx, 1) = it.second[1];
+        }
+        else
+        {
+            edgeFaces(idx, 1) = -1;
+            boundaryEdges.insert(idx);
+        }
+        idx++;
+    }
+    for (Index i = 0; i < nfaces; i++)
+    {
+        for (int j = 0; j < 3; j++)
+        {
+            Index v0 = F(i, j);
+            Index v1 = F(i, (j + 1) % 3);
+            std::pair<Index, Index> e;
+            e.first = std::min(v0, v1);
+            e.second = std::max(v0, v1);
+            faceEdges(i, j) = edgeidx[e];
+        }
+    }
+    
+    bool *deleted = new bool[nfaces];
+    for (Index i = 0; i < nfaces; i++)
+        deleted[i] = false;
+
+    std::set<Index> deletededges;
+    
+    // loop over faces
+    for (Index face = 0; face < nfaces; face++)
+    {
+        // stop at first undeleted face
+        if (deleted[face])
+            continue;
+        deleted[face] = true;
+        std::deque<Index> processEdges;
+        for (int i = 0; i < 3; i++)
+        {
+            Index e = faceEdges(face, i);
+            if (boundaryEdges.count(e))
+                continue;
+            int ndeleted = 0;
+            if (deleted[edgeFaces(e, 0)])
+                ndeleted++;
+            if (deleted[edgeFaces(e, 1)])
+                ndeleted++;
+            if (ndeleted == 1)
+                processEdges.push_back(e);
+        }
+        // delete all faces adjacent to edges with exactly one adjacent face
+        while (!processEdges.empty())
+        {
+            Index nexte = processEdges.front();
+            processEdges.pop_front();
+            Index todelete = nfaces;
+            if (!deleted[edgeFaces(nexte, 0)])
+                todelete = edgeFaces(nexte, 0);
+            if (!deleted[edgeFaces(nexte, 1)])
+                todelete = edgeFaces(nexte, 1);
+            if (todelete != nfaces)
+            {
+                deletededges.insert(nexte);
+                deleted[todelete] = true;
+                for (int i = 0; i < 3; i++)
+                {
+                    Index e = faceEdges(todelete, i);
+                    if (boundaryEdges.count(e))
+                        continue;
+                    int ndeleted = 0;
+                    if (deleted[edgeFaces(e, 0)])
+                        ndeleted++;
+                    if (deleted[edgeFaces(e, 1)])
+                        ndeleted++;
+                    if (ndeleted == 1)
+                        processEdges.push_back(e);
+                }
+            }
+        }
+    }
+    delete[] deleted;
+
+    // accumulated non-deleted edges
+    std::vector<Index> leftedges;
+    for (Index i = 0; i < nedges; i++)
+    {
+        if (!deletededges.count(i))
+            leftedges.push_back(i);
+    }
+
+    deletededges.clear();
+    // prune spines
+    std::map<Index, std::vector<Index> > spinevertedges;
+    for (Index i : leftedges)
+    {
+        spinevertedges[edgeVerts(i, 0)].push_back(i);
+        spinevertedges[edgeVerts(i, 1)].push_back(i);
+    }
+    
+    std::deque<Index> vertsProcess;
+    std::map<Index, int> spinevertnbs;
+    for (auto it : spinevertedges)
+    {
+        spinevertnbs[it.first] = it.second.size();
+        if (it.second.size() == 1)
+            vertsProcess.push_back(it.first);
+    }
+    while (!vertsProcess.empty())
+    {
+        Index vert = vertsProcess.front();
+        vertsProcess.pop_front();
+        for (Index e : spinevertedges[vert])
+        {
+            if (!deletededges.count(e))
+            {
+                deletededges.insert(e);
+                for (int j = 0; j < 2; j++)
+                {
+                    spinevertnbs[edgeVerts(e, j)]--;
+                    if (spinevertnbs[edgeVerts(e, j)] == 1)
+                    {
+                        vertsProcess.push_back(edgeVerts(e, j));
+                    }
+                }
+            }
+        }
+    }
+    std::vector<Index> loopedges;
+    for (Index i : leftedges)
+        if (!deletededges.count(i))
+            loopedges.push_back(i);
+
+    Index nloopedges = loopedges.size();
+    if (nloopedges == 0)
+        return;
+
+    std::map<Index, std::vector<Index> > loopvertedges;
+    for (Index e : loopedges)
+    {
+        loopvertedges[edgeVerts(e, 0)].push_back(e);
+        loopvertedges[edgeVerts(e, 1)].push_back(e);
+    }
+    
+    std::set<Index> usededges;
+    for (Index e : loopedges)
+    {
+        // make a cycle or chain starting from this edge
+        while (!usededges.count(e))
+        {
+            std::vector<Index> cycleverts;
+            std::vector<Index> cycleedges;
+            cycleverts.push_back(edgeVerts(e, 0));
+            cycleverts.push_back(edgeVerts(e, 1));
+            cycleedges.push_back(e);
+
+            std::map<Index, Index> cycleidx;
+            cycleidx[cycleverts[0]] = 0;
+            cycleidx[cycleverts[1]] = 1;
+
+            Index curvert = edgeVerts(e, 1);
+            Index cure = e;
+            bool foundcycle = false;
+            while (curvert != -1 && !foundcycle)
+            {
+                Index nextvert = -1;
+                Index nexte = -1;
+                for (Index cande : loopvertedges[curvert])
+                {
+                    if (!usededges.count(cande) && cande != cure)
+                    {
+                        int vidx = 0;
+                        if (curvert == edgeVerts(cande, vidx))
+                            vidx = 1;
+                        nextvert = edgeVerts(cande, vidx);
+                        nexte = cande;
+                        break;
+                    }
+                }
+                if (nextvert != -1)
+                {
+                    auto it = cycleidx.find(nextvert);
+                    if (it != cycleidx.end())
+                    {
+                        // we've hit outselves
+                        std::vector<Index> cut;
+                        for (Index i = it->second; i < cycleverts.size(); i++)
+                        {
+                            cut.push_back(cycleverts[i]);
+                        }
+                        cut.push_back(nextvert);
+                        cuts.push_back(cut);
+                        for (Index i = it->second; i < cycleedges.size(); i++)
+                        {
+                            usededges.insert(cycleedges[i]);
+                        }
+                        usededges.insert(nexte);
+                        foundcycle = true;
+                    }
+                    else
+                    {
+                        cycleidx[nextvert] = cycleverts.size();
+                        cycleverts.push_back(nextvert);
+                        cycleedges.push_back(nexte);                        
+                    }
+                }                
+                curvert = nextvert;
+                cure = nexte;
+            }
+            if (!foundcycle)
+            {
+                // we've hit a dead end. reverse and try the other direction
+                std::reverse(cycleverts.begin(), cycleverts.end());
+                std::reverse(cycleedges.begin(), cycleedges.end());
+                curvert = cycleverts.back();
+                cure = cycleedges.back();
+                while (curvert != -1 && !foundcycle)
+                {
+                    Index nextvert = -1;
+                    Index nexte = -1;
+                    for (Index cande : loopvertedges[curvert])
+                    {
+                        if (!usededges.count(cande) && cande != cure)
+                        {
+                            int vidx = 0;
+                            if (curvert == edgeVerts(cande, vidx))
+                                vidx = 1;
+                            nextvert = edgeVerts(cande, vidx);
+                            nexte = cande;
+                            break;
+                        }
+                    }
+                    if (nextvert != -1)
+                    {
+                        auto it = cycleidx.find(nextvert);
+                        if (it != cycleidx.end())
+                        {
+                            // we've hit outselves
+                            std::vector<int> cut;
+                            for (Index i = it->second; i < cycleverts.size(); i++)
+                            {
+                                cut.push_back(cycleverts[i]);
+                            }
+                            cut.push_back(nextvert);
+                            cuts.push_back(cut);
+                            for (Index i = it->second; i < cycleedges.size(); i++)
+                            {
+                                usededges.insert(cycleedges[i]);
+                            }
+                            usededges.insert(nexte);
+                            foundcycle = true;
+                        }
+                        else
+                        {
+                            cycleidx[nextvert] = cycleverts.size();
+                            cycleverts.push_back(nextvert);
+                            cycleedges.push_back(nexte);                        
+                        }
+                    }                
+                    curvert = nextvert;
+                    cure = nexte;
+                }
+                if (!foundcycle)
+                {
+                    // we've found a chain
+                    std::vector<Index> cut;
+                    for (Index i = 0; i < cycleverts.size(); i++)
+                    {
+                        cut.push_back(cycleverts[i]);
+                    }
+                    cuts.push_back(cut);
+                    for (Index i = 0; i < cycleedges.size(); i++)
+                    {
+                        usededges.insert(cycleedges[i]);
+                    }                    
+                }
+            }
+        }
+    }
+  }
+}

include/igl/cut_to_disk.h

@@ -0,0 +1,52 @@
+#ifndef IGL_CUT_TO_DISK_H
+#define IGL_CUT_TO_DISK_H
+#include "igl_inline.h"
+
+#include <Eigen/Core>
+
+#include <vector>
+
+namespace igl
+{
+  // Given a triangle mesh, computes a set of edge cuts sufficient to carve the 
+  // mesh into a topological disk, without disconnecting any connected components.
+  // Nothing else about the cuts (including number, total length, or smoothness)
+  // is guaranteed to be optimal.
+  //
+  // Simply-connected components without boundary (topological spheres) are left
+  // untouched (delete any edge if you really want a disk). 
+  // All other connected components are cut into disks. Meshes with boundary are
+  // supported; boundary edges will be included as cuts.
+  //
+  // The cut mesh itself can be materialized using cut_mesh().
+  //
+  // Implements the triangle-deletion approach described by Gu et al's
+  // "Geometry Images."
+  //
+  // Template Parameters:
+  //   Index  Integrable type large enough to represent the total number of faces
+  //     and edges in the surface represented by F, and all entries of F.
+  //
+  // Inputs:
+  //   F  #F by 3 list of the faces (must be triangles)
+  //
+  // Outputs:
+  //   cuts  List of cuts. Each cut is a sequence of vertex indices (where
+  //     pairs of consecutive vertices share a face), is simple, and is either
+  //     a closed loop (in which the first and last indices are identical) or
+  //     an open curve. Cuts are edge-disjoint.
+  //
+  
+  template <
+    typename DerivedF,
+    typename Index>
+  IGL_INLINE void cut_to_disk(
+    const Eigen::MatrixBase<DerivedF> &F,
+    std::vector<std::vector<Index> > &cuts);    
+};
+
+#ifndef IGL_STATIC_LIBRARY
+#include "cut_to_disk.cpp"
+#endif
+
+#endif

include/igl/dijkstra.cpp

@@ -5,14 +5,16 @@
 // This Source Code Form is subject to the terms of the Mozilla Public License 
 // v. 2.0. If a copy of the MPL was not distributed with this file, You can 
 // obtain one at http://mozilla.org/MPL/2.0/.
-#include <igl/dijkstra.h>
+#include "dijkstra.h"
 
 template <typename IndexType, typename DerivedD, typename DerivedP>
-IGL_INLINE int igl::dijkstra_compute_paths(const IndexType &source,
-                                           const std::set<IndexType> &targets,
-                                           const std::vector<std::vector<IndexType> >& VV,
-                                           Eigen::PlainObjectBase<DerivedD> &min_distance,
-                                           Eigen::PlainObjectBase<DerivedP> &previous)
+IGL_INLINE int igl::dijkstra(
+  const IndexType &source,
+  const std::set<IndexType> &targets,
+  const std::vector<std::vector<IndexType> >& VV,
+  const std::vector<double>& weights,
+  Eigen::PlainObjectBase<DerivedD> &min_distance,
+  Eigen::PlainObjectBase<DerivedP> &previous)
 {
   int numV = VV.size();
   min_distance.setConstant(numV, 1, std::numeric_limits<typename DerivedD::Scalar>::infinity());
@@ -37,7 +39,7 @@ IGL_INLINE int igl::dijkstra_compute_paths(const IndexType &source,
          neighbor_iter++)
     {
       IndexType v = *neighbor_iter;
-      typename DerivedD::Scalar distance_through_u = dist + 1.;
+      typename DerivedD::Scalar distance_through_u = dist + weights[u];
       if (distance_through_u < min_distance[v]) {
         vertex_queue.erase(std::make_pair(min_distance[v], v));
 
@@ -53,10 +55,23 @@ IGL_INLINE int igl::dijkstra_compute_paths(const IndexType &source,
   return -1;
 }
 
+template <typename IndexType, typename DerivedD, typename DerivedP>
+IGL_INLINE int igl::dijkstra(
+  const IndexType &source,
+  const std::set<IndexType> &targets,
+  const std::vector<std::vector<IndexType> >& VV,
+  Eigen::PlainObjectBase<DerivedD> &min_distance,
+  Eigen::PlainObjectBase<DerivedP> &previous)
+{
+  std::vector<double> weights(VV.size(), 1.0);
+  return dijkstra(source, targets, VV, weights, min_distance, previous);
+}
+
 template <typename IndexType, typename DerivedP>
-IGL_INLINE void igl::dijkstra_get_shortest_path_to(const IndexType &vertex,
-                                                   const Eigen::PlainObjectBase<DerivedP> &previous,
-                                                   std::vector<IndexType> &path)
+IGL_INLINE void igl::dijkstra(
+  const IndexType &vertex,
+  const Eigen::MatrixBase<DerivedP> &previous,
+  std::vector<IndexType> &path)
 {
   IndexType source = vertex;
   path.clear();
@@ -66,6 +81,6 @@ IGL_INLINE void igl::dijkstra_get_shortest_path_to(const IndexType &vertex,
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-template int igl::dijkstra_compute_paths<int, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(int const&, std::set<int, std::less<int>, std::allocator<int> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template void igl::dijkstra_get_shortest_path_to<int, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(int const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, std::vector<int, std::allocator<int> >&);
+template int igl::dijkstra<int, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(int const&, std::set<int, std::less<int>, std::allocator<int> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::dijkstra<int, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(int const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, std::vector<int, std::allocator<int> >&);
 #endif

include/igl/dijkstra.h

@@ -16,7 +16,7 @@
 
 namespace igl {
 
-  // Dijstra's algorithm for shortest paths, with multiple targets.
+  // Dijkstra's algorithm for vertex-weighted shortest paths, with multiple targets.
   // Adapted from http://rosettacode.org/wiki/Dijkstra%27s_algorithm .
   //
   // Inputs:
@@ -24,19 +24,43 @@ namespace igl {
   //   targets          target vector set
   //   VV               #V list of lists of incident vertices (adjacency list), e.g.
   //                    as returned by igl::adjacency_list
+  //   weights          #V list of scalar vertex weights
   //
   // Output:
   //   min_distance     #V by 1 list of the minimum distances from source to all vertices
   //   previous         #V by 1 list of the previous visited vertices (for each vertex) - used for backtracking
   //
   template <typename IndexType, typename DerivedD, typename DerivedP>
-  IGL_INLINE int dijkstra_compute_paths(const IndexType &source,
-                                        const std::set<IndexType> &targets,
-                                        const std::vector<std::vector<IndexType> >& VV,
-                                        Eigen::PlainObjectBase<DerivedD> &min_distance,
-                                        Eigen::PlainObjectBase<DerivedP> &previous);
+  IGL_INLINE int dijkstra(
+    const IndexType &source,
+    const std::set<IndexType> &targets,
+    const std::vector<std::vector<IndexType> >& VV,
+    const std::vector<double>& weights,
+    Eigen::PlainObjectBase<DerivedD> &min_distance,
+    Eigen::PlainObjectBase<DerivedP> &previous);
 
-  // Backtracking after Dijstra's algorithm, to find shortest path.
+  // Dijkstra's algorithm for shortest paths, with multiple targets.
+  // Adapted from http://rosettacode.org/wiki/Dijkstra%27s_algorithm .
+  //
+  // Inputs:
+  //   source           index of source vertex
+  //   targets          target vector set
+  //   VV               #V list of lists of incident vertices (adjacency list), e.g.
+  //                    as returned by igl::adjacency_list
+  //
+  // Output:
+  //   min_distance     #V by 1 list of the minimum distances from source to all vertices
+  //   previous         #V by 1 list of the previous visited vertices (for each vertex) - used for backtracking
+  //
+  template <typename IndexType, typename DerivedD, typename DerivedP>
+  IGL_INLINE int dijkstra(
+    const IndexType &source,
+    const std::set<IndexType> &targets,
+    const std::vector<std::vector<IndexType> >& VV,
+    Eigen::PlainObjectBase<DerivedD> &min_distance,
+    Eigen::PlainObjectBase<DerivedP> &previous);
+
+  // Backtracking after Dijkstra's algorithm, to find shortest path.
   //
   // Inputs:
   //   vertex           vertex to which we want the shortest path (from same source as above)
@@ -46,9 +70,10 @@ namespace igl {
   //   path             #P by 1 list of vertex indices in the shortest path from source to vertex
   //
   template <typename IndexType, typename DerivedP>
-  IGL_INLINE void dijkstra_get_shortest_path_to(const IndexType &vertex,
-                                                const Eigen::PlainObjectBase<DerivedP> &previous,
-                                                std::vector<IndexType> &path);
+  IGL_INLINE void dijkstra(
+    const IndexType &vertex,
+    const Eigen::MatrixBase<DerivedP> &previous,
+    std::vector<IndexType> &path);
 };
 
 

include/igl/knn.cpp

@@ -102,4 +102,7 @@ namespace igl {
 
 
 
+#ifdef IGL_STATIC_LIBRARY
+template void igl::knn<Eigen::Matrix<double, -1, -1, 0, -1, -1>, int, int, Eigen::Matrix<int, -1, 8, 0, -1, 8>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 8, 0, -1, 8> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+#endif
 

include/igl/polygon_mesh_to_triangle_mesh.cpp

@@ -73,4 +73,5 @@ template void igl::polygon_mesh_to_triangle_mesh<int, Eigen::Matrix<unsigned int
 template void igl::polygon_mesh_to_triangle_mesh<int, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
 template void igl::polygon_mesh_to_triangle_mesh<int, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::polygon_mesh_to_triangle_mesh<int, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> >&);
+template void igl::polygon_mesh_to_triangle_mesh<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> >&);
 #endif

include/igl/remove_unreferenced.cpp

@@ -119,4 +119,5 @@ template void igl::remove_unreferenced<Eigen::Matrix<double, -1, 3, 0, -1, 3>, E
 template void igl::remove_unreferenced<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, 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::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > 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> >&);
 template void igl::remove_unreferenced<Eigen::Matrix<double, -1, 3, 0, -1, 3>, 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::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<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::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::remove_unreferenced<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::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<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::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::remove_unreferenced<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::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<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::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 #endif