Naive surface nets algorithm: how to connect vertices to form triangles?

Question

I'm trying to implement the naive surface nets algorithm (in c++, opengl) as described here. I understand that I need to interpolate endpoints of edges that cross the surface, and then I need to find the center of mass of the results of those interpolations. What I don't understand: I now have a single vertex for each nontrivial cell, how do I order those vertices into triangles that form a mesh?

Answer 1

Let's say your SDF is sampled at half-integer coordinates, i.e. f[0][0][0] is the value in the middle of the cell that spans the 3d space between (0,0,0) and (1,1,1). If you were to voxelize it with minecraft-like cubed mesh, you'd do something along these lines:

for each integer x,y,z {
    // X-major
    if(f[z][y][x-1]*f[z][y][x] < 0) {
        v0 = EmitVertex(vec3(x,y,z));
        v1 = EmitVertex(vec3(x,y+1,z));
        v2 = EmitVertex(vec3(x,y,z+1));
        v3 = EmitVertex(vec3(x,y+1,z+1));
        if(f[z][y][x] < 0) {
            // flip order so normal points along negative X
            swap(v0, v3);
        }
        EmitTriangle(v0, v1, v2);
        EmitTriangle(v2, v1, v3);
    }

    // do the same for Y and Z directions
    // ...
}

(Note: I omit the intricacies of handling f[z][y][x] == 0.)

For surface nets the code remains essentially the same. You use the values at the cell-centers (e.g. f[-1/0][-1/0][-1/0]) to compute the perturbed location p[0][0][0] of the vertex that sits between those cells. On average p[0][0][0] == vec3(0,0,0), but it can be anywhere between (-0.5,-0.5,-0.5) and (0.5,0.5,0.5). Then you run the same code as above, but instead of EmitVertex(vec3(x,y,z)), which would emit the vertex at an integer location, you use the perturbed location: EmitVertex(p[z][y][x]):

// precompute vertex locations
for each integer x,y,z {
    p[z][y][x] = vec3(x,y,z) + deviation(
        f[z-1][y-1][x-1], f[z-1][y-1][x], f[z-1][y][x-1], f[z-1][y][x],
        f[z][y-1][x-1], f[z][y-1][x], f[z][y][x-1], f[z][y][x]
    );
}

// generate mesh
for each integer x,y,z {
    // X-major
    if(f[z][y][x-1]*f[z][y][x] < 0) {
        v0 = EmitVertex(p[z][y][x]);
        v1 = EmitVertex(p[z][y+1][x]);
        v2 = EmitVertex(p[z+1][y][x]);
        v3 = EmitVertex(p[z+1][y+1][x]);
        // same as before ...
    }

    // do the same for Y and Z directions
    // ...
}