diff --git a/Data/Shaders/Lighting.glsl b/Data/Shaders/Lighting.glsl
new file mode 100644
index 0000000..246bd44
--- /dev/null
+++ b/Data/Shaders/Lighting.glsl
@@ -0,0 +1,32 @@
+/**
+ * Simplified Blinn-Phong shading assuming the ambient and diffuse color are equal and the specular color is white.
+ * Additionaly, Fresnel shading is used to enhance the outlines.
+ * Assumes the following global variables are given: cameraPosition.
+ * The camera position is assumed to be the source of a point light.
+*/
+vec3 blinnPhongShadingSurface(
+        in vec3 baseColor, in vec3 fragmentPositionWorld, in vec3 fragmentNormal) {
+    // Blinn-Phong Shading
+    const vec3 lightColor = vec3(1.0);
+    const vec3 ambientColor = baseColor;
+    const vec3 diffuseColor = ambientColor;
+    vec3 phongColor = vec3(0.0);
+
+    const float kA = 0.4;
+    const vec3 Ia = kA * ambientColor;
+    const float kD = 0.6;
+    const float kS = 0.2;
+    const float s = 30;
+
+    const vec3 n = normalize(fragmentNormal);
+    const vec3 v = normalize(cameraPosition - fragmentPositionWorld);
+    const vec3 l = v;//normalize(lightDirection);
+    const vec3 h = normalize(v + l);
+
+    vec3 Id = kD * clamp(abs(dot(n, l)), 0.0, 1.0) * diffuseColor;
+    vec3 Is = kS * pow(clamp(abs(dot(n, h)), 0.0, 1.0), s) * lightColor;
+
+    phongColor = Ia + Id + Is;
+
+    return phongColor;
+}
diff --git a/Data/Shaders/VPT/Clouds.glsl b/Data/Shaders/VPT/Clouds.glsl
index 6830f10..5ee2926 100644
--- a/Data/Shaders/VPT/Clouds.glsl
+++ b/Data/Shaders/VPT/Clouds.glsl
@@ -39,6 +39,10 @@
 pnanovdb_readaccessor_t accessor;
 #endif
 
+#ifdef USE_ISOSURFACES
+vec3 cameraPosition;
+#endif
+
 #include "VptUtils.glsl"
 #include "VptMomentUtils.glsl"
 #include "DeltaTracking.glsl"
@@ -46,6 +50,7 @@ pnanovdb_readaccessor_t accessor;
 #include "ResidualRatioTracking.glsl"
 #include "DecompositionTracking.glsl"
 #include "NextEventTracking.glsl"
+#include "IsosurfaceRendering.glsl"
 
 void pathTraceSample(int i, bool onlyFirstEvent, out ScatterEvent firstEvent){
     uint frame = frameInfo.frameCount + i;
@@ -65,6 +70,10 @@ void pathTraceSample(int i, bool onlyFirstEvent, out ScatterEvent firstEvent){
     vec3 x, w;
     createCameraRay(screenCoord, x, w);
 
+#ifdef USE_ISOSURFACES
+    cameraPosition = x;
+#endif
+
 #ifdef USE_NANOVDB
     accessor = createAccessor();
 #endif
@@ -93,6 +102,8 @@ void pathTraceSample(int i, bool onlyFirstEvent, out ScatterEvent firstEvent){
     vec3 result = nextEventTracking(x, w, firstEvent, onlyFirstEvent);
 #elif defined(USE_NEXT_EVENT_TRACKING_SPECTRAL)
     vec3 result = nextEventTrackingSpectral(x, w, firstEvent, onlyFirstEvent);
+#elif defined(USE_ISOSURFACE_RENDERING)
+    vec3 result = isosurfaceRendering(x, w, firstEvent);
 #endif
 
     if (!onlyFirstEvent) {
diff --git a/Data/Shaders/VPT/IsosurfaceRendering.glsl b/Data/Shaders/VPT/IsosurfaceRendering.glsl
new file mode 100644
index 0000000..ee89bd6
--- /dev/null
+++ b/Data/Shaders/VPT/IsosurfaceRendering.glsl
@@ -0,0 +1,100 @@
+#ifdef USE_ISOSURFACE_RENDERING
+vec3 isosurfaceRendering(vec3 x, vec3 w, out ScatterEvent firstEvent) {
+    firstEvent = ScatterEvent(false, x, 0.0, w, 0.0, 0.0, 0.0);
+
+    vec3 weights = vec3(1, 1, 1);
+    float lastScalarSign, currentScalarSign;
+    bool isFirstPoint = true;
+
+    ivec3 voxelGridSize = textureSize(gridImage, 0);
+    vec3 boxDelta = parameters.boxMax - parameters.boxMin;
+    vec3 voxelSize3d = boxDelta / voxelGridSize;
+    float t = max(voxelSize3d.x, max(voxelSize3d.y, voxelSize3d.z)) * parameters.isoStepWidth;
+
+    bool foundHit = false;
+    int i = 0;
+    float tMin, tMax;
+    if (rayBoxIntersect(parameters.boxMin, parameters.boxMax, x, w, tMin, tMax)) {
+        x += w * tMin;
+        float d = tMax - tMin;
+        while (t <= d) {
+            vec3 xNew = x + w * t;
+            float scalarValue = sampleCloudDirect(xNew);
+
+            currentScalarSign = sign(scalarValue - parameters.isoValue);
+            if (isFirstPoint) {
+                isFirstPoint = false;
+                lastScalarSign = currentScalarSign;
+            }
+
+            if (lastScalarSign != currentScalarSign) {
+                if (!firstEvent.hasValue) {
+                    firstEvent.x = x;
+                    firstEvent.pdf_x = 0;
+                    firstEvent.w = vec3(0.);
+                    firstEvent.pdf_w = 0;
+                    firstEvent.hasValue = true;
+                    firstEvent.density = parameters.extinction.x;
+                    firstEvent.depth = tMax - d + t;
+                }
+                refineIsoSurfaceHit(xNew, x, currentScalarSign);
+                x = xNew;
+                foundHit = true;
+                break;
+            }
+
+            x = xNew;
+            d -= t;
+        }
+    }
+
+    if (foundHit) {
+        vec3 surfaceNormal;
+        vec3 color = getIsoSurfaceHitDirect(x, w, surfaceNormal);
+        weights *= color;
+        x += surfaceNormal * 1e-4;
+
+        vec3 surfaceTangent;
+        vec3 surfaceBitangent;
+        ComputeDefaultBasis(surfaceNormal, surfaceTangent, surfaceBitangent);
+        mat3 frame = mat3(surfaceTangent, surfaceBitangent, surfaceNormal);
+
+        const int numAoSamples = 4;
+        const float MAX_DIST = 0.05;
+        float weight = 1.0;
+
+        for (int i = 0; i < numAoSamples; i++) {
+            w = frame * sampleHemisphere(vec2(random(), random()));
+
+            if (rayBoxIntersect(parameters.boxMin, parameters.boxMax, x, w, tMin, tMax)) {
+                x += w * tMin;
+                float d = tMax - tMin;
+                d = max(d, MAX_DIST);
+                while (t <= d) {
+                    vec3 xNew = x + w * t;
+                    float scalarValue = sampleCloudDirect(xNew);
+
+                    currentScalarSign = sign(scalarValue - parameters.isoValue);
+                    if (isFirstPoint) {
+                        isFirstPoint = false;
+                        lastScalarSign = currentScalarSign;
+                    }
+
+                    if (lastScalarSign != currentScalarSign) {
+                        weight -= 1.0 / float(numAoSamples);
+                        break;
+                    }
+
+                    x = xNew;
+                    d -= t;
+                }
+            }
+        }
+
+        return weights * weight;
+        //return surfaceNormal;
+    }
+
+    return sampleSkybox(w) + sampleLight(w);
+}
+#endif
diff --git a/Data/Shaders/VPT/VptHeader.glsl b/Data/Shaders/VPT/VptHeader.glsl
index 8dd0aeb..7f4edc2 100644
--- a/Data/Shaders/VPT/VptHeader.glsl
+++ b/Data/Shaders/VPT/VptHeader.glsl
@@ -80,6 +80,7 @@ layout (binding = 3) uniform Parameters {
     // Isosurfaces.
     vec3 isoSurfaceColor;
     float isoValue;
+    float isoStepWidth;
 } parameters;
 
 layout (binding = 4) uniform FrameInfo {
diff --git a/Data/Shaders/VPT/VptUtils.glsl b/Data/Shaders/VPT/VptUtils.glsl
index d1d438b..8a0e3e9 100644
--- a/Data/Shaders/VPT/VptUtils.glsl
+++ b/Data/Shaders/VPT/VptUtils.glsl
@@ -590,6 +590,7 @@ float avgComponent(vec3 v) {
 #ifdef USE_ISOSURFACES
 #include "RayTracingUtilities.glsl"
 
+//#define DIFFERENCES_NEIGHBOR
 #define M_PI 3.14159265358979323846
 vec3 computeGradient(vec3 texCoords) {
 #ifdef DIFFERENCES_NEIGHBOR
@@ -606,9 +607,9 @@ vec3 computeGradient(vec3 texCoords) {
             (textureOffset(gridImage, texCoords, ivec3(0, 0, -1)).r
             - textureOffset(gridImage, texCoords, ivec3(0, 0, 1)).r) * 0.5 / dz;
 #else
-    const float dx = 1e-6;
-    const float dy = 1e-6;
-    const float dz = 1e-6;
+    const float dx = 1e-3;
+    const float dy = 1e-3;
+    const float dz = 1e-3;
     float gradX =
             (texture(gridImage, texCoords - vec3(dx, 0.0, 0.0)).r
             - texture(gridImage, texCoords + vec3(dx, 0.0, 0.0)).r) * 0.5 / dx;
@@ -722,3 +723,20 @@ vec3 getIsoSurfaceHit(vec3 currentPoint, inout vec3 w) {
     return color;
 }
 #endif
+
+
+#ifdef USE_ISOSURFACE_RENDERING
+#include "Lighting.glsl"
+
+// Direct illumination
+vec3 getIsoSurfaceHitDirect(vec3 currentPoint, vec3 w, inout vec3 surfaceNormal) {
+    vec3 texCoords = (currentPoint - parameters.boxMin) / (parameters.boxMax - parameters.boxMin);
+    texCoords = texCoords * (parameters.gridMax - parameters.gridMin) + parameters.gridMin;
+    surfaceNormal = computeGradient(texCoords);
+    if (dot(cameraPosition - currentPoint, surfaceNormal) < 0.0) {
+        surfaceNormal = -surfaceNormal;
+    }
+    vec3 color = blinnPhongShadingSurface(parameters.isoSurfaceColor, currentPoint, surfaceNormal);
+    return color;
+}
+#endif
diff --git a/src/PathTracer/VolumetricPathTracingPass.cpp b/src/PathTracer/VolumetricPathTracingPass.cpp
index 2db79f8..ac61d8c 100644
--- a/src/PathTracer/VolumetricPathTracingPass.cpp
+++ b/src/PathTracer/VolumetricPathTracingPass.cpp
@@ -570,6 +570,9 @@ void VolumetricPathTracingPass::flipYZ(bool flip) {
 }
 
 void VolumetricPathTracingPass::setUseIsosurfaces(bool _useIsosurfaces) {
+    if (vptMode == VptMode::ISOSURFACE_RENDERING && !_useIsosurfaces) {
+        _useIsosurfaces = true;
+    }
     if (useIsosurfaces != _useIsosurfaces) {
         useIsosurfaces = _useIsosurfaces;
         if (gridInterpolationType != GridInterpolationType::TRILINEAR) {
@@ -630,6 +633,9 @@ void VolumetricPathTracingPass::updateVptMode() {
     if (accumulationTimer && !reachedTarget) {
         createNewAccumulationTimer = true;
     }
+    if (vptMode == VptMode::ISOSURFACE_RENDERING) {
+        useIsosurfaces = true;
+    }
     if (vptMode == VptMode::RESIDUAL_RATIO_TRACKING && cloudData && !useSparseGrid) {
         superVoxelGridDecompositionTracking = {};
         superVoxelGridResidualRatioTracking = std::make_shared<SuperVoxelGridResidualRatioTracking>(
@@ -843,6 +849,8 @@ void VolumetricPathTracingPass::loadShader() {
         customPreprocessorDefines.insert({ "USE_NEXT_EVENT_TRACKING", "" });
     } else if (vptMode == VptMode::NEXT_EVENT_TRACKING_SPECTRAL) {
         customPreprocessorDefines.insert({ "USE_NEXT_EVENT_TRACKING_SPECTRAL", "" });
+    } else if (vptMode == VptMode::ISOSURFACE_RENDERING) {
+        customPreprocessorDefines.insert({ "USE_ISOSURFACE_RENDERING", "" });
     }
     if (gridInterpolationType == GridInterpolationType::NEAREST) {
         customPreprocessorDefines.insert({ "GRID_INTERPOLATION_NEAREST", "" });
@@ -1129,6 +1137,7 @@ void VolumetricPathTracingPass::_render() {
         }
         uniformData.isoSurfaceColor = isoSurfaceColor;
         uniformData.isoValue = isoValue;
+        uniformData.isoStepWidth = isoStepWidth;
         uniformBuffer->updateData(
                 sizeof(UniformData), &uniformData, renderer->getVkCommandBuffer());
 
@@ -1543,7 +1552,8 @@ bool VolumetricPathTracingPass::renderGuiPropertyEditorNodes(sgl::PropertyEditor
             setShaderDirty();
         }
 
-        if (propertyEditor.addCheckbox("Use Isosurfaces", &useIsosurfaces)) {
+        if (vptMode != VptMode::ISOSURFACE_RENDERING && propertyEditor.addCheckbox(
+                "Use Isosurfaces", &useIsosurfaces)) {
             if (gridInterpolationType != GridInterpolationType::TRILINEAR) {
                 gridInterpolationType = GridInterpolationType::TRILINEAR;
                 updateGridSampler();
@@ -1561,6 +1571,11 @@ bool VolumetricPathTracingPass::renderGuiPropertyEditorNodes(sgl::PropertyEditor
             reRender = true;
             frameInfo.frameCount = 0;
         }
+        if (vptMode == VptMode::ISOSURFACE_RENDERING && propertyEditor.addSliderFloat(
+                "Step Width", &isoStepWidth, 0.1f, 1.0f)) {
+            reRender = true;
+            frameInfo.frameCount = 0;
+        }
         if (useIsosurfaces && propertyEditor.addCombo(
                 "Isosurface Field", (int*)&isosurfaceType,
                 ISOSURFACE_TYPE_NAMES, IM_ARRAYSIZE(ISOSURFACE_TYPE_NAMES))) {
diff --git a/src/PathTracer/VolumetricPathTracingPass.hpp b/src/PathTracer/VolumetricPathTracingPass.hpp
index ab2f7de..0a48a8a 100644
--- a/src/PathTracer/VolumetricPathTracingPass.hpp
+++ b/src/PathTracer/VolumetricPathTracingPass.hpp
@@ -111,8 +111,8 @@ enum class VptMode {
 };
 const char* const VPT_MODE_NAMES[] = {
         "Delta Tracking", "Delta Tracking (Spectral)", "Ratio Tracking",
-        "Decomposition Tracking", "Residual Ratio Tracking", "Next Event Tracking",
-        "Next Event Tracking (Spectral)"
+        "Decomposition Tracking", "Residual Ratio Tracking", "Next Event Tracking", "Next Event Tracking (Spectral)",
+        "Isosurfaces"
 };
 
 enum class GridInterpolationType {
@@ -310,6 +310,7 @@ class VolumetricPathTracingPass : public sgl::vk::ComputePass {
     // Isosurface data.
     bool useIsosurfaces = false;
     float isoValue = 0.5f;
+    float isoStepWidth = 0.25f;
     glm::vec3 isoSurfaceColor = glm::vec3(0.8f, 0.8f, 0.8f);
     IsosurfaceType isosurfaceType = IsosurfaceType::DENSITY;
     SurfaceBrdf surfaceBrdf = SurfaceBrdf::LAMBERTIAN;
@@ -350,6 +351,7 @@ class VolumetricPathTracingPass : public sgl::vk::ComputePass {
         // Isosurfaces.
         glm::vec3 isoSurfaceColor;
         float isoValue = 0.5f;
+        float isoStepWidth = 0.25f;
     };
     UniformData uniformData{};
     sgl::vk::BufferPtr uniformBuffer;