Subject: Description of adaptive marching cubes, for posterity. SUMMARY I have developed an adaptive version of the marching cubes algorithm for isosurface generation; this algorithm will increase both usability and performance of FreeForm. I describe the motivation and goals for the algorithm, the implementation of the recursive isosurface generator, and the implementation of a post-pass which patches the holes left in the surface by the isosurface generator. I also highlight several areas in which the implementation could be straightforwardly improved. MOTIVATION The simple marching cubes algorithm currently employed by most isosurface generators suffers from several performance and output quality problems, fundamentally caused by the creation of millions of redundant, coplanar triangles. Creating between one and five triangles for every voxel through which the isosurface passes, the algorithm, as implemented in FreeForm, fills most of main memory with vertex and triangle data, leading to unacceptable levels of swapping on detailed models, and slowing the time-consuming graphics rendering pipeline to a crawl. In addition, the output files produced by FreeForm are very "triangle-heavy", leading to complaints from customers who have trouble manipulating the meshes in other programs. My goals in this project are to: 1. Reduce the number of triangles and vertices output by the isosurface generator by an order of magnitude. 2. Do not significantly decrease the quality of the output at a given voxel resolution. 3. Decrease, or at least do not significantly increase, the run time of the isosurface generator. 4. In an interactive editing mode of the modeler, allow the option of generating a higher level of detail in the "current region of interest" than in more static regions of the model. More specifically, FreeForm will be improved in these ways: 1. Reduce the number of triangles and vertices output by the isosurface generator by an order of magnitude. a. This will free up memory and speed up the graphics pipeline. b. This will make exported mesh models more "lightweight". 2. Do not significantly decrease the quality of the output at a given voxel resolution. a. Decreasing triangle and vertex count while keeping quality constant at a given resolution will decrease memory and CPU load. b. The extra memory and cycles can be used to increase the operating voxel resolution, improving the upper bound on maximum model detail. 3. Decrease, or at least do not significantly increase, the run time of the isosurface generator. a. The adaptive algorithm must be able to generate isosurfaces in real time, as the model is edited (although, of course, only the modified region need be updated). b. The graphics pipeline will be faster due to the reduced triangle count; thus, some graphics CPU cycles can be traded for isosurface generation cycles. 4. In an interactive editing mode of the modeler, allow the option of generating a higher level of detail in the "current region of interest" than in more static regions of the model. a. By using information about the region of interest of the user (gleaned from pointer location, interaction history, and camera position information), the isosurface generator can avoid rendering complex details in areas of the model that the user cannot see or is not looking at. OVERVIEW Instead of designing an isosurface generation algorithm completely from scratch, I decided to build upon the strength and simplicity of the tried-and-true marching cubes algorithm. In order to make the algorithm adaptive to the level of detail in a region, however, I re-evaluated the usual assumption that the marching cubes algorithm be applied only to unit cubes composed of neighboring voxels; instead, it is possible to apply the algorithm to any cube with sides of integer length. This technique lends itself well to a recursively adaptive algorithm. First, the voxel space is subdivided into cubic chunks with sides of length 2^N. Then, for each cube, the standard marching cubes algorithm is applied to generate between one and five triangles. This rough rendering is compared against the real voxel data within the cube; if the error is above a threshold, the cube is divided into eight subcubes and the algorithm is recursively applied to each subcube. The result is that flat regions result in few triangles, and detailed regions contain many triangles. This approach gives greatly reduced triangle counts over the traditional marching cubes implementation; however, at the boundaries between cubes of different recursion depths, it leaves holes or T vertices in the triangle mesh. Solving this problem is much more difficult than designing the recursive algorithm itself, but the hole-patching algorithm I developed gives very good results. The patching algorithm speed is very reasonable, at the cost of code and data structure complexity. IMPLEMENTATION Details of the Adaptive Marching Cubes Algorithm Most marching cubes algorithms do not place mesh vertices exactly halfway between grid vertices; instead, they linearly interpolate, using the voxel values, to place the mesh vertex as close as possible to the isosurface. When the cube to march has sides of length greater than one, it is possible to go one better than that, and perform a binary search on the voxels along a grid edge to find exactly where the isosurface intersects the edge. More precisely, the algorithm uses a binary search to find which two voxels along the edge the isosurface passes between. Then, linear interpolation between those two voxels is used to generate the exact vertex position, and the vertex is stored in a hash table of grid edges (edges between adjacent voxels) to vertices, to prevent the creation of more than one vertex at the same point. (This becomes more important later on, in the patching phase, although it is generally good practice in any case) The error function used in the current implementation of the adaptive algorithm is fairly simple, but produces good results on certain data sets. Given a voxel region and a list of triangles, it computes the centroid of each triangle and estimates its distance from the isosurface by subtracting the voxel value at the centroid from 127.5. It then returns the average value of these differences. On certain complex datasets, this error function could be misled into producing an inaccurate estimation of the difference between the real dataset and its rendering to triangles. There are many potential improvements to the function; I will outline one that I believe would produce mathematically perfect error values. Given a list of triangles and a voxel region A[x][y][z], we would begin by rasterizing the triangles to a new voxel region B[x][y][z]. (This could be done using our existing rasterization techniques, but I imagine that there could be significant optimizations from using the knowledge of the specific form of triangles generated by applying marching cubes to a single cube.) Then, the error value would be the vector distance between the two voxel fields A and B, i.e. /--------------------------------------- Error = / Sum Sum Sum (A[x][y][z] - B[x][y][z])^2 \/ x y z (forgive the ascii art). Lastly, the adaptive algorithm uses a "maximum depth" mechanism to increase rendered detail in a region of interest, or "active region": the maximum recursion depth allowed to the adaptive marcher decreases as the inverse of the square of the distance to the "active point." Instead of naively sampling the "maxdepth" at the center of the current cube, however, it attempts to take the maximum "maxdepth" over the whole of the interior of the cube. This prevents a blocky look caused by the renderer bailing out too quickly on large cubes. Details of the Mesh Patching Algorithm The recursive algorithm presented above produces the correct result detail-wise, but the mesh that it produces is not a manifold: it contains T vertices and holes between triangles that were generated from cubes of different sizes. Finding and patching these holes is not a trivial task, code-wise or efficiency-wise. In order to facilitate the finding of "trouble areas" in the mesh, code to add meta-information to the mesh is added to the adaptive mesh generation pass. First, connectivity information is added to the mesh (which previously had consisted of a list of vertices and triangles) by associating a list of triangles with each vertex. Second, a four-bit "adjacency" field is added to each vertex, filled out during mesh generation as follows: Every vertex is on a single edge of the voxel grid; thus, every vertex borders exactly four voxel cubes (or, if you like, four "quadrants" of the three-dimensional space). At the creation of a vertex, the adjacency field is all zeros; every time the vertex is used in the rendering of a cube, the appropriate bit in the adjacency field is set (depending on which quadrant the cube is in). If all is well at the end of the mesh generation pass, all bits in all vertices' adjacency fields should be ones, since that means that the isosurface is rendered on all four sides of the vertex. Due to the way I colored them in my example code, vertices with all four adjacency bits set are called "grey" vertices, and those with one or more bits still zero are called "green" vertices. Once the adaptive mesh generation pass is complete, a patching pass is entered. The holes (and T vertices, which are another form of hole) created by this algorithm were experimentally seen to be always of a particular character: they were composed of two adjacent grey vertices and N green vertices, where N > 0. Such a hole can be patched over by finding the common triangle of the two grey vertices and splitting it into N+1 new triangles, effectively stretching the triangle over the hole and attaching it to the green vertices. This algorithm is imperfect, but gives good results on the datasets I've tried it on. For one thing, I have not mathematically guaranteed the character of normal holes as described above; this needs more research. Also, problems occur when two holes adjoin (i.e. when one green vertex is on the edge of two holes). This can probably be fixed the "right way" (by patching over both holes) with some more work into improving the algorithm; however, an easier workaround is to ensure that the depth of the adaptive recursion never changes by more than one level between adjacent cubes. Problem holes only seem to occur when the depth changes by two or more levels, which is an undesirable situation in any case, since it makes the rendering look "blocky". CONCLUSION I have developed an adaptive isosurface generation algorithm based on the marching cubes algorithm. The results are favorable, in both quality and performance; and I see no obvious impediment to incorporating the algorithm into FreeForm. There are some immediate possible improvements to the algorithm, which I have outlined above; also, I have left several possible optimizations and improvements unimplemented in my code, for the sake of simplicity. At this point, it would seem reasonable to move forward on implementing this algorithm in a demo version of FreeForm. This will allow us to determine the primary issues with regards to the incorporation of the algorithm into the primary FreeForm release by the software development team, and will permit us to experiment with the various optimizations and improvements I have outlined above. In conclusion, I recommend that we continue in this direction, to evaluate the benefits of this optimized algorithm to the company and to the end users.