Marching Tetrahedrons is a commonly-used technique in scientific and medical visualization to construct surfaces in a 3-D scalar field. For example, MRI scans produce tissue density data at the vertices of a regular (cubic) 3-D grid. The problem is to construct an isosurface (i.e., a surface connecting all the points in the data that have the same value — also called the isovalue). The data value is chosen based on the density characteristics of the tissue to be visualized (i.e., internal organs or bones).
Most implementations of this technique use OpenGL in a Unix environment because this combination is viewed to be the only viable platform providing the high performance required to produce quality visualizations containing millions of triangles. Furthermore, OpenGL is portable and is viewed to be much simpler to use than Direct3D. In this article, I demonstrate how to use the new Direct3D 8.0 API to produce high-quality images in a modern Windows environment (i.e., with a high-speed CPU, hardware accelerated video card, and copious memory). The article also uses the simplified Direct3D 8.0 API that mitigates the code complexity that discouraged developers under earlier releases of DirectX.
To simplify the code for the article, I create 3-D contours from two mathematical functions rather than from MRI data. Producing contours of well-known functions also provides an idealized environment for testing and debugging visualization code, because the data is noise-free and does not contain complex features such as holes and folds. The basic techniques for visualizing MRI data and 3-D functions are identical, however.
The first function, v=sin(p)*sin(q)*sin(r), produces four hollow “blobs”
given an isovalue of -0.3, and the second function, v= p*p/0.5+q*q+r*r,
produces a hollow ellipsoid given an isovalue of 9.0. The p,
q, and r values range from -pi to +pi for both
functions. The function value ranges from -1.0 to 1.0 for the
first function and from 0.0 to 39.48 for the second function.
I have chosen isovalues to produce intersecting surfaces; I make the four
blobs opaque blue in color and the ellipsoid translucent gray so that the
portions of the blobs that are inside the ellipsoid show through the ellipsoid
(Figure 1).
In this article, the code in geometry.cpp (Listing
1) performs the mathematics needed to generate the surface geometry, and
render.cpp (Listing 2) contains Windows
and Direct3D specific code to render a simple animation of the isosurface
images in a window. You will need to install the DirectX 8.0a SDK and modify
your compiler’s settings to include DirectX’s “include”
and “lib” directories in the appropriate search paths.
The Marching Tetrahedrons Algorithm
The dataset to be visualized consists of a 3-D volume organized in a cubic grid. Each vertex in the grid defines a 3-D position (i.e., x, y, z coordinates) and a value. The fundamental problem is to construct a set of triangles that connects the vertices on the grid with a specific data value. These planar triangles approximate the actual (potentially curved) surface to be visualized. The finer the grid, the closer the approximation will be to the actual surface, resulting in a better-looking image. However, a finer grid will produce a larger number of triangles with a concomitant increase in required memory and processing time.
In GenerateGeometry() in geometry.cpp (Listing
1), I set up a logical grid consisting of MAX_X_SEGMENTS, MAX_Y_SEGMENTS,
and MAX_Z_SEGMENTS segments along the x, y, and z
axis, respectively. You can control the coarseness of the grid by modifying
the number of segments along each axis. The function also sets the global
function pointer variable g_pfnEval to point to the function that will
produce the vertex value. This indirect call to the evaluation function enables
the flexibility required to extract different surfaces with maximum code reuse.
Next, I set the global variable g_fIsovalue to the desired isovalue
and call ProcessCube() for each cube in the grid. After extracting
the triangles for the opaque surface, I similarly extract the translucent
surface. Note that the same vertex buffer, g_pVertBufD3D, is used to
store the triangle data for both surfaces. The global variables g_iNumSineVerts
and g_iNumVertices contain the number of vertices in the first surface
and the total number of vertices in both surfaces, respectively.
The first step in the process of triangle extraction is to divide the cubes formed by the vertices of the mesh into tetrahedrons (hence the name “Marching Tetrahedrons”). Since a tetrahedron itself is composed of four triangles, it is easier to extract triangular surfaces from tetrahedrons than from cubes. Of the several techniques for decomposing cubes into tetrahedrons, I use the one that minimizes the number of tetrahedrons in order to minimize the number of resulting triangles.
Figure 2 depicts the five tetrahedrons that can
be extracted from each cube. Note that two adjacent cubes cannot use the same
decomposition scheme because the adjoining tetrahedron faces from the two
adjacent cubes will not be coincident. Since it is essential that all adjoining
tetrahedron faces coincide in order to produce a continuous surface, a second
decomposition is constructed and the two schemes are applied to alternating
adjacent cubes (e.g., the triangles on the right face of the “Even Decomposition”
will match the triangles on the left face of the “Odd Decomposition”).
The alternating decomposition is achieved by applying the “even”
or “odd” scheme when the sum of the logical x, y,
and z coordinates of vertex 0 of the cube is even or odd, respectively.
ProcessCube() in geometry.cpp (Listing
1) extracts the tetrahedrons from the cube specified by the x,
y, and z parameters. The (x, y, z) triple specifies the
logical coordinates of vertex 0 of the cube. g_aiCubeVert, an 8x3 integer
array, is used to construct the remaining seven vertices of the cube from
the x, y, and z parameters. Each cube vertex is computed
by adding the three elements of each row in g_aiCubeVert to x,
y, and z, respectively. g_aiCubeTets, a 2x5x4 integer
array, is used to extract the tetrahedrons from the cube. The first dimension
specifies the orientation of the cube (i.e., even or odd); the second dimension
specifies the cube number (i.e., 0 through 4); the third dimension specifies
the four cube vertices that make up the tetrahedron.
Next, ProcessCube() converts the logical (x, y, z) coordinates of a tetrahedron vertex into world (x, y, z) coordinates (recall that the world coordinates range from -pi to +pi and therefore have a total range of 2*pi). The tetrahedron vertices reside in the local variable aVert. ProcessTetrahedron() uses this array to construct the surface triangles.
ProcessTetrahedron() constructs the surface triangles from the four
vertex coordinates of the specified tetrahedron. This process is complicated
by the fact that the isovalues will not necessarily be found exactly at the
tetrahedron vertices. It is likely that one of the two vertices along a tetrahedron
edge is less than the isovalue and the other vertex is greater than the isovalue.
In this case, I must use linear interpolation to find the position of the
isovalue point on the line segment connecting the two vertices. Given two
vertices (x1, y1, z1) and (x2, y2, z2) with values v1
and v2 respectively, the Interpolate() function in geometry.cpp
(Listing 1) computes the x, y, and
z coordinates of the isovalue vertex by using the following formulas:
x = x1+(isovalue-v1)*(x2-x1)/(v2-v1)
y = y1+(isovalue-v1)*(y2-y1)/(v2-v1)
z = z1+(isovalue-v1)*(z2-z1)/(v2-v1)
Figure 3 shows the various combinations of the
four vertex values and the resulting isosurface triangles. (Other combinations
are possible, but they do not produce surface triangles and are ignored.)
Vertices marked “-” are less than the isovalue, vertices marked
“=” are equal to the isovalue, and vertices marked “+”
are greater than the isovalue. There are five distinct cases and cases two,
three, four, and five have several instances resulting from rotating the tetrahedron’s
vertices. Note that in case five, two coplanar triangles are extracted from
four interpolated points.
ProcessTetrahedron() first computes the vertex values and determines the triangle extraction strategy to be used. The integer arrays aiEQ, aiGT, and aiLT contain the indices of vertices whose values are equal to, greater than, and less than the isovalue, respectively. The integer variables iNumEQ, iNumGT, and iNumLT track the number of vertices in the arrays described above. As each triangle vertex is computed, it is added to the Direct3D vertex buffer. AddVertex() inserts the vertex into the vertex buffer and updates g_iNumVertices. The function also calls ComputeNormal() to compute the normal vector for the inserted vertex. Direct3D uses the direction of normal vectors to determine the vertex’s orientation with respect to the incident light and the viewer when computing the color for each pixel in a triangle. In this application, the direction of a vertex’s normal vector is given by the gradient of the isosurface at the vertex. Therefore, I compute the normal vector by evaluating the partial derivatives of the surface function at the vertex along the x, y, and z axis. Also note that I normalize the normal vectors (i.e., make them have unit length,) so Direct3D does not need to normalize them in subsequent shading computations.
Initialization
WinMain() in render.c (Listing 2)
calls SetupWindow() to create the display window and the Direct3D objects
and then calls GenerateGeometery() to extract the surface triangles.
After constructing the surface geometry, the function implements a simple
message dispatching loop to render the next frame of the animation that calls
the Display3D() function when no messages are pending. When the display
window is closed, WinMain() calls CleanupWindows() to destroy
the allocated Direct3D and Windows objects.
SetupWindow() in render.cpp registers the window class for the application and creates the window. Next, I call Direct3DCreate8() to obtain the COM interface to Direct3D. I then use the Direct3D interface to create a device that I will use to render the triangles. Since this application runs in windowed mode, the pixel format of the rendering device must be compatible with the pixel format of the current video mode. Therefore, I call GetAdapterDisplayMode() to obtain the current video mode format.
The PresParamsD3D variable specifies the requested device’s characteristics. The PresParamsD3D.Windowed member is set to indicate that the device will operate in a window. The PresParamsD3D.BackBufferCount parameter count is set to indicate that there is one backbuffer on which the animation frames will be rendered. This technique is commonly used in animation where the next frame is rendered onto a non-visible buffer (the backbuffer) while the previous frame is being displayed. As soon as the new frame is rendered, the content of the backbuffer is swapped with the content of the visible buffer (in practice, only buffer pointers may be swapped in order to maximize performance). The PresParamsD3D.SwapEffect member is set to indicate that the contents of the previous frame (visible buffer) be discarded. This option may conserve memory bandwidth on some video hardware. The PresParamsD3D.EnableAutoDepthStencil and PresParamsD3D.AutoDepthStencilFormat members are set to create a 16-bit per pixel z-buffer. The z-buffer (or depth buffer) is used to prevent the pixels of far triangles overwriting pixels of near triangles. This effect is essential because objects closer to the eye occlude farther objects.
Next, I call CreateDevice() to create a device with the characteristics specified above. In addition, the D3DDEVTYPE_HAL parameter specifies that the device should use the HAL (hardware abstraction layer) for rendering. This specification will cause CreateDevice() to fail on video adapters that do not provide hardware acceleration (or whose device drivers do not support the Direct3D 8.0 HAL architecture). The D3DCREATE_HARDWARE_VERTEXPROCESSING parameter indicates to Direct3D that, once generated, the vertices will require processing only by the hardware.
After the device is created, I set up the viewing transformation (i.e., the position of the eye). This transformation is created by calling D3DXMatrixLookAtLH() and SetTransform(). The eye is positioned in a standard left-handed axis system (positive x going right, positive y going up, and positive z going back into the screen) at point (0, 0, -9). Next, I call D3DXMatrixPerspectiveFovLH() and SetTransform() to set up the projection transform. This transform is used to render the 3-D objects on a 2-D screen with foreshortening to simulate the appearance of depth. (Far objects appear smaller than objects closer to the eye.) The field of view is set to pi/3 radians (i.e., 30 degrees). The aspect ratio of the projection is 1.0 because the height and width of the window are the same. The near and far z values of the viewing volume is set to ensure that the entire image is visible without clipping any part of the visualization. Making the view volume arbitrarily large, however, will degrade the discrimination ability of the z-buffer because all possible pixel depths must be scaled down to fit into the 16-bits per pixel z-buffer representation.
Next I set up the scene lighting. I call SetLight() to establish a directional light source for the scene and LightEnable() to turn on the light. Then I call SetRenderState(D3DRS_SPECULARENABLE, TRUE) to enable specular highlighting of triangles. I also create a small amount of gray ambient light by calling SetRenderState(D3DRS_AMBIENT, 0x00202020). This gives a small amount of color to those triangles that are facing away from the light source, as opposed to being completely dark (black).
Direct3D can automatically prevent the drawing of backward-facing triangles. However, the depth buffer also provides the same functionality more reliably, albeit more slowly. Therefore, I call SetRenderState(D3DRS_CULLMODE, D3DCULL_NONE) to disable the back-face culling operation.
Next, I create the vertex buffer that will contain the triangle vertices
for the surfaces. For this application, each vertex buffer entry, defined
by structure VERTEX in marchtet.h (Listing
3), consists of two instances of class D3DXVECTOR3 to store the
position and normal vector of the vertex. Note that I use a “user allocated”
vertex buffer as opposed to a Direct3D vertex buffer (IDirect3DVertexBuffer8).
I then describe the vertex buffer format to Direct3D by calling SetVertexShader(D3DFVF_XYZ
| D3DFVF_NORMAL).
My final setup action is to construct a rotation matrix with approximately 10 degrees of rotation around the x and y axis. I multiply world transform matrix with this rotation matrix after rendering each frame so that the surfaces appear to have been rotated around the origin (0, 0, 0) when the next frame is rendered and displayed.
Rendering the Surface Triangles
The actual rendering occurs in Draw3D(), defined in render.c. In this function, I call Clear() to set the contents of the back-buffer to the background gray color and to set the z-buffer entries to the farthest depth. To produce the transparency effect properly, I first draw the opaque surface. I call SetMaterial() to specify opaque material. The material specification determines the color, shininess, and transparency of the rendered object. Next I call DrawTriangles() to draw the opaque triangles. The first and second parameters specify the starting index of the triangles and the number of triangles to render from the buffer, respectively.
DrawPrimitive() in render.c (Listing
2) draws the requested triangles in a loop. A loop is required because
the maximum number of triangles that can be rendered in a single rendering
call is limited and is reported in the MaxPrimitiveCount member of
the D3DCAPS8 structure. I query the device capabilities into global
variable g_CapsD3D in SetupWindow() after obtaining the Direct3D
device. The actual rendering is performed by calling DrawPrimitiveUP().
The first parameter of the function specifies that the vertex buffer contains
a list of triangles (three vertices per triangle). The second parameter specifies
the number of triangles to render. The third and fourth parameters specify
the memory location of the first vertex and the size of each vertex, respectively.
After drawing the opaque surface, I prepare Direct3D for rendering the translucent surface. The idea is to blend the translucent pixels with the pixels already in the backbuffer. The color of a pixel is specified by its red, green, blue, and alpha values. The alpha value determines the level of opacity of the material (0 being completely transparent and 1.0 being completely opaque). I first call SetRenderState() to enable alpha blending. The next two SetRenderState() calls specify that the color of the source pixel (pixel to be drawn) will be multiplied by the source material’s alpha value and added to the destination pixel (i.e., the corresponding pixel already in the buffer). The fourth SetRenderState() call disables z-buffer updates because the far translucent triangles must show through the near translucent triangles. However, z-buffer testing is still performed because far transparent triangles will be hidden by closer opaque triangles. There are other techniques that produce more realistic images but require a costly sort operation every frame because the triangles must be rendered in a far-to-near order relative to the eye.
Next, I draw the translucent surface by calling DrawTriangles() for the remaining vertices in the vertex buffer. After drawing the translucent surface, I re-enable z-buffer updates and disable alpha blending in preparation for rendering the next frame. Next, I call Present() to display the backbuffer onto the window. And finally, I apply the rotation matrix to the world transform in order to rotate the surfaces for the next frame.
Conclusion
I achieve nearly five frames per second when executing this code on a 800MHz Duron processor with an nVidia GeForce2 MX AGPx4 video card containing 32MB video memory under Windows ME. The initial geometry generation code also only takes a few seconds to generate the 1,663,968 vertices (554,656 triangles). Production applications perform several optimizations to improve performance, reduce memory requirements, and improve the image quality. For example, the image geometry can be simplified by reducing the level of detail for distant surfaces. Also note that because vertices are typically shared by several adjacent triangles, it is more memory efficient to only store unique vertices and to represent each triangle by a trio of integer indices into the vertex buffer. Direct3D supports this representation via the indexed triangle primitive.
This application clearly demonstrates that Windows and Direct3D 8.0 is a viable environment for performing high-quality scientific visualizations. Furthermore, the code complexity in this application is roughly comparable to that of equivalent OpenGL code. This shows that Microsoft is actively addressing programmers’ concerns regarding the unnecessary complexity involved with setting up applications under previous versions of Direct3D.
Yogi Dandass is a researcher at Mississippi State University’s High Performance Computing Laboratory and is a lecturer at the Department of Computer Science. He designs software libraries in C++ for numerical computing, artificial intelligence, and real-time communication middleware. He can be reached at yogi@cs.msstate.edu.
