diff --git a/include/polyscope/render/opengl/shaders/volume_mesh_shaders.h b/include/polyscope/render/opengl/shaders/volume_mesh_shaders.h
new file mode 100644
index 00000000..e0b89487
--- /dev/null
+++ b/include/polyscope/render/opengl/shaders/volume_mesh_shaders.h
@@ -0,0 +1,22 @@
+#pragma once
+
+#include "polyscope/render/opengl/gl_shaders.h"
+
+namespace polyscope {
+namespace render {
+namespace backend_openGL3_glfw {
+
+// High level pipeline
+extern const ShaderStageSpecification SLICE_TETS_VERT_SHADER;
+extern const ShaderStageSpecification SLICE_TETS_GEOM_SHADER;
+extern const ShaderStageSpecification SLICE_TETS_FRAG_SHADER;
+extern const ShaderReplacementRule SLICE_TETS_BASECOLOR_SHADE;
+extern const ShaderReplacementRule SLICE_TETS_MESH_WIREFRAME;
+extern const ShaderReplacementRule SLICE_TETS_PROPAGATE_VALUE;
+extern const ShaderReplacementRule SLICE_TETS_PROPAGATE_VECTOR;
+extern const ShaderReplacementRule SLICE_TETS_VECTOR_COLOR;
+
+
+} // namespace backend_openGL3_glfw
+} // namespace render
+} // namespace polyscope
diff --git a/include/polyscope/slice_plane.h b/include/polyscope/slice_plane.h
index 5f567f81..12509e95 100644
--- a/include/polyscope/slice_plane.h
+++ b/include/polyscope/slice_plane.h
@@ -8,6 +8,7 @@
 
 namespace polyscope {
 
+
 class SlicePlane {
 
 public:
@@ -23,17 +24,21 @@ class SlicePlane {
 
   void buildGUI();
   void draw();
+  void drawGeometry();
+  void resetVolumeSliceProgram();
+  void ensureVolumeInspectValid();
 
   void setSceneObjectUniforms(render::ShaderProgram& p,
                               bool alwaysPass = false); // if alwaysPass, fake values are given so the plane does
                                                         // nothing (regardless of this plane's active setting)
+  void setSliceGeomUniforms(render::ShaderProgram& p);
 
   const std::string name;
   const std::string postfix;
-  
+
   // Set the position and orientation of the plane
   // planePosition is any 3D position which the plane touches (the center of the plane)
-  // planeNormal is a vector giving the normal direction of the plane, objects 
+  // planeNormal is a vector giving the normal direction of the plane, objects
   // in this negative side of the plane will be culled
   void setPose(glm::vec3 planePosition, glm::vec3 planeNormal);
 
@@ -44,45 +49,57 @@ class SlicePlane {
 
   bool getDrawPlane();
   void setDrawPlane(bool newVal);
-  
+
   bool getDrawWidget();
   void setDrawWidget(bool newVal);
 
   glm::mat4 getTransform();
   void setTransform(glm::mat4 newTransform);
-  
+
   void setColor(glm::vec3 newVal);
   glm::vec3 getColor();
-  
+
   void setGridLineColor(glm::vec3 newVal);
   glm::vec3 getGridLineColor();
 
   void setTransparency(double newVal);
   double getTransparency();
 
+  void setVolumeMeshToInspect(std::string meshName);
+  std::string getVolumeMeshToInspect();
+
 protected:
   // = State
-  PersistentValue<bool> active;    // is it actually slicing?
-  PersistentValue<bool> drawPlane; // do we draw the plane onscreen?
+  PersistentValue<bool> active;     // is it actually slicing?
+  PersistentValue<bool> drawPlane;  // do we draw the plane onscreen?
   PersistentValue<bool> drawWidget; // do we draw the widget onscreen?
   PersistentValue<glm::mat4> objectTransform;
   PersistentValue<glm::vec3> color;
   PersistentValue<glm::vec3> gridLineColor;
   PersistentValue<float> transparency;
 
+  // DON'T make these persistent, because it is unintitive to re-add a scene slice plane and have it immediately start
+  // slicing
+  bool shouldInspectMesh;
+  std::string inspectedMeshName;
+
+  std::shared_ptr<render::ShaderProgram> volumeInspectProgram;
+
   // Widget that wraps the transform
   TransformationGizmo transformGizmo;
 
   std::shared_ptr<render::ShaderProgram> planeProgram;
 
   // Helpers
+  void setSliceAttributes(render::ShaderProgram& p);
+  void createVolumeSliceProgram();
   void prepare();
   glm::vec3 getCenter();
   glm::vec3 getNormal();
   void updateWidgetEnabled();
 };
 
-SlicePlane* addSceneSlicePlane(bool initiallyVisible=false);
+SlicePlane* addSceneSlicePlane(bool initiallyVisible = false);
 void removeLastSceneSlicePlane();
 void buildSlicePlaneGUI();
 
diff --git a/include/polyscope/volume_mesh.h b/include/polyscope/volume_mesh.h
index 17ed3949..8aea696f 100644
--- a/include/polyscope/volume_mesh.h
+++ b/include/polyscope/volume_mesh.h
@@ -6,7 +6,6 @@
 
 #include "polyscope/affine_remapper.h"
 #include "polyscope/color_management.h"
-#include "polyscope/polyscope.h"
 #include "polyscope/render/engine.h"
 #include "polyscope/standardize_data_array.h"
 #include "polyscope/structure.h"
@@ -76,7 +75,7 @@ class VolumeMesh : public QuantityStructure<VolumeMesh> {
   
 	// = Vectors (expect vector array, inner type must be indexable with correct dimension (3 for extrinsic, 2 for intrinsic) 
 	template <class T> VolumeMeshVertexVectorQuantity* addVertexVectorQuantity(std::string name, const T& vectors, VectorType vectorType = VectorType::STANDARD); 
-	template <class T> VolumeMeshCellVectorQuantity* addCellVectorQuantity(std::string name, const T& vectors, VectorType vectorType = VectorType::STANDARD); 
+	template <class T> VolumeMeshCellVectorQuantity* addCellVectorQuantity(std::string name, const T& vectors, VectorType vectorType = VectorType::STANDARD);
 
   // clang-format on
 
@@ -99,7 +98,6 @@ class VolumeMesh : public QuantityStructure<VolumeMesh> {
   size_t faceDataSize;
   size_t cellDataSize;
 
-
   // === Manage the mesh itself
 
   // Core data
@@ -117,15 +115,23 @@ class VolumeMesh : public QuantityStructure<VolumeMesh> {
   std::vector<double> cellAreas;
   std::vector<double> faceAreas;
   std::vector<double> vertexAreas;
-  std::vector<char> faceIsInterior; // a flat array whos order matches the iteration order of the mesh
+  std::vector<char> faceIsInterior; // a flat array whose order matches the iteration order of the mesh
 
   // = Mesh helpers
   VolumeCellType cellType(size_t i) const;
   void computeCounts();       // call to populate counts and indices
   void computeGeometryData(); // call to populate normals/areas/lengths
-  std::vector<std::string> addVolumeMeshRules(std::vector<std::string> initRules, bool withSurfaceShade = true);
+  std::vector<std::string> addVolumeMeshRules(std::vector<std::string> initRules, bool withSurfaceShade = true,
+                                              bool isSlice = false);
   glm::vec3 cellCenter(size_t iC);
 
+  // Manage a separate tetrahedral representation used for volumetric visualizations
+  // (for a pure-tet mesh this will be the same as the cells array)
+  std::vector<std::array<int64_t, 4>> tets;
+  size_t nTets();
+  void computeTets();    // fills tet buffer
+  void ensureHaveTets(); //  ensure the tet buffer is filled (but don't rebuild if already done)
+
   // === Member variables ===
   static const std::string structureTypeName;
 
@@ -152,11 +158,22 @@ class VolumeMesh : public QuantityStructure<VolumeMesh> {
   VolumeMesh* setEdgeWidth(double newVal);
   double getEdgeWidth();
 
+  VolumeMeshVertexScalarQuantity* getLevelSetQuantity();
+  void setLevelSetQuantity(VolumeMeshVertexScalarQuantity* _levelSet);
+
   // Rendering helpers used by quantities
   void setVolumeMeshUniforms(render::ShaderProgram& p);
   void fillGeometryBuffers(render::ShaderProgram& p);
+  void fillSliceGeometryBuffers(render::ShaderProgram& p);
   static const std::vector<std::vector<std::array<size_t, 3>>>& cellStencil(VolumeCellType type);
 
+  // Slice plane listeners
+  std::vector<polyscope::SlicePlane*> volumeSlicePlaneListeners;
+  void addSlicePlaneListener(polyscope::SlicePlane* sp);
+  void removeSlicePlaneListener(polyscope::SlicePlane* sp);
+  void refreshVolumeMeshListeners();
+
+
 private:
   // Visualization settings
   PersistentValue<glm::vec3> color;
@@ -165,6 +182,11 @@ class VolumeMesh : public QuantityStructure<VolumeMesh> {
   PersistentValue<std::string> material;
   PersistentValue<float> edgeWidth;
 
+  // Level sets
+  // TODO: not currently really supported
+  float activeLevelSetValue;
+  VolumeMeshVertexScalarQuantity* activeLevelSetQuantity;
+
   // Do setup work related to drawing, including allocating openGL data
   void prepare();
   void preparePick();
@@ -200,6 +222,8 @@ class VolumeMesh : public QuantityStructure<VolumeMesh> {
   // clang-format off
   static const std::vector<std::vector<std::array<size_t, 3>>> stencilTet;
   static const std::vector<std::vector<std::array<size_t, 3>>> stencilHex;
+  static const std::array<std::array<size_t, 8>, 8> rotationMap;
+  static const std::array<std::array<std::array<size_t, 4>, 6>, 4> diagonalMap;
 
   // clang-format off
 
@@ -214,8 +238,8 @@ class VolumeMesh : public QuantityStructure<VolumeMesh> {
 
   // === Helper implementations
 
-  void setVertexTangentBasisXImpl(const std::vector<glm::vec3>& vectors);
-  void setFaceTangentBasisXImpl(const std::vector<glm::vec3>& vectors);
+  //void setVertexTangentBasisXImpl(const std::vector<glm::vec3>& vectors);
+  //void setFaceTangentBasisXImpl(const std::vector<glm::vec3>& vectors);
   // clang-format on
 };
 
diff --git a/include/polyscope/volume_mesh_color_quantity.h b/include/polyscope/volume_mesh_color_quantity.h
index 765ce653..95bb905c 100644
--- a/include/polyscope/volume_mesh_color_quantity.h
+++ b/include/polyscope/volume_mesh_color_quantity.h
@@ -24,6 +24,7 @@ class VolumeMeshColorQuantity : public VolumeMeshQuantity {
   // UI internals
   const std::string definedOn;
   std::shared_ptr<render::ShaderProgram> program;
+  std::shared_ptr<render::ShaderProgram> sliceProgram;
 
   // Helpers
   virtual void createProgram() = 0;
@@ -38,8 +39,12 @@ class VolumeMeshVertexColorQuantity : public VolumeMeshColorQuantity {
   VolumeMeshVertexColorQuantity(std::string name, std::vector<glm::vec3> values_, VolumeMesh& mesh_);
 
   virtual void createProgram() override;
+  virtual std::shared_ptr<render::ShaderProgram> createSliceProgram() override;
+  void fillSliceColorBuffers(render::ShaderProgram& p);
   void fillColorBuffers(render::ShaderProgram& p);
 
+  virtual void drawSlice(polyscope::SlicePlane *sp) override;
+
   void buildVertexInfoGUI(size_t vInd) override;
 
   // === Members
diff --git a/include/polyscope/volume_mesh_quantity.h b/include/polyscope/volume_mesh_quantity.h
index 6a5b5970..8d4ad993 100644
--- a/include/polyscope/volume_mesh_quantity.h
+++ b/include/polyscope/volume_mesh_quantity.h
@@ -15,6 +15,9 @@ class VolumeMeshQuantity : public Quantity<VolumeMesh> {
 public:
   VolumeMeshQuantity(std::string name, VolumeMesh& parentStructure, bool dominates = false);
   ~VolumeMeshQuantity() {};
+  // virtual std::shared_ptr<render::ShaderProgram> tryCreateSliceProgram(){ return nullptr; }; 
+  virtual std::shared_ptr<render::ShaderProgram> createSliceProgram(){ return nullptr; };
+  virtual void drawSlice(polyscope::SlicePlane *sp){};
 
 public:
   // Build GUI info about this element
diff --git a/include/polyscope/volume_mesh_scalar_quantity.h b/include/polyscope/volume_mesh_scalar_quantity.h
index 432212f2..42e34976 100644
--- a/include/polyscope/volume_mesh_scalar_quantity.h
+++ b/include/polyscope/volume_mesh_scalar_quantity.h
@@ -23,6 +23,7 @@ class VolumeMeshScalarQuantity : public VolumeMeshQuantity, public ScalarQuantit
 protected:
   const std::string definedOn;
   std::shared_ptr<render::ShaderProgram> program;
+  std::shared_ptr<render::ShaderProgram> sliceProgram;
 
   // Helpers
   virtual void createProgram() = 0;
@@ -38,10 +39,31 @@ class VolumeMeshVertexScalarQuantity : public VolumeMeshScalarQuantity {
                               DataType dataType_ = DataType::STANDARD);
 
   virtual void createProgram() override;
+  virtual std::shared_ptr<render::ShaderProgram> createSliceProgram() override;
+  virtual void draw() override;
+  virtual void drawSlice(polyscope::SlicePlane *sp) override;
+
+
+  void setLevelSetValue(float f);
+  void setEnabledLevelSet(bool v);
+  void setLevelSetVisibleQuantity(std::string name);
+  void setLevelSetUniforms(render::ShaderProgram &p);
+  void fillLevelSetData(render::ShaderProgram &p);
+  std::shared_ptr<render::ShaderProgram> levelSetProgram;
 
   void fillColorBuffers(render::ShaderProgram& p);
 
+  void fillSliceColorBuffers(render::ShaderProgram& p);
+
+  virtual void buildCustomUI() override;
   void buildVertexInfoGUI(size_t vInd) override;
+  virtual void refresh() override;
+
+  float levelSetValue;
+  bool isDrawingLevelSet;
+  VolumeMeshVertexScalarQuantity* showQuantity;
+
+
 };
 
 
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index ed098b89..5e5ab626 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -28,6 +28,7 @@ if("${POLYSCOPE_BACKEND_OPENGL3_GLFW}")
     render/opengl/shaders/gizmo_shaders.cpp  
     render/opengl/shaders/histogram_shaders.cpp  
     render/opengl/shaders/surface_mesh_shaders.cpp  
+    render/opengl/shaders/volume_mesh_shaders.cpp  
     render/opengl/shaders/vector_shaders.cpp  
     render/opengl/shaders/sphere_shaders.cpp  
     render/opengl/shaders/ribbon_shaders.cpp  
@@ -85,6 +86,7 @@ if("${POLYSCOPE_BACKEND_OPENGL_MOCK}")
     render/opengl/shaders/gizmo_shaders.cpp  
     render/opengl/shaders/histogram_shaders.cpp  
     render/opengl/shaders/surface_mesh_shaders.cpp  
+    render/opengl/shaders/volume_mesh_shaders.cpp  
     render/opengl/shaders/vector_shaders.cpp  
     render/opengl/shaders/sphere_shaders.cpp  
     render/opengl/shaders/ribbon_shaders.cpp  
diff --git a/src/polyscope.cpp b/src/polyscope.cpp
index cdeee586..d270e360 100644
--- a/src/polyscope.cpp
+++ b/src/polyscope.cpp
@@ -207,9 +207,13 @@ void drawStructures() {
       s.second->draw();
     }
   }
+ 
+  // Also render any slice plane geometry
+  for (SlicePlane* s : state::slicePlanes) {
+    s->drawGeometry();
+  }
 }
 
-
 namespace {
 
 float dragDistSinceLastRelease = 0.0;
@@ -331,6 +335,7 @@ void processInputEvents() {
   }
 }
 
+
 void renderSlicePlanes() {
   for (SlicePlane* s : state::slicePlanes) {
     s->draw();
@@ -791,6 +796,8 @@ bool registerStructure(Structure* s, bool replaceIfPresent) {
 
 Structure* getStructure(std::string type, std::string name) {
 
+  if (type == "" || name == "") return nullptr;
+
   // If there are no structures of that type it is an automatic fail
   if (state::structures.find(type) == state::structures.end()) {
     error("No structures of type " + type + " registered");
diff --git a/src/render/mock_opengl/mock_gl_engine.cpp b/src/render/mock_opengl/mock_gl_engine.cpp
index c1be5848..384f6077 100644
--- a/src/render/mock_opengl/mock_gl_engine.cpp
+++ b/src/render/mock_opengl/mock_gl_engine.cpp
@@ -22,6 +22,7 @@
 #include "polyscope/render/opengl/shaders/surface_mesh_shaders.h"
 #include "polyscope/render/opengl/shaders/texture_draw_shaders.h"
 #include "polyscope/render/opengl/shaders/vector_shaders.h"
+#include "polyscope/render/opengl/shaders/volume_mesh_shaders.h"
 
 
 #include "stb_image.h"
@@ -1365,6 +1366,7 @@ void MockGLEngine::populateDefaultShadersAndRules() {
 
   // == Load general base shaders
   registeredShaderPrograms.insert({"MESH", {{FLEX_MESH_VERT_SHADER, FLEX_MESH_FRAG_SHADER}, DrawMode::Triangles}});
+  registeredShaderPrograms.insert({"SLICE_TETS", {{SLICE_TETS_VERT_SHADER, SLICE_TETS_GEOM_SHADER, SLICE_TETS_FRAG_SHADER}, DrawMode::Triangles}});
   registeredShaderPrograms.insert({"RAYCAST_SPHERE", {{FLEX_SPHERE_VERT_SHADER, FLEX_SPHERE_GEOM_SHADER, FLEX_SPHERE_FRAG_SHADER}, DrawMode::Points}});
   registeredShaderPrograms.insert({"POINT_QUAD", {{FLEX_POINTQUAD_VERT_SHADER, FLEX_POINTQUAD_GEOM_SHADER, FLEX_POINTQUAD_FRAG_SHADER}, DrawMode::Points}});
   registeredShaderPrograms.insert({"RAYCAST_VECTOR", {{FLEX_VECTOR_VERT_SHADER, FLEX_VECTOR_GEOM_SHADER, FLEX_VECTOR_FRAG_SHADER}, DrawMode::Points}});
@@ -1455,6 +1457,14 @@ void MockGLEngine::populateDefaultShadersAndRules() {
   registeredShaderRules.insert({"CYLINDER_PROPAGATE_PICK", CYLINDER_PROPAGATE_PICK});
   registeredShaderRules.insert({"CYLINDER_CULLPOS_FROM_MID", CYLINDER_CULLPOS_FROM_MID});
 
+  // marching tets things
+  registeredShaderRules.insert({"SLICE_TETS_BASECOLOR_SHADE", SLICE_TETS_BASECOLOR_SHADE});
+  registeredShaderRules.insert({"SLICE_TETS_PROPAGATE_VALUE", SLICE_TETS_PROPAGATE_VALUE});
+  registeredShaderRules.insert({"SLICE_TETS_PROPAGATE_VECTOR", SLICE_TETS_PROPAGATE_VECTOR});
+  registeredShaderRules.insert({"SLICE_TETS_VECTOR_COLOR", SLICE_TETS_VECTOR_COLOR});
+  registeredShaderRules.insert({"SLICE_TETS_MESH_WIREFRAME", SLICE_TETS_MESH_WIREFRAME});
+
+
   // clang-format on
 };
 
diff --git a/src/render/opengl/gl_engine.cpp b/src/render/opengl/gl_engine.cpp
index b0d88729..ec21d712 100644
--- a/src/render/opengl/gl_engine.cpp
+++ b/src/render/opengl/gl_engine.cpp
@@ -22,6 +22,7 @@
 #include "polyscope/render/opengl/shaders/surface_mesh_shaders.h"
 #include "polyscope/render/opengl/shaders/texture_draw_shaders.h"
 #include "polyscope/render/opengl/shaders/vector_shaders.h"
+#include "polyscope/render/opengl/shaders/volume_mesh_shaders.h"
 
 #include "stb_image.h"
 
@@ -2054,6 +2055,7 @@ void GLEngine::populateDefaultShadersAndRules() {
 
   // == Load general base shaders
   registeredShaderPrograms.insert({"MESH", {{FLEX_MESH_VERT_SHADER, FLEX_MESH_FRAG_SHADER}, DrawMode::Triangles}});
+  registeredShaderPrograms.insert({"SLICE_TETS", {{SLICE_TETS_VERT_SHADER, SLICE_TETS_GEOM_SHADER, SLICE_TETS_FRAG_SHADER}, DrawMode::Points}});
   registeredShaderPrograms.insert({"RAYCAST_SPHERE", {{FLEX_SPHERE_VERT_SHADER, FLEX_SPHERE_GEOM_SHADER, FLEX_SPHERE_FRAG_SHADER}, DrawMode::Points}});
   registeredShaderPrograms.insert({"POINT_QUAD", {{FLEX_POINTQUAD_VERT_SHADER, FLEX_POINTQUAD_GEOM_SHADER, FLEX_POINTQUAD_FRAG_SHADER}, DrawMode::Points}});
   registeredShaderPrograms.insert({"RAYCAST_VECTOR", {{FLEX_VECTOR_VERT_SHADER, FLEX_VECTOR_GEOM_SHADER, FLEX_VECTOR_FRAG_SHADER}, DrawMode::Points}});
@@ -2143,6 +2145,13 @@ void GLEngine::populateDefaultShadersAndRules() {
   registeredShaderRules.insert({"CYLINDER_PROPAGATE_PICK", CYLINDER_PROPAGATE_PICK});
   registeredShaderRules.insert({"CYLINDER_CULLPOS_FROM_MID", CYLINDER_CULLPOS_FROM_MID});
 
+  // marching tets things
+  registeredShaderRules.insert({"SLICE_TETS_BASECOLOR_SHADE", SLICE_TETS_BASECOLOR_SHADE});
+  registeredShaderRules.insert({"SLICE_TETS_PROPAGATE_VALUE", SLICE_TETS_PROPAGATE_VALUE});
+  registeredShaderRules.insert({"SLICE_TETS_PROPAGATE_VECTOR", SLICE_TETS_PROPAGATE_VECTOR});
+  registeredShaderRules.insert({"SLICE_TETS_VECTOR_COLOR", SLICE_TETS_VECTOR_COLOR});
+  registeredShaderRules.insert({"SLICE_TETS_MESH_WIREFRAME", SLICE_TETS_MESH_WIREFRAME});
+
   // clang-format on
 };
 
diff --git a/src/render/opengl/shaders/rules.cpp b/src/render/opengl/shaders/rules.cpp
index 7fc57a32..10809d85 100644
--- a/src/render/opengl/shaders/rules.cpp
+++ b/src/render/opengl/shaders/rules.cpp
@@ -300,9 +300,10 @@ const ShaderReplacementRule CULL_POS_FROM_VIEW (
     { /* replacement sources */
       {"GLOBAL_FRAGMENT_FILTER_PREP", R"(
         vec3 cullPos = viewPos;
-      )"}
+      )"},
+    },
+    /* uniforms */ {
     },
-    /* uniforms */ {},
     /* attributes */ {},
     /* textures */ {}
 );
diff --git a/src/render/opengl/shaders/volume_mesh_shaders.cpp b/src/render/opengl/shaders/volume_mesh_shaders.cpp
new file mode 100644
index 00000000..0fd675f6
--- /dev/null
+++ b/src/render/opengl/shaders/volume_mesh_shaders.cpp
@@ -0,0 +1,405 @@
+#include "polyscope/render/opengl/shaders/volume_mesh_shaders.h"
+
+namespace polyscope {
+namespace render {
+namespace backend_openGL3_glfw {
+
+const ShaderStageSpecification SLICE_TETS_VERT_SHADER = {
+
+    ShaderStageType::Vertex,
+
+    // uniforms
+    {},
+
+    // attributes
+    {
+        {"a_point_1", DataType::Vector3Float},
+        {"a_slice_1", DataType::Vector3Float},
+        {"a_point_2", DataType::Vector3Float},
+        {"a_slice_2", DataType::Vector3Float},
+        {"a_point_3", DataType::Vector3Float},
+        {"a_slice_3", DataType::Vector3Float},
+        {"a_point_4", DataType::Vector3Float},
+        {"a_slice_4", DataType::Vector3Float},
+    },
+
+    {}, // textures
+
+    // source
+    R"(
+        ${ GLSL_VERSION }$
+        ${ VERT_DECLARATIONS }$
+
+        in vec3 a_point_1;
+        in vec3 a_point_2;
+        in vec3 a_point_3;
+        in vec3 a_point_4;
+        in vec3 a_slice_1;
+        in vec3 a_slice_2;
+        in vec3 a_slice_3;
+        in vec3 a_slice_4;
+        out vec3 point_1;
+        out vec3 point_2;
+        out vec3 point_3;
+        out vec3 point_4;
+        out vec3 slice_1;
+        out vec3 slice_2;
+        out vec3 slice_3;
+        out vec3 slice_4;
+
+        void main()
+        {
+            point_1 = a_point_1;
+            point_2 = a_point_2;
+            point_3 = a_point_3;
+            point_4 = a_point_4;
+            slice_1 = a_slice_1;
+            slice_2 = a_slice_2;
+            slice_3 = a_slice_3;
+            slice_4 = a_slice_4;
+            ${ VERT_ASSIGNMENTS }$
+        }
+)"};
+
+
+const ShaderStageSpecification SLICE_TETS_GEOM_SHADER = {
+
+    ShaderStageType::Geometry,
+
+    // uniforms
+    {
+        {"u_modelView", DataType::Matrix44Float},
+        {"u_projMatrix", DataType::Matrix44Float},
+        {"u_slicePoint", DataType::Float},
+        {"u_sliceVector", DataType::Vector3Float},
+    },
+
+    // attributes
+    {},
+
+    {}, // textures
+
+    // source
+    R"(
+        ${ GLSL_VERSION }$
+
+        layout(points) in;
+        layout(triangle_strip, max_vertices=4) out;
+        uniform mat4 u_modelView;
+        uniform mat4 u_projMatrix;
+        uniform float u_slicePoint;
+        uniform vec3 u_sliceVector;
+        in vec3 point_1[];
+        in vec3 point_2[];
+        in vec3 point_3[];
+        in vec3 point_4[];
+        in vec3 slice_1[];
+        in vec3 slice_2[];
+        in vec3 slice_3[];
+        in vec3 slice_4[];
+        out vec3 a_barycoordToFrag;
+        out vec3 a_normalToFrag;
+        ${ GEOM_DECLARATIONS }$
+
+        void main() {
+
+            vec3 s[4] = vec3[](slice_1[0], slice_2[0], slice_3[0], slice_4[0]);
+            vec3 p[4] = vec3[](point_1[0], point_2[0], point_3[0], point_4[0]);
+            ${ GEOM_INIT_DECLARATIONS }$
+            int ordering[4] = int[](0, 1, 2, 3);
+            
+            float d[4]; 
+            for (int i = 0; i < 4; i++ ) d[i] = dot(u_sliceVector, s[i]) - u_slicePoint;
+
+            vec3 q[4];
+            int n = 0;
+            for( int i = 0; i < 4; i++ ) {
+                for( int j = i+1; j < 4; j++ ) {
+                    if( d[i]*d[j] < 0. ) {
+                        float t = (0-d[i])/(d[j]-d[i]);
+                        ${ GEOM_INTERPOLATE }$
+                        q[n] = ( (1.-t)*p[i] + t*p[j] );
+                        n++;
+                    }
+                }
+            }
+
+            if(n == 4){
+                vec3 cross13 = cross(q[1] - q[0], q[3] - q[0]);
+                vec3 cross23 = cross(q[2] - q[0], q[3] - q[0]);
+                if(dot(cross13, cross23) > 0){
+                    if(dot(cross23, cross23) < dot(cross13, cross13)){
+                        ordering[2] = 3;
+                        ordering[3] = 2;
+                    }else{
+                        ordering[1] = 3;
+                        ordering[3] = 1;
+                    }
+                }
+            }
+
+            vec3 cross12 = cross(q[1] - q[0], q[2] - q[0]);
+            if(dot(cross12, u_sliceVector) < 0){
+                int temp = ordering[1];
+                ordering[1] = ordering[2];
+                ordering[2] = temp;
+                cross12 *= -1;
+            }
+            // Offset slice so that tet edges don't clip through
+            vec3 offset = u_sliceVector * 1e-4;
+            // Emit the vertices as a triangle strip
+            mat4 toScreen = u_projMatrix * u_modelView;
+            for (int i = 0; i < n; i++){
+                a_normalToFrag = cross12;
+                a_barycoordToFrag = vec3(0, 0, 0);
+                a_barycoordToFrag[i % 3] = 1.0;
+                ${ GEOM_ASSIGNMENTS }$
+                gl_Position = toScreen * vec4(q[ordering[i]] - offset, 1.0); 
+                EmitVertex();
+            }
+            EndPrimitive();
+        
+        }
+
+)"};
+
+const ShaderStageSpecification SLICE_TETS_FRAG_SHADER = {
+
+    ShaderStageType::Fragment,
+
+    // uniforms
+    {},
+
+    {}, // attributes
+
+    // textures
+    {},
+
+    // source
+    R"(
+        ${ GLSL_VERSION }$
+        in vec3 a_normalToFrag;
+        in vec3 a_barycoordToFrag;
+        layout(location = 0) out vec4 outputF;
+
+        ${ FRAG_DECLARATIONS }$
+        ${ SLICE_TETS_FRAG_DECLARATIONS }$ 
+
+        void main()
+        {
+           float depth = gl_FragCoord.z;
+           ${ GLOBAL_FRAGMENT_FILTER_PREP }$
+           ${ GLOBAL_FRAGMENT_FILTER }$
+        
+          
+           // Shading
+           ${ GENERATE_SHADE_VALUE }$
+           ${ GENERATE_SHADE_COLOR }$
+           
+           // Handle the wireframe
+           ${ APPLY_WIREFRAME }$
+
+           // Lighting
+           vec3 shadeNormal = a_normalToFrag;
+           ${ PERTURB_SHADE_NORMAL }$
+           ${ GENERATE_LIT_COLOR }$
+
+           // Set alpha
+           float alphaOut = 1.0;
+           ${ GENERATE_ALPHA }$
+           
+           // silly dummy usage to ensure normal and barycoords are always used; otherwise we get errors
+           float dummyVal = a_normalToFrag.x + a_barycoordToFrag.x;
+           alphaOut = alphaOut + dummyVal * (1e-12);
+
+           ${ PERTURB_LIT_COLOR }$
+
+           // Write output
+           outputF = vec4(litColor, alphaOut);
+        }
+)"};
+
+const ShaderReplacementRule SLICE_TETS_BASECOLOR_SHADE(
+    /* rule name */ "SLICE_TETS_BASECOLOR_SHADE",
+    {/* replacement sources */
+     {"FRAG_DECLARATIONS", R"(
+          uniform vec3 u_baseColor1;
+          in float a_faceColorTypeToFrag;
+        )"},
+     {"GENERATE_SHADE_COLOR", R"(
+          vec3 albedoColor = u_baseColor1;
+        )"}},
+    /* uniforms */
+    {
+        {"u_baseColor1", DataType::Vector3Float},
+    },
+    /* attributes */ {},
+    /* textures */ {});
+
+const ShaderReplacementRule SLICE_TETS_MESH_WIREFRAME(
+    /* rule name */ "SLICE_TETS_MESH_WIREFRAME",
+    {
+        /* replacement sources */
+        {"GEOM_DECLARATIONS", R"(
+          out vec3 a_edgeIsRealToFrag;
+        )"},
+        {"GEOM_ASSIGNMENTS", R"(
+          vec3 edgeRealV = vec3(1, 1, 1);
+          a_edgeIsRealToFrag = edgeRealV;
+        )"},
+        {"FRAG_DECLARATIONS", R"(
+          in vec3 a_edgeIsRealToFrag;
+
+          uniform float u_edgeWidth;
+          uniform vec3 u_edgeColor;
+      
+          float getEdgeFactor(vec3 UVW, vec3 edgeReal, float width) {
+            // The Nick Sharp Edge Function. There are many like it, but this one is mine.
+            float slopeWidth = 1.;
+            
+            vec3 fw = fwidth(UVW);
+            vec3 realUVW = max(UVW, 1.0 - edgeReal.yzx);
+            vec3 baryWidth = slopeWidth * fw;
+
+            vec3 end = width*fw;
+            vec3 dist = smoothstep(end - baryWidth, end, realUVW);
+
+            float e = 1.0 - min(min(dist.x, dist.y), dist.z);
+            return e;
+          }
+        )"},
+        {"APPLY_WIREFRAME", R"(
+          float edgeFactor = getEdgeFactor(a_barycoordToFrag, a_edgeIsRealToFrag, u_edgeWidth);
+          albedoColor = mix(albedoColor, u_edgeColor, edgeFactor);
+      )"},
+    },
+    /* uniforms */
+    {
+        {"u_edgeColor", DataType::Vector3Float},
+        {"u_edgeWidth", DataType::Float},
+    },
+    /* attributes */ {},
+    /* textures */ {});
+const ShaderReplacementRule SLICE_TETS_PROPAGATE_VECTOR(
+    /* rule name */ "SLICE_TETS_PROPAGATE_VECTOR",
+    {
+        /* replacement sources */
+        {"VERT_DECLARATIONS", R"(
+          in vec3 a_value_1;
+          in vec3 a_value_2;
+          in vec3 a_value_3;
+          in vec3 a_value_4;
+          out vec3 value_1;
+          out vec3 value_2;
+          out vec3 value_3;
+          out vec3 value_4;
+        )"},
+        {"VERT_ASSIGNMENTS", R"(
+          value_1 = a_value_1;
+          value_2 = a_value_2;
+          value_3 = a_value_3;
+          value_4 = a_value_4;
+        )"},
+        {"GEOM_DECLARATIONS", R"(
+          in vec3 value_1[];
+          in vec3 value_2[];
+          in vec3 value_3[];
+          in vec3 value_4[];
+          out vec3 a_valueToFrag;
+      )"},
+        {"GEOM_INIT_DECLARATIONS", R"(
+          vec3 v[4] = vec3[](value_1[0], value_2[0], value_3[0], value_4[0]);
+          vec3 out_v[4] = vec3[](vec3(0), vec3(0), vec3(0), vec3(0));
+      )"},
+        {"GEOM_INTERPOLATE", R"(
+          out_v[n] = ( (1.-t)*v[i] + t*v[j] );
+      )"},
+        {"GEOM_ASSIGNMENTS", R"(
+          a_valueToFrag = out_v[ordering[i]];
+      )"},
+        {"FRAG_DECLARATIONS", R"(
+          in vec3 a_valueToFrag;
+        )"},
+        {"GENERATE_SHADE_VALUE", R"(
+          vec3 shadeValue = a_valueToFrag;
+        )"},
+    },
+    /* uniforms */ {},
+    /* attributes */
+    {
+        {"a_value_1", DataType::Vector3Float},
+        {"a_value_2", DataType::Vector3Float},
+        {"a_value_3", DataType::Vector3Float},
+        {"a_value_4", DataType::Vector3Float},
+    },
+    /* textures */ {});
+
+const ShaderReplacementRule SLICE_TETS_VECTOR_COLOR(
+    /* rule name */ "SLICE_TETS_VECTOR_COLOR",
+    {/* replacement sources */
+     {"GENERATE_SHADE_COLOR", R"(
+          vec3 albedoColor = shadeValue;
+        )"}},
+    /* uniforms */
+    {
+    },
+    /* attributes */ {},
+    /* textures */ {});
+
+const ShaderReplacementRule SLICE_TETS_PROPAGATE_VALUE(
+    /* rule name */ "SLICE_TETS_PROPAGATE_VALUE",
+    {
+        /* replacement sources */
+        {"VERT_DECLARATIONS", R"(
+          in float a_value_1;
+          in float a_value_2;
+          in float a_value_3;
+          in float a_value_4;
+          out float value_1;
+          out float value_2;
+          out float value_3;
+          out float value_4;
+        )"},
+        {"VERT_ASSIGNMENTS", R"(
+          value_1 = a_value_1;
+          value_2 = a_value_2;
+          value_3 = a_value_3;
+          value_4 = a_value_4;
+        )"},
+        {"GEOM_DECLARATIONS", R"(
+          in float value_1[];
+          in float value_2[];
+          in float value_3[];
+          in float value_4[];
+          out float a_valueToFrag;
+      )"},
+        {"GEOM_INIT_DECLARATIONS", R"(
+          float v[4] = float[](value_1[0], value_2[0], value_3[0], value_4[0]);
+          float out_v[4] = float[](0, 0, 0, 0);
+      )"},
+        {"GEOM_INTERPOLATE", R"(
+          out_v[n] = ( (1.-t)*v[i] + t*v[j] );
+      )"},
+        {"GEOM_ASSIGNMENTS", R"(
+          a_valueToFrag = out_v[ordering[i]];
+      )"},
+        {"FRAG_DECLARATIONS", R"(
+          in float a_valueToFrag;
+        )"},
+        {"GENERATE_SHADE_VALUE", R"(
+          float shadeValue = a_valueToFrag;
+        )"},
+    },
+    /* uniforms */ {},
+    /* attributes */
+    {
+        {"a_value_1", DataType::Float},
+        {"a_value_2", DataType::Float},
+        {"a_value_3", DataType::Float},
+        {"a_value_4", DataType::Float},
+    },
+    /* textures */ {});
+
+} // namespace backend_openGL3_glfw
+} // namespace render
+}; // namespace polyscope
\ No newline at end of file
diff --git a/src/slice_plane.cpp b/src/slice_plane.cpp
index 7431c1d0..7d0b17dc 100644
--- a/src/slice_plane.cpp
+++ b/src/slice_plane.cpp
@@ -1,6 +1,7 @@
 #include "polyscope/slice_plane.h"
 
 #include "polyscope/polyscope.h"
+#include "polyscope/volume_mesh.h"
 
 namespace polyscope {
 
@@ -28,6 +29,9 @@ SlicePlane* addSceneSlicePlane(bool initiallyVisible) {
     sceneSlicePlanes.back()->setDrawPlane(false);
     sceneSlicePlanes.back()->setDrawWidget(false);
   }
+  for (size_t i = 0; i < sceneSlicePlanes.size(); i++) {
+    sceneSlicePlanes[i]->resetVolumeSliceProgram();
+  }
   return sceneSlicePlanes.back();
 }
 
@@ -35,6 +39,9 @@ void removeLastSceneSlicePlane() {
   if (sceneSlicePlanes.empty()) return;
   delete sceneSlicePlanes.back();
   sceneSlicePlanes.pop_back();
+  for (size_t i = 0; i < sceneSlicePlanes.size(); i++) {
+    sceneSlicePlanes[i]->resetVolumeSliceProgram();
+  }
 }
 
 void buildSlicePlaneGUI() {
@@ -60,14 +67,21 @@ void buildSlicePlaneGUI() {
   }
 }
 
+void SlicePlane::setSliceGeomUniforms(render::ShaderProgram& p) {
+  glm::vec3 norm = getNormal();
+  p.setUniform("u_sliceVector", norm);
+  p.setUniform("u_slicePoint", glm::dot(getCenter(), norm));
+}
+
+
 SlicePlane::SlicePlane(std::string name_)
     : name(name_), postfix(std::to_string(state::slicePlanes.size())), active("SlicePlane#" + name + "#active", true),
       drawPlane("SlicePlane#" + name + "#drawPlane", true), drawWidget("SlicePlane#" + name + "#drawWidget", true),
       objectTransform("SlicePlane#" + name + "#object_transform", glm::mat4(1.0)),
       color("SlicePlane#" + name + "#color", getNextUniqueColor()),
       gridLineColor("SlicePlane#" + name + "#gridLineColor", glm::vec3{.97, .97, .97}),
-      transparency("SlicePlane#" + name + "#transparency", 0.5),
-      transformGizmo("SlicePlane#" + name + "#transform_gizmo", objectTransform.get(), &objectTransform) {
+      transparency("SlicePlane#" + name + "#transparency", 0.5), shouldInspectMesh(false), inspectedMeshName(""),
+      transformGizmo("SlicePlane#" + name + "#transformGizmo", objectTransform.get(), &objectTransform) {
   state::slicePlanes.push_back(this);
   render::engine->addSlicePlane(postfix);
   transformGizmo.enabled = true;
@@ -75,6 +89,8 @@ SlicePlane::SlicePlane(std::string name_)
 }
 
 SlicePlane::~SlicePlane() {
+  ensureVolumeInspectValid();
+  setVolumeMeshToInspect(""); // disable any slicing
   render::engine->removeSlicePlane(postfix);
   auto pos = std::find(state::slicePlanes.begin(), state::slicePlanes.end(), this);
   if (pos == state::slicePlanes.end()) return;
@@ -105,28 +121,140 @@ void SlicePlane::prepare() {
   planeProgram->setAttribute("a_position", positions);
 }
 
+void SlicePlane::setVolumeMeshToInspect(std::string meshname) {
+  VolumeMesh* oldMeshToInspect = polyscope::getVolumeMesh(inspectedMeshName);
+  if (oldMeshToInspect != nullptr) {
+    oldMeshToInspect->removeSlicePlaneListener(this);
+  }
+  inspectedMeshName = meshname;
+  VolumeMesh* meshToInspect = polyscope::getVolumeMesh(inspectedMeshName);
+  if (meshToInspect == nullptr) {
+    inspectedMeshName = "";
+    shouldInspectMesh = false;
+    volumeInspectProgram.reset();
+    return;
+  }
+  drawPlane = false;
+  meshToInspect->addSlicePlaneListener(this);
+  meshToInspect->setCullWholeElements(false);
+  meshToInspect->ensureHaveTets(); // do this as early as possible because it is expensive
+  shouldInspectMesh = true;
+  volumeInspectProgram.reset();
+}
+
+std::string SlicePlane::getVolumeMeshToInspect() { return inspectedMeshName; }
+
+void SlicePlane::ensureVolumeInspectValid() {
+  if (!shouldInspectMesh) return;
+
+  // This method exists to save us in any cases where we might be inspecting a volume mesh when that mesh is deleted. We
+  // can't just call setVolumeMeshToInspect(""), because that tries to look up the volume mesh.
+
+  if (!hasVolumeMesh(inspectedMeshName)) {
+    inspectedMeshName = "";
+    shouldInspectMesh = false;
+    volumeInspectProgram = nullptr;
+  }
+}
+
+void SlicePlane::createVolumeSliceProgram() {
+  VolumeMesh* meshToInspect = polyscope::getVolumeMesh(inspectedMeshName);
+  volumeInspectProgram = render::engine->requestShader(
+      "SLICE_TETS", meshToInspect->addVolumeMeshRules({"SLICE_TETS_BASECOLOR_SHADE"}, true, true));
+  meshToInspect->fillSliceGeometryBuffers(*volumeInspectProgram);
+  render::engine->setMaterial(*volumeInspectProgram, meshToInspect->getMaterial());
+}
+
+void SlicePlane::resetVolumeSliceProgram() { volumeInspectProgram.reset(); }
+
+void SlicePlane::setSliceAttributes(render::ShaderProgram& p) {
+  VolumeMesh* meshToInspect = polyscope::getVolumeMesh(inspectedMeshName);
+  std::vector<glm::vec3> point1;
+  std::vector<glm::vec3> point2;
+  std::vector<glm::vec3> point3;
+  std::vector<glm::vec3> point4;
+  size_t cellCount = meshToInspect->nCells();
+  point1.resize(cellCount);
+  point2.resize(cellCount);
+  point3.resize(cellCount);
+  point4.resize(cellCount);
+  for (size_t iC = 0; iC < cellCount; iC++) {
+    const std::array<int64_t, 8>& cell = meshToInspect->cells[iC];
+    point1[iC] = meshToInspect->vertices[cell[0]];
+    point2[iC] = meshToInspect->vertices[cell[1]];
+    point3[iC] = meshToInspect->vertices[cell[2]];
+    point4[iC] = meshToInspect->vertices[cell[3]];
+  }
+  glm::vec3 normal = glm::vec3(-1, 0, 0);
+
+  p.setAttribute("a_slice_1", point1);
+  p.setAttribute("a_slice_2", point2);
+  p.setAttribute("a_slice_3", point3);
+  p.setAttribute("a_slice_4", point4);
+}
+
+void SlicePlane::drawGeometry() {
+  if (!active.get()) return;
+
+  ensureVolumeInspectValid();
+
+  if (shouldInspectMesh) {
+    VolumeMesh* vMesh = polyscope::getVolumeMesh(inspectedMeshName);
+
+    // guard against situations where the volume mesh we are slicing has been deleted
+    if (vMesh == nullptr) {
+      setVolumeMeshToInspect("");
+      return;
+    }
+
+    if (vMesh->wantsCullPosition()) return;
+
+    if (volumeInspectProgram == nullptr) {
+      createVolumeSliceProgram();
+    }
+
+
+    if (vMesh->dominantQuantity == nullptr) {
+      vMesh->setStructureUniforms(*volumeInspectProgram);
+      setSceneObjectUniforms(*volumeInspectProgram, true);
+      setSliceGeomUniforms(*volumeInspectProgram);
+      vMesh->setVolumeMeshUniforms(*volumeInspectProgram);
+      volumeInspectProgram->setUniform("u_baseColor1", vMesh->getColor());
+      volumeInspectProgram->draw();
+    }
+
+    for (auto it = vMesh->quantities.begin(); it != vMesh->quantities.end(); it++) {
+      if (!it->second->isEnabled()) continue;
+      it->second->drawSlice(this);
+    }
+  }
+}
+
+
 void SlicePlane::draw() {
-  if (!drawPlane.get() || !active.get()) return;
-
-  // Set uniforms
-  glm::mat4 viewMat = view::getCameraViewMatrix();
-  planeProgram->setUniform("u_viewMatrix", glm::value_ptr(viewMat));
-  glm::mat4 projMat = view::getCameraPerspectiveMatrix();
-  planeProgram->setUniform("u_projMatrix", glm::value_ptr(projMat));
-
-  planeProgram->setUniform("u_objectMatrix", glm::value_ptr(objectTransform.get()));
-  planeProgram->setUniform("u_lengthScale", state::lengthScale);
-  planeProgram->setUniform("u_color", getColor());
-  planeProgram->setUniform("u_gridLineColor", getGridLineColor());
-  planeProgram->setUniform("u_transparency", transparency.get());
-
-  // glm::vec3 center{objectTransform.get()[3][0], objectTransform.get()[3][1], objectTransform.get()[3][2]};
-  // planeProgram->setUniform("u_center", center);
-
-  render::engine->setDepthMode(DepthMode::Less);
-  render::engine->setBackfaceCull(false);
-  render::engine->applyTransparencySettings();
-  planeProgram->draw();
+  if (!active.get()) return;
+
+  if (drawPlane.get()) {
+    // Set uniforms
+    glm::mat4 viewMat = view::getCameraViewMatrix();
+    planeProgram->setUniform("u_viewMatrix", glm::value_ptr(viewMat));
+    glm::mat4 projMat = view::getCameraPerspectiveMatrix();
+    planeProgram->setUniform("u_projMatrix", glm::value_ptr(projMat));
+
+    planeProgram->setUniform("u_objectMatrix", glm::value_ptr(objectTransform.get()));
+    planeProgram->setUniform("u_lengthScale", state::lengthScale);
+    planeProgram->setUniform("u_color", color.get());
+    planeProgram->setUniform("u_gridLineColor", getGridLineColor());
+    planeProgram->setUniform("u_transparency", transparency.get());
+
+    // glm::vec3 center{objectTransform.get()[3][0], objectTransform.get()[3][1], objectTransform.get()[3][2]};
+    // planeProgram->setUniform("u_center", center);
+
+    render::engine->setDepthMode(DepthMode::Less);
+    render::engine->setBackfaceCull(false);
+    render::engine->applyTransparencySettings();
+    planeProgram->draw();
+  }
 }
 
 void SlicePlane::buildGUI() {
@@ -135,7 +263,7 @@ void SlicePlane::buildGUI() {
   if (ImGui::Checkbox(name.c_str(), &active.get())) {
     setActive(getActive());
   }
-  
+
   ImGui::SameLine();
 
   { // Color transparency box
@@ -161,12 +289,45 @@ void SlicePlane::buildGUI() {
   if (ImGui::Checkbox("draw widget", &drawWidget.get())) {
     setDrawWidget(getDrawWidget());
   }
+
+  bool haveVolumeMeshes = state::structures.find("Volume Mesh") != state::structures.end();
+
+  if (haveVolumeMeshes) {
+
+    if (ImGui::Button("Inspect")) {
+      ImGui::OpenPopup("InspectPopup");
+    }
+    if (ImGui::BeginPopup("InspectPopup")) {
+
+      //  Loop over volume meshes and offer them to be inspected
+      std::map<std::string, Structure*>::iterator it;
+      for (it = state::structures["Volume Mesh"].begin(); it != state::structures["Volume Mesh"].end(); it++) {
+        std::string vMeshName = it->first;
+        if (ImGui::MenuItem(vMeshName.c_str(), NULL, inspectedMeshName == vMeshName)) {
+          setVolumeMeshToInspect(vMeshName);
+        }
+      }
+
+      // "None" option
+      if (ImGui::MenuItem("None", NULL, inspectedMeshName == "")) {
+        setVolumeMeshToInspect("");
+      }
+
+      ImGui::EndPopup();
+    }
+  }
+
+
   ImGui::Unindent(16.);
 
   ImGui::PopID();
 }
 
 void SlicePlane::setSceneObjectUniforms(render::ShaderProgram& p, bool alwaysPass) {
+  if (!p.hasUniform("u_slicePlaneNormal_" + postfix)) {
+    return;
+  }
+
   glm::vec3 normal, center;
 
   if (alwaysPass) {
diff --git a/src/structure.cpp b/src/structure.cpp
index 0f4e6773..4a20600b 100644
--- a/src/structure.cpp
+++ b/src/structure.cpp
@@ -232,8 +232,10 @@ void Structure::setStructureUniforms(render::ShaderProgram& p) {
   glm::mat4 viewMat = getModelView();
   p.setUniform("u_modelView", glm::value_ptr(viewMat));
 
-  glm::mat4 projMat = view::getCameraPerspectiveMatrix();
-  p.setUniform("u_projMatrix", glm::value_ptr(projMat));
+  if (p.hasUniform("u_projMatrix")) {
+    glm::mat4 projMat = view::getCameraPerspectiveMatrix();
+    p.setUniform("u_projMatrix", glm::value_ptr(projMat));
+  }
 
   if (render::engine->transparencyEnabled()) {
     if (p.hasUniform("u_transparency")) {
diff --git a/src/volume_mesh.cpp b/src/volume_mesh.cpp
index e2656e10..6dcbebd7 100644
--- a/src/volume_mesh.cpp
+++ b/src/volume_mesh.cpp
@@ -6,15 +6,15 @@
 #include "polyscope/pick.h"
 #include "polyscope/polyscope.h"
 #include "polyscope/render/engine.h"
+#include "polyscope/volume_mesh_quantity.h"
 
 #include "imgui.h"
 
+#include <algorithm>
+#include <numeric>
 #include <unordered_map>
 #include <utility>
 
-using std::cout;
-using std::endl;
-
 namespace polyscope {
 
 // Initialize statics
@@ -29,6 +29,57 @@ const std::vector<std::vector<std::array<size_t, 3>>> VolumeMesh::stencilTet =
    {{1,2,3}},
  };
 
+// Indirection to place vertex 0 always in the bottom left corner
+const std::array<std::array<size_t, 8>, 8> VolumeMesh::rotationMap = 
+ {{
+   {0, 1, 2, 3, 4, 5, 7, 6}, 
+   {1, 0, 4, 5, 2, 3, 6, 7}, 
+   {2, 1, 5, 6, 3, 0, 7, 4}, 
+   {3, 0, 1, 2, 7, 4, 6, 5}, 
+   {4, 0, 3, 7, 5, 1, 6, 2}, 
+   {5, 1, 0, 4, 7, 2, 6, 3}, 
+   {7, 3, 2, 6, 4, 0, 5, 1},
+   {6, 2, 1, 5, 7, 3, 4, 0} 
+ }};
+
+// Map indirected cube to tets
+const std::array<std::array<std::array<size_t, 4>, 6>, 4> VolumeMesh::diagonalMap = 
+ {{
+    {{
+      {0, 1, 2, 5},
+      {0, 2, 6, 5},
+      {0, 2, 3, 6},
+      {0, 5, 6, 4},
+      {2, 6, 5, 7},
+      {0, 0, 0, 0}
+    }},
+    {{
+      {0, 5, 6, 4},
+      {0, 1, 6, 5},
+      {1, 7, 6, 5},
+      {0, 6, 2, 3},
+      {0, 6, 1, 2},
+      {1, 6, 7, 2}
+    }},
+    {{
+      {0, 4, 5, 7},
+      {0, 3, 6, 7},
+      {0, 6, 4, 7},
+      {0, 1, 2, 5},
+      {0, 3, 7, 2},
+      {0, 7, 5, 2}
+    }},
+    {{
+      {0, 2, 3, 7},
+      {0, 3, 6, 7},
+      {0, 6, 4, 7},
+      {0, 5, 7, 4},
+      {1, 5, 7, 0},
+      {1, 7, 2, 0}
+    }}
+ }};
+
+
 const std::vector<std::vector<std::array<size_t, 3>>> VolumeMesh::stencilHex = 
   // numbered like in this diagram, except with 6/7 swapped
   // https://vtk.org/wp-content/uploads/2015/04/file-formats.pdf
@@ -48,7 +99,7 @@ VolumeMesh::VolumeMesh(std::string name, const std::vector<glm::vec3>& vertexPos
       color(uniquePrefix() + "color", getNextUniqueColor()),
       interiorColor(uniquePrefix() + "interiorColor", color.get()),
       edgeColor(uniquePrefix() + "edgeColor", glm::vec3{0., 0., 0.}), material(uniquePrefix() + "material", "clay"),
-      edgeWidth(uniquePrefix() + "edgeWidth", 0.) {
+      edgeWidth(uniquePrefix() + "edgeWidth", 0.), activeLevelSetQuantity(nullptr) {
   cullWholeElements.setPassive(true);
 
   // set the interior color to be a desaturated version of the normal one
@@ -56,12 +107,190 @@ VolumeMesh::VolumeMesh(std::string name, const std::vector<glm::vec3>& vertexPos
   desatColorHSV.y *= 0.3;
   interiorColor.setPassive(HSVtoRGB(desatColorHSV));
 
-
   updateObjectSpaceBounds();
   computeCounts();
   computeGeometryData();
 }
 
+void VolumeMesh::computeTets() {
+  // Algorithm from
+  // https://www.researchgate.net/profile/Julien-Dompierre/publication/221561839_How_to_Subdivide_Pyramids_Prisms_and_Hexahedra_into_Tetrahedra/links/0912f509c0b7294059000000/How-to-Subdivide-Pyramids-Prisms-and-Hexahedra-into-Tetrahedra.pdf?origin=publication_detail
+  // It's a bit hard to look at but it works
+  // Uses vertex numberings to ensure consistent diagonals between faces, and keeps tet counts to 5 or 6 per hex
+  size_t tetCount = 0;
+  // Get number of tets first
+  for (size_t iC = 0; iC < nCells(); iC++) {
+    switch (cellType(iC)) {
+    case VolumeCellType::HEX: {
+      std::array<size_t, 8> sortedNumbering;
+      std::iota(sortedNumbering.begin(), sortedNumbering.end(), 0);
+      std::sort(sortedNumbering.begin(), sortedNumbering.end(),
+                [this, iC](size_t a, size_t b) -> bool { return cells[iC][a] < cells[iC][b]; });
+      std::array<size_t, 8> rotatedNumbering;
+      std::copy(rotationMap[sortedNumbering[0]].begin(), rotationMap[sortedNumbering[0]].end(),
+                rotatedNumbering.begin());
+      size_t n = 0;
+      size_t diagCount = 0;
+      auto checkDiagonal = [this, rotatedNumbering, iC](size_t a1, size_t a2, size_t b1, size_t b2) {
+        return (cells[iC][rotatedNumbering[a1]] < cells[iC][rotatedNumbering[b1]] &&
+                cells[iC][rotatedNumbering[a1]] < cells[iC][rotatedNumbering[b2]]) ||
+               (cells[iC][rotatedNumbering[a2]] < cells[iC][rotatedNumbering[b1]] &&
+                cells[iC][rotatedNumbering[a2]] < cells[iC][rotatedNumbering[b2]]);
+      };
+      if (checkDiagonal(1, 7, 2, 5)) {
+        n += 4;
+        diagCount++;
+      }
+      if (checkDiagonal(3, 7, 2, 6)) {
+        n += 2;
+        diagCount++;
+      }
+      if (checkDiagonal(4, 7, 5, 6)) {
+        n += 1;
+        diagCount++;
+      }
+      if (diagCount == 0) {
+        tetCount += 5;
+      } else {
+        tetCount += 6;
+      }
+      break;
+    }
+    case VolumeCellType::TET:
+      tetCount += 1;
+      break;
+    }
+  }
+  // Mark each edge as real or not (in the original mesh)
+  std::vector<std::array<bool, 6>> realEdges;
+  // Each hex can make up to 6 tets
+  tets.resize(tetCount);
+  realEdges.resize(tetCount);
+  size_t tetIdx = 0;
+  for (size_t iC = 0; iC < nCells(); iC++) {
+    switch (cellType(iC)) {
+    case VolumeCellType::HEX: {
+      std::array<size_t, 8> sortedNumbering;
+      std::iota(sortedNumbering.begin(), sortedNumbering.end(), 0);
+      std::sort(sortedNumbering.begin(), sortedNumbering.end(),
+                [this, iC](size_t a, size_t b) -> bool { return cells[iC][a] < cells[iC][b]; });
+      std::array<size_t, 8> rotatedNumbering;
+      std::copy(rotationMap[sortedNumbering[0]].begin(), rotationMap[sortedNumbering[0]].end(),
+                rotatedNumbering.begin());
+      size_t n = 0;
+      size_t diagCount = 0;
+      // Diagonal exists on the pair of vertices which contain the minimum vertex number
+      auto checkDiagonal = [this, rotatedNumbering, iC](size_t a1, size_t a2, size_t b1, size_t b2) {
+        return (cells[iC][rotatedNumbering[a1]] < cells[iC][rotatedNumbering[b1]] &&
+                cells[iC][rotatedNumbering[a1]] < cells[iC][rotatedNumbering[b2]]) ||
+               (cells[iC][rotatedNumbering[a2]] < cells[iC][rotatedNumbering[b1]] &&
+                cells[iC][rotatedNumbering[a2]] < cells[iC][rotatedNumbering[b2]]);
+      };
+      // Minimum vertex will always have 3 diagonals, check other three faces
+      if (checkDiagonal(1, 7, 2, 5)) {
+        n += 4;
+        diagCount++;
+      }
+      if (checkDiagonal(3, 7, 2, 6)) {
+        n += 2;
+        diagCount++;
+      }
+      if (checkDiagonal(4, 7, 5, 6)) {
+        n += 1;
+        diagCount++;
+      }
+      // Rotate by 120 or 240 degrees depending on diagonal positions
+      if (n == 1 || n == 6) {
+        size_t temp = rotatedNumbering[1];
+        rotatedNumbering[1] = rotatedNumbering[4];
+        rotatedNumbering[4] = rotatedNumbering[3];
+        rotatedNumbering[3] = temp;
+        temp = rotatedNumbering[5];
+        rotatedNumbering[5] = rotatedNumbering[6];
+        rotatedNumbering[6] = rotatedNumbering[2];
+        rotatedNumbering[2] = temp;
+      } else if (n == 2 || n == 5) {
+        size_t temp = rotatedNumbering[1];
+        rotatedNumbering[1] = rotatedNumbering[3];
+        rotatedNumbering[3] = rotatedNumbering[4];
+        rotatedNumbering[4] = temp;
+        temp = rotatedNumbering[5];
+        rotatedNumbering[5] = rotatedNumbering[2];
+        rotatedNumbering[2] = rotatedNumbering[6];
+        rotatedNumbering[6] = temp;
+      }
+
+      // Map final tets according to diagonalMap and the number of diagonals not incident to V_0
+      std::array<std::array<size_t, 4>, 6> tetMap = diagonalMap[diagCount];
+      for (size_t k = 0; k < (diagCount == 0 ? 5 : 6); k++) {
+        for (size_t i = 0; i < 4; i++) {
+          tets[tetIdx][i] = cells[iC][rotatedNumbering[tetMap[k][i]]];
+        }
+        tetIdx++;
+      }
+      break;
+    }
+    case VolumeCellType::TET:
+      for (size_t i = 0; i < 4; i++) {
+        tets[tetIdx][i] = cells[iC][i];
+      }
+      tetIdx++;
+      break;
+    }
+  }
+}
+
+void VolumeMesh::ensureHaveTets() {
+  if (tets.empty()) {
+    computeTets();
+  }
+}
+
+size_t VolumeMesh::nTets() {
+  ensureHaveTets();
+  return tets.size();
+}
+
+void VolumeMesh::addSlicePlaneListener(polyscope::SlicePlane* sp) { volumeSlicePlaneListeners.push_back(sp); }
+
+void VolumeMesh::removeSlicePlaneListener(polyscope::SlicePlane* sp) {
+  for (size_t i = 0; i < volumeSlicePlaneListeners.size(); i++) {
+    if (volumeSlicePlaneListeners[i] == sp) {
+      volumeSlicePlaneListeners.erase(volumeSlicePlaneListeners.begin() + i);
+      break;
+    }
+  }
+}
+
+void VolumeMesh::fillSliceGeometryBuffers(render::ShaderProgram& program) {
+
+  ensureHaveTets();
+
+  std::vector<glm::vec3> point1;
+  std::vector<glm::vec3> point2;
+  std::vector<glm::vec3> point3;
+  std::vector<glm::vec3> point4;
+  size_t tetCount = tets.size();
+  point1.resize(tetCount);
+  point2.resize(tetCount);
+  point3.resize(tetCount);
+  point4.resize(tetCount);
+  for (size_t tetIdx = 0; tetIdx < tets.size(); tetIdx++) {
+    point1[tetIdx] = vertices[tets[tetIdx][0]];
+    point2[tetIdx] = vertices[tets[tetIdx][1]];
+    point3[tetIdx] = vertices[tets[tetIdx][2]];
+    point4[tetIdx] = vertices[tets[tetIdx][3]];
+  }
+
+  program.setAttribute("a_point_1", point1);
+  program.setAttribute("a_point_2", point2);
+  program.setAttribute("a_point_3", point3);
+  program.setAttribute("a_point_4", point4);
+  program.setAttribute("a_slice_1", point1);
+  program.setAttribute("a_slice_2", point2);
+  program.setAttribute("a_slice_3", point3);
+  program.setAttribute("a_slice_4", point4);
+}
 
 void VolumeMesh::computeCounts() {
 
@@ -141,6 +370,15 @@ void VolumeMesh::computeCounts() {
 
 void VolumeMesh::computeGeometryData() {}
 
+VolumeMeshVertexScalarQuantity* VolumeMesh::getLevelSetQuantity() { return activeLevelSetQuantity; }
+
+void VolumeMesh::setLevelSetQuantity(VolumeMeshVertexScalarQuantity* quantity) {
+  if (activeLevelSetQuantity != nullptr && activeLevelSetQuantity != quantity) {
+    activeLevelSetQuantity->isDrawingLevelSet = false;
+  }
+  activeLevelSetQuantity = quantity;
+}
+
 
 void VolumeMesh::draw() {
   if (!isEnabled()) {
@@ -162,12 +400,21 @@ void VolumeMesh::draw() {
     // Set uniforms
     setStructureUniforms(*program);
     setVolumeMeshUniforms(*program);
+    glm::mat4 viewMat = getModelView();
+    glm::mat4 projMat = view::getCameraPerspectiveMatrix();
     program->setUniform("u_baseColor1", getColor());
     program->setUniform("u_baseColor2", getInteriorColor());
 
     program->draw();
   }
 
+  if (activeLevelSetQuantity != nullptr && activeLevelSetQuantity->isEnabled()) {
+    // Draw the quantities
+    activeLevelSetQuantity->draw();
+
+    return;
+  }
+
   // Draw the quantities
   for (auto& x : quantities) {
     x.second->draw();
@@ -184,6 +431,7 @@ void VolumeMesh::drawPick() {
   }
 
   // Set uniforms
+  setVolumeMeshUniforms(*pickProgram);
   setStructureUniforms(*pickProgram);
 
   pickProgram->draw();
@@ -191,7 +439,6 @@ void VolumeMesh::drawPick() {
 
 void VolumeMesh::prepare() {
   program = render::engine->requestShader("MESH", addVolumeMeshRules({"MESH_PROPAGATE_TYPE_AND_BASECOLOR2_SHADE"}));
-
   // Populate draw buffers
   fillGeometryBuffers(*program);
   render::engine->setMaterial(*program, getMaterial());
@@ -200,7 +447,7 @@ void VolumeMesh::prepare() {
 void VolumeMesh::preparePick() {
 
   // Create a new program
-  pickProgram = render::engine->requestShader("MESH", addVolumeMeshRules({"MESH_PROPAGATE_PICK"}, false),
+  pickProgram = render::engine->requestShader("MESH", addVolumeMeshRules({"MESH_PROPAGATE_PICK"}),
                                               render::ShaderReplacementDefaults::Pick);
 
   // == Sort out element counts and index ranges
@@ -222,10 +469,12 @@ void VolumeMesh::preparePick() {
   std::vector<glm::vec3> positions;
   std::vector<glm::vec3> normals;
   std::vector<glm::vec3> bcoord;
+  std::vector<glm::vec3> edgeReal;
   std::vector<std::array<glm::vec3, 3>> vertexColors, edgeColors, halfedgeColors;
   std::vector<glm::vec3> faceColor;
   std::vector<glm::vec3> barycenters;
 
+  bool wantsEdge = (getEdgeWidth() > 0);
   bool wantsBarycenters = wantsCullPosition();
 
   // Reserve space
@@ -239,6 +488,9 @@ void VolumeMesh::preparePick() {
   if (wantsBarycenters) {
     barycenters.resize(3 * nFacesTriangulation());
   }
+  if (wantsEdge) {
+    edgeReal.resize(3 * nFacesTriangulation());
+  }
 
 
   size_t iFront = 0;
@@ -311,6 +563,18 @@ void VolumeMesh::preparePick() {
             barycenters[iData + k] = barycenter;
           }
         }
+        if (wantsEdge) {
+          glm::vec3 edgeRealV{0., 1., 0.};
+          if (j == 0) {
+            edgeRealV.x = 1.;
+          }
+          if (j + 1 == face.size()) {
+            edgeRealV.z = 1.;
+          }
+          for (int k = 0; k < 3; k++) {
+            edgeReal[iData + k] = edgeRealV;
+          }
+        }
       }
 
       iF++;
@@ -326,24 +590,28 @@ void VolumeMesh::preparePick() {
   pickProgram->setAttribute<glm::vec3, 3>("a_edgeColors", edgeColors);
   pickProgram->setAttribute<glm::vec3, 3>("a_halfedgeColors", halfedgeColors);
   pickProgram->setAttribute("a_faceColor", faceColor);
-  if (wantsCullPosition()) {
+  if (wantsBarycenters) {
     pickProgram->setAttribute("a_cullPos", barycenters);
   }
+  if (wantsEdge) {
+    pickProgram->setAttribute("a_edgeIsReal", edgeReal);
+  }
 }
 
-std::vector<std::string> VolumeMesh::addVolumeMeshRules(std::vector<std::string> initRules, bool withSurfaceShade) {
+std::vector<std::string> VolumeMesh::addVolumeMeshRules(std::vector<std::string> initRules, bool withSurfaceShade,
+                                                        bool isSlice) {
 
   initRules = addStructureRules(initRules);
 
   if (withSurfaceShade) {
     if (getEdgeWidth() > 0) {
-      initRules.push_back("MESH_WIREFRAME");
+      initRules.push_back(isSlice ? "SLICE_TETS_MESH_WIREFRAME" : "MESH_WIREFRAME");
     }
   }
 
   initRules.push_back("MESH_BACKFACE_NORMAL_FLIP");
 
-  if (wantsCullPosition()) {
+  if (wantsCullPosition() && !isSlice) {
     initRules.push_back("MESH_PROPAGATE_CULLPOS");
   }
 
@@ -647,17 +915,24 @@ void VolumeMesh::buildCustomOptionsUI() {
 }
 
 
+void VolumeMesh::refreshVolumeMeshListeners() {
+  for (size_t i = 0; i < volumeSlicePlaneListeners.size(); i++) {
+    volumeSlicePlaneListeners[i]->resetVolumeSliceProgram();
+  }
+}
+
 void VolumeMesh::refresh() {
   computeGeometryData();
   program.reset();
   pickProgram.reset();
+  refreshVolumeMeshListeners();
   requestRedraw();
   QuantityStructure<VolumeMesh>::refresh(); // call base class version, which refreshes quantities
 }
 
-void VolumeMesh::geometryChanged() { 
+void VolumeMesh::geometryChanged() {
   // TODO overkill
-  refresh(); 
+  refresh();
 }
 
 VolumeCellType VolumeMesh::cellType(size_t i) const {
@@ -729,7 +1004,7 @@ VolumeMesh* VolumeMesh::setEdgeWidth(double newVal) {
 double VolumeMesh::getEdgeWidth() { return edgeWidth.get(); }
 
 
-// === Quantity adders
+// === Quantity adder}
 
 VolumeMeshVertexColorQuantity* VolumeMesh::addVertexColorQuantityImpl(std::string name,
                                                                       const std::vector<glm::vec3>& colors) {
diff --git a/src/volume_mesh_color_quantity.cpp b/src/volume_mesh_color_quantity.cpp
index 8a0e7d4b..961188db 100644
--- a/src/volume_mesh_color_quantity.cpp
+++ b/src/volume_mesh_color_quantity.cpp
@@ -32,7 +32,61 @@ VolumeMeshVertexColorQuantity::VolumeMeshVertexColorQuantity(std::string name, s
                                                              VolumeMesh& mesh_)
     : VolumeMeshColorQuantity(name, mesh_, "vertex"), values(std::move(values_))
 
-{}
+{
+  parent.refreshVolumeMeshListeners(); // just in case this quantity is being drawn
+}
+
+void VolumeMeshVertexColorQuantity::drawSlice(polyscope::SlicePlane* sp) {
+  if (!isEnabled()) return;
+
+  if (sliceProgram == nullptr) {
+    sliceProgram = createSliceProgram();
+  }
+  parent.setStructureUniforms(*sliceProgram);
+  // Ignore current slice plane
+  sp->setSceneObjectUniforms(*sliceProgram, true);
+  sp->setSliceGeomUniforms(*sliceProgram);
+  parent.setVolumeMeshUniforms(*sliceProgram);
+  sliceProgram->draw();
+}
+
+std::shared_ptr<render::ShaderProgram> VolumeMeshVertexColorQuantity::createSliceProgram() {
+  std::shared_ptr<render::ShaderProgram> p = render::engine->requestShader(
+      "SLICE_TETS", parent.addVolumeMeshRules({"SLICE_TETS_PROPAGATE_VECTOR", "SLICE_TETS_VECTOR_COLOR"}, true, true));
+
+  // Fill color buffers
+  parent.fillSliceGeometryBuffers(*p);
+  fillSliceColorBuffers(*p);
+  render::engine->setMaterial(*p, parent.getMaterial());
+  return p;
+}
+
+void VolumeMeshVertexColorQuantity::fillSliceColorBuffers(render::ShaderProgram& p) {
+  size_t tetCount = parent.nTets();
+  std::vector<glm::vec3> colorval_1;
+  std::vector<glm::vec3> colorval_2;
+  std::vector<glm::vec3> colorval_3;
+  std::vector<glm::vec3> colorval_4;
+
+  colorval_1.resize(tetCount);
+  colorval_2.resize(tetCount);
+  colorval_3.resize(tetCount);
+  colorval_4.resize(tetCount);
+
+  auto vertices = parent.vertices;
+  for (size_t iT = 0; iT < parent.tets.size(); iT++) {
+    colorval_1[iT] = values[parent.tets[iT][0]];
+    colorval_2[iT] = values[parent.tets[iT][1]];
+    colorval_3[iT] = values[parent.tets[iT][2]];
+    colorval_4[iT] = values[parent.tets[iT][3]];
+  }
+
+  // Store data in buffers
+  p.setAttribute("a_value_1", colorval_1);
+  p.setAttribute("a_value_2", colorval_2);
+  p.setAttribute("a_value_3", colorval_3);
+  p.setAttribute("a_value_4", colorval_4);
+}
 
 void VolumeMeshVertexColorQuantity::createProgram() {
   // Create the program to draw this quantity
@@ -97,6 +151,9 @@ std::string VolumeMeshColorQuantity::niceName() { return name + " (" + definedOn
 
 void VolumeMeshColorQuantity::refresh() {
   program.reset();
+  if (sliceProgram) {
+    sliceProgram.reset();
+  }
   Quantity::refresh();
 }
 
diff --git a/src/volume_mesh_scalar_quantity.cpp b/src/volume_mesh_scalar_quantity.cpp
index 6b66b6f7..05f5124c 100644
--- a/src/volume_mesh_scalar_quantity.cpp
+++ b/src/volume_mesh_scalar_quantity.cpp
@@ -47,6 +47,7 @@ void VolumeMeshScalarQuantity::buildCustomUI() {
 
 void VolumeMeshScalarQuantity::refresh() {
   program.reset();
+  sliceProgram.reset();
   Quantity::refresh();
 }
 
@@ -58,10 +59,158 @@ std::string VolumeMeshScalarQuantity::niceName() { return name + " (" + definedO
 
 VolumeMeshVertexScalarQuantity::VolumeMeshVertexScalarQuantity(std::string name, const std::vector<double>& values_,
                                                                VolumeMesh& mesh_, DataType dataType_)
-    : VolumeMeshScalarQuantity(name, mesh_, "vertex", values_, dataType_)
+    : VolumeMeshScalarQuantity(name, mesh_, "vertex", values_, dataType_), levelSetValue(0), isDrawingLevelSet(false),
+      showQuantity(this)
 
 {
   hist.buildHistogram(values, parent.vertexAreas); // rebuild to incorporate weights
+  parent.refreshVolumeMeshListeners();             // just in case this quantity is being drawn
+}
+void VolumeMeshVertexScalarQuantity::fillLevelSetData(render::ShaderProgram& p) {
+  std::vector<glm::vec3> point1;
+  std::vector<glm::vec3> point2;
+  std::vector<glm::vec3> point3;
+  std::vector<glm::vec3> point4;
+  std::vector<glm::vec3> slice1;
+  std::vector<glm::vec3> slice2;
+  std::vector<glm::vec3> slice3;
+  std::vector<glm::vec3> slice4;
+  size_t tetCount = parent.nTets();
+  slice1.resize(tetCount);
+  slice2.resize(tetCount);
+  slice3.resize(tetCount);
+  slice4.resize(tetCount);
+  point1.resize(tetCount);
+  point2.resize(tetCount);
+  point3.resize(tetCount);
+  point4.resize(tetCount);
+  auto vertices = parent.vertices;
+  for (size_t i = 0; i < parent.nTets(); i++) {
+    point1[i] = vertices[parent.tets[i][0]];
+    point2[i] = vertices[parent.tets[i][1]];
+    point3[i] = vertices[parent.tets[i][2]];
+    point4[i] = vertices[parent.tets[i][3]];
+    slice1[i] = glm::vec3(values[parent.tets[i][0]], 0, 0);
+    slice2[i] = glm::vec3(values[parent.tets[i][1]], 0, 0);
+    slice3[i] = glm::vec3(values[parent.tets[i][2]], 0, 0);
+    slice4[i] = glm::vec3(values[parent.tets[i][3]], 0, 0);
+  }
+  p.setAttribute("a_point_1", point1);
+  p.setAttribute("a_point_2", point2);
+  p.setAttribute("a_point_3", point3);
+  p.setAttribute("a_point_4", point4);
+  p.setAttribute("a_slice_1", slice1);
+  p.setAttribute("a_slice_2", slice2);
+  p.setAttribute("a_slice_3", slice3);
+  p.setAttribute("a_slice_4", slice4);
+}
+
+void VolumeMeshVertexScalarQuantity::setLevelSetUniforms(render::ShaderProgram& p) {
+  p.setUniform("u_sliceVector", glm::vec3(1, 0, 0));
+  p.setUniform("u_slicePoint", levelSetValue);
+}
+
+void VolumeMeshVertexScalarQuantity::draw() {
+  if (!isEnabled()) return;
+
+  auto programToDraw = program;
+  if (isDrawingLevelSet) {
+    if (levelSetProgram == nullptr) {
+      levelSetProgram = createSliceProgram();
+      fillLevelSetData(*levelSetProgram);
+    }
+    setLevelSetUniforms(*levelSetProgram);
+    programToDraw = levelSetProgram;
+  } else if (program == nullptr) {
+    createProgram();
+    programToDraw = program;
+  }
+
+  // Set uniforms
+  parent.setStructureUniforms(*programToDraw);
+  parent.setVolumeMeshUniforms(*programToDraw);
+  setScalarUniforms(*programToDraw);
+
+  programToDraw->draw();
+}
+
+void VolumeMeshVertexScalarQuantity::setLevelSetValue(float f) { levelSetValue = f; }
+
+void VolumeMeshVertexScalarQuantity::setEnabledLevelSet(bool v) {
+  if (v) {
+    isDrawingLevelSet = true;
+    setEnabled(true);
+    parent.setLevelSetQuantity(this);
+  } else {
+    isDrawingLevelSet = false;
+    parent.setLevelSetQuantity(nullptr);
+  }
+}
+
+void VolumeMeshVertexScalarQuantity::drawSlice(polyscope::SlicePlane* sp) {
+  if (!isEnabled()) return;
+
+  if (sliceProgram == nullptr) {
+    sliceProgram = createSliceProgram();
+  }
+  parent.setStructureUniforms(*sliceProgram);
+  // Ignore current slice plane
+  sp->setSceneObjectUniforms(*sliceProgram, true);
+  sp->setSliceGeomUniforms(*sliceProgram);
+  parent.setVolumeMeshUniforms(*sliceProgram);
+  setScalarUniforms(*sliceProgram);
+  sliceProgram->draw();
+}
+
+void VolumeMeshVertexScalarQuantity::setLevelSetVisibleQuantity(std::string name) {
+  auto pair = parent.quantities.find(name);
+  if (pair == parent.quantities.end()) {
+    return;
+  }
+  VolumeMeshQuantity* vmq = pair->second.get();
+  VolumeMeshVertexScalarQuantity* q = dynamic_cast<VolumeMeshVertexScalarQuantity*>(vmq);
+  if (q == nullptr) {
+    return;
+  }
+  levelSetProgram = render::engine->requestShader(
+      "SLICE_TETS", parent.addVolumeMeshRules(addScalarRules({"SLICE_TETS_PROPAGATE_VALUE"}), true, true));
+
+  // Fill color buffers
+  parent.fillSliceGeometryBuffers(*levelSetProgram);
+  q->fillSliceColorBuffers(*levelSetProgram);
+  render::engine->setMaterial(*levelSetProgram, parent.getMaterial());
+  fillLevelSetData(*levelSetProgram);
+  setLevelSetUniforms(*levelSetProgram);
+  showQuantity = q;
+}
+
+void VolumeMeshVertexScalarQuantity::buildCustomUI() {
+  VolumeMeshScalarQuantity::buildCustomUI();
+  /* TODO disabled for now
+  if (ImGui::Checkbox("Level Set", &isDrawingLevelSet)) {
+    setEnabledLevelSet(isDrawingLevelSet);
+  }
+  if (isDrawingLevelSet) {
+    ImGui::DragFloat("", &levelSetValue, 0.01f, (float)hist.colormapRange.first, (float)hist.colormapRange.second);
+    if (ImGui::BeginMenu("Show Quantity")) {
+      std::map<std::string, std::unique_ptr<polyscope::VolumeMeshQuantity>>::iterator it;
+      for (it = parent.quantities.begin(); it != parent.quantities.end(); it++) {
+        std::string quantityName = it->first;
+        VolumeMeshQuantity* vmq = it->second.get();
+        VolumeMeshVertexScalarQuantity* vmvsq = dynamic_cast<VolumeMeshVertexScalarQuantity*>(vmq);
+        if (vmvsq != nullptr && ImGui::MenuItem(quantityName.c_str(), NULL, showQuantity == it->second.get())) {
+          setLevelSetVisibleQuantity(quantityName);
+        }
+      }
+      ImGui::EndMenu();
+    }
+  }
+  */
+}
+
+void VolumeMeshVertexScalarQuantity::refresh() {
+  VolumeMeshScalarQuantity::refresh();
+  levelSetProgram.reset();
 }
 
 void VolumeMeshVertexScalarQuantity::createProgram() {
@@ -74,6 +223,17 @@ void VolumeMeshVertexScalarQuantity::createProgram() {
   render::engine->setMaterial(*program, parent.getMaterial());
 }
 
+std::shared_ptr<render::ShaderProgram> VolumeMeshVertexScalarQuantity::createSliceProgram() {
+  std::shared_ptr<render::ShaderProgram> p = render::engine->requestShader(
+      "SLICE_TETS", parent.addVolumeMeshRules(addScalarRules({"SLICE_TETS_PROPAGATE_VALUE"}), true, true));
+
+  // Fill color buffers
+  parent.fillSliceGeometryBuffers(*p);
+  fillSliceColorBuffers(*p);
+  render::engine->setMaterial(*p, parent.getMaterial());
+  return p;
+}
+
 
 void VolumeMeshVertexScalarQuantity::fillColorBuffers(render::ShaderProgram& p) {
   std::vector<double> colorval;
@@ -114,6 +274,34 @@ void VolumeMeshVertexScalarQuantity::fillColorBuffers(render::ShaderProgram& p)
   p.setTextureFromColormap("t_colormap", cMap.get());
 }
 
+void VolumeMeshVertexScalarQuantity::fillSliceColorBuffers(render::ShaderProgram& p) {
+  size_t tetCount = parent.nTets();
+  std::vector<double> colorval_1;
+  std::vector<double> colorval_2;
+  std::vector<double> colorval_3;
+  std::vector<double> colorval_4;
+
+  colorval_1.resize(tetCount);
+  colorval_2.resize(tetCount);
+  colorval_3.resize(tetCount);
+  colorval_4.resize(tetCount);
+
+  auto vertices = parent.vertices;
+  for (size_t iT = 0; iT < parent.tets.size(); iT++) {
+    colorval_1[iT] = values[parent.tets[iT][0]];
+    colorval_2[iT] = values[parent.tets[iT][1]];
+    colorval_3[iT] = values[parent.tets[iT][2]];
+    colorval_4[iT] = values[parent.tets[iT][3]];
+  }
+
+  // Store data in buffers
+  p.setAttribute("a_value_1", colorval_1);
+  p.setAttribute("a_value_2", colorval_2);
+  p.setAttribute("a_value_3", colorval_3);
+  p.setAttribute("a_value_4", colorval_4);
+  p.setTextureFromColormap("t_colormap", cMap.get());
+}
+
 void VolumeMeshVertexScalarQuantity::buildVertexInfoGUI(size_t vInd) {
   ImGui::TextUnformatted(name.c_str());
   ImGui::NextColumn();
@@ -121,6 +309,7 @@ void VolumeMeshVertexScalarQuantity::buildVertexInfoGUI(size_t vInd) {
   ImGui::NextColumn();
 }
 
+
 // ========================================================
 // ==========            Cell Scalar             ==========
 // ========================================================
diff --git a/test/src/basics_test.cpp b/test/src/basics_test.cpp
index ffbf815c..c92f65f3 100644
--- a/test/src/basics_test.cpp
+++ b/test/src/basics_test.cpp
@@ -910,6 +910,28 @@ TEST_F(PolyscopeTest, VolumeMeshCellVector) {
   polyscope::removeAllStructures();
 }
 
+TEST_F(PolyscopeTest, VolumeMeshInspect) {
+  std::vector<glm::vec3> verts;
+  std::vector<std::array<int, 8>> cells;
+  std::tie(verts, cells) = getVolumeMeshData();
+  polyscope::VolumeMesh* psVol = polyscope::registerVolumeMesh("vol", verts, cells);
+
+  // plain old inspecting
+  polyscope::SlicePlane* p = polyscope::addSceneSlicePlane();
+  p->setVolumeMeshToInspect("vol");
+  polyscope::show(3);
+
+  // with a scalar quantity
+  std::vector<float> vals(verts.size(), 0.44);
+  auto q1 = psVol->addVertexScalarQuantity("vals", vals);
+  q1->setEnabled(true);
+  polyscope::show(3);
+  polyscope::removeAllStructures();
+
+  polyscope::removeLastSceneSlicePlane();
+}
+
+
 // ============================================================
 // =============== Ground plane tests
 // ============================================================