diff --git a/include/igl/opengl/MeshGL.cpp b/include/igl/opengl/MeshGL.cpp
index d5c732114..621f10fea 100644
--- a/include/igl/opengl/MeshGL.cpp
+++ b/include/igl/opengl/MeshGL.cpp
@@ -13,7 +13,8 @@
 #include <iostream>
 
 IGL_INLINE igl::opengl::MeshGL::MeshGL():
-  tex_filter(GL_LINEAR)
+  tex_filter(GL_LINEAR),
+  tex_wrap(GL_REPEAT)
 {
 }
 
@@ -93,8 +94,8 @@ IGL_INLINE void igl::opengl::MeshGL::bind_mesh()
   glBindTexture(GL_TEXTURE_2D, vbo_tex);
   if (dirty & MeshGL::DIRTY_TEXTURE)
   {
-    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
-    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_REPEAT);
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, tex_wrap);
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, tex_wrap);
     glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, tex_filter);
     glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, tex_filter);
     glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
diff --git a/include/igl/opengl/MeshGL.h b/include/igl/opengl/MeshGL.h
index 6b89d57f0..e5ecc76c5 100644
--- a/include/igl/opengl/MeshGL.h
+++ b/include/igl/opengl/MeshGL.h
@@ -84,6 +84,7 @@ public:
   int tex_u;
   int tex_v;
   GLint tex_filter;
+  GLint tex_wrap;
   Eigen::Matrix<char,Eigen::Dynamic,1> tex;
 
   Eigen::Matrix<unsigned, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> F_vbo;
diff --git a/include/igl/opengl/ViewerData.cpp b/include/igl/opengl/ViewerData.cpp
index 805430ea9..404de3713 100644
--- a/include/igl/opengl/ViewerData.cpp
+++ b/include/igl/opengl/ViewerData.cpp
@@ -141,11 +141,12 @@ IGL_INLINE void igl::opengl::ViewerData::set_colors(const Eigen::MatrixXd &C)
 {
   using namespace std;
   using namespace Eigen;
+  // This Gouraud coloring should be deprecated in favor of Phong coloring in
+  // set-data
   if(C.rows()>0 && C.cols() == 1)
   {
-    Eigen::MatrixXd C3;
-    igl::parula(C,true,C3);
-    return set_colors(C3);
+    assert(false && "deprecated: call set_data directly instead");
+    return set_data(C);
   }
   // Ambient color should be darker color
   const auto ambient = [](const MatrixXd & C)->MatrixXd
@@ -263,9 +264,25 @@ IGL_INLINE void igl::opengl::ViewerData::set_texture(
   dirty |= MeshGL::DIRTY_TEXTURE;
 }
 
-IGL_INLINE void igl::opengl::ViewerData::set_data(const Eigen::VectorXd & D)
+IGL_INLINE void igl::opengl::ViewerData::set_data(
+  const double caxis_min,
+  const double caxis_max,
+  const Eigen::VectorXd & D)
 {
-  set_uv((D/D.maxCoeff()).replicate(1,2));
+  if(!show_texture)
+  {
+    Eigen::MatrixXd CM;
+    igl::parula(Eigen::VectorXd::LinSpaced(21,0,1).eval(),false,CM);
+    set_colormap(CM);
+  } 
+  set_uv(((D.array()-caxis_min)/(caxis_max-caxis_min)).replicate(1,2));
+}
+
+IGL_INLINE void igl::opengl::ViewerData::set_data( const Eigen::VectorXd & D)
+{
+  const double caxis_min = D.minCoeff();
+  const double caxis_max = D.maxCoeff();
+  return set_data(caxis_min,caxis_max,D);
 }
 
 IGL_INLINE void igl::opengl::ViewerData::set_colormap(const Eigen::MatrixXd & CM)
@@ -282,6 +299,7 @@ IGL_INLINE void igl::opengl::ViewerData::set_colormap(const Eigen::MatrixXd & CM
   set_texture(R,G,B);
   show_texture = true;
   meshgl.tex_filter = GL_NEAREST;
+  meshgl.tex_wrap = GL_CLAMP_TO_EDGE;
 }
 
 IGL_INLINE void igl::opengl::ViewerData::set_points(
@@ -431,6 +449,7 @@ IGL_INLINE void igl::opengl::ViewerData::clear()
   labels_strings.clear();
 
   face_based = false;
+  show_texture = false;
 }
 
 IGL_INLINE void igl::opengl::ViewerData::compute_normals()
diff --git a/include/igl/opengl/ViewerData.h b/include/igl/opengl/ViewerData.h
index 3b293f151..7d28724da 100644
--- a/include/igl/opengl/ViewerData.h
+++ b/include/igl/opengl/ViewerData.h
@@ -98,8 +98,16 @@ public:
   // Set pseudo-colorable scalar data associated with the mesh.
   //
   // Inputs:
+  //   caxis_min  caxis minimum bound
+  //   caxis_max  caxis maximum bound
   //   D  #V by 3 list of colors
+  //
   // To-do: support #F by 1 per-face data
+  IGL_INLINE void set_data(
+    const double caxis_min, 
+    const double caxis_max, 
+    const Eigen::VectorXd & D);
+  // Use min(D) and max(D) to set caxis.
   IGL_INLINE void set_data(const Eigen::VectorXd & D);
   // Not to be confused with set_colors, this creates a _texture_ that will be
   // referenced to pseudocolor accordint to the scalar field passed to set_data.
diff --git a/tutorial/104_Colors/main.cpp b/tutorial/104_Colors/main.cpp
index c79d1e460..66ee7ff8b 100755
--- a/tutorial/104_Colors/main.cpp
+++ b/tutorial/104_Colors/main.cpp
@@ -1,6 +1,5 @@
 #include <igl/readOFF.h>
 #include <igl/opengl/glfw/Viewer.h>
-#include <igl/jet.h>
 #include "tutorial_shared_path.h"
 
 Eigen::MatrixXd V;
@@ -16,11 +15,10 @@ int main(int argc, char *argv[])
   igl::opengl::glfw::Viewer viewer;
   viewer.data().set_mesh(V, F);
 
-  // Use the z coordinate as a scalar field over the surface
-  Eigen::VectorXd Z = V.col(2);
-
-  // Compute per-vertex colors
-  igl::jet(Z,true,C);
+  // Use the (normalized) vertex positions as colors
+  C = 
+    (V.rowwise()            - V.colwise().minCoeff()).array().rowwise()/
+    (V.colwise().maxCoeff() - V.colwise().minCoeff()).array();
 
   // Add per-vertex colors
   viewer.data().set_colors(C);
diff --git a/tutorial/202_GaussianCurvature/main.cpp b/tutorial/202_GaussianCurvature/main.cpp
index 734888d7f..bbf882aa3 100644
--- a/tutorial/202_GaussianCurvature/main.cpp
+++ b/tutorial/202_GaussianCurvature/main.cpp
@@ -3,7 +3,6 @@
 #include <igl/invert_diag.h>
 #include <igl/readOFF.h>
 #include <igl/opengl/glfw/Viewer.h>
-#include <igl/jet.h>
 #include "tutorial_shared_path.h"
 
 int main(int argc, char *argv[])
@@ -24,13 +23,9 @@ int main(int argc, char *argv[])
   // Divide by area to get integral average
   K = (Minv*K).eval();
 
-  // Compute pseudocolor
-  MatrixXd C;
-  igl::jet(K,true,C);
-
   // Plot the mesh with pseudocolors
   igl::opengl::glfw::Viewer viewer;
   viewer.data().set_mesh(V, F);
-  viewer.data().set_colors(C);
+  viewer.data().set_data(K);
   viewer.launch();
 }
diff --git a/tutorial/203_CurvatureDirections/main.cpp b/tutorial/203_CurvatureDirections/main.cpp
index 5dc859ffa..73914aa86 100755
--- a/tutorial/203_CurvatureDirections/main.cpp
+++ b/tutorial/203_CurvatureDirections/main.cpp
@@ -46,11 +46,7 @@ int main(int argc, char *argv[])
   igl::opengl::glfw::Viewer viewer;
   viewer.data().set_mesh(V, F);
 
-
-  // Compute pseudocolor
-  MatrixXd C;
-  igl::parula(H,true,C);
-  viewer.data().set_colors(C);
+  viewer.data().set_data(H);
 
   // Average edge length for sizing
   const double avg = igl::avg_edge_length(V,F);
diff --git a/tutorial/204_Gradient/main.cpp b/tutorial/204_Gradient/main.cpp
index 45c9cf443..888d93a75 100755
--- a/tutorial/204_Gradient/main.cpp
+++ b/tutorial/204_Gradient/main.cpp
@@ -35,12 +35,7 @@ int main(int argc, char *argv[])
   igl::opengl::glfw::Viewer viewer;
   viewer.data().set_mesh(V, F);
 
-  // Compute pseudocolor for original function
-  MatrixXd C;
-  igl::jet(U,true,C);
-  // // Or for gradient magnitude
-  //igl::jet(GU_mag,true,C);
-  viewer.data().set_colors(C);
+  viewer.data().set_data(U);
 
   // Average edge length divided by average gradient (for scaling)
   const double max_size = igl::avg_edge_length(V,F) / GU_mag.mean();
diff --git a/tutorial/206_GeodesicDistance/main.cpp b/tutorial/206_GeodesicDistance/main.cpp
index c88cb663b..8041583de 100755
--- a/tutorial/206_GeodesicDistance/main.cpp
+++ b/tutorial/206_GeodesicDistance/main.cpp
@@ -1,8 +1,9 @@
 #include <igl/readOBJ.h>
 #include <igl/opengl/glfw/Viewer.h>
 #include <igl/exact_geodesic.h>
-#include <igl/colormap.h>
 #include <igl/unproject_onto_mesh.h>
+#include <igl/parula.h>
+#include <igl/isolines_map.h>
 #include <igl/PI.h>
 #include <iostream>
 #include "tutorial_shared_path.h"
@@ -29,15 +30,12 @@ int main(int argc, char *argv[])
     Eigen::VectorXd d;
     std::cout<<"Computing geodesic distance to vertex "<<vid<<"..."<<std::endl;
     igl::exact_geodesic(V,F,VS,FS,VT,FT,d);
-    const double strip_size = 0.05;
-    // The function should be 1 on each integer coordinate
-    d = (d/strip_size*igl::PI).array().sin().abs().eval();
-    // Compute per-vertex colors
-    Eigen::MatrixXd C;
-    igl::colormap(igl::COLOR_MAP_TYPE_INFERNO,d,false,C);
     // Plot the mesh
-    viewer.data().set_mesh(V, F);
-    viewer.data().set_colors(C);
+    Eigen::MatrixXd CM;
+    igl::parula(Eigen::VectorXd::LinSpaced(21,0,1).eval(),false,CM);
+    igl::isolines_map(Eigen::MatrixXd(CM),CM);
+    viewer.data().set_colormap(CM);
+    viewer.data().set_data(d);
   };
 
   // Plot a distance when a vertex is picked
@@ -68,6 +66,7 @@ int main(int argc, char *argv[])
     return false;
   };
   viewer.data().set_mesh(V,F);
+  viewer.data().show_lines = false;
 
   cout << "Click on mesh to define new source.\n" << std::endl;
   update_distance(0);
diff --git a/tutorial/303_LaplaceEquation/main.cpp b/tutorial/303_LaplaceEquation/main.cpp
index 3a7aa1614..6a1b1a91e 100644
--- a/tutorial/303_LaplaceEquation/main.cpp
+++ b/tutorial/303_LaplaceEquation/main.cpp
@@ -60,14 +60,10 @@ int main(int argc, char *argv[])
   igl::min_quad_with_fixed_precompute((-L).eval(),b,Aeq,true,mqwf);
   igl::min_quad_with_fixed_solve(mqwf,B,bc,Beq,Z);
 
-  // Pseudo-color based on solution
-  MatrixXd C;
-  igl::jet(Z,true,C);
-
   // Plot the mesh with pseudocolors
   igl::opengl::glfw::Viewer viewer;
   viewer.data().set_mesh(V, F);
   viewer.data().show_lines = false;
-  viewer.data().set_colors(C);
+  viewer.data().set_data(Z);
   viewer.launch();
 }
diff --git a/tutorial/304_LinearEqualityConstraints/main.cpp b/tutorial/304_LinearEqualityConstraints/main.cpp
index 8b1489d8c..17d78a2b2 100644
--- a/tutorial/304_LinearEqualityConstraints/main.cpp
+++ b/tutorial/304_LinearEqualityConstraints/main.cpp
@@ -61,30 +61,23 @@ int main(int argc, char *argv[])
   }
 
 
-  // Pseudo-color based on solution
-  struct Data{
-    MatrixXd C,C_const;
-  } data;
   // Use same color axes
   const double min_z = std::min(Z.minCoeff(),Z_const.minCoeff());
   const double max_z = std::max(Z.maxCoeff(),Z_const.maxCoeff());
-  igl::jet(      Z,min_z,max_z,data.C);
-  igl::jet(Z_const,min_z,max_z,data.C_const);
 
   // Plot the mesh with pseudocolors
   igl::opengl::glfw::Viewer viewer;
   viewer.data().set_mesh(V, F);
   viewer.data().show_lines = false;
-  viewer.data().set_colors(data.C);
+  viewer.data().set_data(min_z,max_z,Z);
 
   viewer.callback_key_down = 
-    [](igl::opengl::glfw::Viewer& viewer,unsigned char key,int mod)->bool
+    [&Z,&Z_const,&min_z,&max_z](igl::opengl::glfw::Viewer& viewer,unsigned char key,int mod)->bool
     {
       if(key == ' ')
       {
-        Data & data = *static_cast<Data*>(viewer.callback_key_down_data);
         static bool toggle = true;
-        viewer.data().set_colors(toggle?data.C_const:data.C);
+        viewer.data().set_data(min_z,max_z,toggle?Z_const:Z);
         toggle = !toggle;
         return true;
       }else
@@ -92,7 +85,6 @@ int main(int argc, char *argv[])
         return false;
       }
     };
-  viewer.callback_key_down_data = &data;
   cout<<
     "Press [space] to toggle between unconstrained and constrained."<<endl;
   viewer.launch();
diff --git a/tutorial/305_QuadraticProgramming/main.cpp b/tutorial/305_QuadraticProgramming/main.cpp
index 3b142326d..8a79634e0 100644
--- a/tutorial/305_QuadraticProgramming/main.cpp
+++ b/tutorial/305_QuadraticProgramming/main.cpp
@@ -20,10 +20,7 @@ void solve(igl::opengl::glfw::Viewer &viewer)
   igl::active_set_params as;
   as.max_iter = 8;
   igl::active_set(Q,B,b,bc,Aeq,Beq,Aieq,Bieq,lx,ux,as,Z);
-  // Pseudo-color based on solution
-  Eigen::MatrixXd C;
-  igl::jet(Z,0,1,C);
-  viewer.data().set_colors(C);
+  viewer.data().set_data(Z);
 }
 
 bool key_down(igl::opengl::glfw::Viewer &viewer, unsigned char key, int mod)
diff --git a/tutorial/306_EigenDecomposition/main.cpp b/tutorial/306_EigenDecomposition/main.cpp
index a2ad06f96..5fdce75e9 100644
--- a/tutorial/306_EigenDecomposition/main.cpp
+++ b/tutorial/306_EigenDecomposition/main.cpp
@@ -52,20 +52,19 @@ int main(int argc, char * argv[])
         // Rescale eigen vectors for visualization
         VectorXd Z =
           bbd*0.5*U.col(c);
-        Eigen::MatrixXd C;
-        igl::parula(U.col(c).eval(),false,C);
-        c = (c+1)%U.cols();
         if(twod)
         {
           V.col(2) = Z;
+          viewer.data().set_mesh(V,F);
+          viewer.data().compute_normals();
         }
-        viewer.data().set_mesh(V,F);
-        viewer.data().compute_normals();
-        viewer.data().set_colors(C);
+        viewer.data().set_data(U.col(c).eval());
+        c = (c+1)%U.cols();
         return true;
       }
     }
   };
+  viewer.data().set_mesh(V,F);
   viewer.callback_key_down(viewer,' ',0);
   viewer.data().show_lines = false;
   std::cout<<
diff --git a/tutorial/403_BoundedBiharmonicWeights/main.cpp b/tutorial/403_BoundedBiharmonicWeights/main.cpp
index 12c157fc4..7f606bf48 100755
--- a/tutorial/403_BoundedBiharmonicWeights/main.cpp
+++ b/tutorial/403_BoundedBiharmonicWeights/main.cpp
@@ -85,13 +85,6 @@ bool pre_draw(igl::opengl::glfw::Viewer & viewer)
   return false;
 }
 
-void set_color(igl::opengl::glfw::Viewer &viewer)
-{
-  Eigen::MatrixXd C;
-  igl::jet(W.col(selected).eval(),true,C);
-  viewer.data().set_colors(C);
-}
-
 bool key_down(igl::opengl::glfw::Viewer &viewer, unsigned char key, int mods)
 {
   switch(key)
@@ -102,12 +95,12 @@ bool key_down(igl::opengl::glfw::Viewer &viewer, unsigned char key, int mods)
     case '.':
       selected++;
       selected = std::min(std::max(selected,0),(int)W.cols()-1);
-      set_color(viewer);
+      viewer.data().set_data(W.col(selected));
       break;
     case ',':
       selected--;
       selected = std::min(std::max(selected,0),(int)W.cols()-1);
-      set_color(viewer);
+      viewer.data().set_data(W.col(selected));
       break;
   }
   return true;
@@ -162,7 +155,7 @@ int main(int argc, char *argv[])
   // Plot the mesh with pseudocolors
   igl::opengl::glfw::Viewer viewer;
   viewer.data().set_mesh(U, F);
-  set_color(viewer);
+  viewer.data().set_data(W.col(selected));
   viewer.data().set_edges(C,BE,sea_green);
   viewer.data().show_lines = false;
   viewer.data().show_overlay_depth = false;
diff --git a/tutorial/602_Matlab/main.cpp b/tutorial/602_Matlab/main.cpp
index 644e256b5..10cedb9cb 100755
--- a/tutorial/602_Matlab/main.cpp
+++ b/tutorial/602_Matlab/main.cpp
@@ -1,11 +1,14 @@
-#include <igl/jet.h>
+#include "tutorial_shared_path.h"
+
 #include <igl/readOFF.h>
 #include <igl/cotmatrix.h>
 #include <igl/matlab/matlabinterface.h>
 #include <igl/opengl/glfw/Viewer.h>
 #include <iostream>
 
-#include "tutorial_shared_path.h"
+// On mac you may need to issue something like:
+//
+//     PATH=$PATH:/Applications/MATLAB_R2019a.app/bin/ ./tutorial/602_Matlab_bin
 
 // Base mesh
 Eigen::MatrixXd V;
@@ -17,29 +20,13 @@ Engine* engine;
 // Eigenvectors of the laplacian
 Eigen::MatrixXd EV;
 
-void plotEV(igl::opengl::glfw::Viewer& viewer, int id)
-{
-  Eigen::VectorXd v = EV.col(id);
-  v = v.array() - v.minCoeff();
-  v = v.array() / v.maxCoeff();
-
-  // Map to colors using jet colorramp
-  Eigen::MatrixXd C(V.rows(),3);
-  for (unsigned i=0; i<V.rows(); ++i)
-  {
-    double r,g,b;
-    igl::jet(v(i),r,g,b);
-    C.row(i) << r,g,b;
-  }
-
-  viewer.data().set_colors(C);
-}
-
 // This function is called every time a keyboard button is pressed
 bool key_down(igl::opengl::glfw::Viewer& viewer, unsigned char key, int modifier)
 {
   if (key >= '1' && key <= '9')
-    plotEV(viewer,(key - '1') + 1);
+  {
+    viewer.data().set_data(EV.col((key - '1') + 1));
+  }
 
   return false;
 }
@@ -77,7 +64,7 @@ int main(int argc, char *argv[])
   viewer.data().set_mesh(V, F);
 
   // Plot the first non-trivial eigenvector
-  plotEV(viewer,1);
+  viewer.data().set_data(EV.col(1));
 
   // Launch the viewer
   viewer.launch();
diff --git a/tutorial/704_SignedDistance/main.cpp b/tutorial/704_SignedDistance/main.cpp
index e05cd9536..3624c69e5 100755
--- a/tutorial/704_SignedDistance/main.cpp
+++ b/tutorial/704_SignedDistance/main.cpp
@@ -78,14 +78,8 @@ void update_visualization(igl::opengl::glfw::Viewer & viewer)
     // Bunny is a watertight mesh so use pseudonormal for signing
     signed_distance_pseudonormal(V_vis,V,F,tree,FN,VN,EN,EMAP,S_vis,I,C,N);
   }
-  // push to [0,1] range
-  S_vis.array() = 0.5*(S_vis.array()/max_distance)+0.5;
-  MatrixXd C_vis;
-  // color without normalizing
-  igl::parula(S_vis,false,C_vis);
 
-
-  const auto & append_mesh = [&C_vis,&F_vis,&V_vis](
+  const auto & append_mesh = [&F_vis,&V_vis](
     const Eigen::MatrixXd & V,
     const Eigen::MatrixXi & F,
     const RowVector3d & color)
@@ -94,8 +88,6 @@ void update_visualization(igl::opengl::glfw::Viewer & viewer)
     F_vis.bottomRows(F.rows()) = F.array()+V_vis.rows();
     V_vis.conservativeResize(V_vis.rows()+V.rows(),3);
     V_vis.bottomRows(V.rows()) = V;
-    C_vis.conservativeResize(C_vis.rows()+V.rows(),3);
-    C_vis.bottomRows(V.rows()).rowwise() = color;
   };
   if(overlay)
   {
@@ -103,7 +95,7 @@ void update_visualization(igl::opengl::glfw::Viewer & viewer)
   }
   viewer.data().clear();
   viewer.data().set_mesh(V_vis,F_vis);
-  viewer.data().set_colors(C_vis);
+  viewer.data().set_data(S_vis);
   viewer.core().lighting_factor = overlay;
 }
 
diff --git a/tutorial/712_DataSmoothing/main.cpp b/tutorial/712_DataSmoothing/main.cpp
index 2f863de59..fe9ed0e24 100755
--- a/tutorial/712_DataSmoothing/main.cpp
+++ b/tutorial/712_DataSmoothing/main.cpp
@@ -2,8 +2,8 @@
 #include <igl/hessian_energy.h>
 #include <igl/massmatrix.h>
 #include <igl/cotmatrix.h>
-#include <igl/jet.h>
-#include <igl/edges.h>
+#include <igl/isolines_map.h>
+#include <igl/parula.h>
 #include <igl/vertex_components.h>
 #include <igl/remove_unreferenced.h>
 #include <igl/opengl/glfw/Viewer.h>
@@ -18,7 +18,6 @@
 
 #include "tutorial_shared_path.h"
 
-#include <igl/isolines.h>
 
 
 int main(int argc, char * argv[])
@@ -27,12 +26,11 @@ int main(int argc, char * argv[])
 
     //Read our mesh
     Eigen::MatrixXd V;
-    Eigen::MatrixXi F, E;
+    Eigen::MatrixXi F;
     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()
@@ -88,14 +86,7 @@ int main(int argc, char * argv[])
             default:
                 return false;
         }
-        Eigen::MatrixXd isoV;
-        Eigen::MatrixXi isoE;
-        if(key!='2')
-            igl::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.data().set_data(*z);
         return true;
     };
     std::cout << R"(Usage:
@@ -105,6 +96,11 @@ int main(int argc, char * argv[])
 4  Biharmonic smoothing (natural Hessian boundary)
 
 )";
+    Eigen::MatrixXd CM;
+    igl::parula(Eigen::VectorXd::LinSpaced(21,0,1).eval(),false,CM);
+    igl::isolines_map(Eigen::MatrixXd(CM),CM);
+    viewer.data().set_colormap(CM);
+    viewer.data().set_data(znoisy);
     viewer.launch();
 
     return 0;
diff --git a/tutorial/716_HeatGeodesics/main.cpp b/tutorial/716_HeatGeodesics/main.cpp
index 830dcf9eb..c2947d5fd 100755
--- a/tutorial/716_HeatGeodesics/main.cpp
+++ b/tutorial/716_HeatGeodesics/main.cpp
@@ -4,6 +4,7 @@
 #include <igl/heat_geodesics.h>
 #include <igl/unproject_onto_mesh.h>
 #include <igl/avg_edge_length.h>
+#include <igl/isolines_map.h>
 #include <igl/opengl/glfw/Viewer.h>
 #include <igl/opengl/create_shader_program.h>
 #include <igl/opengl/destroy_shader_program.h>
@@ -21,6 +22,7 @@ void set_colormap(igl::opengl::glfw::Viewer & viewer)
     CM(i,1) = std::max(std::min(2.0*t-1.0,1.0),0.0);
     CM(i,2) = std::max(std::min(6.0*t-5.0,1.0),0.0);
   }
+  igl::isolines_map(Eigen::MatrixXd(CM),CM);
   viewer.data().set_colormap(CM);
 }
 
@@ -56,19 +58,32 @@ int main(int argc, char *argv[])
     if(igl::unproject_onto_mesh(Eigen::Vector2f(x,y), viewer.core().view,
       viewer.core().proj, viewer.core().viewport, V, F, fid, bc))
     {
-      // 3d position of hit
-      const Eigen::RowVector3d m3 =
-        V.row(F(fid,0))*bc(0) + V.row(F(fid,1))*bc(1) + V.row(F(fid,2))*bc(2);
-      int cid = 0;
-      Eigen::Vector3d(
-          (V.row(F(fid,0))-m3).squaredNorm(),
-          (V.row(F(fid,1))-m3).squaredNorm(),
-          (V.row(F(fid,2))-m3).squaredNorm()).minCoeff(&cid);
-      const int vid = F(fid,cid);
-      C.row(vid)<<1,0,0;
-      Eigen::VectorXd D = Eigen::VectorXd::Zero(data.Grad.cols());
-	  D(vid) = 1;
-      igl::heat_geodesics_solve(data,(Eigen::VectorXi(1,1)<<vid).finished(),D);
+      Eigen::VectorXd D;
+      // if big mesh, just use closest vertex. Otherwise, blend distances to
+      // vertices of face using barycentric coordinates.
+      if(F.rows()>100000)
+      {
+        // 3d position of hit
+        const Eigen::RowVector3d m3 =
+          V.row(F(fid,0))*bc(0) + V.row(F(fid,1))*bc(1) + V.row(F(fid,2))*bc(2);
+        int cid = 0;
+        Eigen::Vector3d(
+            (V.row(F(fid,0))-m3).squaredNorm(),
+            (V.row(F(fid,1))-m3).squaredNorm(),
+            (V.row(F(fid,2))-m3).squaredNorm()).minCoeff(&cid);
+        const int vid = F(fid,cid);
+        igl::heat_geodesics_solve(data,(Eigen::VectorXi(1,1)<<vid).finished(),D);
+      }else
+      {
+        D = Eigen::VectorXd::Zero(V.rows());
+        for(int cid = 0;cid<3;cid++)
+        {
+          const int vid = F(fid,cid);
+          Eigen::VectorXd Dc;
+          igl::heat_geodesics_solve(data,(Eigen::VectorXi(1,1)<<vid).finished(),Dc);
+          D += Dc*bc(cid);
+        }
+      }
       viewer.data().set_data(D);
       return true;
     }