How could I polygonise a bool[,,]

Question

If anyone cares, I'm using WPF....

Beginning of story:

I got a series of gray-scale images(slices of a model). The user inputs a "range" of gray-scale values to build a 3D model upon. Thus I created a 3D bool array to make it easier to build the 3D model. This array represents a box of pixels, indicating whether each pixel should be built/generated/filled or not.

---

The tie:

Using a bool[,,], I want to retrieve a List<Point3D[]> where each Point3D[] has a length of 3 and represents a triangle in 3D space.

---

Additional info:

The generated model is going to be 3D printed. In the bool[,,], true indicates the presence of matter, while false indicates the absence of matter. I need to avoid cubic models, where each pixel is replaced with a cube(That would be useless according to my purpose). The model should be as smooth as possible and as accurate as possible.

---

What I tried to do:

1- I implemented the marching cubes algorithm, but it simply doesn't seem to be created to accept a "range" of values.

2- I kept struggling for a couple of days trying to make my own algorithm, but I partially failed. (It's real complicated. If you want to know more about it, please inform)

---

Some expected solutions that I have no idea how to implement:

1- Modify the Marching cubes algorithm in order to make it work using a bool[,,]

2- Modify the Marching cubes algorithm in order to make it work using a work using a "range" of isolevel values

3- Showing how to implement another algorithm that is suitable for my purpose (Maybe the Ray-Casting algorithm) in WPF.

4- Requesting the source of the algorithm I tried to implement then showing me how to solve it.(It was made to polygonise a bool[,,] in the first place)

5- Some other magical solution.

Thanks in advance.

---

EDIT:

By saying

Using a bool[,,], I want to retrieve a List<Point3D[]> where each Point3D[] has a length of 3 and represents a triangle in 3D space.

I mean that I want to retrieve a group of Triangles. Each triangle should be represented as 3 Point3Ds. If you don't know what a Point3Dis, it is a struct containing 3 doubles (X,Y and Z) that's used to represent a position in 3D space.

The problem with the marching cubes algorithm is that it's kinda vague. I mean, what do you understand by doing that??

cubeindex = 0; 
if (grid.val[0] < isolevel) cubeindex |= 1; 
if (grid.val[1] < isolevel) cubeindex |= 2; 
if (grid.val[2] < isolevel) cubeindex |= 4; 
if (grid.val[3] < isolevel) cubeindex |= 8; 
if (grid.val[4] < isolevel) cubeindex |= 16; 
if (grid.val[5] < isolevel) cubeindex |= 32; 
if (grid.val[6] < isolevel) cubeindex |= 64; 
if (grid.val[7] < isolevel) cubeindex |= 128; 

Thus, I simply have no idea where to start modifying it.

---

Another EDIT:

After some experimenting with the marching cubes algorithm, I noticed that polygonising a bool[,,] is not going to result in the smooth 3D model I'm looking forward to, so my first expected solution and my fourth expected solution are both out of the play. After a lot of research, I knew about the Ray-Casting method of volume rendering. It seems really cool, but I'm not really sure if it fits to my needs. I can't really understand how it is implemented although I know how it works.

---

One more EDIT: TriTable.LookupTable and EdgeTable.LookupTable are the tables present here

Here's the MarchingCubes class:

public class MarchingCubes
{
    public static Point3D VertexInterp(double isolevel, Point3D p1, Point3D p2, double valp1, double valp2)
    {
        double mu;
        Point3D p = new Point3D();

        if (Math.Abs(isolevel - valp1) < 0.00001)
            return (p1);

        if (Math.Abs(isolevel - valp2) < 0.00001)
            return (p2);

        if (Math.Abs(valp1 - valp2) < 0.00001)
            return (p1);

        mu = (isolevel - valp1) / (valp2 - valp1);

        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);
    }

    public static void Polygonise (ref List<Point3D[]> Triangles, int isolevel, GridCell gridcell)
    {
        int cubeindex = 0;
        if (grid.val[0] < isolevel) cubeindex |= 1;
        if (grid.val[1] < isolevel) cubeindex |= 2;
        if (grid.val[2] < isolevel) cubeindex |= 4;
        if (grid.val[3] < isolevel) cubeindex |= 8;
        if (grid.val[4] < isolevel) cubeindex |= 16;
        if (grid.val[5] < isolevel) cubeindex |= 32;
        if (grid.val[6] < isolevel) cubeindex |= 64;
        if (grid.val[7] < isolevel) cubeindex |= 128;

        // Cube is entirely in/out of the surface 
        if (EdgeTable.LookupTable[cubeindex] == 0)
            return;

        Point3D[] vertlist = new Point3D[12];

        // Find the vertices where the surface intersects the cube 
        if ((EdgeTable.LookupTable[cubeindex] & 1) > 0)
            vertlist[0] = VertexInterp(isolevel, grid.p[0], grid.p[1], grid.val[0], grid.val[1]);

        if ((EdgeTable.LookupTable[cubeindex] & 2) > 0)
            vertlist[1] = VertexInterp(isolevel, grid.p[1], grid.p[2], grid.val[1], grid.val[2]);

        if ((EdgeTable.LookupTable[cubeindex] & 4) > 0)
            vertlist[2] = VertexInterp(isolevel, grid.p[2], grid.p[3], grid.val[2], grid.val[3]);

        if ((EdgeTable.LookupTable[cubeindex] & 8) > 0)
            vertlist[3] = VertexInterp(isolevel, grid.p[3], grid.p[0], grid.val[3], grid.val[0]);

        if ((EdgeTable.LookupTable[cubeindex] & 16) > 0)
            vertlist[4] = VertexInterp(isolevel, grid.p[4], grid.p[5], grid.val[4], grid.val[5]);

        if ((EdgeTable.LookupTable[cubeindex] & 32) > 0)
            vertlist[5] = VertexInterp(isolevel, grid.p[5], grid.p[6], grid.val[5], grid.val[6]);

        if ((EdgeTable.LookupTable[cubeindex] & 64) > 0)
            vertlist[6] = VertexInterp(isolevel, grid.p[6], grid.p[7], grid.val[6], grid.val[7]);

        if ((EdgeTable.LookupTable[cubeindex] & 128) > 0)
            vertlist[7] = VertexInterp(isolevel, grid.p[7], grid.p[4], grid.val[7], grid.val[4]);

        if ((EdgeTable.LookupTable[cubeindex] & 256) > 0)
            vertlist[8] = VertexInterp(isolevel, grid.p[0], grid.p[4], grid.val[0], grid.val[4]);

        if ((EdgeTable.LookupTable[cubeindex] & 512) > 0)
            vertlist[9] = VertexInterp(isolevel, grid.p[1], grid.p[5], grid.val[1], grid.val[5]);

        if ((EdgeTable.LookupTable[cubeindex] & 1024) > 0)
            vertlist[10] = VertexInterp(isolevel, grid.p[2], grid.p[6], grid.val[2], grid.val[6]);

        if ((EdgeTable.LookupTable[cubeindex] & 2048) > 0)
            vertlist[11] = VertexInterp(isolevel, grid.p[3], grid.p[7], grid.val[3], grid.val[7]);

        // Create the triangle 
        for (int i = 0; TriTable.LookupTable[cubeindex, i] != -1; i += 3)
        {
            Point3D[] aTriangle = new Point3D[3];

            aTriangle[0] = vertlist[TriTable.LookupTable[cubeindex, i]];
            aTriangle[1] = vertlist[TriTable.LookupTable[cubeindex, i + 1]];
            aTriangle[2] = vertlist[TriTable.LookupTable[cubeindex, i + 2]];

            theTriangleList.Add(aTriangle);
        }
    }
}

Here's the GridCell class:

public class GridCell
{
    public Point3D[] p = new Point3D[8];
    public Int32[] val = new Int32[8];
}

Ok, that's what I'm doing:

List<int[][]> Slices = new List<int[][]> (); //Each slice is a 2D jagged array. This is a list of slices, each slice is a value between about -1000 to 2300

public static void Main ()
{
    List <Point3D[]> Triangles = new List<Point3D[]> ();
    int RowCount;//That is my row count
    int ColumnCount;//That is my column count
    //Here I fill the List with my values
    //Now I'll polygonise each GridCell
    for (int i = 0; i < Slices.Count - 1; i++)
    {
        int[][] Slice1 = Slices[i];
        int[][] Slice2 = Slices[i + 1];
        for (int j = 0; j < RowCount - 1; j++)
        {
            for (int k = 0; k < ColumnCount - 1; k++)
            {
                GridCell currentCell = GetCurrentCell (Slice1, Slice2, j, k);
                Polygonise (ref Triangles, int isoLevel, GridCell currentCell);
            }
        }
    }
    //I got the "Triangles" :D
}

//What this simply does is that it returns in 3D space from RI and CI
public static GridCell GetCurrentCell (int[][] CTSliceFront, int[][] CTSliceBack, int RI, int CI)//RI is RowIndex and CI is ColumnIndex
{
    //Those are preset indicating X,Y or Z coordinates of points/edges
    double X_Left_Front;
    double X_Right_Front;
    double X_Left_Back;
    double X_Right_Back;
    double Y_Top_Front;
    double Y_Botton_Front;
    double Y_Top_Back;
    double Y_Botton_Back;
    double Z_Front;
    double Z_Back;

    GridCell currentCell = new GridCell();
    currentCell.p[0] = new Point3D(X_Left_Back, Y_Botton_Back, Z_Back);
    currentCell.p[1] = new Point3D(X_Right_Back, Y_Botton_Back, Z_Back);
    currentCell.p[2] = new Point3D(X_Right_Front, Y_Botton_Front, Z_Front);
    currentCell.p[3] = new Point3D(X_Left_Front, Y_Botton_Front, Z_Front);
    currentCell.p[4] = new Point3D(X_Left_Back, Y_Top_Back, Z_Back);
    currentCell.p[5] = new Point3D(X_Right_Back, Y_Top_Back, Z_Back);
    currentCell.p[6] = new Point3D(X_Right_Front, Y_Top_Front, Z_Front);
    currentCell.p[7] = new Point3D(X_Left_Front, Y_Top_Front, Z_Front);

    currentCell.val[] = CTSliceBack[theRowIndex + 1][theColumnIndex];
    currentCell.val[] = CTSliceBack[theRowIndex + 1][theColumnIndex + 1];
    currentCell.val[] = CTSliceFront[theRowIndex + 1][theColumnIndex + 1];
    currentCell.val[] = CTSliceFront[theRowIndex][theColumnIndex];
    currentCell.val[] = CTSliceBack[theRowIndex][theColumnIndex + 1];
    currentCell.val[] = CTSliceFront[theRowIndex][theColumnIndex + 1];
    currentCell.val[] = CTSliceFront[theRowIndex][theColumnIndex];

    return currentCell;
}

I made that right from Stack-Overflow, so please point out any typos. That is working, This results in a shell of a model. I want to do is replace the isolevel variable in the Polygonise function with 2 integers: isolevelMin and isolevelMax. Not only 1 int, but a range.

---

EDIT: Guys, please help. 50 rep is almost 1/2 my reputation. I hope I can change my expectations. Now I just want any idea how I can implement what I want, any snippet on how to do it or if anyone gives out some keywords that I can google which I used to solve my problem he'll get the rep. I just need some light - it's all dark in here.

---

EDITs are awesome: After even more, more and much more research, I found out that the volumetric ray-marching algorithm does not actually generate surfaces. Since I need to print the generated 3D model, I only have one way out now. I gotta make the marching cubes algorithm generate a hollow model from the maximum and minimum isolevel values.

Answer 1

One of the assumptions of Marching Cube algorithm is that every cube geometry is a separate entity and doesn't depend on other cubes (which is not completely true but let's stick to this for now).

Given that, you can split entire mesh into cubes and solve each one separately. You may also notice, that although the final mesh depends on exact values of your 3D scalar field, the topology of the mesh depends only on fact if some given cube corner is inside of a mesh or outside.

For example: Let's say your isolevel is set to 0 and your 3D field has positive and negative numbers. If you now change some value from 0.1 to 1.0 then you'll change only proportions of some triangles but not their number and/or topology. On the other hand, if you change some value from 0.1 to -0.1 then your triangles may switch to another arrangement and their number may change.

Thanks to this we can use some clever optimization - for each of 8 corners you check if the corner is inside or outside of a mesh, you use those 8 bits to build a number - which becomes an index to your precomputed table of possible cube solutions, let's call them Variants.

Each Cube Variant is a list of triangles for this specific configuration of corners. Triangles in Marching Cubes always span between cube edges, so we use edge indices to describe it. Those indices are the exact data you see in those 'magic' tables in Paul Bourke code :).

But this is not the final mesh yet. The triangles you get here are based just on a fact if some corner is inside or outside of a mesh, but how far is it from the surface? Does it have any influence on the result? Here comes your grid values once again - when you've selected your Variant, got the list of triangles and 3 edges for each of those triangles, then you get the values on ends of the edge and interpolate them to find 0 value (we've set it as isolevel, remember?). Those final, interpolated vertices should be used for your mesh.

That became quite a long introduction, I hope you understand what I wrote.

My suggestion: The simplest change is with corners, instead of doing this:

if (grid.val[0] < isolevel) cubeindex |= 1;

Try to change it this way:

if (grid.val[0] > isolevelMin && grid.val[0] < isolevelMax) cubeindex |= 1;
// And the same for other lines...

Final interpolation is a bit trickier and I have no way to test it but let's try:

• Modify your VertexInterp to pass two isolevels - min and max
• Inside VertexInterp check which isolevel is between valp1 and valp2 values
• Interpolate as in the current code, but using selected isolevel (min or max)

And let me know if it worked :) or if you have any questions.