How to construct the marching cubes tables algorithmically?

Question

I want to implement the marching cubes algorithm from scratch, but I'm stuck at the polygon generation phase (building edge loops with correct orientation and triangulating them). Obviously, I should consider all 2^8 = 256 cases, find out zero-crossing edges (build 12-bit edge masks) and figure out the triangulation in each case.

Could you please give me some ideas on how to correctly construct and triangulate edge loops?

I'd like to avoid considering each MC configuration (out of 15) separately as it seems somewhat tedious and messy to me.

Answer 1

I also had exactly the same thoughts as you. For my final university project I studied different methods of voxel mesh smoothing. The best method I found was Surface Nets. It produces a result that looks very similar to March Cubes but without all that lookup table hassle. You can also choose how smooth you want the object by performing more passes of the algorithm.

The basic procedure is to do what Marching Cubes does to start with:

• Find all the points of voxel volume that intersect the current voxel cube
• Interpolate to get the edge of that volume

You then take these interpolated points and find their 'centre of mass', this gives you a point to place a vertex at. From this vertex position you connect outwards to other vertices in neighbouring voxels.

I currently do not have access to any of my source code to help expand on this solution but I will edit this answer as soon as possible to make it clearer and to support the article I've referenced.

---

Edit:

There's going to be a lot of code splurged below, followed by an open source licence, feel free to ask me more questions!

So here is my GetCentreOfMass function, I'll try to break down in comments:

Get Centre Of Mass Method

/// <summary>
/// Gets the centre of mass of a voxel for use in the SurfaceNets mesh generation     
/// algorithm.
/// </summary>
/// <param name="voxels">The voxels.</param>
/// <returns></returns>
/// <exception cref="Project.VoxelEngine.Visualisation.VoxelNotOnSurfaceException">Voxels not on surface\nvoxels before:  + voxelStringBefore + \nedge positions:  + edgePositionsString + \nvoxels after:  + voxelStringAfter + , edgePosition count:  + edgePositions.Count</exception>
public Vector3 GetCentreOfMass(Voxel[] voxels)
{
    string voxelStringBefore = voxels.ElementsToString(", ");
    string edgePositionsString = "";
    List<Vector3> edgePositions = new List<Vector3>();

    // For each edge definition: here an edge definition is a relationship between two of 
    // the vertices of a voxel, so it will contain 2 numbers which are an indices that             
    // point to voxels relative to the current voxel we want to find the centre of mass 
    // for. The index 0 would be this current voxel, while the index 1 means the voxel in 
    // front of the current one in the x direction, 2 would mean forward in the y 
    // direction etc
    for (int i = 0; i < edgeDefinitions.Length; i++)
    {
        Voxel voxelA = voxels[(int)edgeDefinitions[i].x];
        Voxel voxelB = voxels[(int)edgeDefinitions[i].y];


        // Evaluate whether or not the surface crosses the edge,
        // here we want to work out whether or not the density of the voxels at either end
        // of the current edge definition changes from positive to 0 across the edge, this
        // is similar to part of the marching cubes process

        if(voxelA.GetDensity() == 0 && voxelB.GetDensity() != 0 ||
            voxelA.GetDensity() != 0 && voxelB.GetDensity() == 0)
        {

            // Interpolate the position that the surface crosses this edge, we    
            // mathematically find a fast and rough guess for where the 'physical' surface
            // of the voxel volume lines on this edge definition. This will be where one
            // of the faces of our mesh will cut across between two voxels.
            // Once we've found a value we add it the edge positions list, a collection
            // of positions where the surface intersects with the boundaries of this 
            // current voxel and its neighbours.

            int difference = voxelB.GetDensity() - voxelA.GetDensity();
            float lerp = (1f / Voxel.maxDensity) * Mathf.Abs(difference);
            if (difference < 0)
            {
                // In this part VoxelGridData.neightbourVectors is a structure that
                // I have in a static class that gives me the physical vector of a  
                // neighbour based on the index we have from the edge definition (bit
                // convoluted I know)
                edgePositions.Add(Vector3.Lerp(
                    VoxelGridData.neighbourVectors[(int)edgeDefinitions[i].x],
                    VoxelGridData.neighbourVectors[(int)edgeDefinitions[i].y],
                    lerp));
            }
            else
            {
                edgePositions.Add(Vector3.Lerp(
                    VoxelGridData.neighbourVectors[(int)edgeDefinitions[i].y],
                    VoxelGridData.neighbourVectors[(int)edgeDefinitions[i].x],
                    lerp));
            }

            // Debugging stuff :P
            edgePositionsString += "\nEdge Definition: " + i + ", Edge Position " + (edgePositions.Count-1) + ": " + edgePositions[edgePositions.Count-1] + ", maxDensity: " + Voxel.maxDensity + ", difference: " + difference + ", lerp: " + lerp;
        }
    }

    // If we found some intersections, i.e. the voxel actually sits on the surface, if 
    // it doesn't my algorithms that call this method have actually failed so I throw an
    // exception
    if (edgePositions.Count > 0)
    {
        // Find the centre of mass of edge intersections with the surface using simple                
        // maths!
        Vector3 centreOfMass = new Vector3();
        centreOfMass.x = edgePositions.Sum(i => i.x) / edgePositions.Count;
        centreOfMass.y = edgePositions.Sum(i => i.y) / edgePositions.Count;
        centreOfMass.z = edgePositions.Sum(i => i.z) / edgePositions.Count;
        return centreOfMass * this.voxelData.GetScale();
    }
    else
    {
        string voxelStringAfter = voxels.ElementsToString(", ");
        throw new VoxelNotOnSurfaceException("Voxels not on surface\nvoxels before: " + voxelStringBefore + "\nedge positions: " + edgePositionsString + "\nvoxels after: " + voxelStringAfter + ", edgePosition count: " + edgePositions.Count);
    }
}

This is pretty hefty and it's only getting the centre of mass but it does the job well. To fully build the mesh I store this centre of mass and an edge mask (for each voxel 8 bits: 1 for an edge on the surface and 0 for an edge not on the surface) together in a structure called a NetVertex. I then use these to create the mesh, the vertices are easy:

Mesh Vertex Generation Method

foreach (KeyValuePair<Vector3, NetVertex> netVertex in this.netVertices)
{
    vertices[i] = netVertex.Value.GetVertexPosition(); // Set the actual vertex from the netVertex
    netVertex.Value.SetVertexIndex(i); // Store which actual vertex relates to this NetVertex
    i += 1;
}

Triangles of the mesh are somewhat more complicated:

Mesh Triangle Generation Method

// For every net vertex
foreach (KeyValuePair<Vector3, NetVertex> netVertex in this.netVertices)
{
    // Get the edgeMask of the net vertex.
    int[] edgeMask = netVertex.Value.GetEdgeMask();
    // For each edge that can have a triangle crossing it
    for (int j = 0; j < triangleEdges.Length; j++)
    {
        // If the edgeMask says that the surface crosses this edge
        if (edgeMask[triangleEdges[j]] != 0)
        {
            // Get the vectors of the two new points of the triangle
            Vector2 triangleVectorIndices = triangleVectors[(j % 3)];
            Vector3 origin = netVertex.Key;

            // Find out the triangle facing direction
            int sign = 0;
            if (j <= 2) sign = -1;
            else if (j >= 3) sign = 1;

            // Get the two net vertex positions of the other points of the triangles
            Vector3 netVertexPosition1 = origin + (neighbourVectors[(int)triangleVectorIndices.x] * sign);
            Vector3 netVertexPosition2 = origin + (neighbourVectors[(int)triangleVectorIndices.y] * sign);

            // If these net vertices exist
            if (this.netVertices.ContainsKey(netVertexPosition1) && this.netVertices.ContainsKey(netVertexPosition2))
            {
                // Add the three vertex indicies to build the triangle, taking care of facing direction
                triangles.Add(netVertex.Value.GetVertexIndex());
                if (edgeMask[triangleEdges[j]] == 1)
                {
                    triangles.Add(this.netVertices[netVertexPosition1].GetVertexIndex());
                    triangles.Add(this.netVertices[netVertexPosition2].GetVertexIndex());
                }
                else if (edgeMask[triangleEdges[j]] == -1)
                {
                    triangles.Add(this.netVertices[netVertexPosition2].GetVertexIndex());
                    triangles.Add(this.netVertices[netVertexPosition1].GetVertexIndex());
                }
            }
        }
    }
}

---

All the source code I've added to this post follows the standard MIT open source licence:

The MIT License (MIT)

Copyright (c) 2015 Matthew Torr

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.

Answer 2

Disclaimer: Not tested yet. Edit: Tested and works, code on demand.

The configuration, or rather the mesh, can be calculated in a slightly simpler, in my opinion, way, although this is less efficient, but can be memoized, I think.

Lets call the nodes that have the active/inactive state as "ControlNode", and the 3 nodes that each of those hold "SubNode".
So each Cube holds 8 ControlNodes and 12 SubNodes.

The Cubes need to store this data in such a way that you know which SubNodes are neighboring which ControlNodes. Remember that the ControlNodes themselves have 3 SubNodes where one or more usually land outside the Cube, we don't want those.

Now to get the vertices that we need to build the mesh for this Cube:
1) Create an array for the states of the SubNodes. Yes SubNodes.
2) Iterate over the active ControlNodes, turning "on" the adjacent SubNodes, or turning them "off" if they were already "on".
3) The SubNodes that are "on", after step 2, are the vertices you need to use to build the mesh.
4) Alas, we don't know the direction of the surface, but we can calculate it! The Normal can be calculated by simply calculating the average positions of the selected SubNodes (A) and active ControlNodes (B).
Vector AB should be the Normal.
5) Having the vertices, and the normal, you can construct triangles from the data collected.

Extra: If you want to make this work faster, save the configuration as a number, as in the original algorithm, and store the results for each configuration as you find them.

Edit: A simple iteration will not work for cases where the configuration creates 2 separate meshes!
A solution to that is placing the active nodes in a list, and running a "paint bucket" style search, such as BFS, to get a cluster of active ControlNodes that are all "connected" (meaning forming one mesh), removing the affected ControlNodes from the list. And if any ControlNodes remain, then we do the same for that group as well.
Adjacency for the BFS can be determined by proximity, since the distance between neighboring ControlNodes is known.

Not as simple as I first thought, but it should provide a general case solution.
Hope that answers the question.