Motivation

SPH fluid simulator solves the Navier-Stokes equations by animating particles that represent water-droplets density (or other fluid) over time. To render the actual fluid, one first has to reconstruct the fluid surface from density data and then apply a rendering method using the fluid BxDF properties.

This is where Marching Cubes algorithm kicks in – it allows us to reconstruct surface from density as demonstrated in Figure 1. In fact, the algorithm works for any scalar-field that is usually being stored in a 3D grid (texture).

Figure 1a: SPH simulator before surface reconstruction (source).

Figure 1b: SPH simulator after surface reconstruction and rendering (source).

It’s Demo Time

Don’t worry if you don’t know how fluid simulation works. In this post we will implement a 2D version of the Marching Cubes algorithm (namely – Marching Squares).

The density in our demonstration is being generated using Meta-balls animated over time as suggested by Jamie Wong, though as I mentioned, any scalar-field would work. I strongly suggest that you first read his blog post where he covers the theory.

The goal of this post is to visualize the Marching Squares algorithm in practice and provide glsl code reference for others. To that end, the simulation runs entirely on the GPU and the code was not optimized for performance, as I wanted it to be readable by others. I also stayed consistent with the mentioned post coding schemes of the cells’ corners.

Without further ado, here’s the shadertoy:

Figure 2: Marching Squares with Meta-balls shadertoy demo.

Please watch at least one full cycle of the shader (76 seconds) prior to reading on.

Shader Breakdown

Let’s unpack this shader: when the value of the scalar-field, generated using Meta-balls, approaches 1, the shape perimeters turn purple. You can think of the Meta-balls generated data has a “height” pointing towards the screen where black is “far away” and white represents surface peaks (see Figure 3).

Figure 3: Meta-balls scalar-field data

Reconstructing the surface when we actually have an analytic expression seems superfluous, but it’s still useful: it will serve as a quick sanity check for the algorithm’s correctness and quality.

First, I voxelized the function into a low-resolution grid and recorded the value at each of the cells’ center (Figure 4). Voxels with scalar-field value of lower than 1 are considered “empty” and marked as red – otherwise they are “non-empty” (e.g. have fluid density) and are painted green. Note that due to under-sampling some of the “smaller” balls would potentially not be recorded as they travel through voxels. The net result is voxels “popping” in and out – this is expected and can be fixed by increasing the sampling rate (with additional memory cost).

Figure 4: Voxelizing the scalar-field.

Voxelizing the scalar-field on the GPU is definitely an overkill and the data should be provided via texture (e.g. MRI scan). However, I’ve decided to pack the whole simulation into a single fragment shader for simplicity and readability.

Moving on to the actual algorithm implementation. Here we have two flavors: surface reconstruction by using the “classic” 16 states of the corners (in cyan, in Figure 5a) and a reconstruction that uses the corners’ weights to linearly interpolate the voxel-boundary and surface intersection-point (in yellow, in Figure 5b). As you can see in Figures 5a,b, the linear-interpolation yields a “less-blocky” reconstruction and animates better as the Meta-balls pass through voxels.

Figure 5a: Marching Squares without linear-interpolation.

Figure 5b: Marching Squares with linear-interpolation.

In my implementation, the surface is being drawn directly on the GPU by calculating line equations and its SDF within each cell on a per-pixel basis. In practice though, this is the place where the implementation will feed the vertex shader with triangles and potentially merge vertices to compress vertex data buffering. The nice thing about drawing the surface contour in the fragment shader is that we avoid triangle facing issues (i.e. the ambiguity in cases 5 and 10 as mentioned here). Full implementation will have to address these issues.

As we saw in the last section, both reconstructions suffer from under-sampling issues. As expected, by increasing our grid resolution we get superior results (Figures 6a,b). Note that linear-interpolation yields better results at relatively cheap costs.

Figure 6a: Marching Squares without linear-interpolation and higher resolution.

Figure 6b: Marching Squares with linear-interpolation and higher resolution.

Observe also that the under-sampling issues are mitigated and the “small” Meta-balls are being recorded in voxels and their surfaces can be reconstructed. The higher volumetric resolution is coming, of course, with an additional memory and texture bandwidth costs and the best grid size depends on the application and the desired quality.

Finally, I’ll note that the algorithm can be easily extended to extract iso-surfaces with different iso-values as shown in Figure 7. This can be useful to generate something like contour maps (e.g. Figure 8).

Figure 7: Marching Squares with iso-surfaces of 0.5, 1.0, 2.0 and 3.0.

Figure 8: Example contour map.

Closing Thoughts

Now that we explained how the shader works, please try to re-watch a full animation cycle (76 seconds) below and take notes of the generated surfaces.

The next step would be, naturally, to extend the implementation for full 3D, reconstruct surface normals and shade the produced surface. This is left as an exercise for the reader.

I hope you found this visualization to be both helpful and inspiring. If you happen to extend this work, please let me know so I can link to your extensions.