## Sunday, 9 November 2014

### 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;
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];
const int c2 = edgevmap[i];

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, int dir, IndexBuffer& indexBuffer)
{
int minSize = 1000000;  // arbitrary big number
int minIndex = 0;
int indices = { -1, -1, -1, -1 };
bool flip = false;
bool signChange = { false, false, false, false };

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

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);
indexBuffer.push_back(indices);
indexBuffer.push_back(indices);

indexBuffer.push_back(indices);
indexBuffer.push_back(indices);
indexBuffer.push_back(indices);
}
else
{
indexBuffer.push_back(indices);
indexBuffer.push_back(indices);
indexBuffer.push_back(indices);

indexBuffer.push_back(indices);
indexBuffer.push_back(indices);
indexBuffer.push_back(indices);
}
}
}
```

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, GL_STATIC_DRAW);

glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, indexBuffer_);
glBufferData(GL_ELEMENT_ARRAY_BUFFER, sizeof(int) * indices.size(), &indices, 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 = { -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:

#### 42 comments:

1. Very helpful! Though my implementation is very slow :/ Could you do a tutorial on how to effectively implement this, maybe even on a gpu? That would be awesome :)

2. Glad you found it useful :)

Yes, this implementation is very naive so it will be quite slow. The simplest way to speed it up would be to cache some of the density calculations that use the noise function.

I have some vague plans to talk about my current OpenCL implementation, I dunno when I'll get around to writing it though.

3. This comment has been removed by the author.

4. Look for the GenerateVertexIndices (or something similar) function. As each leaf node is processed its data is added to the vertices and the previous size is used for the index. E.g. it'll be something like drawinfo->index = vertexBuffer.size(); vertexBuffer.push_back(...);

1. i found the problem, seems to be a typical syntax error from my side. you used a holder variable d for drawinfo and put the indexes inside, which doesnt change the value of the original drawinfo in c# .... but now after all it seems for my implementation that at least the Y coords of the normals are - when they should be + and the other way around. is this the same in your code or did translating this to c# messed sth up? because i didnt realy changed some of the math parts

2. Hm, hard to say without seeing the code. One possible cause of that would be if you defined 'solid' as +ve and air as -ve which is the opposite to my code.

3. yeah that seems to be the issue, i still dont understand why everybody defines solid as negative in voxel worlds xD for me solid means sth like "there is sth" which is not very intuitive to negative numbers ^^ but great tutorial anyways, helped me alot getting dual contouring to unity engine. now i have to face the chunk transition problem too :-/

4. That's great you've got it working, good job.

Yeah that confused me too, wince with marching cubes or Minecraft style its a sort of boolean, 0 or not 0. The reason inside is negative is we're now using a distance field, so its "how far away is the surface?"

5. This comment has been removed by the author.

6. This comment has been removed by the author.

7. I'm using this implementation and I find it very slow, how can I speed it up? I've made some tests and I have noticed that the heaviest parts are that of the QEF stuff...

8. That's a good question! You can see if you can find a better QEF implementation or come up with a different scheme for positioning the vertices, e.g. if you search for Leonardo Schmitz Dual Contouring on the GPU dissertation he has an interesting method. The other slow part is the noise function. You can try caching some of these values or do what I have done which is use OpenCL.

1. The OpenCL solution seems pretty smart, the unluck is that I am newbie with that tecnology. Do you have very important improvments from the CPU implementation? I believe yes watching your videos... the engine runs pretty good while my implementation freeze at generating only one 32x32x32 octree... Have you got some tips to begin coding with OpenCL? Maybe some implementation that I could look at... As usually thanks for your help! I am really learning stuff here! :)

2. I have the 'voxel discovery' phase in OpenCL currently, which means I can specify a min position for a chunk and the OpenCL code will examine every voxel in that volume and return all the voxels that contain the surface. That is one of the main bottle necks and fits with OpenCL very well. Once I have the info for each active voxel (e.g. how many edge crossing, position and normal for each edge crossing...) I can then construct the octree upwards from the voxels (which are the leaf nodes of the octree).

A good place to start (as always!) is the ProcWorld blog: http://procworld.blogspot.co.uk/2010/11/opencl-first-iteration.html

There's also a Marching Cubes sample in either the AMD or NVIDIA OpenCL samples which was a big help.

This is the first time I've used OpenCL too so it took a lot longer than I had hoped/expected to get my first version working but I think its worth it in the long term. Good luck!

3. So you basically build the octrees voxels in the openCL and then you build the octree using the CPU and that runs fine? This is very cool because doing everything on the CPU isn't going fluidly, it takes to much time to build one single (32 sized) octree.

4. Yes that's right. The QEFs are the next problem :)

For reference my chunks are now 128^3 since the OpenCL impl is so much faster.

5. what do you think about building the entire octree with openCL instead of building only the leaves? would it be faster? or maybe it's not a good idea?

6. That bit is not particularly slow, and as I understand it tree traversals don't map well to GPU coding.

9. Hi guys, can anyone give me a hint on why my implementation (c#) is building vertex and tris and feeds them into the arrays but no mesh is created?

1. Found it it was the OctreeDrawInfo d in GenerateVertexIndices(...).
Just use this:
node.drawInfo.index = meshData.vertices.Count;
meshData.AddVertex(node.drawInfo.position);
meshData.AddNormal(node.drawInfo.averageNormal);

10. Hi Andrea. The could be any number of reasons, you'll need to narrow the scope a bit!

When you say "no mesh is created" do you mean that nothing is drawn? If so you could try just setting the mesh to draw one triangle and make sure you can see that. If that works your DC mesh might be not be visible if the triangles are drawing back to front so you could flip the winding order. I could go on... :) Another thing to try is using a really simple distance field, e.g. I use a flat plane like float density(vec3& pos) { return pos.y - 10.f; } to debug stuff quite often.

1. The problem was OctreeDrawInfo d = node.drawInfo; in void GenerateVertexIndices(OctreeNode node, MeshData meshData). I had to use node.drawInfo (that was due to a blind conversion from C++ to C#), anyway thanks for the plane density, it may come in handy.

11. Hi Andrea. The could be any number of reasons, you'll need to narrow the scope a bit!

When you say "no mesh is created" do you mean that nothing is drawn? If so you could try just setting the mesh to draw one triangle and make sure you can see that. If that works your DC mesh might be not be visible if the triangles are drawing back to front so you could flip the winding order. I could go on... :) Another thing to try is using a really simple distance field, e.g. I use a flat plane like float density(vec3& pos) { return pos.y - 10.f; } to debug stuff quite often.

1. Hi Nick! First of all thanks for your dedication to this article, very useful to understand the mechanics under this algorithm!
I am working together with Andrea and as you noticed we are struggling with the mesh generation part.
We managed to build the octree and extract meshes from it, but apparently something isnt working properly.
While the octree seems fine, the meshes look twice the size they should be, and deformed around the edges. I wrote simple "debug" function to draw the borders of the octree and this is the result:
http://oi57.tinypic.com/jgqxb9.jpg
http://i58.tinypic.com/1eq8na.png
I am almost sure this has something to do with QEF calculations, however that part wasn't the easiest one, and Im afraid I didnt fully understand its logic, or I didn't consider some pointers by mistake passing from c++ to c#.
Do you (or anyone else) have any suggestion on how to fix this? Thanks a lot!

2. Funny when I find the problems by myself, especially when they are stupid! :)
For anyone using a language different from C++, NEVER EVER use external functions to solve QEF, especially for vectors.
It appears that using Unity's built-in class Vector3 and relative functions, instead of Vec3 used in the svd lib, cause some problems with scaling and normalizing.
It is possible to use Vector3 without problems (I kept it), but be sure to use the functions provided in VecUtils.

Anyway back on topic: I read in the comments about optimization. Now, this actually works reasonably fast for a 16 or 32 sized chunk, but a little improvement always helps. You suggested some caching, which data is worth "chaching" somewhere in your opinion?
Thanks again ^^

3. Hi Edoardo, glad you figured it out :)

One possible caching strategy is (assuming you're using a heightmap...) is to store the heightmap values for the chunk (so a 2D array) and use those to calculate the 3D density values. Something like:

float density(vec3 position)
{
return heightmap.value(position.x, position.z) - position.y;
}

That's a good place to start.

More complicated caching systems might not neccesarily provide much speed up, and might compilcate your code. Bear in my mind if you do lots of caching that memory accesses are not free -- you could end up trading CPU calculations for trips to main RAM or even paging, so make sure you can test that your caching strategy is actually effective!

Good luck :)

4. Saving the noise functions for the future could save a lot of time indeed! Right now I just have simple test functions like a plane (as you suggested, pretty useful to debug stuff), or a sphere with some simplex noise added to it to break the regular shape a bit...
Building a "good looking" terrain I guess it will need a much more complex function with different parameters, using a heightmap would speed it up a lot.
I was thinking about some caching for the QEF solve method, but I was still surprised by the overall speed. Let's think about optimization after filling the cracks between the chunks :P
Thanks as always!

12. Hey Nick,

Just wanted to give you a huge thank here for sharing this. I'm a newbie in this voxel tech and i needed something to help me with a first step. I adapted your code to Unreal Engine and threw in a 3D perlin noise:
https://www.youtube.com/watch?v=SB5uk-4JvL8
Obviously lots of work ahead but it's a promising first step, all thanks to you !
Cheers
Cedric

1. Hi. Cool stuff, good job getting it working in Unreal. Glad I could help :)

2. Hey, what was the difficulty curve in getting this to work in Unreal? I'm trying to do something similar, but am unsure where to start.

13. Hey Andrew,

Sorry for the delayed answer, i wasn't noticed so i saw your post only when coming back to read more info on this -amazing- blog :-)

I have (heavily modified then) deleted the code i used back then, only keeping the qef part and trying to rewrite the rest by myself (it currently works but 10 times slower than Nick's code and with nasty glitches !!)

From what i recall it was a fairly easy process, basically change the vec3 into FVector, use TARRAYs for every array little things like that.

The rendering part was made with one UProceduralMeshComponent for the entire Octree, using one call to the CreateMeshSection method (hence the mandatory use of TARRAYs), with a custom UE4 material.

The time it took is difficult to estimate, a few days working a few hours per day, so maybe around 20 hours, more or less a lot :-))

If you don't know how to start, this could give you a starting point:
1/ learn to draw one simple mesh with UProceduralMeshComponent::CreateMeshSection
2/ try to compile successfully Nick's code in UE4
3/ connect 1 and 2 by feeding CreateMeshSection with the Arrays provided by Nick's code.

Cheers

Cedric

14. This comment has been removed by the author.

15. Would changing to an intersection free implementation require a refactoring of the octree? I'm noticing some cracks in my mesh that I would like to be rid of. Also do you have any pointers on where to start looking for a clear grasp of the no-intersections method?

1. Hi. I've not really looked at the no intersection version myself. I did find some code online at some point, and it looked like it was doing quite a lot of work so I didn't bother with it. Instead I use a mesh simplification step after generating the mesh with dual contouring (so I never need to collapse nodes, etc).

2. Alright thanks for the quick answer!

16. Hi Nick! Thanks for sharing you code!
I was particularly excited about the ability of representing sharp edges they claim in the paper.

Fist test is render a cube: perfect! Just 2 polygons per face, perfectly cubic.

Second test is render a rotated cube...
I apply the rotation matrices to worldPosition in the Density_Func but the result is not what I was hoping for:
https://ibb.co/gZkg1Q

Do you think this is a flaw in the dual contouring algorithm or there is something wrong with the implementation?

1. Hi Emanuele,

I'm glad you've found it useful :)

If you are using the DualContouringSample then I think the problem is the QEF impl included in there has some problems. You could try taking the glsl_svd.cpp from my "qef" repo and using that instead, it should fix your problem.

17. Hi Nick thank you for your hard work!
I'm trying to implement it in UE4 but i have a problem with my triangle array (indexbuffer).
In ue4 triangle arrays work by 3 and reference the vertexbuffer array like this
//1 triangle
tri=0; p0
tri=1; p1
tri=2; p2

for (int32 Index = 0; Index != ArrIndices.Num(); Index = Index + 6)
{
Triangles.Add(ArrIndices[Index + 0] + Index/6*4);
Triangles.Add(ArrIndices[Index + 1] + Index/6*4);
Triangles.Add(ArrIndices[Index + 2] + Index/6*4);
Triangles.Add(ArrIndices[Index + 3] + Index/6*4);
Triangles.Add(ArrIndices[Index + 4] + Index/6*4);
Triangles.Add(ArrIndices[Index + 5] + Index/6*4);
}
As you can see I increment ArrIndices(the IndexBuffer) because they are the values inside ContourProcessEdge and to build my tri i must referre all vertices.
It make me crazy :p
I can give you my work to have a look
http://imgur.com/y7l1nb9

1. Hi. I'm not at all familiar with UE but I strongly suspect you are not building the mesh correctly. I would suggest looking for a tutorial which builds a mesh in the same way (without using DC) and then look to see if you can see any problems with your code. (Off the top of my head it seems you are perhaps mixing up the positions with indices, and the "Index/6*4" looks wrong.

2. I compiled your code it run without problem and I looked at the array of indexing triangles
I noticed that the error did not come from my way of making the triangles but from an error I made elsewhere that cause the indexbuffer to stay at value between 0 to 3.
I need to do further investigation. thank you nick.
And yes i was wrong with what i've done i just have to put your indexbuffer but i need to find my error :p

18. How can you implement voxel data, i just can't find it?

19. Very helpful. Thanks for your work!
And i have a question: what is "MAX_CROSSINGS = 6" mean in octree.cpp and why is 6?