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.

Angi

"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": 1. On a legal pad, write out in longhand precisely what you mean by that, your guess being certainly a great deal better than mine. 2. Refine that until it looks like pseudocode. 3. Easy part: Translate that into real code. 4. Ask for help if its borken. – Sharecrop

"I implemented the marching cubes algorithm, but it simply doesn't seem to be created to accept a "range" of values": When a question says "doesn't seem", I die a little. The information in that sentence rounds to zero. If you have a concrete problem implementing the algorithm, post your code and say where it hurts. But you're not describing a concrete implementation problem; you're describing an emotional state. Maybe somebody here can write prescriptions, I dunno. I can't. – Sharecrop

here is a paper that might help that talks about how to smooth your mesh when using the marching cube algorithm – Yolk

@ people casting close vote - Do not close this please. The asker is a new person. He has clearly put in effort into his research and written a well formatted question and just because there isn't enough code (i'm sure its because of complexity rather than inability to code) it's getting close votes. Ed Plunkett is right, we need precise information but I'm sure the OP is happy to furnish it. Please at least let the bounty time die before you cast votes. – Mongolic

@Angi people will not post answers until they have precise code to deal with. You say you've seen pseudo-code that you dont' understand, link that and other source material that you've dealt with in (even if it doesn't work). Setup a github minimal complete sample project with all your code in its current state and link it here so we can see exactly what you've done and what needs fixing. As it stands right now, your bounty rep will get wasted. – Mongolic

in your case i would just a dictionary like Dic<bool, Point3d> so you will have points in one hand and bools in another. – Josselyn

How did you arrive at the size of your bool array? Is your data model to scale of the range of the 3D printer? – Blennioid

The input is a series of gray-scale slices(CT scan). The length of the bool array is the length of 1 slices, the width of the array is the width of 1 slice, and the height of the array is the number of slices. The dimensions of all slices are the same. My data model is to scale of the range of the 3D printer. – Angi

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.

Shealy

Answer 2

disclaimer: it won't be a full answer, as I don't have the time to put one right now, I'm not rushing for the bounty.

I'll start an idea for a way to do ray-tracing, that should produce a list of points visible on the outside of your shape.

basically for raytracing, you trace a ray from each point of each face (and maybe diagonal), this won't be the best way, but that should output something.

the ray is a vector, in this case if you go from the x,y plane on one side to the other it'll be something like this [pseudocode]

bool[,,] points = getpoints();
bool[,,] outsidepoints = new bool[width,height,depth];
//this is traceray from x,y,0 to x,y,depth
for(int x=0;x<width;x++){
  for(int y=0;y<height;y++){
    for(int z=0;z<depth;z++){
      if (points[x,y,z]==true){
        outsidepoints[x,y,z]=true
        z=depth;//we break the last for
      }
    }
  } 
}

do the same for all 6 combinations in the last loop (x=0..width, x=width..0, y=0..height etc...)

that should work well enough if you have no overhangs (e.g. a single convex shape), if you do they might get cut off.

once you have your external points you just need an algorithm to make them into triangles (something like making slices of connections in each plane along an axis, then sewing for the last 2 connections.)

ps: if what you want is just display, you can use raytracing by taking the value of first point that connects with the ray for each point of your display (expect an image every couple of seconds depending on your hardware and picture size.) for this you can use the Line drawing algorithm (check wikipedia, there's a sample, in 2d though) best thing about it is that you can directly get the color from your picture (just set a threshold at which the pixel is void (the ray passes through))

another edit: if you need any precision, pm me, or comment down here

using pass through rays:

bool lastwasempty=true;
for(int x=0;x<width;x++){
  for(int y=0;y<height;y++){
    for(int z=0;z<depth;z++){
      if (points[x,y,z]==true){
        if(lastwasempty){
           outsidepoints[x,y,z]=true
        }
        lastwasempty=false;
      }else{
        lastwasempty=true;
      }
    }
  } 
}

Adroit

The problem is that I need the image to be hollow according to a minimum preset value. I'd give a +1 if you show me how to get the points on the edges inside the sphere not just the outer border. Thanks for the "Line drawing algorithm". – Angi

@Angi you could try to make the rays pass through (don't break at the first hit), then check if the point before was transparent (false), I'll add some code for that in the answer. – Adroit

@Angi , the code will detect if there's an empty space after, as you parse all sides (and top and bottom), before one way is after the other way. to create a printable model, you could probably just scan by line or make a cube centered around each point, which would probably make something ugly, though. So something I thought, if you need a smoother model, you can use the Catmull–Clark subdivision surface algorithm. to make a 3d printable file, you can export it to obj hodge.net.au/sam/blog/wp-content/uploads/obj_format.txt and use a slicer, or just generate g-code. – Adroit