Implementing Dual Contouring

If you've read the Dual Contouring paper and the ProcWorld post From Voxels To Polygons but haven't been able to get an implementation working, then this post is for you. I've put together a Dual Contouring implementation and added it to a Github repo. The rest of this post will explain how it works.

The sample application rendering a mesh

My starting point when I originally implemented DC was the reference implementation released by the paper's authors. Their implementation provides the contouring algorithm itself and a QEF implementation, which solves two of the bigger problems. Unfortunately it only supports constructing octrees from the PLY and SOG model formats, which isn't much use it you want to use DC to render terrains.

I had originally planned to build my sample implementation using the reference implementation without modification. I tried that initially but I had forgotten how much work I had actually put in to refactor the code before being able to use it, so rather than do all that again I've just lifted bits from my own renderer. The main changes are reworking the OctreeNode class to remove the class heirarchy and building a wrapper class for the QEF code. I think both of these changes make things a bit easier to understand.

Constructing the Octree

This is the OctreeNode class:

enum OctreeNodeType
{
 Node_None,
 Node_Internal,
 Node_Psuedo,
 Node_Leaf,
};

class OctreeNode
{
public:

 OctreeNode()
  : type(Node_None)
  , min(0, 0, 0)
  , size(0)
  , drawInfo(nullptr)
 {
  for (int i = 0; i < 8; i++)
  {
   children[i] = nullptr;
  }
 }

 OctreeNode(const OctreeNodeType _type)
  : type(_type)
  , min(0, 0, 0)
  , size(0)
  , drawInfo(nullptr)
 {
  for (int i = 0; i < 8; i++)
  {
   children[i] = nullptr;
  }
 }

 OctreeNodeType type;
 ivec3   min;
 int    size;
 OctreeNode*  children[8];
 OctreeDrawInfo* drawInfo;
};

I think this is a fairly easy to understand node definition. The type of the node can be changed which will come in handy when we get to the octree simplification, this is a divergence from the reference implementation which preferred a class heirarchy. The bounding box for the node is determined by the min and size members, and each node has up to 8 children. The drawInfo member is a pointer as only some nodes will need this data. It contains the information required to draw the node as part of the mesh, as such its completely useless on branch (or Internal) nodes.

struct OctreeDrawInfo 
{
 OctreeDrawInfo()
  : index(-1)
  , corners(0)
 {
 }

 int    index;
 int    corners;
 vec3   position;
 vec3   averageNormal;
 QEF    qef;
};

The draw info consists of an index into a vertex buffer for this node's vertex (used to construct polygons to send to the graphics API), the position of the vertex and its normal and two more interesting values. The corners variable contains information about which edges contain the surface, and the qef variable stores the state for this node's QEF which allows us to calculate the vertex position.

Before we can run the contouring algorithm we of course need to actually construct an octree. Conceptually, do this by sampling the volume as 1x1x1 cubes (i.e. voxels) to determine if the surface of the volume crosses any of the cube's edges. The octree is then constructed containing all the cube's which contain the surface, with the 1x1x1 cubes being the the leaf nodes.

My aim with this sample is to provide a simple implemenation. With that in mind we'll actually construct the octree from the root node down. This approach is quite wasteful as all the intermediate parent branches of all leafs are constructed before the leaf nodes are evaluated to determine if they contain the surface. This means some parent nodes will be constructed just to be discarded again. Ideally we'd avoid all this useless work, but for the sake of simplicity we'll just take the performance hit.

OctreeNode* ConstructOctreeNodes(OctreeNode* node)
{
 if (!node)
 {
  return nullptr;
 }

 if (node->size == 1)
 {
  return ConstructLeaf(node);
 }
 
 const int childSize = node->size / 2;
 bool hasChildren = false;

 for (int i = 0; i < 8; i++)
 {
  OctreeNode* child = new OctreeNode;
  child->size = childSize;
  child->min = node->min + (CHILD_MIN_OFFSETS[i] * childSize);
  child->type = Node_Internal;

  node->children[i] = ConstructOctreeNodes(child);
  hasChildren |= (node->children[i] != nullptr);
 }

 if (!hasChildren)
 {
  delete node;
  return nullptr;
 }

 return node;
}

This is pretty straightforward: the root node is constructed and then the recursive calls to ConstructOctreeNode are kicked off, constructing the octree nodes depth first. The terminating case is when the node's size is 1, which makes it a leaf node (or voxel). ConstructLeaf will return nullptr when the leaf does not contain the surface. If all the leaf nodes belonging to a parent are null then the parent node will be deleted. If none of the parent's sibling nodes contain the surface then their parent will be deleted, and so on.

One interesting thing in this function is how we calculate the child node's min value. This works by taking advantage of the indexing of the child nodes. E.g. we know that the child with index 2 is always going to be in the same place with the same offset from the parent. The CHILD_MIN_OFFSETS array holds the offset values:

const ivec3 CHILD_MIN_OFFSETS[] =
{
 // needs to match the vertMap from Dual Contouring impl
 ivec3( 0, 0, 0 ),
 ivec3( 0, 0, 1 ),
 ivec3( 0, 1, 0 ),
 ivec3( 0, 1, 1 ),
 ivec3( 1, 0, 0 ),
 ivec3( 1, 0, 1 ),
 ivec3( 1, 1, 0 ),
 ivec3( 1, 1, 1 ),
};

Since we know the parent node's min and the size of its child nodes, we can calculate the min for any child.

The real meat of the implementation is in the CreateLeaf function, but before we can create the leaf nodes we need some volume data to sample. The volume will be represented by a Signed Distance Function. Using the SDF we can determine if any 3D point is inside or outside the volume. The simplest possible SDF is a sphere volume:

float Sphere(const vec3& worldPosition, const vec3& origin, float radius)
{
 return length(worldPosition - origin) - radius;
}

Given an origin and radius for the sphere we can determine if any point in 3D space is inside or outside the sphere. Points inside the sphere will have a negative density and points outside the sphere will have a positive density value. Using the corners of the leaf node's bounding box as input to the SDF allow us to determine whether any of the edges of the bounding box contain a sign change, i.e. one end of the edge has a negative density value and the other has a positive density value. If an edge has a sign change then that means the surface of the volume (in this case a sphere) crosses the bounding box edge. If any of the leaf nodes have at least one edge with a sign change then we create a vertex for the leaf node, otherwise the leaf node is discarded.

Dual Contouring relies on Hermite data to determine the position the leaf node vertices. This means we need not only the position of the surface on each edge, but also the normal of the surface at the position. Mathematically the surface normal can be determined by finding the derivative of the density function. We can approximate this with any SDF by using the Finite Difference method. This works by going back to the basic principle of derivatives: finding the rate of change. Using the zero crossing position (i.e. the position of the surface on each edge) as a base we can measure the density values on each axis either side of the surface to determine the rate of change on each axis. These values then determine the surface normal:

vec3 CalculateSurfaceNormal(const vec3& p)
{
 const float H = 0.001f;
 const float dx = Density_Func(p + vec3(H, 0.f, 0.f)) - Density_Func(p - vec3(H, 0.f, 0.f));
 const float dy = Density_Func(p + vec3(0.f, H, 0.f)) - Density_Func(p - vec3(0.f, H, 0.f));
 const float dz = Density_Func(p + vec3(0.f, 0.f, H)) - Density_Func(p - vec3(0.f, 0.f, H));

 return glm::normalize(vec3(dx, dy, dz));
}

Using this information we can now create the leaf nodes for the octree. The first step is to determine whether the leaf node contains the surface:

 int corners = 0;
 for (int i = 0; i < 8; i++)
 {
  const ivec3 cornerPos = leaf->min + CHILD_MIN_OFFSETS[i];
  const float density = Density_Func(vec3(cornerPos));
  const int material = density < 0.f ? MATERIAL_SOLID : MATERIAL_AIR;
  corners |= (material << i);
 }

 if (corners == 0 || corners == 255)
 {
  // voxel is full inside or outside the volume
  delete leaf;
  return nullptr;
 }

If you've implemented Marching Cubes then the corners variable should look familiar. Each of the first 8 bits represents a corner of the node's bounding box. If the density value at the corner position is negative the corner is inside the volume, otherwise the corner is outside the volume. If all the corners are outside, or all the corners inside, the volume the node is discarded. Otherwise we need to determine the position of the leaf node's vertex.

 const int MAX_CROSSINGS = 6;
 vec3 positions[MAX_CROSSINGS], normals[MAX_CROSSINGS];
 int edgeCount = 0;

 for (int i = 0; i < 12 && edgeCount < MAX_CROSSINGS; i++)
 {
  const int c1 = edgevmap[i][0];
  const int c2 = edgevmap[i][1];

  const int m1 = (corners >> c1) & 1;
  const int m2 = (corners >> c2) & 1;

  if ((m1 == MATERIAL_AIR && m2 == MATERIAL_AIR) ||
   (m1 == MATERIAL_SOLID && m2 == MATERIAL_SOLID))
  {
   // no zero crossing on this edge
   continue;
  }

  const vec3 p1 = vec3(leaf->min + CHILD_MIN_OFFSETS[c1]);
  const vec3 p2 = vec3(leaf->min + CHILD_MIN_OFFSETS[c2]);
  const vec3 p = ApproximateZeroCrossingPosition(p1, p2);

  positions[edgeCount] = p;
  normals[edgeCount] = CalculateSurfaceNormal(p);

  edgeCount++;
 }

 The edgevmap array is from the reference implementation. It allows us to look up the indices which determine the end points of each edge. Using these indices we can determine whether the responding bit is set in the corners variable we generated when determining if the leaf node contained the surface. Only some of the edges in each leaf node will be set, so we use the corner values to determine each individual edge contains the surface. If it does then we find the position along the edge where the surface crosses (i.e. the density value is zero) and calculate the surface normal using this position. We keep track of all the positions & normals for the leaf node and use these to determine the vertex's position and normal.

There are lots of ways to determine the position of the zero crossing. I think this is probably the simplest: take N samples along the edge and use the position of the sample with the lowest absolute value. Ideally somewhere along the line we'll find the exact position where the density is zero, but generally we'll find a position that is not quite the zero crossing. The amount of error introduced by this approach is determined by the number of samples, but even with N=8 the results are good enough.

vec3 ApproximateZeroCrossingPosition(const vec3& p0, const vec3& p1)
{
 // approximate the zero crossing by finding the min value along the edge
 float minValue = 100000.f;
 float t = 0.f;
 float currentT = 0.f;
 const int steps = 8;
 const float increment = 1.f / (float)steps;
 while (currentT <= 1.f)
 {
  const vec3 p = p0 + ((p1 - p0) * currentT);
  const float density = glm::abs(Density_Func(p));
  if (density < minValue)
  {
   minValue = density;
   t = currentT;
  }

  currentT += increment;
 }

 return p0 + ((p1 - p0) * t);
}

With all edges evaluated the gathered positions and normals can be fed in to the Quadric Error Function. Reading about the QEFs in the Dual Contouring paper was off-putting: I didn't understand the maths and had no idea how to implement it. Thanks to the reference implementation I didn't need to worry about how to implement it. This sample contains my refactoring of the original implementation which makes it a bit easier to use.

 OctreeDrawInfo* drawInfo = new OctreeDrawInfo;
 drawInfo->qef.initialise(edgeCount, normals, positions);
 drawInfo->position = drawInfo->qef.solve();

 const vec3 min = vec3(leaf->min);
 const vec3 max = vec3(leaf->min + ivec3(leaf->size));
 if (drawInfo->position.x < min.x || drawInfo->position.x > max.x ||
  drawInfo->position.y < min.y || drawInfo->position.y > max.y ||
  drawInfo->position.z < min.z || drawInfo->position.z > max.z)
 {
  drawInfo->position = drawInfo->qef.masspoint;
 }

The QEF is initialised using the positions and normals and then solved. The solve function find the position for the vertex with the smallest error, i.e. the vertex is placed as close to the surface as possible using the supplied inputs. Occasionally the QEF will calculate the position outside of the bounding box, so we need to catch this case an clamp the position to the masspoint of the QEF which is just the average of all the input positions.

All that remains is to fill out the rest of the drawInfo for this leaf node:

 for (int i = 0; i < edgeCount; i++)
 {
  drawInfo->averageNormal += normals[i];
 }
 drawInfo->averageNormal = glm::normalize(drawInfo->averageNormal);

 drawInfo->corners = corners;

 leaf->type = Node_Leaf;
 leaf->drawInfo = drawInfo;

 return leaf;

The normal is calculated as an average of the input normals. The corners variable is stored as this is used by the contouring algorithm, and the node's type is changed to a Leaf. Once all the nodes have been evaluated/constructed the octree construction is complete.

Generating the Mesh

Once the octree is constructed we can generate a polygonal mesh to send to a graphics API. The reference implementation provides all the code for this, which makes things considerably easier. There are two variants of the contouring algorithm provided by: the original and a "no intersections" variant. In the sample is code I've refactored from the simpler, original variant.

struct MeshVertex
{
 MeshVertex(const glm::vec3& _xyz, const glm::vec3& _normal)
  : xyz(_xyz)
  , normal(_normal)
 {
 }

 glm::vec3  xyz, normal;
};

typedef std::vector<MeshVertex> VertexBuffer;
typedef std::vector<int> IndexBuffer;

void GenerateMeshFromOctree(OctreeNode* node, VertexBuffer& vertexBuffer, IndexBuffer& indexBuffer)
{
 if (!node)
 {
  return;
 }

 vertexBuffer.clear();
 indexBuffer.clear();

 GenerateVertexIndices(node, vertexBuffer);
 ContourCellProc(node, indexBuffer);
}

GenerateVertexIndices is very straightfoward: the octree is traversed depth-first with each leaf node adding a vertex to the buffer and being assigned a corresponding index. ContourCellProc kicks off a series of recursive calls to three functions: ContourCellProc, ContourFaceProc and ContourEdgeProc. This is discussed in both the original DC paper and the Secret Sauce paper so I'll be brief. The cellproc is only interested in Internal nodes. It calls faceproc for each of the faces of the node edgeproc for each of the edges, and calls cellproc on each of its children. The faceproc is again only interested in Internal nodes. It calls edgeproc for each of the edges on the face and calls faceproc on its children. If edgeproc encounters an edge where none of the nodes sharing the edge are Internal nodes then it calls ContourProcessEdge which generates the polygons. Otherwise it calls edgeproc on its children. 

All those recursive calls eventually result in this function being called:

void ContourProcessEdge(OctreeNode* node[4], int dir, IndexBuffer& indexBuffer)
{
 int minSize = 1000000;  // arbitrary big number
 int minIndex = 0;
 int indices[4] = { -1, -1, -1, -1 };
 bool flip = false;
 bool signChange[4] = { false, false, false, false };

 for (int i = 0; i < 4; i++)
 {
  const int edge = processEdgeMask[dir][i];
  const int c1 = edgevmap[edge][0];
  const int c2 = edgevmap[edge][1];

  const int m1 = (node[i]->drawInfo->corners >> c1) & 1;
  const int m2 = (node[i]->drawInfo->corners >> c2) & 1;

  if (node[i]->size < minSize)
  {
   minSize = node[i]->size;
   minIndex = i;
   flip = m1 != MATERIAL_AIR; 
  }

  indices[i] = node[i]->drawInfo->index;

  signChange[i] = 
   (m1 == MATERIAL_AIR && m2 != MATERIAL_AIR) ||
   (m1 != MATERIAL_AIR && m2 == MATERIAL_AIR);
 }

 if (signChange[minIndex])
 {
  if (!flip)
  {
   indexBuffer.push_back(indices[0]);
   indexBuffer.push_back(indices[1]);
   indexBuffer.push_back(indices[3]);

   indexBuffer.push_back(indices[0]);
   indexBuffer.push_back(indices[3]);
   indexBuffer.push_back(indices[2]);
  }
  else
  {
   indexBuffer.push_back(indices[0]);
   indexBuffer.push_back(indices[3]);
   indexBuffer.push_back(indices[1]);

   indexBuffer.push_back(indices[0]);
   indexBuffer.push_back(indices[2]);
   indexBuffer.push_back(indices[3]);
  }
 }
}

Firstly each of the four neighbouring nodes are processed. The processEdgeMask array lets us look up the appropriate edge for each node, and the materials on either end of the edge are recovered from the corners variable.  The minSize check is neccessary to account for differently sized leaf nodes which are created when the octree is simplified, we'll get to that later. The indices array is filled out using the indexes into the VertexBuffer generated when GenerateVertexIndices was called. 

Finally we only want to produce polygons when appropriate: the signChange array is used to track whether the effected edge of each node contains the surface. When this is true of the minSize node we can generate two triangles connecting the vertices of the four nodes. 

Once GenerateMeshForOctree has finished executing the VertexBuffer and IndexBuffer will be filled with all the data needed to render the surface. The sample uses OpenGL 3.3 in a very basic setup. Firstly the data is uploaded to the GPU:

void Mesh::uploadData(const VertexBuffer& vertices, const IndexBuffer& indices)
{
 if (vertices.empty() || indices.empty())
 {
  return;
 }

 glBindVertexArray(vertexArrayObj_);

 glBindBuffer(GL_ARRAY_BUFFER, vertexBuffer_);
 glBufferData(GL_ARRAY_BUFFER, sizeof(MeshVertex) * vertices.size(), &vertices[0], GL_STATIC_DRAW);

 glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, indexBuffer_);
 glBufferData(GL_ELEMENT_ARRAY_BUFFER, sizeof(int) * indices.size(), &indices[0], GL_STATIC_DRAW);
 numIndices_ = indices.size();

 printf("Mesh: %d vertices %d triangles\n", vertices.size(), indices.size() / 3);

 glBindVertexArray(0);
}

And is then drawn:

 glBindVertexArray(mesh.vertexArrayObj_);
 glDrawElements(GL_TRIANGLES, mesh.numIndices_, GL_UNSIGNED_INT, (void*)0);

Simplifying

The final point to cover is simplifying the octree which allows us to produce a mesh representing the same volume with fewer vertices and triangles. When we used the QEF to calculate the position of the leaf nodes we were in fact also recording the error value produced by the QEF. In the case of leaf nodes this doesn't have much use since we always need to produce a vertex regardless of the error. However, when simplifying the octree the error produced by the QEF can be used control how aggressive we are in collapsing the child nodes.

The mesh generated after very aggressive simplification

At the core of the process is another recursive function, SimplifyOctree. It contains this loop which handles the recursion:

 QEF qef;
 int signs[8] = { -1, -1, -1, -1, -1, -1, -1, -1 };
 int midsign = -1;
 int edgeCount = 0;
 bool isCollapsible = true;

 for (int i = 0; i < 8; i++)
 {
  node->children[i] = SimplifyOctree(node->children[i], threshold);
  if (node->children[i])
  {
   OctreeNode* child = node->children[i];
   if (child->type == Node_Internal)
   {
    isCollapsible = false;
   }
   else
   {
    qef.add(child->drawInfo->qef);

    midsign = (child->drawInfo->corners >> (7 - i)) & 1; 
    signs[i] = (child->drawInfo->corners >> i) & 1; 

    edgeCount++;
   }
  }
 }

 if (!isCollapsible)
 {
  // at least one child is an internal node, can't collapse
  return node;
 }

For each Internal node in the octree SimplifyOctree is called on the node's children. In the case where a node whose children are all Leaf nodes the QEF values of the Leafs nodes are accumulated to produce a new QEF for the parent node. If all the node's children are Leaf nodes then we can remove all the children and convert the parent node into a Pseudo Leaf Node (i.e. a Leaf node with size > 1). As with the Leaf nodes an OctreeDrawInfo instance is created and the new QEF is used to determine the position of the node's vertex. Once the solve()  function has been called the error value for the new QEF is known. We can use this value to determine whether the proposed new node is to be created or whether it should be discarded.

Once the octree is simplified it can be contoured just as before by calling GenerateMeshFromOctree. This is one of the strengths of the Dual Contouring algorithm, it works on any simplified octree without needing to change the contouring process.

Further Reading

If that all made sense and you want to find out more about Dual Contouring on Distance Functions here are some links:

• Inigo Quilez has a great page about distance functions and their uses (I used this info to construct the scene in the demo)
• The Upvoid dev blog has a nice article about CSG operations using distance fields