basic marching cubes working

.vscode/launch.json

@@ -11,7 +11,20 @@
             "program": "${workspaceRoot}/build/FluidSimulation",
             "args": [],
             "cwd": "${workspaceRoot}/build",
-            "MIMode": "lldb"
+            "MIMode": "gdb",
+            "setupCommands": [
+                {
+                    "text": "-enable-pretty-printing",
+                    "description": "enable pretty printing",
+                    "ignoreFailures": false
+                },
+                {
+                    "text": "source ${workspaceRoot}/Seele/gdb/Seele.py",
+                    "description": "load Seele python script",
+                    "ignoreFailures": false
+                }
+            ],
+            
         }
     ]
 }

res/shaders/FluidGridData.slang

@@ -1,4 +1,4 @@
-static const uint3 gridSize = uint3(64, 64, 64);
+static const uint3 gridSize = uint3(32, 32, 32);
 struct FluidGridData<T>
 {
     RWStructuredBuffer<T> dataGrid;

res/shaders/Render.slang

@@ -1,32 +1,32 @@
+import Common;
 import FluidGridData;
 
 struct Params
 {
-    FluidGridData<float> density;
+    StructuredBuffer<float> vertexBuffer;
+    StructuredBuffer<uint> indexBuffer;
 };
 ParameterBlock<Params> params;
 
 struct VertexOut
 {
     float4 position : SV_Position;
-    float2 uv : UV;
 };
 
 [shader("vertex")]
 VertexOut vertexMain(uint vertexId : SV_VertexID)
 {
     VertexOut output;
-    output.uv = float2((vertexId << 1) & 2, vertexId & 2);
+    uint index = params.indexBuffer[vertexId];
-    output.position = float4(output.uv * 2.0f - 1.0f, 0, 1);
+    float3 vertex = float3(params.vertexBuffer[index * 3 + 0],
+                           params.vertexBuffer[index * 3 + 1],
+                           params.vertexBuffer[index * 3 + 2]) * 5;
+    output.position = mul(pViewParams.viewProjectionMatrix, float4(vertex, 1));
     return output;
 }
 
 [shader("pixel")]
 float4 fragmentMain(VertexOut input) : SV_Target
 {
-    FluidGridData<float> density = params.density;
+    return float4(0, 0, 1, 1);
-    int x = int(input.uv.x * (gridSize.x - 2) + 1);
-    int y = int(input.uv.y * (gridSize.y - 2) + 1);
-    float d = density[x, y, 1];
-    return float4(d, d, d, 1);
 }

res/shaders/SignedDistance.slang

@@ -0,0 +1,65 @@
+import FluidGridData;
+
+struct Params
+{
+    float dtau;
+    // read-write: phi being reinitialized
+    FluidGridData<float> phi;
+    // read-only: snapshot of phi before reinitialization (for sign)
+    FluidGridData<float> phi0;
+};
+ParameterBlock<Params> params;
+
+// Smooth sign function to avoid discontinuity at zero
+float smoothSign(float x, float dx)
+{
+    return x / sqrt(x * x + dx * dx);
+}
+
+[shader("compute")]
+[numthreads(32, 8, 1)]
+void reinitialize(uint3 dispatchThreadID : SV_DispatchThreadID)
+{
+    uint x = dispatchThreadID.x + 1;
+    uint y = dispatchThreadID.y + 1;
+    uint z = dispatchThreadID.z + 1;
+    FluidGridData<float> phi = params.phi;
+    FluidGridData<float> phi0 = params.phi0;
+    if(x >= gridSize.x - 1 || y >= gridSize.y - 1 || z >= gridSize.z - 1) return;
+
+    float dtau = params.dtau;
+    float dx = 1.0f;
+
+    // Current value and original sign
+    float p = phi[x, y, z];
+    float s = smoothSign(phi0[x, y, z], dx);
+
+    // Upwind finite differences for |grad phi|
+    // Forward and backward differences
+    float dxp = phi[x + 1, y, z] - p;
+    float dxm = p - phi[x - 1, y, z];
+    float dyp = phi[x, y + 1, z] - p;
+    float dym = p - phi[x, y - 1, z];
+    float dzp = phi[x, y, z + 1] - p;
+    float dzm = p - phi[x, y, z - 1];
+
+    // Godunov upwind scheme
+    float gradPhi;
+    if (s > 0)
+    {
+        float ax = max(max(dxm, 0.0f), -min(dxp, 0.0f));
+        float ay = max(max(dym, 0.0f), -min(dyp, 0.0f));
+        float az = max(max(dzm, 0.0f), -min(dzp, 0.0f));
+        gradPhi = sqrt(ax * ax + ay * ay + az * az) / dx;
+    }
+    else
+    {
+        float ax = max(-min(dxm, 0.0f), max(dxp, 0.0f));
+        float ay = max(-min(dym, 0.0f), max(dyp, 0.0f));
+        float az = max(-min(dzm, 0.0f), max(dzp, 0.0f));
+        gradPhi = sqrt(ax * ax + ay * ay + az * az) / dx;
+    }
+
+    // Update: push |grad phi| toward 1
+    phi[x, y, z] = p - dtau * s * (gradPhi - 1.0f);
+}

src/Common.h

@@ -1,8 +1,7 @@
 #pragma once
 #include "MinimalEngine.h"
 
-constexpr static uint32 gridSize = 64;
+constexpr static uint32 gridSize = 32;
 constexpr static float viscosity = 0.001f;
 constexpr static float diffusion = 0.001f;
-constexpr static float iso = 0.3f;
 constexpr static float dt = 0.1f;

src/FluidRenderPass.cpp

@@ -13,10 +13,15 @@ FluidRenderPass::FluidRenderPass(Seele::Gfx::PGraphics graphics)
   pipelineLayout = graphics->createPipelineLayout("FluidPipelineLayout");
   descriptorLayout = graphics->createDescriptorLayout("params");
   descriptorLayout->addDescriptorBinding(Gfx::DescriptorBinding{
-      .name = "density",
+      .name = "vertexBuffer",
+      .descriptorType = Gfx::SE_DESCRIPTOR_TYPE_STORAGE_BUFFER,
+  });
+  descriptorLayout->addDescriptorBinding(Gfx::DescriptorBinding{
+      .name = "indexBuffer",
       .descriptorType = Gfx::SE_DESCRIPTOR_TYPE_STORAGE_BUFFER,
   });
   descriptorLayout->create();
+  pipelineLayout->addDescriptorLayout(viewParamsLayout);
   pipelineLayout->addDescriptorLayout(descriptorLayout);
 
   graphics->beginShaderCompilation(ShaderCompilationInfo{
@@ -36,18 +41,20 @@ FluidRenderPass::~FluidRenderPass() {}
 void FluidRenderPass::beginFrame(const Seele::Component::Camera &camera,
                                  const Seele::Component::Transform &transform) {
   updateViewParameters(camera, transform);
-  descriptorSet = descriptorLayout->allocateDescriptorSet();
-  descriptorSet->updateBuffer("density", 0, densityBuffer);
-  descriptorSet->writeChanges();
 }
 
 void FluidRenderPass::render() {
+  viewParamsSet = createViewParamsSet();
+  descriptorSet = descriptorLayout->allocateDescriptorSet();
+  descriptorSet->updateBuffer("vertexBuffer", 0, vertexBuffer);
+  descriptorSet->updateBuffer("indexBuffer", 0, indexBuffer);
+  descriptorSet->writeChanges();
   graphics->beginRenderPass(renderPass);
   Gfx::ORenderCommand renderCommand = graphics->createRenderCommand();
   renderCommand->setViewport(viewport);
   renderCommand->bindPipeline(pipeline);
-  renderCommand->bindDescriptor(descriptorSet);
+  renderCommand->bindDescriptor({viewParamsSet, descriptorSet});
-  renderCommand->draw(3, 1, 0, 0);
+  renderCommand->drawIndirect(surfaceCounts, 0, 1, 0);
   graphics->executeCommands(std::move(renderCommand));
   graphics->endRenderPass();
 }
@@ -65,7 +72,9 @@ void FluidRenderPass::publishOutputs() {
 }
 
 void FluidRenderPass::createRenderPass() {
-  densityBuffer = resources->requestBuffer("DENSITY");
+  vertexBuffer = resources->requestBuffer("SURFACE_VERTICES");
+  indexBuffer = resources->requestBuffer("SURFACE_INDICES");
+  surfaceCounts = resources->requestBuffer("SURFACE_COUNTS");
   colorAttachment = Gfx::RenderTargetAttachment(
       viewport, Gfx::SE_IMAGE_LAYOUT_UNDEFINED,
       Gfx::SE_IMAGE_LAYOUT_PRESENT_SRC_KHR, Gfx::SE_ATTACHMENT_LOAD_OP_CLEAR,
@@ -126,7 +135,7 @@ void FluidRenderPass::createRenderPass() {
           },
       .rasterizationState =
           {
-              .cullMode = Gfx::SE_CULL_MODE_BACK_BIT,
+              .cullMode = Gfx::SE_CULL_MODE_FRONT_BIT,
           },
       .colorBlend =
           {

src/FluidRenderPass.h

@@ -22,6 +22,8 @@ private:
   Seele::Gfx::RenderTargetAttachment depthAttachment;
   Seele::Gfx::OTexture2D depthTexture;
 
+  Seele::Gfx::ODescriptorSet viewParamsSet;
+
   Seele::Gfx::OPipelineLayout pipelineLayout;
   Seele::Gfx::ODescriptorLayout descriptorLayout;
   Seele::Gfx::ODescriptorSet descriptorSet;
@@ -29,6 +31,8 @@ private:
   Seele::Gfx::OFragmentShader fragmentShader;
   Seele::Gfx::PGraphicsPipeline pipeline;
 
-  Seele::Gfx::PShaderBuffer densityBuffer;
+Seele::Gfx::PShaderBuffer vertexBuffer;
+Seele::Gfx::PShaderBuffer indexBuffer;
+Seele::Gfx::PShaderBuffer surfaceCounts;
 };
 DEFINE_REF(FluidRenderPass)

src/FluidSimulationView.cpp

@@ -1,16 +1,22 @@
 #include "FluidSimulationView.h"
 #include "FluidRenderPass.h"
+#include "Graphics/Enums.h"
 #include "SimulationComputePass.h"
+#include "SurfaceExtractPass.h"
+#include <iostream>
 
 using namespace Seele;
 
-FluidSimulationView::FluidSimulationView(Gfx::PGraphics graphics, Seele::PWindow window, const Seele::ViewportCreateInfo &createInfo)
+FluidSimulationView::FluidSimulationView(
+    Gfx::PGraphics graphics, Seele::PWindow window,
+    const Seele::ViewportCreateInfo &createInfo)
     : View(graphics, window, createInfo, "Fluid Simulation View") {
   camera = Component::Camera{};
   transform = Component::Transform{
       .transform = Math::Transform({0.0f, 0.0f, -5.0f}),
   };
   renderGraph.addPass(new SimulationComputePass(graphics));
+  renderGraph.addPass(new SurfaceExtractPass(graphics));
   renderGraph.addPass(new FluidRenderPass(graphics));
   renderGraph.setViewport(viewport);
   renderGraph.createRenderPass();
@@ -29,8 +35,65 @@ void FluidSimulationView::prepareRender() {}
 void FluidSimulationView::render() { renderGraph.render(camera, transform); }
 
 void FluidSimulationView::applyArea(Seele::URect area) {}
-void FluidSimulationView::keyCallback(Seele::KeyCode code, Seele::InputAction action, Seele::KeyModifier modifier) {}
+
-void FluidSimulationView::mouseMoveCallback(double xPos, double yPos) {}
+void FluidSimulationView::keyCallback(Seele::KeyCode code,
-void FluidSimulationView::mouseButtonCallback(Seele::MouseButton button, Seele::InputAction action, Seele::KeyModifier modifier) {}
+                                      Seele::InputAction action,
+                                      Seele::KeyModifierFlags modifier) {
+  float cameraMove = static_cast<float>(Gfx::getCurrentFrameDelta()) * 4;
+  if ((modifier & (uint32)KeyModifier::MOD_SHIFT)) {
+    cameraMove *= 10;
+  }
+  Vector forward = transform.getForward();
+  Vector side = transform.getRight();
+  Vector moveVector = Vector();
+
+  if (code == KeyCode::KEY_J && action != InputAction::RELEASE) {
+    std::cout << transform.getPosition() << std::endl;
+  }
+
+  if (code == KeyCode::KEY_W && action != InputAction::RELEASE) {
+    moveVector += forward * cameraMove;
+  }
+  if (code == KeyCode::KEY_S && action != InputAction::RELEASE) {
+    moveVector += forward * -cameraMove;
+  }
+  if (code == KeyCode::KEY_A && action != InputAction::RELEASE) {
+    moveVector += side * cameraMove;
+  }
+  if (code == KeyCode::KEY_D && action != InputAction::RELEASE) {
+    moveVector += side * -cameraMove;
+  }
+  if (code == KeyCode::KEY_E && action != InputAction::RELEASE) {
+    moveVector += Vector(0, cameraMove, 0);
+  }
+  if (code == KeyCode::KEY_Q && action != InputAction::RELEASE) {
+    moveVector += Vector(0, -cameraMove, 0);
+  }
+  transform.translate(moveVector);
+}
+
+void FluidSimulationView::mouseMoveCallback(double xPos, double yPos) {
+  deltaX = xPos - mouseX;
+  deltaY = yPos - mouseY;
+  mouseX = xPos;
+  mouseY = yPos;
+
+  if (buttonDown) {
+    transform.setRotation(
+        glm::rotate(transform.getRotation(), deltaX / 500, Vector(0, 1, 0)));
+    transform.setRotation(glm::rotate(transform.getRotation(), -deltaY / 500,
+                                      transform.getRight()));
+  }
+}
+
+void FluidSimulationView::mouseButtonCallback(
+    Seele::MouseButton button, Seele::InputAction action,
+    Seele::KeyModifierFlags modifier) {
+  if(button == MouseButton::MOUSE_BUTTON_1) {
+    buttonDown = action != InputAction::RELEASE;
+  }
+}
+
 void FluidSimulationView::scrollCallback(double xOffset, double yOffset) {}
+
 void FluidSimulationView::fileCallback(int count, const char **paths) {}

src/FluidSimulationView.h

@@ -18,9 +18,9 @@ public:
 
 protected:
   virtual void applyArea(Seele::URect area) override;
-  virtual void keyCallback(Seele::KeyCode code, Seele::InputAction action, Seele::KeyModifier modifier) override;
+  virtual void keyCallback(Seele::KeyCode code, Seele::InputAction action, Seele::KeyModifierFlags modifier) override;
   virtual void mouseMoveCallback(double xPos, double yPos) override;
-  virtual void mouseButtonCallback(Seele::MouseButton button, Seele::InputAction action, Seele::KeyModifier modifier) override;
+  virtual void mouseButtonCallback(Seele::MouseButton button, Seele::InputAction action, Seele::KeyModifierFlags modifier) override;
   virtual void scrollCallback(double xOffset, double yOffset) override;
   virtual void fileCallback(int count, const char **paths) override;
 
@@ -28,4 +28,9 @@ private:
   Seele::Component::Camera camera;
   Seele::Component::Transform transform;
   Seele::RenderGraph renderGraph;
+  float deltaX;
+  float deltaY;
+  float mouseX;
+  float mouseY;
+  bool buttonDown = false;
 };

src/SimulationComputePass.cpp

@@ -1,4 +1,5 @@
 #include "SimulationComputePass.h"
+#include "Common.h"
 #include "Graphics/Buffer.h"
 #include "Graphics/Command.h"
 #include "Graphics/Descriptor.h"
@@ -6,10 +7,14 @@
 #include "Graphics/Graphics.h"
 #include "Graphics/Initializer.h"
 #include "Graphics/Shader.h"
-#include "Common.h"
 
 using namespace Seele;
 
+
+constexpr static int reinitInterval = 5;
+constexpr static int reinitIterations = 5;
+constexpr static float reinitDtau = 0.5f;
+
 #define SWAP(a, b)                                                             \
   {                                                                            \
     auto temp = std::move(a);                                                  \
@@ -231,6 +236,44 @@ SimulationComputePass::SimulationComputePass(Seele::Gfx::PGraphics graphics)
             .pipelineLayout = advect.pipelineLayout,
         });
   }
+  {
+    reinitialize.pipelineLayout =
+        graphics->createPipelineLayout("ReinitializePipelineLayout");
+    reinitialize.descriptorLayout =
+        graphics->createDescriptorLayout("params");
+    reinitialize.descriptorLayout->addDescriptorBinding(
+        Gfx::DescriptorBinding{
+            .name = "dtau",
+            .uniformLength = sizeof(float),
+        });
+    reinitialize.descriptorLayout->addDescriptorBinding(
+        Gfx::DescriptorBinding{
+            .name = "phi",
+            .descriptorType = Gfx::SE_DESCRIPTOR_TYPE_STORAGE_BUFFER,
+        });
+    reinitialize.descriptorLayout->addDescriptorBinding(
+        Gfx::DescriptorBinding{
+            .name = "phi0",
+            .descriptorType = Gfx::SE_DESCRIPTOR_TYPE_STORAGE_BUFFER,
+        });
+    reinitialize.descriptorLayout->create();
+    reinitialize.pipelineLayout->addDescriptorLayout(
+        reinitialize.descriptorLayout);
+
+    graphics->beginShaderCompilation(ShaderCompilationInfo{
+        .name = "SignedDistance",
+        .modules = {"SignedDistance"},
+        .entryPoints = {{"reinitialize", "SignedDistance"}},
+        .rootSignature = reinitialize.pipelineLayout,
+    });
+    reinitialize.pipelineLayout->create();
+    reinitialize.shader = graphics->createComputeShader({0});
+    reinitialize.pipeline =
+        graphics->createComputePipeline(Gfx::ComputePipelineCreateInfo{
+            .computeShader = reinitialize.shader,
+            .pipelineLayout = reinitialize.pipelineLayout,
+        });
+  }
 
   uint64_t bufferSize = gridSize * gridSize * gridSize * sizeof(float);
   Array<float> initialDensity(bufferSize / sizeof(float));
@@ -238,10 +281,29 @@ SimulationComputePass::SimulationComputePass(Seele::Gfx::PGraphics graphics)
   Array<float> initialVelocityY(bufferSize / sizeof(float));
   Array<float> initialVelocityZ(bufferSize / sizeof(float));
   Array<float> initialPhi(bufferSize / sizeof(float));
-  for (uint32 x = 1; x < gridSize - 2; ++x) {
+
-    for (uint32 y = 1; y < gridSize - 2; ++y) {
+  // Initialize level set to positive (outside) everywhere
-        initialDensity[x + gridSize * y + gridSize * gridSize] = 1.0f;
+  float halfGrid = static_cast<float>(gridSize) * 0.5f;
-        initialPhi[x + gridSize * y + gridSize * gridSize] = iso - 1.0f;
+  float radius = static_cast<float>(gridSize) * 0.2f;
+  // Place sphere center in the lower portion of the grid
+  float cx = halfGrid;
+  float cy = halfGrid;
+  float cz = static_cast<float>(gridSize) * 0.3f;
+  for (uint32 z = 1; z < gridSize - 1; ++z) {
+    for (uint32 y = 1; y < gridSize - 1; ++y) {
+      for (uint32 x = 1; x < gridSize - 1; ++x) {
+        uint32 idx = x + gridSize * y + gridSize * gridSize * z;
+        float dx = static_cast<float>(x) - cx;
+        float dy = static_cast<float>(y) - cy;
+        float dz = static_cast<float>(z) - cz;
+        float dist = std::sqrt(dx * dx + dy * dy + dz * dz);
+        initialPhi[idx] = dist - radius;
+        if (dist < radius) {
+          initialDensity[idx] = 1.0f;
+          // Give an initial upward velocity to kick the simulation
+          initialVelocityZ[idx] = 5.0f;
+        }
+      }
     }
   }
 
@@ -256,6 +318,7 @@ SimulationComputePass::SimulationComputePass(Seele::Gfx::PGraphics graphics)
       .sourceData =
           {
               .size = bufferSize,
+              .data = (uint8 *)initialPhi.data(),
           },
       .name = "Phi0Buffer",
   });
@@ -359,6 +422,7 @@ void SimulationComputePass::beginFrame(
   computeDivergence.descriptorLayout->reset();
   project.descriptorLayout->reset();
   advect.descriptorLayout->reset();
+  reinitialize.descriptorLayout->reset();
 }
 
 void SimulationComputePass::render() {
@@ -453,19 +517,32 @@ void SimulationComputePass::render() {
   setBoundsPass(density0Buffer, 0);
   graphics->endDebugRegion();
 
-  graphics->beginDebugRegion("AdjectLevelSet");
+  graphics->beginDebugRegion("AdvectLevelSet");
   advectPass(phiBuffer, phi0Buffer, velocityX0Buffer, velocityY0Buffer,
              velocityZ0Buffer, dt);
   setBoundsPass(phiBuffer, 0);
-  SWAP(phiBuffer, phi0Buffer);
+
-  // TODO: reinitialize level set every few frames
+  // Reinitialize level set every N frames to restore |grad phi| = 1
+  if (frameCounter % reinitInterval == 0) {
+    reinitializeLevelSetPass(phiBuffer, phi0Buffer);
+  }
+  ++frameCounter;
+
+  phiBuffer->pipelineBarrier(Gfx::SE_ACCESS_SHADER_WRITE_BIT,
+                             Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT,
+                             Gfx::SE_ACCESS_TRANSFER_READ_BIT,
+                             Gfx::SE_PIPELINE_STAGE_TRANSFER_BIT);
+  graphics->copyBuffer(phiBuffer, phi0Buffer);
+  phi0Buffer->pipelineBarrier(
+      Gfx::SE_ACCESS_TRANSFER_WRITE_BIT, Gfx::SE_PIPELINE_STAGE_TRANSFER_BIT,
+      Gfx::SE_ACCESS_TRANSFER_READ_BIT, Gfx::SE_PIPELINE_STAGE_TRANSFER_BIT);
   graphics->endDebugRegion();
 }
 
 void SimulationComputePass::endFrame() {}
 
 void SimulationComputePass::publishOutputs() {
-  resources->registerBufferOutput("DENSITY", densityBuffer);
+  resources->registerBufferOutput("PHI", phi0Buffer);
 }
 
 void SimulationComputePass::createRenderPass() {}
@@ -626,3 +703,46 @@ void SimulationComputePass::advectPass(Gfx::PShaderBuffer quantity,
                             Gfx::SE_ACCESS_SHADER_READ_BIT,
                             Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT);
 }
+
+void SimulationComputePass::reinitializeLevelSetPass(Seele::Gfx::PShaderBuffer phi, Seele::Gfx::PShaderBuffer phi0) {
+  graphics->beginDebugRegion("ReinitializeLevelSet");
+  // Snapshot current phi into phi0 for the sign function
+  phi->pipelineBarrier(Gfx::SE_ACCESS_SHADER_WRITE_BIT,
+                       Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT,
+                       Gfx::SE_ACCESS_TRANSFER_READ_BIT,
+                             Gfx::SE_PIPELINE_STAGE_TRANSFER_BIT);
+  graphics->copyBuffer(phi, phi0);
+  phi0->pipelineBarrier(Gfx::SE_ACCESS_TRANSFER_WRITE_BIT,
+                              Gfx::SE_PIPELINE_STAGE_TRANSFER_BIT,
+                              Gfx::SE_ACCESS_SHADER_READ_BIT,
+                              Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT);
+  phi->pipelineBarrier(Gfx::SE_ACCESS_TRANSFER_READ_BIT,
+                             Gfx::SE_PIPELINE_STAGE_TRANSFER_BIT,
+                             Gfx::SE_ACCESS_SHADER_READ_BIT | Gfx::SE_ACCESS_SHADER_WRITE_BIT,
+                             Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT);
+
+  float dtau = reinitDtau;
+  for (int i = 0; i < reinitIterations; ++i) {
+    Gfx::ODescriptorSet reinitSet =
+        reinitialize.descriptorLayout->allocateDescriptorSet();
+    reinitSet->updateConstants("dtau", 0, &dtau);
+    reinitSet->updateBuffer("phi", 0, phi);
+    reinitSet->updateBuffer("phi0", 0, phi0);
+    reinitSet->writeChanges();
+    Gfx::OComputeCommand reinitCommand = graphics->createComputeCommand();
+    reinitCommand->bindPipeline(reinitialize.pipeline);
+    reinitCommand->bindDescriptor(reinitSet);
+    reinitCommand->dispatch(gridSize / reinitialize.threadGroupSize.x,
+                            gridSize / reinitialize.threadGroupSize.y,
+                            gridSize / reinitialize.threadGroupSize.z);
+    graphics->executeCommands(std::move(reinitCommand));
+
+    phi->pipelineBarrier(Gfx::SE_ACCESS_SHADER_WRITE_BIT,
+                               Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT,
+                               Gfx::SE_ACCESS_SHADER_READ_BIT,
+                               Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT);
+    setBoundsPass(phi, 0);
+    SWAP(phi, phi0);
+  }
+  graphics->endDebugRegion();
+}

src/SimulationComputePass.h

@@ -22,6 +22,7 @@ private:
     void divergencePass(Seele::Gfx::PShaderBuffer divergence, Seele::Gfx::PShaderBuffer pressure, Seele::Gfx::PShaderBuffer velocityX, Seele::Gfx::PShaderBuffer velocityY, Seele::Gfx::PShaderBuffer velocityZ);
     void projectPass(Seele::Gfx::PShaderBuffer velocityX, Seele::Gfx::PShaderBuffer velocityY, Seele::Gfx::PShaderBuffer velocityZ, Seele::Gfx::PShaderBuffer pressure);
     void advectPass(Seele::Gfx::PShaderBuffer quantity, Seele::Gfx::PShaderBuffer quantity0, Seele::Gfx::PShaderBuffer velocityX, Seele::Gfx::PShaderBuffer velocityY, Seele::Gfx::PShaderBuffer velocityZ, float dt);
+    void reinitializeLevelSetPass(Seele::Gfx::PShaderBuffer phi, Seele::Gfx::PShaderBuffer phi0);
     struct SetBounds
     {
         Seele::Gfx::OPipelineLayout pipelineLayout;
@@ -66,6 +67,14 @@ private:
         Seele::Gfx::PComputePipeline pipeline;
         constexpr static Seele::UVector threadGroupSize = {32, 8, 1};
     } advect;
+    struct Reinitialize
+    {
+        Seele::Gfx::OPipelineLayout pipelineLayout;
+        Seele::Gfx::ODescriptorLayout descriptorLayout;
+        Seele::Gfx::OComputeShader shader;
+        Seele::Gfx::PComputePipeline pipeline;
+        constexpr static Seele::UVector threadGroupSize = {32, 8, 1};
+    } reinitialize;
     Seele::Gfx::OShaderBuffer phiBuffer;
     Seele::Gfx::OShaderBuffer phi0Buffer;
     Seele::Gfx::OShaderBuffer densityBuffer;
@@ -80,5 +89,6 @@ private:
     Seele::Gfx::OShaderBuffer velocityX0Buffer;
     Seele::Gfx::OShaderBuffer velocityY0Buffer;
     Seele::Gfx::OShaderBuffer velocityZ0Buffer;
+    uint32 frameCounter = 0;
 };
 DEFINE_REF(SimulationComputePass)

src/SurfaceExtractPass.cpp

@@ -1,53 +1,340 @@
 #include "SurfaceExtractPass.h"
 #include "Common.h"
-#include "Graphics/Descriptor.h"
 #include "Graphics/Enums.h"
 #include "Graphics/Graphics.h"
-#include "Graphics/Shader.h"
+#include "Math/Vector.h"
+#include <algorithm>
+#include <cmath>
 
 using namespace Seele;
 
+// Marching Cubes edge table: for each cube configuration (256), a 12-bit mask
+// indicating which edges are intersected by the isosurface.
+// clang-format off
+static const uint16 edgeTable[256] = {
+    0x000, 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c,
+    0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,
+    0x190, 0x099, 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c,
+    0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90,
+    0x230, 0x339, 0x033, 0x13a, 0x636, 0x73f, 0x435, 0x53c,
+    0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30,
+    0x3a0, 0x2a9, 0x1a3, 0x0aa, 0x7a6, 0x6af, 0x5a5, 0x4ac,
+    0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0,
+    0x460, 0x569, 0x663, 0x76a, 0x066, 0x16f, 0x265, 0x36c,
+    0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,
+    0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0x0ff, 0x3f5, 0x2fc,
+    0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0,
+    0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x055, 0x15c,
+    0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950,
+    0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0x0cc,
+    0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0,
+    0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc,
+    0x0cc, 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0,
+    0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c,
+    0x15c, 0x055, 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,
+    0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc,
+    0x2fc, 0x3f5, 0x0ff, 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0,
+    0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c,
+    0x36c, 0x265, 0x16f, 0x066, 0x76a, 0x663, 0x569, 0x460,
+    0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac,
+    0x4ac, 0x5a5, 0x6af, 0x7a6, 0x0aa, 0x1a3, 0x2a9, 0x3a0,
+    0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c,
+    0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x033, 0x339, 0x230,
+    0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c,
+    0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x099, 0x190,
+    0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c,
+    0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x000,
+};
+
+// Marching Cubes triangle table: for each cube configuration, up to 5
+// triangles specified as sequences of edge indices, terminated by -1.
+static const int triTable[256][16] = {
+    {-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {0,8,3,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {0,1,9,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {1,8,3,9,8,1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {1,2,10,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {0,8,3,1,2,10,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {9,2,10,0,2,9,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {2,8,3,2,10,8,10,9,8,-1,-1,-1,-1,-1,-1,-1},
+    {3,11,2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {0,11,2,8,11,0,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {1,9,0,2,3,11,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {1,11,2,1,9,11,9,8,11,-1,-1,-1,-1,-1,-1,-1},
+    {3,10,1,11,10,3,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {0,10,1,0,8,10,8,11,10,-1,-1,-1,-1,-1,-1,-1},
+    {3,9,0,3,11,9,11,10,9,-1,-1,-1,-1,-1,-1,-1},
+    {9,8,10,10,8,11,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {4,7,8,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {4,3,0,7,3,4,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {0,1,9,8,4,7,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {4,1,9,4,7,1,7,3,1,-1,-1,-1,-1,-1,-1,-1},
+    {1,2,10,8,4,7,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {3,4,7,3,0,4,1,2,10,-1,-1,-1,-1,-1,-1,-1},
+    {9,2,10,9,0,2,8,4,7,-1,-1,-1,-1,-1,-1,-1},
+    {2,10,9,2,9,7,2,7,3,7,9,4,-1,-1,-1,-1},
+    {8,4,7,3,11,2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {11,4,7,11,2,4,2,0,4,-1,-1,-1,-1,-1,-1,-1},
+    {9,0,1,8,4,7,2,3,11,-1,-1,-1,-1,-1,-1,-1},
+    {4,7,11,9,4,11,9,11,2,9,2,1,-1,-1,-1,-1},
+    {3,10,1,3,11,10,7,8,4,-1,-1,-1,-1,-1,-1,-1},
+    {1,11,10,1,4,11,1,0,4,7,11,4,-1,-1,-1,-1},
+    {4,7,8,9,0,11,9,11,10,11,0,3,-1,-1,-1,-1},
+    {4,7,11,4,11,9,9,11,10,-1,-1,-1,-1,-1,-1,-1},
+    {9,5,4,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {9,5,4,0,8,3,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {0,5,4,1,5,0,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {8,5,4,8,3,5,3,1,5,-1,-1,-1,-1,-1,-1,-1},
+    {1,2,10,9,5,4,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {3,0,8,1,2,10,4,9,5,-1,-1,-1,-1,-1,-1,-1},
+    {5,2,10,5,4,2,4,0,2,-1,-1,-1,-1,-1,-1,-1},
+    {2,10,5,3,2,5,3,5,4,3,4,8,-1,-1,-1,-1},
+    {9,5,4,2,3,11,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {0,11,2,0,8,11,4,9,5,-1,-1,-1,-1,-1,-1,-1},
+    {0,5,4,0,1,5,2,3,11,-1,-1,-1,-1,-1,-1,-1},
+    {2,1,5,2,5,8,2,8,11,4,8,5,-1,-1,-1,-1},
+    {10,3,11,10,1,3,9,5,4,-1,-1,-1,-1,-1,-1,-1},
+    {4,9,5,0,8,1,8,10,1,8,11,10,-1,-1,-1,-1},
+    {5,4,0,5,0,11,5,11,10,11,0,3,-1,-1,-1,-1},
+    {5,4,8,5,8,10,10,8,11,-1,-1,-1,-1,-1,-1,-1},
+    {9,7,8,5,7,9,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {9,3,0,9,5,3,5,7,3,-1,-1,-1,-1,-1,-1,-1},
+    {0,7,8,0,1,7,1,5,7,-1,-1,-1,-1,-1,-1,-1},
+    {1,5,3,3,5,7,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {9,7,8,9,5,7,10,1,2,-1,-1,-1,-1,-1,-1,-1},
+    {10,1,2,9,5,0,5,3,0,5,7,3,-1,-1,-1,-1},
+    {8,0,2,8,2,5,8,5,7,10,5,2,-1,-1,-1,-1},
+    {2,10,5,2,5,3,3,5,7,-1,-1,-1,-1,-1,-1,-1},
+    {7,9,5,7,8,9,3,11,2,-1,-1,-1,-1,-1,-1,-1},
+    {9,5,7,9,7,2,9,2,0,2,7,11,-1,-1,-1,-1},
+    {2,3,11,0,1,8,1,7,8,1,5,7,-1,-1,-1,-1},
+    {11,2,1,11,1,7,7,1,5,-1,-1,-1,-1,-1,-1,-1},
+    {9,5,8,8,5,7,10,1,3,10,3,11,-1,-1,-1,-1},
+    {5,7,0,5,0,9,7,11,0,1,0,10,11,10,0,-1},
+    {11,10,0,11,0,3,10,5,0,8,0,7,5,7,0,-1},
+    {11,10,5,7,11,5,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {10,6,5,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {0,8,3,5,10,6,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {9,0,1,5,10,6,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {1,8,3,1,9,8,5,10,6,-1,-1,-1,-1,-1,-1},
+    {1,6,5,2,6,1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {1,6,5,1,2,6,3,0,8,-1,-1,-1,-1,-1,-1},
+    {9,6,5,9,0,6,0,2,6,-1,-1,-1,-1,-1,-1},
+    {5,9,8,5,8,2,5,2,6,3,2,8,-1,-1,-1,-1},
+    {2,3,11,10,6,5,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {11,0,8,11,2,0,10,6,5,-1,-1,-1,-1,-1,-1},
+    {0,1,9,2,3,11,5,10,6,-1,-1,-1,-1,-1,-1},
+    {5,10,6,1,9,2,9,11,2,9,8,11,-1,-1,-1,-1},
+    {6,3,11,6,5,3,5,1,3,-1,-1,-1,-1,-1,-1,-1},
+    {0,8,11,0,11,5,0,5,1,5,11,6,-1,-1,-1,-1},
+    {3,11,6,0,3,6,0,6,5,0,5,9,-1,-1,-1,-1},
+    {6,5,9,6,9,11,11,9,8,-1,-1,-1,-1,-1,-1,-1},
+    {5,10,6,4,7,8,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {4,3,0,4,7,3,6,5,10,-1,-1,-1,-1,-1,-1,-1},
+    {1,9,0,5,10,6,8,4,7,-1,-1,-1,-1,-1,-1,-1},
+    {10,6,5,1,9,7,1,7,3,7,9,4,-1,-1,-1,-1},
+    {6,1,2,6,5,1,4,7,8,-1,-1,-1,-1,-1,-1,-1},
+    {1,2,5,5,2,6,3,0,4,3,4,7,-1,-1,-1,-1},
+    {8,4,7,9,0,5,0,6,5,0,2,6,-1,-1,-1,-1},
+    {7,3,9,7,9,4,3,2,9,5,9,6,2,6,9,-1},
+    {3,11,2,7,8,4,10,6,5,-1,-1,-1,-1,-1,-1,-1},
+    {5,10,6,4,7,2,4,2,0,2,7,11,-1,-1,-1,-1},
+    {0,1,9,4,7,8,2,3,11,5,10,6,-1,-1,-1,-1},
+    {9,2,1,9,11,2,9,4,11,7,11,4,5,10,6,-1},
+    {8,4,7,3,11,5,3,5,1,5,11,6,-1,-1,-1,-1},
+    {5,1,11,5,11,6,1,0,11,7,11,4,0,4,11,-1},
+    {0,5,9,0,6,5,0,3,6,11,6,3,8,4,7,-1},
+    {6,5,9,6,9,11,4,7,9,7,11,9,-1,-1,-1,-1},
+    {10,4,9,6,4,10,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {4,10,6,4,9,10,0,8,3,-1,-1,-1,-1,-1,-1},
+    {10,0,1,10,6,0,6,4,0,-1,-1,-1,-1,-1,-1},
+    {8,3,1,8,1,6,8,6,4,6,1,10,-1,-1,-1,-1},
+    {1,4,9,1,2,4,2,6,4,-1,-1,-1,-1,-1,-1,-1},
+    {3,0,8,1,2,9,2,4,9,2,6,4,-1,-1,-1,-1},
+    {0,2,4,4,2,6,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {8,3,2,8,2,4,4,2,6,-1,-1,-1,-1,-1,-1,-1},
+    {10,4,9,10,6,4,11,2,3,-1,-1,-1,-1,-1,-1},
+    {0,8,2,2,8,11,4,9,10,4,10,6,-1,-1,-1,-1},
+    {3,11,2,0,1,6,0,6,4,6,1,10,-1,-1,-1,-1},
+    {6,4,1,6,1,10,4,8,1,2,1,11,8,11,1,-1},
+    {9,6,4,9,3,6,9,1,3,11,6,3,-1,-1,-1,-1},
+    {8,11,1,8,1,0,11,6,1,9,1,4,6,4,1,-1},
+    {3,11,6,3,6,0,0,6,4,-1,-1,-1,-1,-1,-1,-1},
+    {6,4,8,11,6,8,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {7,10,6,7,8,10,8,9,10,-1,-1,-1,-1,-1,-1,-1},
+    {0,7,3,0,10,7,0,9,10,6,7,10,-1,-1,-1,-1},
+    {10,6,7,1,10,7,1,7,8,1,8,0,-1,-1,-1,-1},
+    {10,6,7,10,7,1,1,7,3,-1,-1,-1,-1,-1,-1,-1},
+    {1,2,6,1,6,8,1,8,9,8,6,7,-1,-1,-1,-1},
+    {2,6,9,2,9,1,6,7,9,0,9,3,7,3,9,-1},
+    {7,8,0,7,0,6,6,0,2,-1,-1,-1,-1,-1,-1,-1},
+    {7,3,2,6,7,2,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {2,3,11,10,6,8,10,8,9,8,6,7,-1,-1,-1,-1},
+    {2,0,7,2,7,11,0,9,7,6,7,10,9,10,7,-1},
+    {1,8,0,1,7,8,1,10,7,6,7,10,2,3,11,-1},
+    {11,2,1,11,1,7,10,6,1,6,7,1,-1,-1,-1,-1},
+    {8,9,6,8,6,7,9,1,6,11,6,3,1,3,6,-1},
+    {0,9,1,11,6,7,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {7,8,0,7,0,6,3,11,0,11,6,0,-1,-1,-1,-1},
+    {7,11,6,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {7,6,11,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {3,0,8,11,7,6,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {0,1,9,11,7,6,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {8,1,9,8,3,1,11,7,6,-1,-1,-1,-1,-1,-1},
+    {10,1,2,6,11,7,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {1,2,10,3,0,8,6,11,7,-1,-1,-1,-1,-1,-1},
+    {2,9,0,2,10,9,6,11,7,-1,-1,-1,-1,-1,-1},
+    {6,11,7,2,10,3,10,8,3,10,9,8,-1,-1,-1,-1},
+    {7,2,3,6,2,7,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {7,0,8,7,6,0,6,2,0,-1,-1,-1,-1,-1,-1,-1},
+    {2,7,6,2,3,7,0,1,9,-1,-1,-1,-1,-1,-1,-1},
+    {1,6,2,1,8,6,1,9,8,8,7,6,-1,-1,-1,-1},
+    {10,7,6,10,1,7,1,3,7,-1,-1,-1,-1,-1,-1,-1},
+    {10,7,6,1,7,10,1,8,7,1,0,8,-1,-1,-1,-1},
+    {0,3,7,0,7,10,0,10,9,6,10,7,-1,-1,-1,-1},
+    {7,6,10,7,10,8,8,10,9,-1,-1,-1,-1,-1,-1,-1},
+    {6,8,4,11,8,6,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {3,6,11,3,0,6,0,4,6,-1,-1,-1,-1,-1,-1,-1},
+    {8,6,11,8,4,6,9,0,1,-1,-1,-1,-1,-1,-1,-1},
+    {9,4,6,9,6,3,9,3,1,11,3,6,-1,-1,-1,-1},
+    {6,8,4,6,11,8,2,10,1,-1,-1,-1,-1,-1,-1,-1},
+    {1,2,10,3,0,11,0,6,11,0,4,6,-1,-1,-1,-1},
+    {4,11,8,4,6,11,0,2,9,2,10,9,-1,-1,-1,-1},
+    {10,9,3,10,3,2,9,4,3,11,3,6,4,6,3,-1},
+    {8,2,3,8,4,2,4,6,2,-1,-1,-1,-1,-1,-1,-1},
+    {0,4,2,4,6,2,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {1,9,0,2,3,4,2,4,6,4,3,8,-1,-1,-1,-1},
+    {1,9,4,1,4,2,2,4,6,-1,-1,-1,-1,-1,-1,-1},
+    {8,1,3,8,6,1,8,4,6,6,10,1,-1,-1,-1,-1},
+    {10,1,0,10,0,6,6,0,4,-1,-1,-1,-1,-1,-1,-1},
+    {4,6,3,4,3,8,6,10,3,0,3,9,10,9,3,-1},
+    {10,9,4,6,10,4,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {4,9,5,7,6,11,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {0,8,3,4,9,5,11,7,6,-1,-1,-1,-1,-1,-1,-1},
+    {5,0,1,5,4,0,7,6,11,-1,-1,-1,-1,-1,-1,-1},
+    {11,7,6,8,3,4,3,5,4,3,1,5,-1,-1,-1,-1},
+    {9,5,4,10,1,2,7,6,11,-1,-1,-1,-1,-1,-1,-1},
+    {6,11,7,1,2,10,0,8,3,4,9,5,-1,-1,-1,-1},
+    {7,6,11,5,4,10,4,2,10,4,0,2,-1,-1,-1,-1},
+    {3,4,8,3,5,4,3,2,5,10,5,2,11,7,6,-1},
+    {7,2,3,7,6,2,5,4,9,-1,-1,-1,-1,-1,-1,-1},
+    {9,5,4,0,8,6,0,6,2,6,8,7,-1,-1,-1,-1},
+    {3,6,2,3,7,6,1,5,0,5,4,0,-1,-1,-1,-1},
+    {6,2,8,6,8,7,2,1,8,4,8,5,1,5,8,-1},
+    {9,5,4,10,1,6,1,7,6,1,3,7,-1,-1,-1,-1},
+    {1,6,10,1,7,6,1,0,7,8,7,0,9,5,4,-1},
+    {4,0,10,4,10,5,0,3,10,6,10,7,3,7,10,-1},
+    {7,6,10,7,10,8,5,4,10,4,8,10,-1,-1,-1,-1},
+    {6,9,5,6,11,9,11,8,9,-1,-1,-1,-1,-1,-1,-1},
+    {3,6,11,0,6,3,0,5,6,0,9,5,-1,-1,-1,-1},
+    {0,11,8,0,5,11,0,1,5,5,6,11,-1,-1,-1,-1},
+    {6,11,3,6,3,5,5,3,1,-1,-1,-1,-1,-1,-1,-1},
+    {1,2,10,9,5,11,9,11,8,11,5,6,-1,-1,-1,-1},
+    {0,11,3,0,6,11,0,9,6,5,6,9,1,2,10,-1},
+    {11,8,5,11,5,6,8,0,5,10,5,2,0,2,5,-1},
+    {6,11,3,6,3,5,2,10,3,10,5,3,-1,-1,-1,-1},
+    {5,8,9,5,2,8,5,6,2,3,8,2,-1,-1,-1,-1},
+    {9,5,6,9,6,0,0,6,2,-1,-1,-1,-1,-1,-1,-1},
+    {1,5,8,1,8,0,5,6,8,3,8,2,6,2,8,-1},
+    {1,5,6,2,1,6,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {1,3,6,1,6,10,3,8,6,5,6,9,8,9,6,-1},
+    {10,1,0,10,0,6,9,5,0,5,6,0,-1,-1,-1,-1},
+    {0,3,8,5,6,10,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {10,5,6,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {11,5,10,7,5,11,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {11,5,10,11,7,5,8,3,0,-1,-1,-1,-1,-1,-1},
+    {5,11,7,5,10,11,1,9,0,-1,-1,-1,-1,-1,-1},
+    {10,7,5,10,11,7,9,8,1,8,3,1,-1,-1,-1,-1},
+    {11,1,2,11,7,1,7,5,1,-1,-1,-1,-1,-1,-1},
+    {0,8,3,1,2,7,1,7,5,7,2,11,-1,-1,-1,-1},
+    {9,7,5,9,2,7,9,0,2,2,11,7,-1,-1,-1,-1},
+    {7,5,2,7,2,11,5,9,2,3,2,8,9,8,2,-1},
+    {2,5,10,2,3,5,3,7,5,-1,-1,-1,-1,-1,-1,-1},
+    {8,2,0,8,5,2,8,7,5,10,2,5,-1,-1,-1,-1},
+    {9,0,1,5,10,3,5,3,7,3,10,2,-1,-1,-1,-1},
+    {9,8,2,9,2,1,8,7,2,10,2,5,7,5,2,-1},
+    {1,3,5,3,7,5,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {0,8,7,0,7,1,1,7,5,-1,-1,-1,-1,-1,-1,-1},
+    {9,0,3,9,3,5,5,3,7,-1,-1,-1,-1,-1,-1,-1},
+    {9,8,7,5,9,7,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {5,8,4,5,10,8,10,11,8,-1,-1,-1,-1,-1,-1,-1},
+    {5,0,4,5,11,0,5,10,11,11,3,0,-1,-1,-1,-1},
+    {0,1,9,8,4,10,8,10,11,10,4,5,-1,-1,-1,-1},
+    {10,11,4,10,4,5,11,3,4,9,4,1,3,1,4,-1},
+    {2,5,1,2,8,5,2,11,8,4,5,8,-1,-1,-1,-1},
+    {0,4,11,0,11,3,4,5,11,2,11,1,5,1,11,-1},
+    {0,2,5,0,5,9,2,11,5,4,5,8,11,8,5,-1},
+    {9,4,5,2,11,3,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {2,5,10,3,5,2,3,4,5,3,8,4,-1,-1,-1,-1},
+    {5,10,2,5,2,4,4,2,0,-1,-1,-1,-1,-1,-1,-1},
+    {3,10,2,3,5,10,3,8,5,4,5,8,0,1,9,-1},
+    {5,10,2,5,2,4,1,9,2,9,4,2,-1,-1,-1,-1},
+    {8,4,5,8,5,3,3,5,1,-1,-1,-1,-1,-1,-1,-1},
+    {0,4,5,1,0,5,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {8,4,5,8,5,3,9,0,5,0,3,5,-1,-1,-1,-1},
+    {9,4,5,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {4,11,7,4,9,11,9,10,11,-1,-1,-1,-1,-1,-1},
+    {0,8,3,4,9,7,9,11,7,9,10,11,-1,-1,-1,-1},
+    {1,10,11,1,11,4,1,4,0,7,4,11,-1,-1,-1,-1},
+    {3,1,4,3,4,8,1,10,4,7,4,11,10,11,4,-1},
+    {4,11,7,9,11,4,9,2,11,9,1,2,-1,-1,-1,-1},
+    {9,7,4,9,11,7,9,1,11,2,11,1,0,8,3,-1},
+    {11,7,4,11,4,2,2,4,0,-1,-1,-1,-1,-1,-1,-1},
+    {11,7,4,11,4,2,8,3,4,3,2,4,-1,-1,-1,-1},
+    {2,9,10,2,7,9,2,3,7,7,4,9,-1,-1,-1,-1},
+    {9,10,7,9,7,4,10,2,7,8,7,0,2,0,7,-1},
+    {3,7,10,3,10,2,7,4,10,1,10,0,4,0,10,-1},
+    {1,10,2,8,7,4,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {4,9,1,4,1,7,7,1,3,-1,-1,-1,-1,-1,-1,-1},
+    {4,9,1,4,1,7,0,8,1,8,7,1,-1,-1,-1,-1},
+    {4,0,3,7,4,3,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {4,8,7,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {9,10,8,10,11,8,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {3,0,9,3,9,11,11,9,10,-1,-1,-1,-1,-1,-1,-1},
+    {0,1,10,0,10,8,8,10,11,-1,-1,-1,-1,-1,-1,-1},
+    {3,1,10,11,3,10,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {1,2,11,1,11,9,9,11,8,-1,-1,-1,-1,-1,-1,-1},
+    {3,0,9,3,9,11,1,2,9,2,11,9,-1,-1,-1,-1},
+    {0,2,11,8,0,11,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {3,2,11,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {2,3,8,2,8,10,10,8,9,-1,-1,-1,-1,-1,-1,-1},
+    {9,10,2,0,9,2,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {2,3,8,2,8,10,0,1,8,1,10,8,-1,-1,-1,-1},
+    {1,10,2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {1,3,8,9,1,8,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {0,9,1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {0,3,8,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+    {-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},
+};
+// clang-format on
+
+// Linearly interpolate the position where the isosurface crosses the edge
+// between two grid points.
+static Vector vertexInterp(float isolevel, const Vector &p1, const Vector &p2,
+                           float v1, float v2) {
+  if (std::abs(v1 - v2) < 1e-10f)
+    return p1;
+  float t = (isolevel - v1) / (v2 - v1);
+  return Vector(p1.x + t * (p2.x - p1.x), p1.y + t * (p2.y - p1.y),
+                p1.z + t * (p2.z - p1.z));
+}
+
 SurfaceExtractPass::SurfaceExtractPass(Seele::Gfx::PGraphics graphics)
     : RenderPass(graphics) {
-  pipelineLayout =
-      graphics->createPipelineLayout("SurfaceExtractPipelineLayout");
-  descriptorLayout = graphics->createDescriptorLayout("params");
-  descriptorLayout->addDescriptorBinding(Gfx::DescriptorBinding{
-      .name = "density",
-      .descriptorType = Gfx::SE_DESCRIPTOR_TYPE_STORAGE_BUFFER,
-  });
-  descriptorLayout->addDescriptorBinding(Gfx::DescriptorBinding{
-      .name = "vertexBuffer",
-      .descriptorType = Gfx::SE_DESCRIPTOR_TYPE_STORAGE_BUFFER,
-  });
-  descriptorLayout->addDescriptorBinding(Gfx::DescriptorBinding{
-      .name = "indexBuffer",
-      .descriptorType = Gfx::SE_DESCRIPTOR_TYPE_STORAGE_BUFFER,
-  });
-  descriptorLayout->addDescriptorBinding(Gfx::DescriptorBinding{
-      .name = "surfaceCounts",
-      .descriptorType = Gfx::SE_DESCRIPTOR_TYPE_STORAGE_BUFFER,
-  });
-  descriptorLayout->create();
   vertexBuffer = graphics->createShaderBuffer(ShaderBufferCreateInfo{
-      .sourceData =
-          {
-              .size = 1024 * 1024, // TODO: make this dynamic
-          },
       .name = "SurfaceExtractVertexBuffer",
   });
   indexBuffer = graphics->createShaderBuffer(ShaderBufferCreateInfo{
-      .sourceData =
-          {
-              .size = 1024 * 1024, // TODO: make this dynamic
-          },
       .name = "SurfaceExtractIndexBuffer",
   });
   surfaceCounts = graphics->createShaderBuffer(ShaderBufferCreateInfo{
       .sourceData =
           {
-              .size = sizeof(uint32) * 3, // TODO: make this dynamic
+              .size = sizeof(uint32) * 4,
           },
+      .numElements = 4,
+      .usage = Gfx::SE_BUFFER_USAGE_INDIRECT_BUFFER_BIT,
+      .name = "SurfaceExtractCountsBuffer",
   });
 }
 
@@ -58,29 +345,148 @@ void SurfaceExtractPass::beginFrame(
     const Seele::Component::Transform &transform) {}
 
 void SurfaceExtractPass::render() {
-  Array<float> densityData(densityBuffer->getNumElements());
+  Array<float> levelSet(gridSize * gridSize * gridSize);
-  densityBuffer->readContents(0, densityData.size() * sizeof(float),
+  phiBuffer->readContents(0, levelSet.size() * sizeof(float),
-                              densityData.data());
+                          levelSet.data());
-  for (int x = 1; x < gridSize - 3; ++x) {
-    for (int y = 1; y < gridSize - 3; ++y) {
-      for (int z = 1; z < gridSize - 3; ++z) {
-        auto getDensity = [&](int x, int y, int z) {
-          return densityData[x + gridSize * (y + gridSize * z)];
-        };
-        float cube[8] = {
-            getDensity(x, y, z),
-            getDensity(x + 1, y, z),
-            getDensity(x, y + 1, z),
-            getDensity(x + 1, y + 1, z),
-            getDensity(x, y, z + 1),
-            getDensity(x + 1, y, z + 1),
-            getDensity(x, y + 1, z + 1),
-            getDensity(x + 1, y + 1, z + 1)
-        };
 
+  Array<Vector> vertices;
+  Array<uint32> indices;
+  vertices.reserve(1024);
+  indices.reserve(1024);
+
+  const float isolevel = 0.0f;
+  const float cellSize = 1.0f / static_cast<float>(gridSize);
+
+  // Helper to sample the level set with bounds checking
+  auto sample = [&](uint32 x, uint32 y, uint32 z) -> float {
+    return levelSet[z * gridSize * gridSize + y * gridSize + x];
+  };
+
+  // Helper to compute world-space position from grid coordinates
+  auto gridPos = [&](uint32 x, uint32 y, uint32 z) -> Vector {
+    return Vector(static_cast<float>(x) * cellSize,
+                  static_cast<float>(y) * cellSize,
+                  static_cast<float>(z) * cellSize);
+  };
+
+  // Marching Cubes: iterate over all cells in the grid
+  for (uint32 z = 1; z < gridSize - 2; ++z) {
+    for (uint32 y = 1; y < gridSize - 2; ++y) {
+      for (uint32 x = 1; x < gridSize - 2; ++x) {
+        // Fetch the 8 corner values of this cell
+        float val[8];
+        val[0] = sample(x, y, z);
+        val[1] = sample(x + 1, y, z);
+        val[2] = sample(x + 1, y + 1, z);
+        val[3] = sample(x, y + 1, z);
+        val[4] = sample(x, y, z + 1);
+        val[5] = sample(x + 1, y, z + 1);
+        val[6] = sample(x + 1, y + 1, z + 1);
+        val[7] = sample(x, y + 1, z + 1);
+
+        // Determine the cube index based on which corners are inside
+        uint32 cubeIndex = 0;
+        if (val[0] < isolevel) cubeIndex |= 1;
+        if (val[1] < isolevel) cubeIndex |= 2;
+        if (val[2] < isolevel) cubeIndex |= 4;
+        if (val[3] < isolevel) cubeIndex |= 8;
+        if (val[4] < isolevel) cubeIndex |= 16;
+        if (val[5] < isolevel) cubeIndex |= 32;
+        if (val[6] < isolevel) cubeIndex |= 64;
+        if (val[7] < isolevel) cubeIndex |= 128;
+
+        // Skip if entirely inside or outside
+        if (edgeTable[cubeIndex] == 0)
+          continue;
+
+        // Compute the 8 corner positions
+        Vector pos[8];
+        pos[0] = gridPos(x, y, z);
+        pos[1] = gridPos(x + 1, y, z);
+        pos[2] = gridPos(x + 1, y + 1, z);
+        pos[3] = gridPos(x, y + 1, z);
+        pos[4] = gridPos(x, y, z + 1);
+        pos[5] = gridPos(x + 1, y, z + 1);
+        pos[6] = gridPos(x + 1, y + 1, z + 1);
+        pos[7] = gridPos(x, y + 1, z + 1);
+
+        // Interpolate vertices along the 12 edges where the surface crosses
+        Vector vertList[12];
+        if (edgeTable[cubeIndex] & 1)
+          vertList[0] = vertexInterp(isolevel, pos[0], pos[1], val[0], val[1]);
+        if (edgeTable[cubeIndex] & 2)
+          vertList[1] = vertexInterp(isolevel, pos[1], pos[2], val[1], val[2]);
+        if (edgeTable[cubeIndex] & 4)
+          vertList[2] = vertexInterp(isolevel, pos[2], pos[3], val[2], val[3]);
+        if (edgeTable[cubeIndex] & 8)
+          vertList[3] = vertexInterp(isolevel, pos[3], pos[0], val[3], val[0]);
+        if (edgeTable[cubeIndex] & 16)
+          vertList[4] = vertexInterp(isolevel, pos[4], pos[5], val[4], val[5]);
+        if (edgeTable[cubeIndex] & 32)
+          vertList[5] = vertexInterp(isolevel, pos[5], pos[6], val[5], val[6]);
+        if (edgeTable[cubeIndex] & 64)
+          vertList[6] = vertexInterp(isolevel, pos[6], pos[7], val[6], val[7]);
+        if (edgeTable[cubeIndex] & 128)
+          vertList[7] = vertexInterp(isolevel, pos[7], pos[4], val[7], val[4]);
+        if (edgeTable[cubeIndex] & 256)
+          vertList[8] = vertexInterp(isolevel, pos[0], pos[4], val[0], val[4]);
+        if (edgeTable[cubeIndex] & 512)
+          vertList[9] = vertexInterp(isolevel, pos[1], pos[5], val[1], val[5]);
+        if (edgeTable[cubeIndex] & 1024)
+          vertList[10] = vertexInterp(isolevel, pos[2], pos[6], val[2], val[6]);
+        if (edgeTable[cubeIndex] & 2048)
+          vertList[11] = vertexInterp(isolevel, pos[3], pos[7], val[3], val[7]);
+
+        // Emit triangles from the triangle table
+        for (int i = 0; triTable[cubeIndex][i] != -1; i += 3) {
+          uint32 base = static_cast<uint32>(vertices.size());
+          vertices.add(vertList[triTable[cubeIndex][i]]);
+          vertices.add(vertList[triTable[cubeIndex][i + 1]]);
+          vertices.add(vertList[triTable[cubeIndex][i + 2]]);
+          indices.add(base);
+          indices.add(base + 1);
+          indices.add(base + 2);
+        }
       }
     }
   }
+
+  uint32 counts[4] = {
+      static_cast<uint32>(indices.size()), // index count
+      1,                                   // instance count
+      0,                                   // first index
+      0,                                   // first instance
+  };
+
+  const uint64 vertexBytes = std::max<uint64>(
+      1, static_cast<uint64>(vertices.size()) * sizeof(Vector));
+  const uint64 indexBytes =
+      std::max<uint64>(1, static_cast<uint64>(indices.size()) * sizeof(uint32));
+
+  vertexBuffer->rotateBuffer(vertexBytes, false);
+  indexBuffer->rotateBuffer(indexBytes, false);
+  surfaceCounts->rotateBuffer(sizeof(counts), false);
+
+  if (!vertices.empty()) {
+    vertexBuffer->updateContents(0, vertices.size() * sizeof(Vector),
+                                 vertices.data());
+    vertexBuffer->pipelineBarrier(Gfx::SE_ACCESS_TRANSFER_WRITE_BIT,
+                                  Gfx::SE_PIPELINE_STAGE_TRANSFER_BIT,
+                                  Gfx::SE_ACCESS_SHADER_READ_BIT,
+                                  Gfx::SE_PIPELINE_STAGE_VERTEX_SHADER_BIT);
+  }
+  if (!indices.empty()) {
+    indexBuffer->updateContents(0, indices.size() * sizeof(uint32),
+                                indices.data());
+    indexBuffer->pipelineBarrier(
+        Gfx::SE_ACCESS_TRANSFER_WRITE_BIT, Gfx::SE_PIPELINE_STAGE_TRANSFER_BIT,
+        Gfx::SE_ACCESS_INDEX_READ_BIT, Gfx::SE_PIPELINE_STAGE_VERTEX_INPUT_BIT);
+  }
+  surfaceCounts->updateContents(0, sizeof(counts), counts);
+  surfaceCounts->pipelineBarrier(Gfx::SE_ACCESS_TRANSFER_WRITE_BIT,
+                                 Gfx::SE_PIPELINE_STAGE_TRANSFER_BIT,
+                                 Gfx::SE_ACCESS_INDIRECT_COMMAND_READ_BIT,
+                                 Gfx::SE_PIPELINE_STAGE_DRAW_INDIRECT_BIT);
 }
 
 void SurfaceExtractPass::endFrame() {}
@@ -92,31 +498,5 @@ void SurfaceExtractPass::publishOutputs() {
 }
 
 void SurfaceExtractPass::createRenderPass() {
-  densityBuffer = resources->requestBuffer("DENSITY");
+  phiBuffer = resources->requestBuffer("PHI");
-}
-
-void SurfaceExtractPass::computeEdgeTable() {
-    for(uint i = 0; i < 256; ++i) {
-        uint32 edgeFlags = 0;
-        for(int j = 0; j < 8; ++j) {
-            if(i & (1 << j)) {
-                edgeFlags |= (1 << j);
-            }
-        }
-        edgeTable[i] = edgeFlags;
-    }
-}
-
-void SurfaceExtractPass::computeTriTable() {
-    for(uint i = 0; i < 256; ++i) {
-        uint32 edgeFlags = edgeTable[i];
-        int triIndex = 0;
-        for(int j = 0; j < 16; ++j) {
-            if(edgeFlags & (1 << j)) {
-                triTable[i][triIndex++] = j;
-            } else {
-                triTable[i][triIndex] = -1;
-            }
-        }
-    }
 }

src/SurfaceExtractPass.h

@@ -17,20 +17,9 @@ public:
   virtual void createRenderPass() override;
 
 private:
-  void computeEdgeTable();
-  void computeTriTable();
-  Seele::Gfx::OPipelineLayout pipelineLayout;
-  Seele::Gfx::ODescriptorLayout descriptorLayout;
-  Seele::Gfx::ODescriptorSet descriptorSet;
-  Seele::Gfx::OComputeShader computeShader;
-  Seele::Gfx::PComputePipeline pipeline;
-
   Seele::Gfx::OShaderBuffer vertexBuffer;
   Seele::Gfx::OShaderBuffer indexBuffer;
   Seele::Gfx::OShaderBuffer surfaceCounts;
-  Seele::Gfx::PShaderBuffer densityBuffer;
+  Seele::Gfx::PShaderBuffer phiBuffer;
-
-  static uint32_t edgeTable[256];
-  static uint32_t triTable[256][16];
 };
 DEFINE_REF(SurfaceExtractPass)