Marching Cubes
I want to write a short beginners article about the marching cubes algorithm. At first I must say that the marching cubes algorithm is patented, and thus may not be used in commercial programs. So don't blame me if you get arrested or whatever. This article is only meant for learning purposes, not for implementation.
Another thing is that I am no expert at marching cubes. Despite that, this article might be interesting for people who already implemented it, as it might feature things that are even new to you.
I hope you read the great article by Paradox in Hugi 21. This is some valuable information even though I will use a different way for displaying the metaballs than described there.
I will now suppose that you have already built up a 3D array of numbers, let me say between 0 and 1, where numbers below 0.5 mean that the point is outside your blobs and a number greater than 0.5 means solid. The first thing that you have to make clear is that with a 3D array you only have an approximation of what you really want to display. If you wanted to display the metaballs correctly you would need all possible points, which are infinitely many. The thing I hear most often now is to "choose" the points that are "on the surface" of your metaballs. As you only have a finite number of points, the points in your grid won't be equal to 0.5. So you will maybe say that all points between 0.4 and 0.6 lie on the surface of the grid.
It should be clear that this approach is far from being perfect. The first disadvantage is that you might miss some points on the surface resulting in holes in your blobs. To reduce these things, you will increase the range to [0.3|0.7] or whatever, which will result in blobs with several surfaces, which will lead to lower frame rates. A second disadvantage is that the edges of your blobs will always lie on the grid, which results in a very ugly movement of the metaballs.
I suggest another solution to find the edges of our metaballs. Let me solve it in 1D first with a simple example:
0.1-----0.2-----0.6-----0.8-----0.4
The numbers shown here lie in the grid, and the - represent the borders in the grid. Where would you set the edges by hand? Let us analyse the above picture. At first I must mention that I see the edge points in the grid as infinitely small. This is the reason why the possibility that the surface of the blob goes through the grid points is zero. So I will in the future suppose that that this is not the case. Even though you should create your implementation in a way that the program will at least not crash...
Now let us have a look at the top picture again. Between 0.1 and 0.2 there is only air. Probably no surface of the blob. There must be at least one surface edge between 0.2 and 0.6. We suppose that there is only one, and as the surface is at 0.5, it will be closer to 0.6. Between 0.6 and 0.8, it is solid. Even if there is any surface, we will not see it. and between 0.8 and 0.4 there will at least be one surface edge. Possibly close to 0.4.
The idea should be clear now. I take every edge in my 3D array, and test whether the borders of it are solid or not. If exactly one is solid and one not, there is a surface edge between them. To find a good and fast solution, I would suggest to use linear interpolation for the position of the point. You may of course also take the middle between the two grid points, or use an interpolation that takes the 1/x^2 character of metaballs into account. It is up to you. Let me show the resulting points of the above example (with linear interpolation):
0.1-----0.2---x-0.6-----0.8---x-0.4
Nice, isn't it? And this is of course the same in 2D and 3D. Let me draw some 2D thing, I will also use later, with just a G where the grid is solid, and an o where it is not:
o------o------o------o------o
| | | | |
| | | | |
| x | | |
| | x | |
o----x-G------G----x-o------o
| | | | |
| | x | |
| | | x |
| | | | |
o--x---G---x--o--x---G--x---o
Fantastic. You will not miss any required point and you also will not get points that are not on the surface of the blob. Believe me. Or try it yourself. So the first step in your metaball implementation will be creating an array of vectors using this algorithm. In addition to this you will also have to store somewhere where the vector of a special border is stored (there is a maximum of one vector per border). The simplest thing would be building up an array for every border in your 3D grid, and store the appropriate vertex number in it. But is this also fast? I am still unsure.
The problem we still have is to find the triangles that create the surface. Let me at first solve the problem in 2D. It's much the same in 3D. We now do not have to find triangles but rather lines. The good thing now is that for finding those, we only need information from one "cell" where I mean 4 gridpoints, building one rectangle as cell. And we even do not need the exact numbers in the edgepoints nor the positions of the surface-edges on the borders of the cell. Let me illustrate this with an example:
G------G
| |
| |
| |
| |
G------o
We do know that in this cell, we will have two edge points, one on the right and one on the bottom border. The exact position of them is unclear, but also unnecessary. G is solid. o is air. We will need exactly one line, from bottom to right. We of course now need the information which vectors are on those borders. We also do know that the "outside" of our line will be towards bottom-right. This might be important if you use backface culling for your triangles.
It should be clear that there are 16 possibilities of cells in 2D. Every edge may be solid or not, that means you have 2*2*2*2 = 16 possible cells. And you will have a maximum of two lines per cell. (In 3D it is 256 possibilities and 4 triangles if I am not wrong). You might go through all possibilities with if..then..else conditional branches, but conditional branching is expansive on modern hardware. So I'd create a bitvector out of your cell, with 0 meaning air and 1 meaning solid, and a fixed position of every gridpoint in the bitvector. So the above cell would have a 1110 bitvector if you store the edges top-down-left-right. Than you use a table to look up the required lines (triangles). This table will perfectly fit into any cache, and you do not need any conditional branching at all. This information must now be converted so that you get the correct vectors for this cell. I am sure you will find a good implementation for this.
Now everything is done, and you will get the following picture (I left away the grid):
o o o o o
x--
/ \---x---\
o x G G x o o
| -/
/ x-/
| / --x-
/ / / \
o x G x o x G x o
The algorithm is especially attractive, as you will not get any redundant data, nor will you get any holes (if you take care with 0.5 at gridpoints). And the other thing is that by using linear interpolation, the viewer will hardly see the grid, as the surfaces move smoothly within the grid, and not only along its borders.
Some thoughts about speed improvement
It should be clear that I only presented the idea of marching cubes. There are many things to improve the speed. I will only give you two hints:
1. If the metaballs do not cover much space, it might be clever to start calculating in the middle of the metaballs and work your way outwards, until you cross the surface. You must then of course avoid double-working on your data.
2. You might at first work on a very small grid, and at "hot" regions, where the gridpoints are around 0.4 and 0.6 you use subdivision. This may of course be done in several steps. But I think that using two different grid depths is far enough.
Some thoughts about vertex normals
So far I have not found any method to calculate good looking normals. I know they exist, but the three articles I read about it didn't deliver the best results in my implementation. I would be happy if I found an article in next Hugi about this. :)