Pavel BAXA+, Václav SKALA+*, Roman MOUCEK+
+ Department of Computer Science and Engineering
University of West Bohemia
Univerzitní 22, Box 314
30614 Plzen
CZECH REPUBLIC
baxa@kiv.zcu.cz, skala@kiv.zcu.cz, moucek@hermes.zcu.cz

* Media Technology Research Centre
University of Bath
Bath BA2 7AY
UNITED KINGDOM

Volume visualization algorithms can be divided into two basic categories: Direct volume rendering and surface fitting.

Direct volume rendering methods are characterized by direct mapping scene elements to the display unit. The most commonly used algorithm belonging to this category is ray casting. It is based on casting a ray through each pixel of a 2D display device to the voxel scene. As the ray travels through the scene, it intersects a sequence of cells, each having defined density (or opacity), until it leaves the scene or until the sum of the densities reaches some limit value. The sum of densities defines the final color of the corresponding pixel. The main disadvantage of the direct volume rendering methods is that each time the image is to be rendered, the entire dataset must be traversed [Yagel et al. 92].

Surface fitting algorithms make use of a completely different idea. With a common assumption that all the voxels inside an object have different value (e.g. density) range than those outside, the surface boundary of the object can be defined by an isosurface. An isosurface is a surface associated with a constant threshold value which separates voxels with values greater than or equal to the threshold from voxels with values less than it.

Marching cubes (MC) presented by [Lorensen-Cline87] is an example of a surface fitting algorithm. The basic notion is that for each eight neighbouring voxels in a dataset we define a cubic cell, vertices of which are represented by the eight voxels. If one or more voxels of the cell have values less than the user-defined threshold, and one or more of them have values greater than this value, we know the cell must contribute some component of the isosurface. There exist 256 different arrangements of voxels, since each of the eight cell vertices can lie on both sides of the isosurface. Due to geometrical symmetry (mirroring, rotations), this amount can be reduced to 15 basic configurations, see Fig. 1. By determining which edges of the cube cell are intersected by the isosurface, we can create triangular patches which divide the cube into regions within the isosurface and regions outside. By connecting the patches from all cubes on the isosurface boundary, we get a surface representation.

One obvious problem with marching cubes is ambiguity in creating triangle patches in case four or more intersection points are to create two or more neighbouring triangles. There is also a possibility of gaps in the isosurface.

Marching tetrahedra (MT) algorithm presented in [Bloomenthal88] is a modification of the marching cubes. The basic idea is that each cell is subdivided into several tetrahedra first, and then intersections of the isosurface and tetrahedra edges are computed. This partialy solves the ambiguity problems, because maximum number of intersections between the isosurface and tetrahedra edges is four, that is, at most two triangles will be created.

One of the most important features that any visualization algorithm must have is accuracy. The visualization process should represent the original data without any loss of information. If the visualization algorithm uses any kind of approximation, the error of such an approximation has to be minimized. This is extremely important when medical data are to be visualized. We decided to determine how the approximations used in MT algorithm affect the accuracy of visualization, and find a way how to minimize the overall error of the results. The error is evaluated as a volume difference between objects and their representation by MT algorithm. A special attention is devoted to small objects of voxel size.

Next section of the paper describes several commonly used tetrahedronization schemes and presents some new ones, which were supposed to give good results. Then, the theoretical error estimation is shown in two-dimensional case for simplicity. The idea is afterwards extended to 3D space. Finally, the estimated values are compared with experimentally achieved errors and the test results are presented.

One of the essential elements of the MT algorithm is a tetrahedronization scheme, i.e. a way how the cell is subdivided into tetrahedra. The common schemes [Kolcun94] divide the cell into five or six tetrahedra, see Fig. 2.

The commonly used 5-6 tetrahedra schemes have several disadvantages. The tetrahedra are uneven, mutualy different and asymmetric, result of which is that the generated boundary representation is not visually acceptable.

We propose three tetrahedronization schemes, that fulfil the following requirements:

• The subdivision is symmetrical.
• There is no information loss in comparison with the original voxel scene. This is very important in medicine, where no changes of original data are permitted.
• Tetrahedra are at the most two shapes: the less different shapes the tetrahedra have, the more simple is a possible hardware implementation.

The new subdivision schemes are shown in Fig. 3. In addition to the values in the original eight cell vertices, some new vertices are used. In Fig. 3a, the centers of each cell side are given a value that is an average value of the nearest 4 vertices values. In Fig. 3b, also the value of the cell’s center of gravity is computed as an average value from the known 8 vertices as a result of trilinear interpolation in the cell. Subdivision into 48 tetrahedra can be achieved from the scheme in Fig. 3b by adding centers of the cell edges and thus splitting each tetrahedron in two.

Such interpolations have no influence on the eight known values. This fact is important especially when this approach is used for medical data visualization purposes. Since we have no information about values between the cell vertices, linear interpolations seem to be suitable enough. Thus, instead of five or six tetrahedra schemes used in the original algorithm, we use 24 and 48-tetrahedral subdivision.

Since the problem of an error estimation in three-dimensional space is extremely difficult to present in a simple, comprehensive way, we decided to reduce it to 2D space for the presentation purposes, so that it can be easily understood from the drawings.

First, it is necessary to define assumptions needed for the estimation. Let the visualized object be defined as a circle with radius r. The scene is a continuous function h defined as a distance from the circle boundary.

The discretized scene is a matrix of points having coordinates i, j and value h_i,j,

The circle center has coordinates [0,0]. The MT algorithm is used to find the boundary of the circle, i.e. isosurface with h=0. Then, the real circle surface area S_R=p r² is compared with the surface area S_D of the approximating polygon, where the relative error e is evaluated as

The most simple case to explain can be seen in Fig. 4a. The circle radius is less than 1, an object of one-voxel size is approximated. Only one quadrant is shown for simplicity. The MT algorithm uses a simple triangularization scheme where each 4-voxel cell is divided by one of its diagonals into two triangles. The aim is to find a point x_p such that h(x_p ,0) = 0. Using linear interpolation on the edge h₀₀-h₁₀, the value of x_p is computed as follows:

Now, it is possible to compute the surface area of both the real body (circle) and the approximated one. Also we can estimate an error e.

Fig. 4b shows similar situation, but the other diagonal is used for division this time. It is necessary to find a point u_p on the diagonal that approximates the intersection point between the circle and the diagonal, while x_p stays unchanged.

The x- and and y-coordinate of u_p is computed using linear interpolation on the edge h₀₀-h₁₁:

It is obvious that u_p lies on the circle boundary. The error in this case is

Thus, the absolute value of the error alternates between 10% and 36%, depending on used diagonal. We presumed that such asymmetry can be avoided when some different triangularization scheme is used. In 2D case, this reduces to scheme shown in Fig 4c. Both cell diagonals are used, so that the cell is subdivided into four triangles. New point C lies in the center of the cell and its value h_C is computed as mean value of the four voxels.

The coordinates of u_p has changed, since it is computed using linear interpolation on the line segment h₀₀-h_C.

This error value is less than the previous, still the average error is surprisingly greater than in the previous case.

Let us suppose now that the circle radius lie in the interval . In Fig. 5, the same cell is shown. The circle intersect the cell edges in points x_r, y_r . In this case,

It can be shown that x_p - x_r < 0 for any r from the interval . If the interpolation on the diagonal is computed, as in Fig. 5b, the intersection point u_p has obviously coordinates

In order to determine the surface area of the approximated circle, also the neighbouring cells must be taken into account, as presented in Fig. 6. Situation in Fig. 6a) again shows the most simple case, where the diagonal has no influence on the approximation.

In this case, the values in the four voxel vertices evaluate as follows:

It can be shown that

.

Even in such a simple case that was shown here, the resulting error cannot be easily estimated, since it varies dependingly on the value of r. Still, there exists a common rule saying that a difference between real circle surface area S_R and surface area of its approximation (gray-shaded in figures) is biggest in cases shown in Fig. 5a, 6a, and smallest in Fig. 5b, 6b, respectively. It means that the finer subdivision of the square cell gives worse results. This effect, which may be surprising a bit, is caused by the added point h_C. This point was given a value that had been linearly interpolated from the cell vertices. Additional interpolation resulted in worse object approximations (from surface-area point of view) in 2D case.

Now let us imagine analogous problem in three-dimensional space. Square cells are replaced with cubes (cubic cells), a volume error of an approximation is determined instead of a surface area error. The same way we divided square cells into two triangles with no additional point or four triangles when h_C had been added, now the cubic cells can be subdivided into 5 or 6 tetrahedra without adding points, 12 (cell‘s center of gravity is added), 24 (centers of cell sides added) or even 48 tetrahedra (edge centers added).

Since the particular situations in three-dimensional space are very difficult to imagine, let us have a look only at the easiest one. The radius of the visualized sphere is and cube division into 5 tetrahedra is used. Since this scheme is asymmetrical, both mutually symmetrical cases have to be considered - see Fig. 7a, b. The schemes consisting of more tetrahedra could not be graphically presented here because of their complexity.

It can be proved that the approximation error is 68.2% in case a), 29.7% in case b) and finally 54.4% for selected 24 tetrahedra subdivision. So again we achieved similar theoretical estimation to 2D case. The finer subdivisions give worse results, when computing relative volume error of approximation.

Regarding the volume error measuring, we generated 200 voxel scenes with sphere radius increasing from 0.0 to 20.0 units by 0.1 and the voxel values represented by both integer and real numbers. For each scene, the real sphere volume was computed, as well as the volume of the approximations. The relative error was computed as

.

The results of the experiments are presented in the Chart 1.

The results of the test tend to the same behaviour that was estimated theoretically in previous sections. Concerning only the volume accuracy, it proves to be unnecessary to implement tetrahedronization schemes containing 24 or even 48 tetrahedra. The idea to add new points to the existing eight cell vertices does not seem to be good enough. The new points bring some error, since the linear interpolations are used to compute their values.

Still, the images generated by MT algorithm are visually obviously better for 24 and 48 schemes, as shown in next section. The reason of this has to be sought in features of human vision. A human eye is not able to percept volume accuracy, but is sensitive to edges and changing gradient of the surface. We calculated that some different criterion should explain why better results were achieved with 24 and 48 tetrahedra subdivision.

We also implemented a similar test, but this time we measured the surface area of the sphere and its approximations instead of volume. The results that we achieved can be seen in Chart 2. They de facto correspond with previous test, i.e. the finer subdivision does not seem to be advantageous when surface area precision is requested.

The implemented MT algorithm was tested not only for volume and surface area accuracy, but also from the visual point of view. Fig. 8 and Fig. 9 show samples of 256x256x109 MRI scanned head visualized using different tetrahedronization schemes. The 24 scheme that we used gives results looking significantly better than commonly used schemes, especially when small details are rendered (see nose, ear and lips in Fig. 8 or detailed view of an ear in Fig. 9).

We presented a new approach to the estimation of errors in the construction of isosurfaces by surface-fitting algorithms. We proposed several tetrahedronization schemes that make use of added points to achieve better visual results. An error estimation was made for existing 5, 6 and 12 tetrahedra per a cube and compared with new 24 and 48 tetrahedra schemes. As a measure of the error we chose a relative volume error. Both theoretical estimation and experiments proved that it is useless to implement 24 or 48 tetrahedra schemes, when only the volume or surface area criterion is taken into account. In the future work, a new criterion will be tested, that will focus on measuring gradient changes of the surface, thus getting near to the human vision system. Besides, an adaptive MT algorithm is to be proposed, which will approximate the objects more accurately both in volume and surface area.

This work was supported by The Ministry of Education of the Czech Republic – projects PG 97193, PG 97185 and VS 97155.

Authors would like to thank Ms. Ivana Kolingerová and Ms. Michelle Dobson for their suggestions and for making the final corrections to the English text.