Generating mesh from SDFs with Dual Contouring

Oct 20, 2021 · 8 minute read

Dual Contouring is a method to generate a mesh from isosurfaces. While relatively easy to understand, it can surely be easier if it includes some interactive visualizations.

What is isosurface? 🔗

Isourface is a set of all the 3d points where we have constant value (e.g. pressure, temperature, density). As we are going to extract the surface out of the SDF we can use the distance as the value and find the isosurface where the values is zero.

Defining the scene 🔗

For the scene I’m going to use a simple SDF scene that is shown below, it includes a sphere and a cube to test if we can handle sharp edges.

You can get a list of SDFs primitives at the Inigo Quilez’s website .

Creating a grid 🔗

First thing we need todo is to create a grid that contains the entire scene.

The grid size is also the resolution so the bigger the grid is the better the resulting mesh will be, but it will have more vertices and take more time to calculate.

You will probably need the AABB of all the primitives, but I will not cover this here, I will just use a fixed size grid as I already know that my primitives will be inside the grid.

The grid will store the vertices index so it will be an 3d array of int:

Array3D<int> voxels;

Each voxel in the grid will may contain a vertice or may not, if not I like to set them to -1.

Efficiently handling edges 🔗

The edges have a high amount of overlap between the voxels, so if we calculate for each voxel we will be doing a lot of unecessary work.

When computing each voxels’s vertex we need to access all the edges of a single voxel, and we need the intersection point of the edges with the surface and the normal.

So we define the Edge like this:

struct Edge{
    vec3 intersection;
    vec3 normal;
    bool crossed;
    bool sign;
};

To efficiently store each edge we need to keep 3 unique edges per voxel, and to access all the 12 edges for the voxel, we get the 3 from this voxel and the remaining 9 from the neighbours, I will show how later.

struct TripleEdge{
    Edge triplets[3];
};

Now to store the edges you need some sort of 3d addressable array.

Array3D<TripleEdge> edges;

Initializing the edges 🔗

Initializing the edges is the first step needed for the algorithm.

To do that you need a 3d loop of all the voxels then for each voxel you need to popule the 3 unique edges.

First detect if the edge cross the surface, for that we just need to see if one of the vertices is inside an primitive and the other is outside, we can know this if we compare the 2 vertices of the edge, if they have a different sign of surface distance then it’s crossing the surface, then we store this inside the edge crossed.

To store the sign of the edge, just choose one vertice of the edge and set true if is inside.

We also need the intersection, for that I will be doing a simple binary search along the edge to find it. A simple fixed step can also be used.

The normal is just the delta of the SDF function.

Initializing the voxels 🔗

After we have all the edges information we need to initialize the voxels.

To map all the 12 edges to a single voxel we use this table:

struct EdgeRemap{ int x, y, z, e; };
EdgeRemap EDGE_REMAP_TABLE[12] = {
    {0,0,0, 0}, {0,0,0, 1}, {0,0,0, 2},
    {1,0,0, 1}, {1,0,0, 2},
    {0,1,0, 0}, {0,1,0, 2},
    {0,0,1, 0}, {0,0,1, 1},
    {1,0,1, 1},
    {1,1,0, 2}, {0,1,1, 0},
};

Now to access the edge we need something like this:

// edge is a number between 0-12
Edge get_edge(int x, int y, int z, int edge) {
    EdgeRemap remap = EDGE_REMAP_TABLE[edge];
    return edges.get(
        x+remap.x,
        y+remap.y,
        z+remap.z
    ).triplets[remap.e];
}

Then we just loop again trough all the voxels and check if at least one of the 12 edges of the voxel has crossed, if so just create a vertice for it.

The vertex position will be the intersection of the planes of all the crossing edges’s position and normal, this point will keep the sharp edges sharp.

This position can be found by using QEF but I can also be approximated by projecting the point to all the planes.

if(any_edge_crossed){
    voxels.at(x, y, z) = vertices.size();
    vertices.push_back(vertex_position);
}else{
    voxels.at(x, y, z) = -1;
}

You should be able to get the result below:

Connecting the triangles 🔗

Now we need to find the triangles, for that the dual countouring works by finding the quad perpendicular of the crossing edges.

For that just loop normally trough all the voxels and connect the 4 voxels that have the edge crossing.

Example pseudo code:

for(int axis = 0; axis < 3; axis++){
    Edge e = edges.get(x, y, z).triplets[axis];
    if(edge.crossed){
        if(axis == 0){
            // Connect 4 vertices in the perpendicular to X
        }else if(axis == 1){
            // Connect 4 vertices in the perpendicular to Y
        } else {
            // Connect 4 vertices in the perpendicular to Z
        }
    }
}

It should give this result:

References 🔗