diff --git a/include/igl/marching_tets.cpp b/include/igl/marching_tets.cpp
new file mode 100644
index 000000000..02fabb014
--- /dev/null
+++ b/include/igl/marching_tets.cpp
@@ -0,0 +1,194 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2018 Francis Williams <francis@fwilliams.info>
+//
+// 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 "marching_tets.h"
+
+#include <unordered_map>
+#include <vector>
+#include <utility>
+#include <cstdint>
+#include <iostream>
+
+template <typename DerivedTV,
+          typename DerivedTT,
+          typename DerivedS,
+          typename DerivedSV,
+          typename DerivedSF,
+          typename DerivedJ,
+          typename BCType>
+void igl::marching_tets(
+    const Eigen::PlainObjectBase<DerivedTV>& TV,
+    const Eigen::PlainObjectBase<DerivedTT>& TT,
+    const Eigen::PlainObjectBase<DerivedS>& isovals,
+    double isovalue,
+    Eigen::PlainObjectBase<DerivedSV>& outV,
+    Eigen::PlainObjectBase<DerivedSF>& outF,
+    Eigen::PlainObjectBase<DerivedJ>& J,
+    Eigen::SparseMatrix<BCType>& BC)
+{
+  using namespace std;
+
+  // We're hashing edges to deduplicate using 64 bit ints. The upper and lower
+  // 32 bits of a key are the indices of vertices in the mesh. The implication is
+  // that you can only have 2^32 vertices which I have deemed sufficient for
+  // anything reasonable.
+  const auto make_edge_key = [](const pair<int32_t, int32_t>& p) -> int64_t
+  {
+    std::int64_t ret = 0;
+    ret |= p.first;
+    ret |= static_cast<std::int64_t>(p.second) << 32;
+    return ret;
+  };
+
+  const int mt_cell_lookup[16][4] =
+  {
+    { -1, -1, -1, -1 },
+    {  0,  2,  1, -1 },
+    {  0,  3,  4, -1 },
+    {  2,  1,  3,  4 },
+    {  5,  3,  1, -1 },
+    {  0,  2,  5,  3 },
+    {  0,  1,  5,  4 },
+    {  2,  5,  4, -1 },
+    {  4,  5,  2, -1 },
+    {  0,  4,  5,  1 },
+    {  0,  3,  5,  2 },
+    {  1,  3,  5, -1 },
+    {  4,  3,  1,  2 },
+    {  0,  4,  3, -1 },
+    {  0,  1,  2, -1 },
+    { -1, -1, -1, -1 },
+  };
+
+  const int mt_edge_lookup[6][2] =
+  {
+    {0, 1},
+    {0, 2},
+    {0, 3},
+    {1, 2},
+    {1, 3},
+    {2, 3},
+  };
+
+  // Store the faces and the tet they are in
+  vector<pair<Eigen::RowVector3i, int>> faces;
+
+  // Store the edges in the tet mesh which we add vertices on
+  // so we can deduplicate
+  vector<pair<int, int>> edge_table;
+
+
+  assert(TT.cols() == 4 && TT.rows() >= 1);
+  assert(TV.cols() == 3 && TV.rows() >= 4);
+  assert(isovals.cols() == 1);
+
+  // For each tet
+  for (int i = 0; i < TT.rows(); i++)
+  {
+    uint8_t key = 0;
+    for (int v = 0; v < 4; v++)
+    {
+      const int vid = TT(i, v);
+      const uint8_t flag = isovals[vid] > isovalue;
+      key |= flag << v;
+    }
+
+    // This will contain the index in TV of each vertex in the tet
+    int v_ids[4] = {-1, -1, -1, -1};
+
+    // Insert any vertices if the tet intersects the level surface
+    for (int e = 0; e < 4 && mt_cell_lookup[key][e] != -1; e++)
+    {
+      const int tv1_idx = TT(i, mt_edge_lookup[mt_cell_lookup[key][e]][0]);
+      const int tv2_idx = TT(i, mt_edge_lookup[mt_cell_lookup[key][e]][1]);
+      const int vertex_id = edge_table.size();
+      edge_table.push_back(make_pair(min(tv1_idx, tv2_idx), max(tv1_idx, tv2_idx)));
+      v_ids[e] = vertex_id;
+    }
+
+    // Insert the corresponding faces
+    if (v_ids[0] != -1)
+    {
+      bool is_quad = mt_cell_lookup[key][3] != -1;
+      if (is_quad)
+      {
+        const Eigen::RowVector3i f1(v_ids[0], v_ids[1], v_ids[3]);
+        const Eigen::RowVector3i f2(v_ids[1], v_ids[2], v_ids[3]);
+        faces.push_back(make_pair(f1, i));
+        faces.push_back(make_pair(f2, i));
+      }
+      else
+      {
+        const Eigen::RowVector3i f(v_ids[0], v_ids[1], v_ids[2]);
+        faces.push_back(make_pair(f, i));
+      }
+
+    }
+  }
+
+  int num_unique = 0;
+  outV.resize(edge_table.size(), 3);
+  outF.resize(faces.size(), 3);
+  J.resize(faces.size());
+
+  // Sparse matrix triplets for BC
+  vector<Eigen::Triplet<BCType>> bc_triplets;
+  bc_triplets.reserve(edge_table.size());
+
+  // Deduplicate vertices
+  unordered_map<int64_t, int> emap;
+  emap.max_load_factor(0.5);
+  emap.reserve(edge_table.size());
+
+  for (int f = 0; f < faces.size(); f++)
+  {
+    for (int v = 0; v < 3; v++)
+    {
+      const int vi = faces[f].first[v];
+      const int ti = faces[f].second;
+      const pair<int32_t, int32_t> edge = edge_table[vi];
+      const int64_t key = make_edge_key(edge);
+      auto it = emap.find(key);
+      if (it == emap.end()) // New unique vertex, insert it
+      {
+        // Typedef to make sure we handle floats properly
+        typedef Eigen::Matrix<typename DerivedTV::Scalar, 1, 3, Eigen::RowMajor, 1, 3> RowVector;
+        const RowVector v1 = TV.row(edge.first);
+        const RowVector v2 = TV.row(edge.second);
+        const double a = fabs(isovals[edge.first] - isovalue);
+        const double b = fabs(isovals[edge.second] - isovalue);
+        const double w = a / (a+b);
+
+        // Create a casted copy in case BCType is a float and we need to downcast
+        const BCType bc_w = static_cast<BCType>(w);
+        bc_triplets.push_back(Eigen::Triplet<BCType>(num_unique, edge.first, 1-bc_w));
+        bc_triplets.push_back(Eigen::Triplet<BCType>(num_unique, edge.second, bc_w));
+
+        // Create a casted copy in case DerivedTV::Scalar is a float and we need to downcast
+        const typename DerivedTV::Scalar v_w = static_cast<typename DerivedTV::Scalar>(w);
+        outV.row(num_unique) = (1-v_w)*v1 + v_w*v2;
+        outF(f, v) = num_unique;
+        J[f] = ti;
+
+        emap.emplace(key, num_unique);
+        num_unique += 1;
+      } else {
+        outF(f, v) = it->second;
+      }
+    }
+  }
+  outV.conservativeResize(num_unique, 3);
+  J.conservativeResize(num_unique, 1);
+  BC.resize(num_unique, TV.rows());
+  BC.setFromTriplets(bc_triplets.begin(), bc_triplets.end());
+}
+
+
+#ifdef IGL_STATIC_LIBRARY
+template void igl::marching_tets<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<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, double>(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> > const&, double, 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::SparseMatrix<double, 0, int>&);
+#endif // IGL_STATIC_LIBRARY
diff --git a/include/igl/marching_tets.h b/include/igl/marching_tets.h
new file mode 100644
index 000000000..9f3fa5435
--- /dev/null
+++ b/include/igl/marching_tets.h
@@ -0,0 +1,196 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2018 Francis Williams <francis@fwilliams.info>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef IGL_MARCHING_TETS_H
+#define IGL_MARCHING_TETS_H
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+
+namespace igl {
+  // marching_tets( TV, TT, S, isovalue, SV, SF, J, BC)
+  //
+  // performs the marching tetrahedra algorithm on a tet mesh defined by TV and
+  // TT with scalar values defined at each vertex in TV. The output is a
+  // triangle mesh approximating the isosurface coresponding to the value
+  // isovalue.
+  //
+  // Input:
+  //  TV  #tet_vertices x 3 array -- The vertices of the tetrahedral mesh
+  //  TT  #tets x 4 array --  The indexes of each tet in the tetrahedral mesh
+  //  S  #tet_vertices x 1 array -- The values defined on each tet vertex
+  //  isovalue  scalar -- The isovalue of the level set we want to compute
+  //
+  // Output:
+  //  SV  #SV x 3 array -- The vertices of the output level surface mesh
+  //  SF  #SF x 3 array -- The face indexes of the output level surface mesh
+  //  J   #SF list of indices into TT revealing which tet each face comes from
+  //  BC  #SV x #TV list of barycentric coordinates so that SV = BC*TV
+  template <typename DerivedTV,
+            typename DerivedTT,
+            typename DerivedS,
+            typename DerivedSV,
+            typename DerivedSF,
+            typename DerivedJ,
+            typename BCType>
+  IGL_INLINE void marching_tets(
+      const Eigen::PlainObjectBase<DerivedTV>& TV,
+      const Eigen::PlainObjectBase<DerivedTT>& TT,
+      const Eigen::PlainObjectBase<DerivedS>& S,
+      double isovalue,
+      Eigen::PlainObjectBase<DerivedSV>& SV,
+      Eigen::PlainObjectBase<DerivedSF>& SF,
+      Eigen::PlainObjectBase<DerivedJ>& J,
+      Eigen::SparseMatrix<BCType>& BC);
+
+  // marching_tets( TV, TT, S, SV, SF, J, BC)
+  //
+  // Performs the marching tetrahedra algorithm on a tet mesh defined by TV and
+  // TT with scalar values defined at each vertex in TV. The output is a
+  // triangle mesh approximating the isosurface coresponding to an isovalue of 0.
+  //
+  // Input:
+  //  TV  #tet_vertices x 3 array -- The vertices of the tetrahedral mesh
+  //  TT  #tets x 4 array --  The indexes of each tet in the tetrahedral mesh
+  //  S  #tet_vertices x 1 array -- The values defined on each tet vertex
+  //  isovalue  scalar -- The isovalue of the level set we want to compute
+  //
+  // Output:
+  //  SV  #SV x 3 array -- The vertices of the output level surface mesh
+  //  SF  #SF x 3 array -- The face indexes of the output level surface mesh
+  //  J   #SF list of indices into TT revealing which tet each face comes from
+  //  BC  #SV x #TV list of barycentric coordinates so that SV = BC*TV
+  template <typename DerivedTV,
+            typename DerivedTT,
+            typename DerivedS,
+            typename DerivedSV,
+            typename DerivedSF,
+            typename DerivedJ,
+            typename BCType>
+  IGL_INLINE void marching_tets(
+      const Eigen::PlainObjectBase<DerivedTV>& TV,
+      const Eigen::PlainObjectBase<DerivedTT>& TT,
+      const Eigen::PlainObjectBase<DerivedS>& S,
+      Eigen::PlainObjectBase<DerivedSV>& SV,
+      Eigen::PlainObjectBase<DerivedSF>& SF,
+      Eigen::PlainObjectBase<DerivedJ>& J,
+      Eigen::SparseMatrix<BCType>& BC) {
+    return igl::marching_tets(TV, TT, S, 0.0, SV, SF, J, BC);
+  }
+
+  // marching_tets( TV, TT, S, isovalue, SV, SF, J)
+  //
+  // performs the marching tetrahedra algorithm on a tet mesh defined by TV and
+  // TT with scalar values defined at each vertex in TV. The output is a
+  // triangle mesh approximating the isosurface coresponding to the value
+  // isovalue.
+  //
+  // Input:
+  //  TV  #tet_vertices x 3 array -- The vertices of the tetrahedral mesh
+  //  TT  #tets x 4 array --  The indexes of each tet in the tetrahedral mesh
+  //  S  #tet_vertices x 1 array -- The values defined on each tet vertex
+  //  isovalue  scalar -- The isovalue of the level set we want to compute
+  //
+  // Output:
+  //  SV  #SV x 3 array -- The vertices of the output level surface mesh
+  //  SF  #SF x 3 array -- The face indexes of the output level surface mesh
+  //  J   #SF list of indices into TT revealing which tet each face comes from
+  template <typename DerivedTV,
+            typename DerivedTT,
+            typename DerivedS,
+            typename DerivedSV,
+            typename DerivedSF,
+            typename DerivedJ>
+  IGL_INLINE void marching_tets(
+      const Eigen::PlainObjectBase<DerivedTV>& TV,
+      const Eigen::PlainObjectBase<DerivedTT>& TT,
+      const Eigen::PlainObjectBase<DerivedS>& S,
+      double isovalue,
+      Eigen::PlainObjectBase<DerivedSV>& SV,
+      Eigen::PlainObjectBase<DerivedSF>& SF,
+      Eigen::PlainObjectBase<DerivedJ>& J) {
+    Eigen::SparseMatrix<double> _BC;
+    return igl::marching_tets(TV, TT, S, isovalue, SV, SF, J, _BC);
+  }
+
+  // marching_tets( TV, TT, S, isovalue, SV, SF, BC)
+  //
+  // performs the marching tetrahedra algorithm on a tet mesh defined by TV and
+  // TT with scalar values defined at each vertex in TV. The output is a
+  // triangle mesh approximating the isosurface coresponding to the value
+  // isovalue.
+  //
+  // Input:
+  //  TV  #tet_vertices x 3 array -- The vertices of the tetrahedral mesh
+  //  TT  #tets x 4 array --  The indexes of each tet in the tetrahedral mesh
+  //  S  #tet_vertices x 1 array -- The values defined on each tet vertex
+  //  isovalue  scalar -- The isovalue of the level set we want to compute
+  //
+  // Output:
+  //  SV  #SV x 3 array -- The vertices of the output level surface mesh
+  //  SF  #SF x 3 array -- The face indexes of the output level surface mesh
+  //  BC  #SV x #TV list of barycentric coordinates so that SV = BC*TV
+  template <typename DerivedTV,
+            typename DerivedTT,
+            typename DerivedS,
+            typename DerivedSV,
+            typename DerivedSF,
+            typename BCType>
+  IGL_INLINE void marching_tets(
+      const Eigen::PlainObjectBase<DerivedTV>& TV,
+      const Eigen::PlainObjectBase<DerivedTT>& TT,
+      const Eigen::PlainObjectBase<DerivedS>& S,
+      double isovalue,
+      Eigen::PlainObjectBase<DerivedSV>& SV,
+      Eigen::PlainObjectBase<DerivedSF>& SF,
+      Eigen::SparseMatrix<BCType>& BC) {
+    Eigen::VectorXi _J;
+    return igl::marching_tets(TV, TT, S, isovalue, SV, SF, _J, BC);
+  }
+
+  // marching_tets( TV, TT, S, isovalue, SV, SF)
+  //
+  // performs the marching tetrahedra algorithm on a tet mesh defined by TV and
+  // TT with scalar values defined at each vertex in TV. The output is a
+  // triangle mesh approximating the isosurface coresponding to the value
+  // isovalue.
+  //
+  // Input:
+  //  TV  #tet_vertices x 3 array -- The vertices of the tetrahedral mesh
+  //  TT  #tets x 4 array --  The indexes of each tet in the tetrahedral mesh
+  //  S  #tet_vertices x 1 array -- The values defined on each tet vertex
+  //  isovalue  scalar -- The isovalue of the level set we want to compute
+  //
+  // Output:
+  //  SV  #SV x 3 array -- The vertices of the output level surface mesh
+  //  SF  #SF x 3 array -- The face indexes of the output level surface mesh
+  template <typename DerivedTV,
+            typename DerivedTT,
+            typename DerivedS,
+            typename DerivedSV,
+            typename DerivedSF>
+  IGL_INLINE void marching_tets(
+      const Eigen::PlainObjectBase<DerivedTV>& TV,
+      const Eigen::PlainObjectBase<DerivedTT>& TT,
+      const Eigen::PlainObjectBase<DerivedS>& S,
+      double isovalue,
+      Eigen::PlainObjectBase<DerivedSV>& SV,
+      Eigen::PlainObjectBase<DerivedSF>& SF) {
+    Eigen::VectorXi _J;
+    Eigen::SparseMatrix<double> _BC;
+    return igl::marching_tets(TV, TT, S, isovalue, SV, SF, _J, _BC);
+  }
+
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "marching_tets.cpp"
+#endif
+
+#endif // IGL_MARCHING_TETS_H
diff --git a/include/igl/slice_tets.cpp b/include/igl/slice_tets.cpp
deleted file mode 100644
index dbe042fa1..000000000
--- a/include/igl/slice_tets.cpp
+++ /dev/null
@@ -1,356 +0,0 @@
-// This file is part of libigl, a simple c++ geometry processing library.
-// 
-// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
-// obtain one at http://mozilla.org/MPL/2.0/.
-#include "slice_tets.h"
-#include "LinSpaced.h"
-#include "sort.h"
-#include "edges.h"
-#include "slice.h"
-#include "cat.h"
-#include "ismember.h"
-#include "unique_rows.h"
-#include <cassert>
-#include <algorithm>
-#include <vector>
-
-template <
-  typename DerivedV, 
-  typename DerivedT, 
-  typename DerivedS,
-  typename DerivedSV,
-  typename DerivedSF,
-  typename DerivedJ,
-  typename BCType>
-IGL_INLINE void igl::slice_tets(
-  const Eigen::MatrixBase<DerivedV>& V,
-  const Eigen::MatrixBase<DerivedT>& T,
-  const Eigen::MatrixBase<DerivedS> & S,
-  Eigen::PlainObjectBase<DerivedSV>& SV,
-  Eigen::PlainObjectBase<DerivedSF>& SF,
-  Eigen::PlainObjectBase<DerivedJ>& J,
-  Eigen::SparseMatrix<BCType> & BC)
-{
-  Eigen::MatrixXi sE;
-  Eigen::Matrix<typename DerivedSV::Scalar,Eigen::Dynamic,1> lambda;
-  igl::slice_tets(V,T,S,SV,SF,J,sE,lambda);
-  const int ns = SV.rows();
-  std::vector<Eigen::Triplet<BCType> > BCIJV(ns*2);
-  for(int i = 0;i<ns;i++)
-  {
-    BCIJV[2*i+0] = Eigen::Triplet<BCType>(i,sE(i,0),    lambda(i));
-    BCIJV[2*i+1] = Eigen::Triplet<BCType>(i,sE(i,1),1.0-lambda(i));
-  }
-  BC.resize(SV.rows(),V.rows());
-  BC.setFromTriplets(BCIJV.begin(),BCIJV.end());
-}
-
-template <
-  typename DerivedV, 
-  typename DerivedT, 
-  typename DerivedS,
-  typename DerivedSV,
-  typename DerivedSF,
-  typename DerivedJ>
-IGL_INLINE void igl::slice_tets(
-  const Eigen::MatrixBase<DerivedV>& V,
-  const Eigen::MatrixBase<DerivedT>& T,
-  const Eigen::MatrixBase<DerivedS> & S,
-  Eigen::PlainObjectBase<DerivedSV>& SV,
-  Eigen::PlainObjectBase<DerivedSF>& SF,
-  Eigen::PlainObjectBase<DerivedJ>& J)
-{
-  Eigen::MatrixXi sE;
-  Eigen::Matrix<typename DerivedSV::Scalar,Eigen::Dynamic,1> lambda;
-  igl::slice_tets(V,T,S,SV,SF,J,sE,lambda);
-}
-
-template <
-  typename DerivedV, 
-  typename DerivedT, 
-  typename DerivedS,
-  typename DerivedSV,
-  typename DerivedSF,
-  typename DerivedJ,
-  typename DerivedsE,
-  typename Derivedlambda>
-IGL_INLINE void igl::slice_tets(
-  const Eigen::MatrixBase<DerivedV>& V,
-  const Eigen::MatrixBase<DerivedT>& T,
-  const Eigen::MatrixBase<DerivedS> & S,
-  Eigen::PlainObjectBase<DerivedSV>& SV,
-  Eigen::PlainObjectBase<DerivedSF>& SF,
-  Eigen::PlainObjectBase<DerivedJ>& J,
-  Eigen::PlainObjectBase<DerivedsE>& sE,
-  Eigen::PlainObjectBase<Derivedlambda>& lambda)
-{
-
-  using namespace Eigen;
-  using namespace std;
-  assert(V.cols() == 3 && "V should be #V by 3");
-  assert(T.cols() == 4 && "T should be #T by 4");
-
-  static const Eigen::Matrix<int,12,4> flipped_order = 
-    (Eigen::Matrix<int,12,4>(12,4)<<
-      3,2,0,1,
-      3,1,2,0,
-      3,0,1,2,
-      2,3,1,0,
-      2,1,0,3,
-      2,0,3,1,
-      1,3,0,2,
-      1,2,3,0,
-      1,0,2,3,
-      0,3,2,1,
-      0,2,1,3,
-      0,1,3,2
-    ).finished();
-
-  // number of tets
-  const size_t m = T.rows();
-
-  typedef typename DerivedS::Scalar Scalar;
-  typedef typename DerivedT::Scalar Index;
-  typedef Matrix<Scalar,Dynamic,1> VectorXS;
-  typedef Matrix<Scalar,Dynamic,4> MatrixX4S;
-  typedef Matrix<Scalar,Dynamic,3> MatrixX3S;
-  typedef Matrix<Scalar,Dynamic,2> MatrixX2S;
-  typedef Matrix<Index,Dynamic,4> MatrixX4I;
-  typedef Matrix<Index,Dynamic,3> MatrixX3I;
-  typedef Matrix<Index,Dynamic,2> MatrixX2I;
-  typedef Matrix<Index,Dynamic,1> VectorXI;
-  typedef Array<bool,Dynamic,1> ArrayXb;
-  
-  MatrixX4S IT(m,4);
-  for(size_t t = 0;t<m;t++)
-  {
-    for(size_t c = 0;c<4;c++)
-    {
-      IT(t,c) = S(T(t,c));
-    }
-  }
-
-  // Essentially, just a glorified slice(X,1)
-  // 
-  // Inputs:
-  //   T  #T by 4 list of tet indices into V
-  //   IT  #IT by 4 list of isosurface values at each tet
-  //   I  #I list of bools whether to grab data corresponding to each tet
-  const auto & extract_rows = [](
-    const MatrixBase<DerivedT> & T,
-    const MatrixX4S & IT,
-    const ArrayXb & I,
-    MatrixX4I  & TI,
-    MatrixX4S & ITI,
-    VectorXI & JI)
-  {
-    const Index num_I = std::count(I.data(),I.data()+I.size(),true);
-    TI.resize(num_I,4);
-    ITI.resize(num_I,4);
-    JI.resize(num_I,1);
-    {
-      size_t k = 0;
-      for(size_t t = 0;t<(size_t)T.rows();t++)
-      {
-        if(I(t))
-        {
-          TI.row(k) = T.row(t);
-          ITI.row(k) = IT.row(t);
-          JI(k) = t;
-          k++;
-        }
-      }
-      assert(k == num_I);
-    }
-  };
-
-  ArrayXb I13 = (IT.array()<0).rowwise().count()==1;
-  ArrayXb I31 = (IT.array()>0).rowwise().count()==1;
-  ArrayXb I22 = (IT.array()<0).rowwise().count()==2;
-  MatrixX4I T13,T31,T22;
-  MatrixX4S IT13,IT31,IT22;
-  VectorXI J13,J31,J22;
-  extract_rows(T,IT,I13,T13,IT13,J13);
-  extract_rows(T,IT,I31,T31,IT31,J31);
-  extract_rows(T,IT,I22,T22,IT22,J22);
-
-  const auto & apply_sort4 = [] (
-     const MatrixX4I & T, 
-     const MatrixX4I & sJ, 
-     MatrixX4I & sT)
-  {
-    sT.resize(T.rows(),4);
-    for(size_t t = 0;t<(size_t)T.rows();t++)
-    {
-      for(size_t c = 0;c<4;c++)
-      {
-        sT(t,c) = T(t,sJ(t,c));
-      }
-    }
-  };
-
-  const auto & apply_sort2 = [] (
-     const MatrixX2I & E, 
-     const MatrixX2I & sJ, 
-     Eigen::PlainObjectBase<DerivedsE>& sE)
-  {
-    sE.resize(E.rows(),2);
-    for(size_t t = 0;t<(size_t)E.rows();t++)
-    {
-      for(size_t c = 0;c<2;c++)
-      {
-        sE(t,c) = E(t,sJ(t,c));
-      }
-    }
-  };
-
-  const auto & one_below = [&apply_sort4](
-    const MatrixX4I & T,
-    const MatrixX4S & IT,
-    MatrixX2I & U,
-    MatrixX3I & SF)
-  {
-    // Number of tets
-    const size_t m = T.rows();
-    if(m == 0)
-    {
-      U.resize(0,2);
-      SF.resize(0,3);
-      return;
-    }
-    MatrixX4S sIT;
-    MatrixX4I sJ;
-    sort(IT,2,true,sIT,sJ);
-    MatrixX4I sT;
-    apply_sort4(T,sJ,sT);
-    U.resize(3*m,2);
-    U<<
-      sT.col(0),sT.col(1),
-      sT.col(0),sT.col(2),
-      sT.col(0),sT.col(3);
-    SF.resize(m,3);
-    for(size_t c = 0;c<3;c++)
-    {
-      SF.col(c) = 
-        igl::LinSpaced<
-        Eigen::Matrix<typename DerivedSF::Scalar,Eigen::Dynamic,1> >
-        (m,0+c*m,(m-1)+c*m);
-    }
-    ArrayXb flip;
-    {
-      VectorXi _;
-      ismember_rows(sJ,flipped_order,flip,_);
-    }
-    for(int i = 0;i<m;i++)
-    {
-      if(flip(i))
-      {
-        SF.row(i) = SF.row(i).reverse().eval();
-      }
-    }
-  };
-
-  const auto & two_below = [&apply_sort4](
-    const MatrixX4I & T,
-    const MatrixX4S & IT,
-    MatrixX2I & U,
-    MatrixX3I & SF)
-  {
-    // Number of tets
-    const size_t m = T.rows();
-    if(m == 0)
-    {
-      U.resize(0,2);
-      SF.resize(0,3);
-      return;
-    }
-    MatrixX4S sIT;
-    MatrixX4I sJ;
-    sort(IT,2,true,sIT,sJ);
-    MatrixX4I sT;
-    apply_sort4(T,sJ,sT);
-    U.resize(4*m,2);
-    U<<
-      sT.col(0),sT.col(2),
-      sT.col(0),sT.col(3),
-      sT.col(1),sT.col(2),
-      sT.col(1),sT.col(3);
-    SF.resize(2*m,3);
-    SF.block(0,0,m,1) = igl::LinSpaced<VectorXI >(m,0+0*m,(m-1)+0*m);
-    SF.block(0,1,m,1) = igl::LinSpaced<VectorXI >(m,0+1*m,(m-1)+1*m);
-    SF.block(0,2,m,1) = igl::LinSpaced<VectorXI >(m,0+3*m,(m-1)+3*m);
-    SF.block(m,0,m,1) = igl::LinSpaced<VectorXI >(m,0+0*m,(m-1)+0*m);
-    SF.block(m,1,m,1) = igl::LinSpaced<VectorXI >(m,0+3*m,(m-1)+3*m);
-    SF.block(m,2,m,1) = igl::LinSpaced<VectorXI >(m,0+2*m,(m-1)+2*m);
-    ArrayXb flip;
-    {
-      VectorXi _;
-      ismember_rows(sJ,flipped_order,flip,_);
-    }
-    for(int i = 0;i<m;i++)
-    {
-      if(flip(i))
-      {
-        SF.row(i  ) = SF.row(i  ).reverse().eval();
-        SF.row(i+m) = SF.row(i+m).reverse().eval();
-      }
-    }
-  };
-
-  MatrixX3I SF13,SF31,SF22;
-  MatrixX2I U13,U31,U22;
-  one_below(T13, IT13,U13,SF13);
-  one_below(T31,-IT31,U31,SF31);
-  two_below(T22, IT22,U22,SF22);
-  // https://forum.kde.org/viewtopic.php?f=74&t=107974
-  const MatrixX2I U = 
-    (MatrixX2I(U13.rows()+ U31.rows()+ U22.rows(),2)<<U13,U31,U22).finished();
-  MatrixX2I sU;
-  {
-    MatrixX2I _;
-    sort(U,2,true,sU,_);
-  }
-  MatrixX2I E;
-  VectorXI uI,uJ;
-  unique_rows(sU,E,uI,uJ);
-  MatrixX2S IE(E.rows(),2);
-  for(size_t t = 0;t<E.rows();t++)
-  {
-    for(size_t c = 0;c<2;c++)
-    {
-      IE(t,c) = S(E(t,c));
-    }
-  }
-  MatrixX2S sIE;
-  MatrixX2I sJ;
-  sort(IE,2,true,sIE,sJ);
-  apply_sort2(E,sJ,sE);
-  lambda = sIE.col(1).array() / (sIE.col(1)-sIE.col(0)).array();
-  SV.resize(sE.rows(),V.cols());
-  for(int e = 0;e<sE.rows();e++)
-  {
-    SV.row(e) = V.row(sE(e,0)).template cast<Scalar>()*lambda(e) + 
-                V.row(sE(e,1)).template cast<Scalar>()*(1.0-lambda(e));
-  }
-  SF.resize( SF13.rows()+SF31.rows()+SF22.rows(),3);
-  SF<<
-    SF13,
-    U13.rows()+           SF31.rowwise().reverse().array(),
-    U13.rows()+U31.rows()+SF22.array();
-
-  std::for_each(
-    SF.data(),
-    SF.data()+SF.size(),
-    [&uJ](typename DerivedSF::Scalar & i){i=uJ(i);});
-
-  J.resize(SF.rows());
-  J<<J13,J31,J22,J22;
-}
-
-#ifdef IGL_STATIC_LIBRARY
-// Explicit template instantiation
-template void igl::slice_tets<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<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::SparseMatrix<double, 0, int>&);
-#endif
diff --git a/include/igl/slice_tets.h b/include/igl/slice_tets.h
deleted file mode 100644
index 349c5d790..000000000
--- a/include/igl/slice_tets.h
+++ /dev/null
@@ -1,96 +0,0 @@
-// This file is part of libigl, a simple c++ geometry processing library.
-// 
-// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
-// obtain one at http://mozilla.org/MPL/2.0/.
-#ifndef IGL_SLICE_TETS_H
-#define IGL_SLICE_TETS_H
-#include "igl_inline.h"
-
-#include <Eigen/Dense>
-#include <Eigen/Sparse>
-
-#include <vector>
-
-namespace igl
-{
-  // SLICE_TETS Slice through a tet mesh (V,T) along a given plane (via its
-  // implicit equation).
-  //
-  // Inputs:
-  //   V  #V by 3 list of tet mesh vertices
-  //   T  #T by 4 list of tet indices into V 
-  ////   plane  list of 4 coefficients in the plane equation: [x y z 1]'*plane = 0
-  //   S  #V list of values so that S = 0 is the desired isosurface
-  // Outputs:
-  //   SV  #SV by 3 list of triangle mesh vertices along slice
-  //   SF  #SF by 3 list of triangles indices into SV
-  //   J  #SF list of indices into T revealing from which tet each faces comes
-  //   BC  #SU by #V list of barycentric coordinates (or more generally: linear
-  //     interpolation coordinates) so that SV = BC*V
-  // 
-  template <
-    typename DerivedV, 
-    typename DerivedT, 
-    typename DerivedS,
-    typename DerivedSV,
-    typename DerivedSF,
-    typename DerivedJ,
-    typename BCType>
-  IGL_INLINE void slice_tets(
-    const Eigen::MatrixBase<DerivedV>& V,
-    const Eigen::MatrixBase<DerivedT>& T,
-    const Eigen::MatrixBase<DerivedS> & S,
-    Eigen::PlainObjectBase<DerivedSV>& SV,
-    Eigen::PlainObjectBase<DerivedSF>& SF,
-    Eigen::PlainObjectBase<DerivedJ>& J,
-    Eigen::SparseMatrix<BCType> & BC);
-  template <
-    typename DerivedV, 
-    typename DerivedT, 
-    typename DerivedS,
-    typename DerivedSV,
-    typename DerivedSF,
-    typename DerivedJ>
-  IGL_INLINE void slice_tets(
-    const Eigen::MatrixBase<DerivedV>& V,
-    const Eigen::MatrixBase<DerivedT>& T,
-    const Eigen::MatrixBase<DerivedS> & S,
-    Eigen::PlainObjectBase<DerivedSV>& SV,
-    Eigen::PlainObjectBase<DerivedSF>& SF,
-    Eigen::PlainObjectBase<DerivedJ>& J);
-  // Outputs:
-  //   sE  #SV by 2 list of sorted edge indices into V
-  //   lambda  #SV by 1 list of parameters along each edge in sE so that:
-  //     SV(i,:) = V(sE(i,1),:)*lambda(i) + V(sE(i,2),:)*(1-lambda(i));
-  template <
-    typename DerivedV, 
-    typename DerivedT, 
-    typename DerivedS,
-    typename DerivedSV,
-    typename DerivedSF,
-    typename DerivedJ,
-    typename DerivedsE,
-    typename Derivedlambda
-    >
-  IGL_INLINE void slice_tets(
-    const Eigen::MatrixBase<DerivedV>& V,
-    const Eigen::MatrixBase<DerivedT>& T,
-    const Eigen::MatrixBase<DerivedS> & S,
-    Eigen::PlainObjectBase<DerivedSV>& SV,
-    Eigen::PlainObjectBase<DerivedSF>& SF,
-    Eigen::PlainObjectBase<DerivedJ>& J,
-    Eigen::PlainObjectBase<DerivedsE>& sE,
-    Eigen::PlainObjectBase<Derivedlambda>& lambda);
-
-}
-
-#ifndef IGL_STATIC_LIBRARY
-#  include "slice_tets.cpp"
-#endif
-
-#endif
-
-
diff --git a/include/igl/writeDMAT.cpp b/include/igl/writeDMAT.cpp
index c8f6ae0f6..656dae693 100644
--- a/include/igl/writeDMAT.cpp
+++ b/include/igl/writeDMAT.cpp
@@ -89,4 +89,5 @@ template bool igl::writeDMAT<Eigen::Matrix<double, 1, 3, 1, 1, 3> >(std::string,
 template bool igl::writeDMAT<Eigen::Matrix<float, 1, 3, 1, 1, 3> >(std::string, Eigen::MatrixBase<Eigen::Matrix<float, 1, 3, 1, 1, 3> > 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::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, bool);
 template bool igl::writeDMAT<Eigen::Matrix<int, -1, 2, 0, -1, 2> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::MatrixBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, bool);
+template bool igl::writeDMAT<Eigen::Matrix<double, -1, 1, 0, -1, 1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, bool);
 #endif
diff --git a/python/py_doc.cpp b/python/py_doc.cpp
index c16997f82..ef9dffea2 100644
--- a/python/py_doc.cpp
+++ b/python/py_doc.cpp
@@ -1275,7 +1275,7 @@ const char *__doc_igl_slice_mask = R"igl_Qu8mg5v7(// Act like the matlab X(row_m
   //   Y  #trues-in-R by #trues-in-C matrix
   //
   // See also: slice_mask)igl_Qu8mg5v7";
-const char *__doc_igl_slice_tets = R"igl_Qu8mg5v7(// SLICE_TETS Slice through a tet mesh (V,T) along a given plane (via its
+const char *__doc_igl_marching_tets = R"igl_Qu8mg5v7(// SLICE_TETS Slice through a tet mesh (V,T) along a given plane (via its
   // implicit equation).
   //
   // Inputs:
diff --git a/python/py_doc.h b/python/py_doc.h
index 4446fa90f..c5165e216 100644
--- a/python/py_doc.h
+++ b/python/py_doc.h
@@ -108,7 +108,7 @@ extern const char *__doc_igl_signed_distance_winding_number;
 extern const char *__doc_igl_slice;
 extern const char *__doc_igl_slice_into;
 extern const char *__doc_igl_slice_mask;
-extern const char *__doc_igl_slice_tets;
+extern const char *__doc_igl_marching_tets;
 extern const char *__doc_igl_sortrows;
 extern const char *__doc_igl_streamlines_init;
 extern const char *__doc_igl_streamlines_next;
@@ -127,4 +127,4 @@ extern const char *__doc_igl_winding_number_2;
 extern const char *__doc_igl_writeMESH;
 extern const char *__doc_igl_writeOBJ;
 extern const char *__doc_igl_writePLY;
-extern const char *__doc_igl_readPLY;
\ No newline at end of file
+extern const char *__doc_igl_readPLY;
diff --git a/python/py_igl.cpp b/python/py_igl.cpp
index 2f6f71bcc..551babcbc 100644
--- a/python/py_igl.cpp
+++ b/python/py_igl.cpp
@@ -92,7 +92,7 @@
 #include <igl/slice.h>
 #include <igl/slice_into.h>
 #include <igl/slice_mask.h>
-#include <igl/slice_tets.h>
+#include <igl/marching_tets.h>
 #include <igl/sortrows.h>
 #include <igl/triangle_triangle_adjacency.h>
 #include <igl/unique.h>
@@ -191,7 +191,7 @@ void python_export_igl(py::module &m)
 #include "py_igl/py_slice.cpp"
 #include "py_igl/py_slice_into.cpp"
 #include "py_igl/py_slice_mask.cpp"
-#include "py_igl/py_slice_tets.cpp"
+#include "py_igl/py_marching_tets.cpp"
 #include "py_igl/py_sortrows.cpp"
 #include "py_igl/py_triangle_triangle_adjacency.cpp"
 #include "py_igl/py_unique.cpp"
diff --git a/python/py_igl/py_slice_tets.cpp b/python/py_igl/py_marching_tets.cpp
similarity index 88%
rename from python/py_igl/py_slice_tets.cpp
rename to python/py_igl/py_marching_tets.cpp
index 7cdeeab6b..ec02e1dee 100644
--- a/python/py_igl/py_slice_tets.cpp
+++ b/python/py_igl/py_marching_tets.cpp
@@ -5,7 +5,7 @@
 // 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/.
-m.def("slice_tets", []
+m.def("marching_tets", []
 (
   const Eigen::MatrixXd& V,
   const Eigen::MatrixXi& T,
@@ -21,8 +21,8 @@ m.def("slice_tets", []
   if (plane.size() != 0)
     planev = plane;
   Eigen::VectorXi Jv;
-  igl::slice_tets(V, T, planev, U, G, Jv, BC);
+  igl::marching_tets(V, T, planev, U, G, Jv, BC);
   J = Jv;
-}, __doc_igl_slice_tets,
+}, __doc_igl_marching_tets,
 py::arg("V"), py::arg("T"), py::arg("plane"), py::arg("U"), py::arg("G"), py::arg("J"), py::arg("BC"));
 
diff --git a/python/python_shared.cpp b/python/python_shared.cpp
index 75ef5bfac..894f0092d 100644
--- a/python/python_shared.cpp
+++ b/python/python_shared.cpp
@@ -147,7 +147,7 @@ PYBIND11_PLUGIN(pyigl) {
            slice
            slice_into
            slice_mask
-           slice_tets
+           marching_tets
            sortrows
            streamlines
            triangle_triangle_adjacency
diff --git a/python/tutorial/702_WindingNumber.py b/python/tutorial/702_WindingNumber.py
index c76f70d54..4387afdf8 100755
--- a/python/tutorial/702_WindingNumber.py
+++ b/python/tutorial/702_WindingNumber.py
@@ -35,7 +35,7 @@ def update(viewer):
     F_vis = igl.eigen.MatrixXi()
     J = igl.eigen.MatrixXi()
     bary = igl.eigen.SparseMatrixd()
-    igl.slice_tets(V, T, plane, V_vis, F_vis, J, bary)
+    igl.marching_tets(V, T, plane, V_vis, F_vis, J, bary)
     W_vis = igl.eigen.MatrixXd()
     igl.slice(W, J, W_vis)
     C_vis = igl.eigen.MatrixXd()
diff --git a/python/tutorial/704_SignedDistance.py b/python/tutorial/704_SignedDistance.py
index 22a83b00d..cbc9d829d 100644
--- a/python/tutorial/704_SignedDistance.py
+++ b/python/tutorial/704_SignedDistance.py
@@ -55,7 +55,7 @@ def update_visualization(viewer):
     # Extract triangle mesh slice through volume mesh and subdivide nasty triangles
     J = igl.eigen.MatrixXi()
     bary = igl.eigen.SparseMatrixd()
-    igl.slice_tets(V, T, plane, V_vis, F_vis, J, bary)
+    igl.marching_tets(V, T, plane, V_vis, F_vis, J, bary)
     max_l = 0.03
     while True:
         l = igl.eigen.MatrixXd()
diff --git a/tutorial/702_WindingNumber/main.cpp b/tutorial/702_WindingNumber/main.cpp
index f9ff70e36..d65bdc0c8 100755
--- a/tutorial/702_WindingNumber/main.cpp
+++ b/tutorial/702_WindingNumber/main.cpp
@@ -3,7 +3,7 @@
 #include <igl/parula.h>
 #include <igl/readMESH.h>
 #include <igl/slice.h>
-#include <igl/slice_tets.h>
+#include <igl/marching_tets.h>
 #include <igl/winding_number.h>
 #include <igl/opengl/glfw/Viewer.h>
 #include <Eigen/Sparse>
@@ -40,7 +40,7 @@ void update_visualization(igl::opengl::glfw::Viewer & viewer)
         V.col(1)*plane(1) + 
         V.col(2)*plane(2)).array()
       + plane(3);
-    igl::slice_tets(V,T,IV,V_vis,F_vis,J,bary);
+    igl::marching_tets(V,T,IV,V_vis,F_vis,J,bary);
   }
   VectorXd W_vis;
   igl::slice(W,J,W_vis);
diff --git a/tutorial/704_SignedDistance/main.cpp b/tutorial/704_SignedDistance/main.cpp
index 37e00b3dc..789f4cdb7 100755
--- a/tutorial/704_SignedDistance/main.cpp
+++ b/tutorial/704_SignedDistance/main.cpp
@@ -8,7 +8,7 @@
 #include <igl/readMESH.h>
 #include <igl/signed_distance.h>
 #include <igl/slice_mask.h>
-#include <igl/slice_tets.h>
+#include <igl/marching_tets.h>
 #include <igl/upsample.h>
 #include <igl/opengl/glfw/Viewer.h>
 #include <igl/writeOBJ.h>
@@ -48,7 +48,7 @@ void update_visualization(igl::opengl::glfw::Viewer & viewer)
          V.col(1)*plane(1) + 
          V.col(2)*plane(2)).array()
         + plane(3);
-      igl::slice_tets(V,T,IV,V_vis,F_vis,J,bary);
+      igl::marching_tets(V,T,IV,V_vis,F_vis,J,bary);
       igl::writeOBJ("vis.obj",V_vis,F_vis);
     }
     while(true)
diff --git a/tutorial/714_MarchingTets/CMakeLists.txt b/tutorial/714_MarchingTets/CMakeLists.txt
new file mode 100644
index 000000000..391af8e91
--- /dev/null
+++ b/tutorial/714_MarchingTets/CMakeLists.txt
@@ -0,0 +1,5 @@
+get_filename_component(PROJECT_NAME ${CMAKE_CURRENT_SOURCE_DIR} NAME)
+project(${PROJECT_NAME})
+
+add_executable(${PROJECT_NAME}_bin main.cpp)
+target_link_libraries(${PROJECT_NAME}_bin igl::core igl::opengl igl::opengl_glfw igl::tetgen tutorials)
diff --git a/tutorial/714_MarchingTets/main.cpp b/tutorial/714_MarchingTets/main.cpp
new file mode 100644
index 000000000..bcbd9176a
--- /dev/null
+++ b/tutorial/714_MarchingTets/main.cpp
@@ -0,0 +1,48 @@
+#include <igl/opengl/glfw/Viewer.h>
+#include <igl/copyleft/tetgen/tetrahedralize.h>
+#include <igl/readOBJ.h>
+#include <igl/marching_tets.h>
+#include <Eigen/Core>
+
+#include "tutorial_shared_path.h"
+
+
+int main(int argc, char * argv[])
+{
+
+  // Load a surface mesh which is a cube
+  Eigen::MatrixXd surfaceV;
+  Eigen::MatrixXi surfaceF;
+  igl::readOBJ(TUTORIAL_SHARED_PATH "/cube.obj", surfaceV, surfaceF);
+
+  // Find the centroid of the loaded mesh
+  Eigen::RowVector3d surfaceCenter = surfaceV.colwise().sum() / surfaceV.rows();
+
+  // Center the mesh about the origin
+  surfaceV.rowwise() -= surfaceCenter;
+
+  // Tetrahedralize the surface mesh
+  Eigen::MatrixXd TV; // Tet mesh vertices
+  Eigen::MatrixXi TF; // Tet mesh boundary face indices
+  Eigen::MatrixXi TT; // Tet mesh tetrahedron indices
+  igl::copyleft::tetgen::tetrahedralize(surfaceV, surfaceF, "pq1.414a0.0001", TV, TT, TF);
+
+  // Compute a scalar at each tet vertex which is the distance from the vertex to the origin
+  Eigen::VectorXd S = TV.rowwise().norm();
+
+  // Compute a mesh (stored in SV, SF) representing the iso-level-set for the isovalue 0.5
+  Eigen::MatrixXd SV;
+  Eigen::MatrixXi SF;
+  igl::marching_tets(TV, TT, S, 0.45, SV, SF);
+
+  // Draw the mesh stored in (SV, SF)
+  igl::opengl::glfw::Viewer viewer;
+  viewer.data().set_mesh(SV, SF);
+  viewer.callback_key_down =
+    [&](igl::opengl::glfw::Viewer & viewer, unsigned char key, int mod)->bool
+    {
+      viewer.data().set_face_based(true);
+      return true;
+    };
+  viewer.launch();
+}
diff --git a/tutorial/CMakeLists.txt b/tutorial/CMakeLists.txt
index 947fb86cc..b7b256c25 100644
--- a/tutorial/CMakeLists.txt
+++ b/tutorial/CMakeLists.txt
@@ -149,6 +149,9 @@ if(TUTORIALS_CHAPTER7)
   add_subdirectory("711_Subdivision")
   add_subdirectory("712_DataSmoothing")
   add_subdirectory("713_ShapeUp")
+  if(LIBIGL_WITH_TETGEN)
+    add_subdirectory("714_MarchingTets")
+  endif()
 endif()