Signed Distance Function(SDF) to Mesh
by Ethan Alexander Shulman (November 30, 2021)

In this article I'm going to describe an algorithm for generating triangle meshes from signed or unsigned distance functions. The algorithm works by generating a voxel mesh and snapping the vertices to the nearest point on the surface of the shape defined by the distance function. I'll be writing all the example code in C# with the Unity engine but it can be applied to other engines or languages.


A sphere and cuboid generated using the algorithm.



Voxels
The algorithm starts with voxel mesh generation, loop through an array of voxels and create triangles for the exposed faces of filled in voxels. We will use Unity's new NativeArray's to store the vertices and indices we generate, they are allocated in the Awake() function. You can see the mesh generation in the Generate() function.
using UnityEngine; using Unity.Collections; public class VoxelChunk : MonoBehaviour { public int width = 16, height = 16, length = 16; [HideInInspector] public bool[] voxels;//Our voxel data, true where voxels exist. Index from voxel position = x+y*width+z*width*height. [HideInInspector] public Mesh mesh; //Buffer arrays for storing generated vertices and indices private NativeArray<Vector3> vertices; private NativeArray<ushort> indices; private void Awake() { //Create Mesh object and attach to MeshFilter. MeshFilter meshFilter = GetComponent<MeshFilter>(); if (meshFilter != null) { mesh = new Mesh(); meshFilter.sharedMesh = mesh; } //Allocate buffers. voxels = new bool[width*height*length]; vertices = new NativeArray<Vector3>(voxels.Length*12, Allocator.Persistent); indices = new NativeArray<ushort>(voxels.Length*12, Allocator.Persistent); } private void OnDestroy() { //Free NativeArrays. if (vertices.IsCreated) { vertices.Dispose(); indices.Dispose(); } } public void Generate() { //Variables for array indexing. int nv = 0, ni = 0, stride = width*height, id = 0, len = length, het = height, wdt = width, l1 = len-1, h1 = het-1, w1 = wdt-1; //Voxel scaling/offset values, mesh has normalized coordinates between -0.5 and 0.5. float ws = 1f/(float)wdt, hs = 1f/(float)het, ls = 1f/(float)len, vx = ws*0.5f, vy = hs*0.5f, vz = ls*0.5f; //Loop through all voxels in volume. for (int z = 0; z < len; z++) { for (int y = 0; y < het; y++) { for (int x = 0; x < wdt; x++) { //Check if voxel exists at point in volume. if (voxels[id]) { //Voxel center position. float px = (x+0.5f)*ws-0.5f, py = (y+0.5f)*hs-0.5f, pz = (z+0.5f)*ls-0.5f; //Add bottom face if at border between filled and empty. if (y == 0 || !voxels[id-width]) { vertices[nv++] = new Vector3(px-vx, py-vy, pz-vz); vertices[nv++] = new Vector3(px+vx, py-vy, pz-vz); vertices[nv++] = new Vector3(px-vx, py-vy, pz+vz); vertices[nv++] = new Vector3(px+vx, py-vy, pz+vz); indices[ni++] = (ushort)(nv-4); indices[ni++] = (ushort)(nv-3); indices[ni++] = (ushort)(nv-2); indices[ni++] = (ushort)(nv-1); indices[ni++] = (ushort)(nv-2); indices[ni++] = (ushort)(nv-3); } //Add left face if at border between filled and empty. if (x == 0 || !voxels[id-1]) { vertices[nv++] = new Vector3(px-vx, py-vy, pz-vz); vertices[nv++] = new Vector3(px-vx, py+vy, pz-vz); vertices[nv++] = new Vector3(px-vx, py-vy, pz+vz); vertices[nv++] = new Vector3(px-vx, py+vy, pz+vz); indices[ni++] = (ushort)(nv-2); indices[ni++] = (ushort)(nv-3); indices[ni++] = (ushort)(nv-4); indices[ni++] = (ushort)(nv-3); indices[ni++] = (ushort)(nv-2); indices[ni++] = (ushort)(nv-1); } //Add back face if at border between filled and empty. if (z == 0 || !voxels[id-stride]) { vertices[nv++] = new Vector3(px-vx, py-vy, pz-vz); vertices[nv++] = new Vector3(px-vx, py+vy, pz-vz); vertices[nv++] = new Vector3(px+vx, py-vy, pz-vz); vertices[nv++] = new Vector3(px+vx, py+vy, pz-vz); indices[ni++] = (ushort)(nv-4); indices[ni++] = (ushort)(nv-3); indices[ni++] = (ushort)(nv-2); indices[ni++] = (ushort)(nv-1); indices[ni++] = (ushort)(nv-2); indices[ni++] = (ushort)(nv-3); } //Add right face if at border between filled and empty. if (x == w1 || !voxels[id+1]) { vertices[nv++] = new Vector3(px+vx, py-vy, pz-vz); vertices[nv++] = new Vector3(px+vx, py+vy, pz-vz); vertices[nv++] = new Vector3(px+vx, py-vy, pz+vz); vertices[nv++] = new Vector3(px+vx, py+vy, pz+vz); indices[ni++] = (ushort)(nv-4); indices[ni++] = (ushort)(nv-3); indices[ni++] = (ushort)(nv-2); indices[ni++] = (ushort)(nv-1); indices[ni++] = (ushort)(nv-2); indices[ni++] = (ushort)(nv-3); } //Add front face if at border between filled and empty. if (z == l1 || !voxels[id+stride]) { vertices[nv++] = new Vector3(px-vx, py-vy, pz+vz); vertices[nv++] = new Vector3(px-vx, py+vy, pz+vz); vertices[nv++] = new Vector3(px+vx, py-vy, pz+vz); vertices[nv++] = new Vector3(px+vx, py+vy, pz+vz); indices[ni++] = (ushort)(nv-2); indices[ni++] = (ushort)(nv-3); indices[ni++] = (ushort)(nv-4); indices[ni++] = (ushort)(nv-3); indices[ni++] = (ushort)(nv-2); indices[ni++] = (ushort)(nv-1); } //Add bottom face if at border between filled and empty. if (y == h1 || !voxels[id+width]) { vertices[nv++] = new Vector3(px-vx, py+vy, pz-vz); vertices[nv++] = new Vector3(px+vx, py+vy, pz-vz); vertices[nv++] = new Vector3(px-vx, py+vy, pz+vz); vertices[nv++] = new Vector3(px+vx, py+vy, pz+vz); indices[ni++] = (ushort)(nv-2); indices[ni++] = (ushort)(nv-3); indices[ni++] = (ushort)(nv-4); indices[ni++] = (ushort)(nv-3); indices[ni++] = (ushort)(nv-2); indices[ni++] = (ushort)(nv-1); } } id++; } } } if (ni > 65535) Debug.LogError("Total vertices must be under 65535 due to 16bit index buffer limits."); //Set Unity Mesh to generated vertices/indices. mesh.SetVertices<Vector3>(vertices, 0, nv); mesh.SetIndices<ushort>(indices, 0, ni, MeshTopology.Triangles, 0); //Generate normals from triangles. mesh.RecalculateNormals(); mesh.UploadMeshData(true); } }


We will test the generation function by randomly filling the voxels and calling Generate().
using UnityEngine; public class ProceduralTest : MonoBehaviour { private void Start() { VoxelChunk chunk = GetComponent<VoxelChunk>(); for (int i = 0; i < chunk.voxels.Length; i++) chunk.voxels[i] = Random.value>0.5f;//50/50 chance of voxel being filled chunk.Generate(); } }


Attach these 2 components VoxelChunk and ProceduralTest to a GameObject with a MeshFilter and MeshRenderer in Unity.

Press play in the Unity editor and you should see a mesh that looks something like this.



Signed Distance Function
Next we will define our signed distance function, I'm going to use a sphere which is very simple. The distance from the sphere center minus radius. We will test it by filling voxels where distance is smaller or equal to 0 which means the point is inside the sphere. Put the following code inside ProceduralTest.
private void Start() { VoxelChunk chunk = GetComponent<VoxelChunk>(); int id = 0; for (int z = 0; z < chunk.length; z++) { float nz = z/(float)chunk.length-0.5f;//normalized coordinates in range -0.5 to 0.5 for (int y = 0; y < chunk.height; y++) { float ny = y/(float)chunk.height-0.5f; for (int x = 0; x < chunk.width; x++) { float nx = x/(float)chunk.width-0.5f; chunk.voxels[id++] = sphereDistance(new Vector3(nx,ny,nz)) <= 0f;//fill in voxels where distance <= 0 } } } chunk.Generate(); } private float sphereDistance(Vector3 pos) { return pos.magnitude-0.4f; }


Entering play mode in Unity you should see a sphere made up of voxels.



Now were going to make the generator work directly from the distance function, create a copy of the VoxelChunk class and name it SDFChunk. Instead of a bool array holding voxels, we have a distance function reference using a delegate(C# version of function pointers).
public class SDFChunk : MonoBehaviour { //Delegate for our distance function. public delegate float DistanceFunction(Vector3 pos); public DistanceFunction distance;//Distance function reference used to describe the scene.


And instead of checking if a voxel bool is true we check if the distance function returns smaller or equal to 0.
//Voxel center position. float px = (x+0.5f)*ws-0.5f, py = (y+0.5f)*hs-0.5f, pz = (z+0.5f)*ls-0.5f; if (distance(new Vector3(px, py, pz)) <= 0f) {//If distance <= 0 then inside shape. //Add bottom face if at border between filled and empty. if (y == 0 || distance(new Vector3(px, py-hs, pz)) > 0f) {


Now switching out VoxelChunk for SDFChunk on the GameObject and in ProceduralTest. Executing it in the Unity editor should generate the same voxelized sphere as before.
private void Start() { SDFChunk chunk = GetComponent<SDFChunk>(); chunk.distance = sphereDistance; chunk.Generate(); } private float sphereDistance(Vector3 pos) { return pos.magnitude-0.4f; }



Snapping to Surface
We have voxels that meet our shape up to the precision of the voxels but we can do better, we can deform the voxels to match the surface. The difference gradient between 2 points in the distance function give us the direction/normal to the surface and scaling that by the distance gives us the offset we need for a vertex to match the surface. You can see the gradient function below.
private const float GRAD_EPSILON = 1e-3f;//Too small introduces precision errors, too high decreases accuracy. public Vector3 gradient(Vector3 p) { float center = distance(p); //Returns offset from point to surface. return new Vector3(distance(p+new Vector3(GRAD_EPSILON, 0f, 0f))-center, distance(p+new Vector3(0f, GRAD_EPSILON, 0f))-center, distance(p+new Vector3(0f, 0f, GRAD_EPSILON))-center).normalized*-center; }


In our code before we prepared the voxel cube vertices for each face we needed. We want all the corners to share the same gradient deformation so we will define them once and add the gradient for each corner. The code now looks like this.
//corner vertices Vector3 c000 = new Vector3(px-vx, py-vy, pz-vz), c100 = new Vector3(px+vx, c000.y, c000.z), c001 = new Vector3(c000.x, c000.y, pz+vz), c101 = new Vector3(c100.x, c000.y, c001.z), c010 = new Vector3(c000.x, py+vy, c000.z), c011 = new Vector3(c000.x, c010.y, c001.z), c110 = new Vector3(c100.x, c010.y, c000.z), c111 = new Vector3(c100.x, c010.y, c001.z); //snap corners to surface using distance gradient c000 += gradient(c000); c100 += gradient(c100); c001 += gradient(c001); c101 += gradient(c101); c010 += gradient(c010); c011 += gradient(c011); c110 += gradient(c110); c111 += gradient(c111); //Add bottom face if at border between filled and empty. if (y == 0 || distance(new Vector3(px,py-hs,pz)) > 0f) { vertices[nv++] = c000; vertices[nv++] = c100; vertices[nv++] = c001; vertices[nv++] = c101; indices[ni++] = (ushort)(nv-4); indices[ni++] = (ushort)(nv-3); indices[ni++] = (ushort)(nv-2); indices[ni++] = (ushort)(nv-1); indices[ni++] = (ushort)(nv-2); indices[ni++] = (ushort)(nv-3); } //Add left face if at border between filled and empty. if (x == 0 || distance(new Vector3(px-ws, py, pz)) > 0f) { vertices[nv++] = c000; vertices[nv++] = c010; vertices[nv++] = c001; vertices[nv++] = c011; indices[ni++] = (ushort)(nv-2); indices[ni++] = (ushort)(nv-3); indices[ni++] = (ushort)(nv-4); indices[ni++] = (ushort)(nv-3); indices[ni++] = (ushort)(nv-2); indices[ni++] = (ushort)(nv-1); } //Add back face if at border between filled and empty. if (z == 0 || distance(new Vector3(px, py, pz-ls)) > 0f) { vertices[nv++] = c000; vertices[nv++] = c010; vertices[nv++] = c100; vertices[nv++] = c110; indices[ni++] = (ushort)(nv-4); indices[ni++] = (ushort)(nv-3); indices[ni++] = (ushort)(nv-2); indices[ni++] = (ushort)(nv-1); indices[ni++] = (ushort)(nv-2); indices[ni++] = (ushort)(nv-3); } //Add right face if at border between filled and empty. if (x == w1 || distance(new Vector3(px+ws, py, pz)) > 0f) { vertices[nv++] = c100; vertices[nv++] = c110; vertices[nv++] = c101; vertices[nv++] = c111; indices[ni++] = (ushort)(nv-4); indices[ni++] = (ushort)(nv-3); indices[ni++] = (ushort)(nv-2); indices[ni++] = (ushort)(nv-1); indices[ni++] = (ushort)(nv-2); indices[ni++] = (ushort)(nv-3); } //Add front face if at border between filled and empty. if (z == l1 || distance(new Vector3(px, py, pz+ls)) > 0f) { vertices[nv++] = c001; vertices[nv++] = c011; vertices[nv++] = c101; vertices[nv++] = c111; indices[ni++] = (ushort)(nv-2); indices[ni++] = (ushort)(nv-3); indices[ni++] = (ushort)(nv-4); indices[ni++] = (ushort)(nv-3); indices[ni++] = (ushort)(nv-2); indices[ni++] = (ushort)(nv-1); } //Add bottom face if at border between filled and empty. if (y == h1 || distance(new Vector3(px,py+hs,pz)) > 0f) { vertices[nv++] = c010; vertices[nv++] = c110; vertices[nv++] = c011; vertices[nv++] = c111; indices[ni++] = (ushort)(nv-2); indices[ni++] = (ushort)(nv-3); indices[ni++] = (ushort)(nv-4); indices[ni++] = (ushort)(nv-3); indices[ni++] = (ushort)(nv-2); indices[ni++] = (ushort)(nv-1); }


And that's it! Pressing Play in the Unity editor should generate a sphere. You can find the completed SDFChunk code here.



More
Now this just covered the naive algorithm but there's still plenty of optimizations to be done and issues to solve. First of all there's no UV coordinates so no texturing but if you notice in the above image the triangles aren't evenly spaced. This causes texturing to be distorted so it's suggested you use a world space shader for texturing. The meshes generated are far from optimal, occasionally generating small smushed triangles which aren't efficient for rasterization. As well as no greedy meshing flat surfaces have many repeating triangles that can be reduced.

I've made improved versions of the above voxel and SDF generation classes that support tile atlasing.
ProceduralVoxelChunk.cs
ProceduralSDFChunk.cs

An amazing source of distance functions is this article by Inigo Quilez.


Thanks for Reading!
If you have any questions feel free to reach out to me at xaloezcontact@gmail.com.

XALOEZ.com This website has no ads or tracking, support development via Patreon.