I've been trying to implement the Marching Cubes algorithm with limited success. I'm using a 3D array of boolean values as the point data, but I get stuck when I need to interpolate points, I don't really understand what that means or how I should do it. I can loop through the cubes and get the correct edgetable value, but from there on I get confused.

Here is my code so far:


for (int x=0; x < bX; x++) //x axis
{
for (int y=0; y < bY; y++) //y axis
{
for (int z=0; z < bZ; z++) //z axis
{
bool cubeVerts[8];
cubeVerts[0] = points[x][y][z]; // Back top left vert.
cubeVerts[1] = points[x+1][y][z]; // Back top right vert.
cubeVerts[2] = points[x][y+1][z]; // Back bottom left vert.
cubeVerts[3] = points[x+1][y+1][z]; // Back bottom right vert.
cubeVerts[4] = points[x][y][z+1]; // Front top left vert.
cubeVerts[5] = points[x+1][y][z+1]; // Front top right vert.
cubeVerts[6] = points[x][y+1][z+1]; // Front bottom left vert.
cubeVerts[7] = points[x+1][y+1][z+1]; // Front bottom right vert.

int cubeIndex = int(0);
for (int n=0; n < 8; n++)
{
if(cubeVerts[n])
cubeIndex |= (1 << n);
}

if ((cubeIndex != 0) && (cubeIndex != 255))
{
// Do interpolating stuff.
}
}
}
}



Can anybody help me or point me to a resource that explains the process simply? Thanks.

normally marching squares is applied to a real valued function, in which case we can easily interpolate to find a more accurate root along the edges; for instance one vertex f(x,y,z) might be -3 and at the other f(x,y,z) might be 6 so we interpolate and take a point 1/3rd along from the first vertex (and this can be repeated sampling f further).

If the function is simply and 'in-out' valued function the best you can do is to do a binary search for a certain number of iterations and take that as the root of the function where f changes from being 'in' to being 'out'.

All this assumes that the function is in fact continuous and can be sampled at any point in space. If you literally only have values for the vertices of the grid then you can do no interpolation and you are doomed to blocky looking surface extractions



-- edit:

looking at your code you havn't actually got this stage yet. The value you compute based on the in-out at the vertices of the grid cell are normally used to lookup another integer value in a map which defines where we place vertices in the out-polygons; either vertices of the grid cell, or points interpolated (as above) between two vertices along a grid line.

I got it working, thanks for your help luca. This was probably the most helpful resource to understanding the algorithm.




Things still to improve:

1. Interpolate using the Perlin noise value at each point, instead of simply 1 or 0 for whether a point exists there.
2. Use vertex normals instead of face normals.

Here's my code for anyone who is doing something similar.

Typedefs:

struct Vertices
{
D3DXVECTOR3 pos;
float val;
};

struct Triangles
{
D3DXVECTOR3 pos[3];
D3DXVECTOR3 norm[3];
};

int m_NumTriangles;
Triangles m_Triangles[3*bX*bY*bZ];



Helper functions:
D3DXVECTOR3 VertexInterp(D3DXVECTOR3 p1, D3DXVECTOR3 p2, float v1, float v2, float isoLevel)
{
double mu;
D3DXVECTOR3 p;

if (abs(isoLevel-v1) < 0.00001)
return(p1);
if (abs(isoLevel-v2) < 0.00001)
return(p2);
if (abs(v1-v2) < 0.00001)
return(p1);

mu = (isoLevel - v1) / (v2 - v1);
p.x = p1.x + mu * (p2.x - p1.x);
p.y = p1.y + mu * (p2.y - p1.y);
p.z = p1.z + mu * (p2.z - p1.z);

return(p);
}

D3DXVECTOR3 TriangleNormal(D3DXVECTOR3* vVertex1, D3DXVECTOR3* vVertex2, D3DXVECTOR3* vVertex3)
{
D3DXVECTOR3 vNormal;
D3DXVECTOR3 v1;
D3DXVECTOR3 v2;

D3DXVec3Subtract(&v1, vVertex2, vVertex1);
D3DXVec3Subtract(&v2, vVertex3, vVertex1);

D3DXVec3Cross(&vNormal, &v1, &v2);

D3DXVec3Normalize(&vNormal, &vNormal);

return vNormal;
}


Marching Cubes algorithm:
float isoLevel = 0.5f;
m_NumTriangles = 0;
int i;

for (i = 0; i < 3*bX*bY*bZ; i++)
{
for (int j = 0; j < 3; j++)
{
m_Triangles.pos[j].x = 0.0f;
m_Triangles.pos[j].y = 0.0f;
m_Triangles.pos[j].z = 0.0f;
m_Triangles.norm[j].x = 0.0f;
m_Triangles.norm[j].y = 0.0f;
m_Triangles.norm[j].z = 0.0f;
}
}

for (int x=0; x < bX-1; x++) //x axis
{
for (int y=0; y < bY-1; y++) //y axis
{
for (int z=0; z < bZ-1; z++) //z axis
{
Vertices verts[8];
verts[0].val = points[x][y][z]; // Back Bottom Left
verts[0].pos = D3DXVECTOR3(x, y, z); // Back Bottom Left

verts[1].val = points[x+1][y][z]; // Back Bottom Right
verts[1].pos = D3DXVECTOR3(x+1, y, z); // Back Bottom Left

verts[2].val = points[x+1][y][z+1]; // Front Bottom Right
verts[2].pos = D3DXVECTOR3(x+1, y, z+1); // Back Bottom Left

verts[3].val = points[x][y][z+1]; // Front Bottom Left
verts[3].pos = D3DXVECTOR3(x, y, z+1); // Back Bottom Left

verts[4].val = points[x][y+1][z]; // Back Top Left
verts[4].pos = D3DXVECTOR3(x, y+1, z); // Back Bottom Left

verts[5].val = points[x+1][y+1][z]; // Back Top Right
verts[5].pos = D3DXVECTOR3(x+1, y+1, z); // Back Bottom Left

verts[6].val = points[x+1][y+1][z+1]; // Front Top Right
verts[6].pos = D3DXVECTOR3(x+1, y+1, z+1); // Back Bottom Left

verts[7].val = points[x][y+1][z+1]; // Front Top Left
verts[7].pos = D3DXVECTOR3(x, y+1, z+1); // Back Bottom Left

/*
Determine the index into the edge table which
tells us which vertices are inside of the surface
*/
int cubeIndex = int(0);
for (int n=0; n < 8; n++)
{
if(verts[n].val > 0.5f)
cubeIndex |= (1 << n);
}

/* Cube is not entirely in/out of the surface */
if ((cubeIndex != 0) && (cubeIndex != 255))
{
D3DXVECTOR3 vertList[12];

/* Find the vertices where the surface intersects the cube */
if (edgeTable[cubeIndex] & 1)
vertList[0] = VertexInterp(verts[0].pos, verts[1].pos, verts[0].val, verts[1].val, isoLevel);

if (edgeTable[cubeIndex] & 2)
vertList[1] = VertexInterp(verts[1].pos, verts[2].pos, verts[1].val, verts[2].val, isoLevel);

if (edgeTable[cubeIndex] & 4)
vertList[2] = VertexInterp(verts[2].pos, verts[3].pos, verts[2].val, verts[3].val, isoLevel);

if (edgeTable[cubeIndex] & 8)
vertList[3] = VertexInterp(verts[3].pos, verts[0].pos, verts[3].val, verts[0].val, isoLevel);

if (edgeTable[cubeIndex] & 16)
vertList[4] = VertexInterp(verts[4].pos, verts[5].pos, verts[4].val, verts[5].val, isoLevel);

if (edgeTable[cubeIndex] & 32)
vertList[5] = VertexInterp(verts[5].pos, verts[6].pos, verts[5].val, verts[6].val, isoLevel);

if (edgeTable[cubeIndex] & 64)
vertList[6] = VertexInterp(verts[6].pos, verts[7].pos, verts[6].val, verts[7].val, isoLevel);

if (edgeTable[cubeIndex] & 128)
vertList[7] = VertexInterp(verts[7].pos, verts[4].pos, verts[7].val, verts[4].val, isoLevel);

if (edgeTable[cubeIndex] & 256)
vertList[8] = VertexInterp(verts[0].pos, verts[4].pos, verts[0].val, verts[4].val, isoLevel);

if (edgeTable[cubeIndex] & 512)
vertList[9] = VertexInterp(verts[1].pos, verts[5].pos, verts[1].val, verts[5].val, isoLevel);

if (edgeTable[cubeIndex] & 1024)
vertList[10] = VertexInterp(verts[2].pos, verts[6].pos, verts[2].val, verts[6].val, isoLevel);

if (edgeTable[cubeIndex] & 2048)
vertList[11] = VertexInterp(verts[3].pos, verts[7].pos, verts[3].val, verts[7].val, isoLevel);

/* Create the triangle */
for (i = 0; triTable[cubeIndex] != -1; i += 3)
{
D3DXVECTOR3 *v1 = new D3DXVECTOR3(vertList[triTable[cubeIndex][i ]]);
D3DXVECTOR3 *v2 = new D3DXVECTOR3(vertList[triTable[cubeIndex][i+1]]);
D3DXVECTOR3 *v3 = new D3DXVECTOR3(vertList[triTable[cubeIndex][i+2]]);

m_Triangles[m_NumTriangles].pos[0] = *v1;
m_Triangles[m_NumTriangles].pos[1] = *v2;
m_Triangles[m_NumTriangles].pos[2] = *v3;

m_Triangles[m_NumTriangles].norm[0] = TriangleNormal(v1, v2, v3);
m_Triangles[m_NumTriangles].norm[1] = TriangleNormal(v1, v2, v3);
m_Triangles[m_NumTriangles].norm[2] = TriangleNormal(v1, v2, v3);

m_NumTriangles++;
}
}
}
}
}


Rendering:
m_VertexCount = m_NumTriangles*3;
m_IndexCount = m_NumTriangles*3;

// Create the vertex array.
vertices = new VertexType[m_VertexCount];
if(!vertices){ return false; }

m_Indices = new unsigned long[m_IndexCount];
if(!m_Indices){ return false; }

// load the vertex array and index array with data.
int triCounter = 0;
for(i = 0; i < m_VertexCount; i += 3)
{
vertices.position = m_Triangles[triCounter].pos[0]*m_Scale;
vertices.colour = D3DXVECTOR4(1.0f, 1.0f, 1.0f, 1.0f);
vertices.normal = m_Triangles[triCounter].norm[0];

vertices[i+1].position = m_Triangles[triCounter].pos[1]*m_Scale;
vertices[i+1].colour = D3DXVECTOR4(1.0f, 1.0f, 1.0f, 1.0f);
vertices[i+1].normal = m_Triangles[triCounter].norm[1];

vertices[i+2].position = m_Triangles[triCounter].pos[2]*m_Scale;
vertices[i+2].colour = D3DXVECTOR4(1.0f, 1.0f, 1.0f, 1.0f);
vertices[i+2].normal = m_Triangles[triCounter].norm[2];
triCounter++;
}

for(i = 0; i < m_IndexCount; i += 3)
{
m_Indices[i+0] = i+0;
m_Indices[i+1] = i+1;
m_Indices[i+2] = i+2;
}

// The rest is standard D3D rendering stuff.



Feel free to critique, I'm still quite inexperienced. Also, if anyone could advise on the best way to store the Triangles it would be much appreciated, because if I increase the size of the model much more the compiler says the array is too big.