diff --git a/tutorial/712_DataSmoothing/CMakeLists.txt b/tutorial/712_DataSmoothing/CMakeLists.txt
new file mode 100644
index 000000000..cd71d8881
--- /dev/null
+++ b/tutorial/712_DataSmoothing/CMakeLists.txt
@@ -0,0 +1,8 @@
+cmake_minimum_required(VERSION 2.8.12)
+project(712_DataSmoothing)
+
+add_executable(${PROJECT_NAME}_bin
+  main.cpp)
+target_include_directories(${PROJECT_NAME}_bin PRIVATE ${LIBIGL_INCLUDE_DIRS})
+target_compile_definitions(${PROJECT_NAME}_bin PRIVATE ${LIBIGL_DEFINITIONS})
+target_link_libraries(${PROJECT_NAME}_bin ${LIBIGL_LIBRARIES} ${LIBIGL_VIEWER_EXTRA_LIBRARIES} ${LIBIGL_OPENGL_EXTRA_LIBRARIES} ${LIBIGL_OPENGL_GLFW_EXTRA_LIBRARIES})
diff --git a/tutorial/712_DataSmoothing/isolines.h b/tutorial/712_DataSmoothing/isolines.h
new file mode 100644
index 000000000..e7c807749
--- /dev/null
+++ b/tutorial/712_DataSmoothing/isolines.h
@@ -0,0 +1,90 @@
+
+#include <vector>
+#include <limits>
+#include <stdlib.h>
+
+#include <igl/remove_unreferenced.h>
+
+
+static void isolines(const Eigen::MatrixXd& V,
+    const Eigen::MatrixXi& F,
+    const Eigen::VectorXd& z,
+    const int grads,
+    Eigen::MatrixXd& isoV,
+    Eigen::MatrixXi& isoE)
+{
+    const double min = z.minCoeff(), max = z.maxCoeff();
+
+    //Following http://www.alecjacobson.com/weblog/?p=2529
+    Eigen::VectorXd iso(grads+1);
+    for(int i=0; i<iso.size(); ++i)
+        iso(i) = double(i)/double(grads)*(max-min) + min;
+
+    Eigen::MatrixXd t12(F.rows(), iso.size()), t23(F.rows(), iso.size()),
+        t31(F.rows(), iso.size());
+    for(int i=0; i<t12.rows(); ++i) {
+        const double z1=z(F(i,0)), z2=z(F(i,1)), z3=z(F(i,2));
+        const double s12 = z2-z1;
+        const double s23 = z3-z2;
+        const double s31 = z1-z3;
+        for(int j=0; j<t12.cols(); ++j) {
+            t12(i,j) = (iso(j)-z1) / s12;
+            t23(i,j) = (iso(j)-z2) / s23;
+            t31(i,j) = (iso(j)-z3) / s31;
+            if(t12(i,j)<0 || t12(i,j)>1)
+                t12(i,j) = std::numeric_limits<double>::quiet_NaN();
+            if(t23(i,j)<0 || t23(i,j)>1)
+                t23(i,j) = std::numeric_limits<double>::quiet_NaN();
+            if(t31(i,j)<0 || t31(i,j)>1)
+                t31(i,j) = std::numeric_limits<double>::quiet_NaN();
+        }
+    }
+
+    std::vector<int> F12, F23, F31, I12, I23, I31;
+    for(int i=0; i<t12.rows(); ++i) {
+        for(int j=0; j<t12.cols(); ++j) {
+            if(std::isfinite(t23(i,j)) && std::isfinite(t31(i,j))) {
+                F12.push_back(i);
+                I12.push_back(j);
+            }
+            if(std::isfinite(t31(i,j)) && std::isfinite(t12(i,j))) {
+                F23.push_back(i);
+                I23.push_back(j);
+            }
+            if(std::isfinite(t12(i,j)) && std::isfinite(t23(i,j))) {
+                F31.push_back(i);
+                I31.push_back(j);
+            }
+        }
+    }
+
+    const int K = F12.size()+F23.size()+F31.size();
+    isoV.resize(2*K, 3);
+    int b = 0;
+    for(int i=0; i<F12.size(); ++i) {
+        isoV.row(b+i) = (1.-t23(F12[i],I12[i]))*V.row(F(F12[i],1))
+            + t23(F12[i],I12[i])*V.row(F(F12[i],2));
+        isoV.row(K+b+i) = (1.-t31(F12[i],I12[i]))*V.row(F(F12[i],2))
+            + t31(F12[i],I12[i])*V.row(F(F12[i],0));
+    }
+    b += F12.size();
+    for(int i=0; i<F23.size(); ++i) {
+        isoV.row(b+i) = (1.-t31(F23[i],I23[i]))*V.row(F(F23[i],2))
+            + t31(F23[i],I23[i])*V.row(F(F23[i],0));
+        isoV.row(K+b+i) = (1.-t12(F23[i],I23[i]))*V.row(F(F23[i],0))
+            + t12(F23[i],I23[i])*V.row(F(F23[i],1));
+    }
+    b += F23.size();
+    for(int i=0; i<F31.size(); ++i) {
+        isoV.row(b+i) = (1.-t12(F31[i],I31[i]))*V.row(F(F31[i],0))
+            + t12(F31[i],I31[i])*V.row(F(F31[i],1));
+        isoV.row(K+b+i) = (1.-t23(F31[i],I31[i]))*V.row(F(F31[i],1))
+            + t23(F31[i],I31[i])*V.row(F(F31[i],2));
+    }
+
+    isoE.resize(K,2);
+    for(int i=0; i<K; ++i)
+        isoE.row(i) << i, K+i;
+
+}
+
diff --git a/tutorial/712_DataSmoothing/main.cpp b/tutorial/712_DataSmoothing/main.cpp
new file mode 100755
index 000000000..58aeef323
--- /dev/null
+++ b/tutorial/712_DataSmoothing/main.cpp
@@ -0,0 +1,101 @@
+#include <igl/read_triangle_mesh.h>
+#include <igl/hessian_energy.h>
+#include <igl/massmatrix.h>
+#include <igl/cotmatrix.h>
+#include <igl/jet.h>
+#include <igl/edges.h>
+#include <igl/components.h>
+#include <igl/remove_unreferenced.h>
+#include <igl/viewer/Viewer.h>
+
+#include <Eigen/Core>
+#include <Eigen/SparseCholesky>
+
+#include <iostream>
+#include <set>
+#include <limits>
+#include <stdlib.h>
+
+#include "tutorial_shared_path.h"
+
+#include "isolines.h"
+
+
+int main(int argc, char * argv[])
+{
+    typedef Eigen::SparseMatrix<double> SparseMat;
+
+    //Read our mesh
+    Eigen::MatrixXd V;
+    Eigen::MatrixXi F, E;
+    if(!igl::read_triangle_mesh(
+        argc>1?argv[1]: TUTORIAL_SHARED_PATH "/beetle.off",V,F)) {
+        std::cout << "Failed to load mesh." << std::endl;
+    }
+    igl::edges(F,E);
+
+    //Constructing an exact function to smooth
+    Eigen::VectorXd zexact = V.block(0,2,V.rows(),1).array()
+        + 0.5*V.block(0,1,V.rows(),1).array()
+        + V.block(0,1,V.rows(),1).array().pow(2)
+        + V.block(0,2,V.rows(),1).array().pow(3);
+
+    //Make the exact function noisy
+    srand(0);
+    const double s = 0.2*(zexact.maxCoeff() - zexact.minCoeff());
+    Eigen::VectorXd znoisy = zexact + s*Eigen::VectorXd::Random(zexact.size());
+
+    //Constructing the squared Laplacian and squared Hessian energy
+    SparseMat L, M;
+    igl::cotmatrix(V, F, L);
+    igl::massmatrix(V, F, igl::MASSMATRIX_TYPE_BARYCENTRIC, M);
+    Eigen::SimplicialLDLT<SparseMat> solver(M);
+    SparseMat MinvL = solver.solve(L);
+    SparseMat QL = L.transpose()*MinvL;
+    SparseMat QH;
+    igl::hessian_energy(V, F, QH);
+
+    //Solve to find Laplacian-smoothed and Hessian-smoothed solutions
+    const double al = 5e-4;
+    Eigen::SimplicialLDLT<SparseMat> lapSolver(al*QL + (1.-al)*M);
+    Eigen::VectorXd zl = lapSolver.solve(al*M*znoisy);
+    const double ah = 3e-6;
+    Eigen::SimplicialLDLT<SparseMat> hessSolver(ah*QH + (1.-ah)*M);
+    Eigen::VectorXd zh = hessSolver.solve(ah*M*znoisy);
+
+    //Viewer that shows all functions: zexact, znoisy, zl, zh
+    igl::viewer::Viewer viewer;
+    viewer.data.set_mesh(V,F);
+    viewer.core.show_lines = false;
+    viewer.callback_key_down =
+      [&](igl::viewer::Viewer & viewer, unsigned char key, int mod)->bool {
+        //Graduate result to show isolines, then compute color matrix
+        const Eigen::VectorXd* z;
+        switch(key) {
+            case '1':
+                z = &zexact;
+                break;
+            case '2':
+                z = &znoisy;
+                break;
+            case '3':
+                z = &zl;
+                break;
+            case '4':
+                z = &zh;
+                break;
+            default:
+                return false;
+        }
+        Eigen::MatrixXd isoV;
+        Eigen::MatrixXi isoE;
+        isolines(V, F, *z, 30, isoV, isoE);
+        viewer.data.set_edges(isoV,isoE,Eigen::RowVector3d(0,0,0));
+        Eigen::MatrixXd colors;
+        igl::jet(*z, true, colors);
+        viewer.data.set_colors(colors);
+    };
+    viewer.launch();
+
+    return 0;
+}