Volume Rendering: Marching Cubes Algorithm
CS 580 Programming Assignment 5
Lavanya Viswanathan
December 1, 1997
Volume Rendering
Volume rendering techniques enable the visualization of sampled scalar
functions of three spatial dimensions. Typical data that needs to be rendered
this way consists of data from X-ray Computer Tomography (CT) scanners,
and PET scanners. These data collect a series of two-dimensional arrays
of densities ordered by depth planes.
Marching Cubes algorithm
The marching cubes algorithm is used to construct contours from volumetric
data by applying an isovalue surface detector (with threshold values for
densities). These contours form boundaries between distinct regions
of the data. For instance, in a data set containing densities from a CT
scan of a human head, one can generate contours that allow us to look only
at the skin portions and other contours that allow us to look at only the
bones instead.
The marching cubes algorithm is implemented as follows:
First a three dimensional grid of cubes is generated and each vertex
of each cube is assigned a density value (read in from the data file).
Then, we examine each cube in turn and determine, for each of its vertices,
which side of the isosurface contour it lies on. There are two trivial
cases: when all the vertices of the cube are on one side of the isosurface.
This gives us a set of edges of the cube that are intersected by the isosurface.
The exact point of intersection of the isosurface is then computed by interpolating
the densities at the vertices to get the desired threshold density. Once
the intersection points for a cube have been determined, we compute a three-dimensional
polygonal model for the contour. This is a complicated process that first
requires that we compute all possible permutations of the cube. The fact
that each cube possesses eight vertices and that there are two states for
each vertex ("inside" or "outside" the isosurface) means that there are
28 = 256 ways in which the contour can intersect the cube. Because
of symmetry, we can reduce these cases to 15 basic cases that are shown
in the diagram below.
Below, I present pseudocode for the implementation of the marching
cubes algorithm.
Read the volume data set from file
(see section on Data Sets for details
on how to do this)
Generate Facet Data, i.e.,
compute the facets for each of the
256 possible permutations of a cube.
Generate Edge Data, i.e.,
compute the intersection edges for
each of the 256 possible permutations of a cube.
Initialize memory for all the data
structures in the program.
NI = number of rows in a slice of data
NJ = number of columns in a slice
of data
NK = number of slices of data
for slices k=0 through (NK-1), do
{
if (k == 0) /* i.e., for
the first slice */
{
Read in data for slice 0
compute XEDGE data and YEDGE data for
initial conditions
by interpolation of densities
(note: these data structures are used
to store the intersections that have already been computed for earlier
slices. Because of the continuity and continguity of slices, we do not
need to recompute these intersections for every cube. This saves a lot
of computation.)
}
Read in data for slice (k+1)
for rows i=0 through (NI-1), do
{
}
}
The diagrams below explain how the data structures in this program are
organized:
-
Vertex numbering for each cube : 8 vertices, numbered from 0 to 7.
-
Edge numbering for each cube : 12 edges, numbered from 0 to 11.
-
Computations performed for each cube :
The above diagram shows the computations performed for each cube. For
nine out of 12 edges, the intersection with the isosurface is found by
using previously stored values (XEDGE, YEDGE, ZEDGE, TOPYEDGE, OLD6PT or
OLD10PT). New intersections are found only for the remaining three edges
(numbered 5, 6 and 10). These edges are highlighted in the diagram with
red crosses. Once these new intersections are computed, the existing stored
values are updated to contain the new values for the next iteration. This
is done in the function LOAD_FACET. New intersection points are computed
through linear interpolation of the densities at the cube vertices.
Now, I present the pseudocode for the function LOAD_FACET:
LOAD_FACET
For each edge in the EDGE_DATA table
for the current cube (index in permutation table = cube_index), do
(note: these are the edges on which
intersections with the contour surface would occur for the current cube)
{
Find the intersection point with the
contour surface
(this is done by considering several
cases, as shown in the diagram above: we store an array of pointers into
a list of intersection points, 'edge_table')
Update the stored values XEDGE, YEDGE,
ZEDGE, TOPYEDGE, OLD6PT and OLD10PT
}
For each facet in the FACET_DATA table
for the current cube
compute the vertices of the facet.
(note: each facet is stored as a triangle
that is rendered later)
Examples
Data sets
Three data sets were used:
Data is stored in these files as bytes. In each of these files, data is
stored in the following order:
Color coding
The table below shows the volume data sets for some Z slices. The data
is visualized using a simple color coding scheme, where the intensity of
the pixel color is proportional to the density of the volumetric data at
that point. Since the data range from 0 to 256, the densities have been
normalized by dividing each value by 256. No thresholding operation was
performed. The first row shows the human torso data set (resolution 512
* 512 * 114); the second row shows the lobster data set; and the third
row shows the hydrogen molecule data set.
Please click on the images to see them in full size.
Z slice = 1
|
Z slice = 30
|
Z slice = 60
|
Z slice = 90
|
Z slice = 114
|
Z slice = 3
|
Z slice = 10
|
Z slice = 17
|
Z slice = 23
|
Z slice = 30
|
Z slice = 20
|
Z slice = 26
|
Z slice = 32
|
Z slice = 38
|
Z slice = 44
|
Rendering with the marching cubes algorithm
I tested my implementation of the marching cubes algorithm on the lobster
data set. Some output images that were generated are shown below: the threshold
for this simulation was 200, so we can see only the much denser skeleta
portions of the lobster. The first figure shows the lobster seen from the
top, the second as seen from underneath and the third slightly rotated.
Please click on the images to see them in full size.
References
Watt, A. and Watt, M. (1996) "Advanced Animation and Rendering Techniques:
Theory and Practice". New York: ACM Press.
Implementation issues
The particle animation tool was programmed in C using the GLUT libraries
and OpenGL. The computer platform used was an SGI machine running on a
REALITY engine. The program executable can be downloaded by clicking
here. The code for this application was divided into the following
files:
-
Makefile: This makefile works for
on the REALITY ENGINE in my lab and should run on other SGI machines too.
-
const.h : Constants declaration for
the simulation.
-
funcs.h : Function declaration for
the simulation.
-
data.h : The header file for the file-reading
and writing functions in 'data.c'.
-
cube.h : The header file containing
cube-related data structures.
-
marching_cubes.c : This file
contains the main program for this simulation.
-
callbacks.c : Contains all the
callback registration functions for this simulation, including the display
and reshape functions, the idle function, keyboard callback and all menu
handling functions. This file also contains all the marching cubes functions.
-
data.c : Functions to read volumetric
data sets from files and memory management for the required data structures.