# JVRB - Journal of Virtual Reality and Broadcasting

##### Sections
Home / Topologically Accurate Dual Isosurfacing Using Ray Intersection

# Topologically Accurate Dual Isosurfacing Using Ray Intersection

1. Jaya Sreevalsan-Nair Institute for Data Analysis and Visualizationand Department of Computer Science
2. Lars Linsen School of Engineering and Science
3. Bernd Hamann Institute for Data Analysis and Visualizationand Department of Computer Science

## Abstract

“Dual contouring” approaches provide an alternative to standard Marching Cubes (MC) method to extract and approximate an isosurface from trivariate data given on a volumetric mesh. These dual approaches solve some of the problems encountered by the MC methods. We present a simple method based on the MC method and the ray intersection technique to compute isosurface points in the cell interior. One of the advantages of our method is that it does not require us to use Hermite interpolation scheme, unlike other dual contouring methods. We perform a complete analysis of all possible configurations to generate a look-up table for all configurations. We use the look-up table to optimize the ray-intersection method to obtain minimum number of points necessarily sufficient for defining topologically correct isosurfaces in all possible configurations. Isosurface points are connected using a simple strategy.

1. published: 2007-08-30

## 1.  Introduction

Given a scalar field discretized on a three-dimensional grid, isosurfaces represent the geometry of a trivariate function being equal to a constant scalar value. The original MC method used for isosurface extraction is a simple algorithm of marching through all rectilinear hexahedral cells of a volumetric grid and computing two-manifold isosurface for each cell independently, in a piecewise fashion []. In each cell, the isosurface points are computed as zero intersections on the edges which are linearly interpolating the scalar values at its endpoints. With the help of a look-up table, a triangular mesh approximation to the isosurface is obtained. The overall 256 possible topological cases depending on the configuration of the sign of vertices in a cell can be condensed to 14 by considering rotational symmetry and mirror cases, that is, cases obtained by interchanging the positive and negative signs.

However, while the look-up table of the original algorithm represents single topological result for each of the 14 cases, further research has proven more than one topological configurations in some cases. Chernyaev [] extended the case-table of 14 “sign-configurations” to obtain 33 “topological configurations” to deal with complex topologies like tunnels. This configuration set was further reduced to 31 distinct topological configurations in [], which are analytically discussed in []. We also use the asymptotic decider test [] and DeVella′s necklace test [] to categorize all the cases.

Recently, there is a lot of work done in another class of algorithms, called “dual” contouring algorithms. The motivation behind these methods has been to generate isosurfaces with better triangles as well as to “preserve features” within each cell. MC method generally misses the details in the cell interior, due to under-sampling of isosurface points by linear interpolation along the edges. Though this shortcoming has been corrected using repeated subdivision of the cells in the MC method, dual approaches eliminate the overhead of extra storage and processing of the subdivided cells. The “SurfaceNets” algorithm [], the “Extended MC” method [], and dual contouring using Hermite interpolation [] represent initial work done in this field. However, these methods do not address the issue of resolving and representing multiple disconnected isosurface components within one cell. Greß et al. [] extended dual contouring [] to use more than one point within each cell, using vertex splits. Dual contouring has been further extended to obtain thin walls in isosurfaces by using the calculated isosurface points from the primal grid to make a dual grid [] and then implement MC on the new grid. Nielson [] improved the MC-surface in a two-step procedure of generating a “MC-Patch” surface and a “MC-dual” surface, which is dual to the MC-patch surface. MC-patch surface is the MC-surface after eliminating edges of triangles within the cell.

We present a dual isosurfacing method based on ray intersection to determine isosurface points [] []. Since we are using a dual approach, our algorithm has the typical advantages of being able to extract special features, generating nice triangulations, and extending isosurface generation to adaptively refined grids without generating discontinuities. Moreover, our algorithm generates topologically correct isosurfaces including tunnel cases and multiple isosurface components per cell, computes exact points on the isosurface with respect to trilinear interpolation, and does not require any normal information.

For the computation of the ray-isosurface intersection, we determine the minimum number of rays necessary for finding sufficient number of points in each cell to accommodate cases of single components, “multiple components” (i.e., disjoint pieces) of the isosurface as well as complex features like tunnels. We generate a look-up table based on all 31 topological configurations of the cells []. Lopes et. al. [] had a similar analysis of 256 cases as ours, however, they stated that finding interior points in the cell did not suffice for good isosurface generation. In our work, we illustrate that a right choice of diagonals can help us find a good sample set of isosurface points to generate the surface.

We provide a brief overview of trilinear interpolation followed by the description of the ray-intersection method. We formulate certain observed properties of the cell diagonals used as rays, and enlist the rules used to optimize the ray-intersections. For polygonization, we cluster isosurface points in the neighborhood of each cell using connectivity information. Our simple strategy for connecting the points also guarantees that multiple isosurface components within each cell are rendered distinctly. We explain the algorithm for computing isosurface points and our simple polygonization strategy.

## 2.  Mathematical Foundation

Let F(x,y,z) be the trivariate scalar field defined over a structured rectilinear hexahedral grid, and f(u,v,w) be its parametric representation over the domain [xmin,xmax] x [ymin, ymax] x [zmin,zmax] ∈ R3 , i. e., .

With (ui,vj,wk) being the parametric representation of the vertices of the cell (or voxel), for i,j,k, ∈ {0,1}, the trilinear model f(u,v,w) is defined as

(1)

where Ui = |u-ui|, Vj = |v - vj|, and Wk = |w - wk|.

A hexahedral cell is represented by vertices with minimum coordinates (u0,v0,w0) = (0,0,0) and maximum coordinates (u1,v1,w1) = (1,1,1), respectively. A point interior to the cell has parameters (u,v,w) ∈ [0,1]3 . An isosurface associated with value I, in parametric form and set form can be written as

T(I) = {(u,v,w) : f(u,v,w) = I, 0 ≤ u,v,w ≤ 1}  (2)

For our ray-intersection approach, the diagonals are used as rays and points of intersection with the trilinear isosurface are computed. The equation of a ray r(t) from vertex V to vertex W is given by

r(t) = V + t(W - V), 0 ≤ t ≤ 1  (3)

where r(t), V and W are position vectors which can be represented as coordinates, or in the parametric form.

Let point p be the intersection of a ray and an implicitly defined isosurface. It should satisfy Equations ( 1 ), ( 2 ) and ( 3 ). These conditions reduce to a cubic equation in variable t, which corresponds to the parametric representation of p on the ray, i.e., p=r(t), 0 ≤ t ≤ 1. The cubic equation can be written as

G3t3 + G2t2 + G1t + G0 = 0.  (4)

The computation of its coefficients G0,G1,G2 and G3 is explained in Appendix A.

We solve Equation ( 4 ) using the analytical Cardan′s method [] in case the equation has points of inflexion or the numerical Newton-Raphson method [] otherwise. The intersection point(s) p is then computed using the real root(s) of Equation ( 4 ) in Equation ( 3 ).

The roots of Equation ( 4 ) depend on the behavior of the discriminant of the cubic or the reduced quadratic equation, summarized in Table 1 . Discriminant D3 for Equation ( 4 ) is evaluated based on the point of inflexion (xN,yN) []. We need to compute

(5)

If G3 = 0, then the cubic equation reduces to a quadratic one, and the quadratic discriminant D2 determines the nature of the roots.

(6)

.

Table 1. Roots of governing cubic equation of trilinear function depending on its degree.

 Degree Discriminant Nature of Roots Cubic D3 > 0 1 real, 2 imaginary (G3 ≠ 0) yN = D3 = 0 3 equal real yN ≠ D3 = 0 1 distinct real, 2 equal real D3 < 0 3 distinct real Quadratic D2 < 0 2 distinct real (G3 = 0, D2 = 0 2 equal real G2 ≠ 0) D2 < 0 2 imaginary

An intersecting ray can produce at most three isopoints, by virtue of the governing cubic equation. However, for most of the cases, either one or two solutions result.

## 3.  Terminology

The sign of a vertex of a cell is positive, when the scalar value at the vertex is greater than or equal to a given isovalue, and negative otherwise. In the following, vertices will be denoted as “(+) vertices” and “(-) vertices,” respectively.

We use cell diagonals to define rays, and depending on the sign of the vertices on the ray, we can have either same-sign ended rays (+/+ or -/-) or different-sign ended rays (+/- or -/+). In the former case, there are either zero or two solutions to the cubic equation for ray-isosurface intersection, and in the latter, there are either one or three solutions. Intuitively, an odd number of intersections (or zeros) causes a sign change, and an even number maintains the same sign at the endpoints. We denote edges with same-sign endpoints as “(+/+) or (-/-) edges,” and different-sign ended edges as “(+/-) edges,” see Figure 2 (a).

When a cell has more (+) vertices than (-) vertices, the cell is said to be a positive cell (“(+) cell”). Similarly, if the cell has more (-) vertices than (+) vertices, it has a negative sign and is called a negative cell (“(-) cell”). If both (+) vertices and (-) vertices are equal in number then, the cell is a neutral cell (“(0) cell”), with no sign. Examples of positive, negative and neutral cells are shownin Figure 2 (b).

A cell face with four same-sign ended diagonals and four (+/-) edges is called an ambiguous face. An ambiguous face separated with respect to a specific sign [] [] means a contour on the face is a rectangular hyperbola and the sign is that of the vertices in the quadrants containing the hyperbola, see Figure 2 (c). The sign, with respect to which an ambiguous face is separated, is determined using the asymptotic decider []. A face separated with respect to positive sign will be denoted as “(+) ambiguous-face,” and that with respect to negative sign will be denoted as “(-) ambiguous-face.” Equations for the asymptotic decider are given in Appendix A.

## 4.  Optimization

We exhaustively use the trilinear function for all possible configurations, and determine the topology of the analytical isosurface. We can observe from the analytical isosurface that there are cases where more than one diagonal of the cell intersect the same isosurface component. Since evaluating the intersection point requires a significant number of computational operations, we optimize our algorithm by finding just one point per disjoint piece of the isosurface within the cell. We analyze each of the 31 cases from the exhaustive MC case-table [] [] [] to find the diagonals used for intersection. A priori knowledge of diagonals for intersection helps us in easy computation of isosurface points. The choice of diagonals for intersection depends on the number of ambiguous faces of the cell, the sign of the ambiguous faces [], and the sign of the diagonals′ vertices. However, the condition of “one-point-per-disjoint-component” is relaxed in the case of complex features inside the cell, like tunnels, as sufficient points are needed to capture the complexity of the manifold.

We categorize the original 14 sign-configurations (numbered 0 to 13) of the MC-case table to five categories, based on their number of ambiguous faces. There can be only cases with zero, one, two, three, and six ambiguous faces, by virtue of the possible combinations of eight connected, signed vertices. We observe that cells with the same number of ambiguous faces show similar choices for optimal number of rays. We define rules to distinguish between different topological configurations, and decide which diagonals are to be used for intersection. The configurations are represented using numerical indices, which were used in earlier works [] [], of the form “x,” “x.y” or “x.y.z,” where “x” is the sign configuration, “y” the topological configuration, and sub-configuration “z”. Nielson [] showed analytically that there can be just three “levels of characterization” for the cell configurations. We use DeVella′s necklace test [] (which leads to a positive result exclusively in the presence of tunnels) in a few cases to distinguish its different topological sub-configurations. DeVella′s necklace test checks whether two vertices on a diagonal are internally connected. Relevant equations of the test are explained in Appendix A.

The following subsections explain the choice of the diagonals for each sign-configuration and their respective topological configurations. Table 2 lists the different types of diagonals in each of the sign configurations, and Table 3 summarizes the rules for the diagonals to be used for the five categories. In this section, for cases of analysis of signed cells, we will use the example of a (+) cell, and similar analysis can be extended to a (-) cell, by interchanging the signs of the diagonals, ambiguous faces and the cell.

Table 2.  Types of diagonals for the 14 sign configurations (for (+) cell, in case of non-(0) cells).

 Case Diagonal types Case Diagonal types 0 4 (-/-) 7 3 (+/-), 1 (+/+) 1 1 (+/-), 3 (+/+) 8 4 (+/-) 2 2 (+/-), 2 (+/+) 9 4 (+/-) 3 2 (+/-), 2 (+/+) 10 2 (+/+), 2 (-/-) 4 1 (-/-), 3 (+/+) 11 4 (+/-) 5 3 (+/-), 1 (+/+) 12 2 (+/-), 1 (+/+), 1 (-/-) 6 1 (+/-), 2 (+/+), 13 4 (+/-), 1 (-/-)

Table 3.  Rules used to choose diagonals in 31 topological configurations based on the number of ambiguous faces (for (+) cell, in case of non-(0) cells). Notations: “AF” implies ambiguous faces, “IP” means isosurface points, T stands for tunnel points, DS represents disjoint surface, and (+/-)c (in case 13) is the (+/-) diagonal through the common vertices of 3 (+) ambiguous faces and 3 (-) ambiguous faces, respectively.

 Cases Diagonals/Rays #IP. 0 AF: 1, 2, 5, For x (≠4),1 (+/-) 1 8, 9, 11, For 4.1, 1 (-/-) 2 4 (4.1, 4.2) For 4.2, 2 (+/+) 4 (T) 1 AF: 3 (3.1,3.2), For 3.1, 2 (+/-) 2 6 (6.1.1, 6.1.2, For 6.1.1, 1 (-/-) 2 6.2) For 6.1.2, 2 (+/+) 4 (T) For x.2, 1 (+/-) 1 2 AF: 10,12 For x.1.1, 1 (+/+) 1 (x.1.1, x.1.2, x.2) For x.1.2, 2 (-/-) 4 (T) (case of 2 (+) AF. For 10.2, 1(+/+) 1 for x.1.z) For 12.2, 1 (+/-) 1 3 AF: 7 (7.1, 7.2, For 7.1, 3 (+/-) 3 7.3, 7.4.1, 7.4.2) For 7.2, 2 (+/-) 2 For 7.3, 1 (+/-) 1 For 7.4.1, 1 (+/+) 2 For 7.4.2, 3 (+/-) 3 (T) 6 AF: 13 (13.1, For 13.1, 4 (+/-) 4 13.2, 13.3, 13.4, For 13.2, 3 (+/-) 3 13.5.1,13.5.2) For 13.3, 2 (+/-) 2 For 13.4, 1 (+/-) 1 For 13.5.1, 1(+/-)c 3 For 13.5.2, 1(+/-)c 1(DS) >& 3(+/-) & 3(T)

Cases without Ambiguous Faces. Sign configurations 0, 1, 2, 4, 5, 8, 9, and 11 have no ambiguous faces, see Figure 3 . We ignore case 0 in our diagonal analysis, as in this case, the cell does not intersect the isosurface. Except for case 4, all the others have just one topological configuration each, and one isosurface component each. It can be seen from the analytical trilinear surface, shown in Figure 3 , that cases 1, 2, 5, 8, 9, and 11 can use a (+/-) diagonal for intersection.

In case 4, two different topologies occur, depending on whether the vertices on the (-/-) diagonal are “internally connected”. In configuration 4.1, the contour forms disjoint surfaces, and in configuration 4.2, they are internally connected to form a tunnel. We use DeVella′s necklace test to distinguish between the two sub-configurations []. For configuration 4.1, the (-/-) diagonal is used for intersection, and for configuration 4.2, 2 of the 3 (+/+) diagonals are used to obtain multiple points for defining the tunnel.

Cases with One Ambiguous Face. Sign configurations 3 and 6 have one ambiguous face each, as shown in Figure 4 . If a (+) cell has a (+) ambiguous face (i.e., the same sign as the cell), a single isosurface component occurs; and if the (+) cell has a (-) ambiguous face, two isosurface components occur. 3.1 and 6.1 are configurations where there are two disjoint isosurface pieces, and 3.2 and 6.2 are configurations with a single isosurface piece. In the case of 6.1, a same-sign ended diagonal intersects the isosurfaces, and the vertices on the diagonal may be internally connected [] to form a tunnel. Thus, there are two sub-configurations 6.1.1 and 6.1.2, based on the absence and presence of a tunnel, respectively, and they can be distinguished using the DeVella′s necklace test. In the case of two isosurface components and no tunnels (cases 3.1 and 6.1.1), for a (+) cell, a (-/-) diagonal is used to obtain two intersection points. In case of 3.1, where there are no (-/-) diagonals, two (+/-) diagonals are used to obtain an intersection point each. For case 6.1.2, both (+/+) diagonals are used to obtain four intersection points on the tunnel. In the case of a single isosurface component, in cases 3.2 and 6.2, a (+/-) diagonal is used to determine an intersection point.

Cases with Two Ambiguous Faces. Sign configurations 10 and 12 have two ambiguous faces, see Figure 5 . In this case, there are two topological configurations based on the signs of the ambiguous faces. If both ambiguous faces are separated with respect to the same sign, then a same-sign ended diagonal of the cell can intersect two disjoint isosurface pieces (cases 10.1.1 and 12.1.1), or the vertices of the diagonal can be internally connected to form tunnels (cases 10.1.2 and 12.1.2). If the ambiguous faces are separated with respect to different signs, then a single isosurface piece occurs (cases 10.2 and 12.2). In the case of a single isosurface component (cases 10.2 and 12.2), a (+/-) diagonal is used to obtain a single intersection point. In the case of disjoint surfaces in a configuration with two (+) ambiguous faces, a (+/+) diagonal is used to obtain two intersection points. Similarly, in the case of two (-) ambiguous faces and disjoint surfaces, a (-/-) diagonal is used. In the case of a tunnel and two (+) ambiguous faces, two (-/-) diagonals (in case 10.1.2) or one (-/-) diagonal and two (+/-) diagonals (in case 12.1.2) are used. A similar extension is done for a tunnel and two (-) ambiguous faces.

Cases with Three Ambiguous Faces. Sign configuration 7 has three ambiguous faces, see Figure 6 . In a (+) cell, there are four possibilities of sign combinations of the three ambiguous faces, represented in set form as (1) (+,+,+), (2) (+,+,-), (3) (+,-,-), and (4) (-,-,-). In configuration 7.1, there are three (+) ambiguous faces, which leads to three disjoint surfaces, and three (+/-) diagonals are used as rays. In configuration 7.2 with two (+) ambiguous faces and one (-) ambiguous face, two disjoint surfaces are formed. In this case, two (+/-) diagonals are used through the (+) vertices in one of the (+) ambiguous faces. These rays are specifically used to ensure that one obtains one intersection point per surface. In configuration 7.3, with one (+) ambiguous face and two (-) ambiguous faces, one surface occurs, and any one of the (+/-) diagonals can be used. In configuration 7.4 with three (-) ambiguous faces, there are two disjoint surfaces. However, in this case the vertices on the (+/+) diagonal can be internally connected to form a tunnel. Thus, there are two topological sub-configurations, 7.4.1 and 7.4.2, with two disjoint surfaces and with a tunnel, respectively. In case 7.4.1, a (+/+) diagonal is used to obtain two intersection points. In case 7.4.2, all three (+/-) diagonals are used to obtain three points on the tunnel.

Cases with Six Ambiguous Faces. Sign configuration 13 has six ambiguous faces and four (+/-) diagonals, see Figure 7 . There are seven possibilities of sign combinations of the six ambiguous faces, which are, (1) (+,+,+,+,+,+), (2) (+,+,+,+,+,-), (3) (+,+,+,+,-,-), (4) (+,+,+,-,-,-), (5) (+,+,-,-,-,-), (6) (+,-,-,-,-,-), and (7) (-,-,-,-,-,-). However, this can be reduced to four cases, as (1) and (7), (2) and (6), and (3) and (5) are pairs of mirror cases.

In configuration 13.1 with all six (+) ambiguous faces, four disjoint surfaces are formed and all the diagonals are used. In configuration 13.2 with five (+) ambiguous faces and one (-) ambiguous face, three disjoint surfaces are formed. Except for a diagonal through one of the two (+) vertices in the (-) ambiguous face, all the other three (+/-) diagonals are used. Specific rays are used in order to generate one intersection point per surface piece.

In configuration 13.3 with four (+) ambiguous faces and two (-) ambiguous faces, there are two disjoint surfaces. We find the (+) vertex that is common to all three (+) ambiguous faces, and the (+/-) diagonal through this vertex is one of the two rays used for isosurface computations. Any one of the remaining three diagonals is used as the second ray. In case of equal number of (+) and (-) ambiguous faces, suppose A and B are two vertices of the cell, such that A is common to all (+) ambiguous faces and B is common to all (-) ambiguous faces. AB is a diagonal of the cell. There are two possibilities of this configuration depending on the signs of A and B.

In configuration 13.4, the sign of the common vertex is different from that of the faces (i.e., A is a (-) vertex and B, a (+) vertex), and in configuration 13.5, it is the same as that of the faces (i.e., A is a (+) vertex and B is a (-) vertex). In configuration 13.4, there is one isosurface piece, and any one of the four (+/-) diagonals can be used.

In configuration 13.5, there can be three disjoint surfaces, which can reduce to two, in case of internal connection. In configuration 13.5.1, there are three disjoint surfaces, and three intersection points are obtained from the (+/-) diagonal through the common vertex of the three same-signed ambiguous faces (i.e., AB). This is the unique case where a (+/-) ray leads to three solutions. In configuration 13.5.2, two disjoint surfaces occur, of which one is a tunnel. Hence, the diagonal through the common vertex of the three same-signed ambiguous faces (AB) is used as a ray. The remaining three diagonals are used as rays to find isosurface points for the tunnel. Thus, all four diagonals are used in case 13.5.2. However the difference in the diagonals used for the two surfaces is to be stored for connectivity during polygonization.

## 5.  Algorithm

Our algorithm involves two major steps: (1) isosurface point computation and (2) triangulating the isosurface points.

We generate a look-up table for the choice of diagonals in each topological configuration. This a priori information saves a lot of computation time. Our algorithm generates a list of all the heterogeneous cells, and evaluates the topological configuration of each cell. Based on its configuration, the diagonals specified in the look-up table are used to find out the isosurface points.

We associate each isosurface point with an edge of the cell, when the edge contains the end-point of the diagonal (used for the respective intersection) nearest to the isosurface point, as shown in Figure 8 (a). Certain configurations require additional associations with vertices to make appropriate connectivity, see Figure 8 (b). The isosurface point is then associated with all the edges containing its associated vertices. An edge is shared by four neighboring cells, and thus, it can get associated with three or four isosurface points (which are in different cells that share the edge). The isosurface points associated with an edge are connected to form one or two distinct non-overlapping triangles. This strategy is similar to the connectivity used in earlier dual isosurfacing works [] [], and it results in a triangular mesh. In cases of tunnels, we treat them as isolated cells with no neighbors and employ polygonization for its isosurface points independently [].

## 6.  Results and Discussion

We have applied our algorithm on a synthetic 2x2x2 dataset especially for the tunnel cases. Figures 9 , 10 , 11 , 12 , and 13 show examples of configurations with zero, one, two, three, and six ambiguous faces. They demonstrate how our algorithm generates isosurfaces which are topologically equivalent to the ones generated by trilinear interpolation. In some cases, MC algorithm fails to produce the same result, which is due to the fact that not all cases of ambiguity, especially the internal ambiguities, are resolved in the MC algorithm.

We have further applied the algorithm on a few of the standard datasets. Figure 1 and 14 show the resulting isosurfaces in comparison to isosurfaces extracted using a state-of-the-art MC method. Superposition of both surfaces show their similarity. Since our algorithm have memory limitations in the polygonization phase, for datasets bigger than 64x64x64, we have to downsample the datasets by a ratio, as shown in the cases of bonsai and skull datasets, as shown in Figure 14 . The computation times required to generate the points on the isosurface are similar to those required by the standard MC algorithm. For these examples, the performance statistics are shown in Table 4 . The results show that time performance depends on the number of tunnel cases.

Table 4.  Performance of the algorithm

 Dataset (Size, Ratio) Computation # Tunnels Time (in s) Neghip (32x32x32, 1:1) 0.96 0 Buckyball (32x32x32, 1:1) 1.67 0 Bonsai (256x256x256, 1:4) 2.36 42 Skull (256x256x256, 1:4) 4.08 67

The polygonization step still needs to be improved. Using efficient data structures like kd-trees [], for example, would improve the performance of this method tremendously by more robust representation of neighborhood. The current polygonization method is O(n) for a grid size n, which is not favorable for grid sizes greater than 64x64x64. Hence an optimized memory management is a necessary step for this method.

We plan to generalize our method to adaptive grids and irregular meshes. Other future directions include developing a meshless surface construction using the isosurface points directly. Moreover, we plan to analytically calculate the translation in the axis, using the analytical trilinear model [], and employ the designated diagonals for intersection.

## 7.  Conclusion

We have presented a dual isosurfacing method based on ray intersection to determine isosurface points. Apart from having the typical advantages of dual approaches (being able to extract special features, generating high quality triangulations, and extending it to isosurface generation in adaptively refined grids without generating discontinuities), our algorithm generates topologically correct isosurfaces, including tunnel cases and multiple isosurface components per cell, computes exact points on the isosurface with respect to trilinear interpolation, and does not require any normal information such as Hermite data.

## 8.  Appendix A

The parametric equation for the isosurface (Equation ( 2 )) in Section 2 can be rewritten as:

Fi(u,v,w,I) = Auvw + Buv + Cvw + Dwv + Eu + Fv + Gw + H,  (7)

where I is the given isovalue and Tijk = f(ui,vj,wk) denote the vertices of a cell, see Figure 2 (d). The values of coefficients are : A = -T000 + T010 - T110 + T100 + T001 - T011 + T111 - T101, B = T000 - T010 - T100 +T110, C = T000 - T010 - T100 + T110, D = T000 - T100 - T001 + T101, E = -T000 + T100, F = -T000 + T010, G = -T000 + T001, and H = T000 - I. Let (ur0 ,vr0 ,wr0 ) and (ur1 ,vr1 ,wr1 ) be the endpoints of the ray (in parametric form). The intersection point is given, in parametric form, by (u*,v*,w*) = (ur0 + tΔur,vr0 +tΔvr,wr0 +tΔwr), where Δur = ur1 - ur0 , Δvr = vr1 - vr0 , Δwr = wr1 - wr0 , and parameter t, 0 ≤ t ≤ 1, computed using the ray equation (Equations ( 3 ) and ( 4 )). The coefficients G0 , G1 , G2 , and G3 , of the cubic equation (Equation ( 4 )) are obtained by using Equation ( 7 ), Fi(u*,v*,w*,I) = 0. G0 = Fi(ur0 ,vr0 ,wr0 ,I), G1 = A(ur0 vr0 Δwr + vr0 wr0 Δur + wr0 ur0 Δvr) +B(vr0 Δur + ur0 Δvr) + C(wr0 Δrr + vr0 Δwr) + D(ur0 Δwr + wr0 Δur) + EΔur + FΔvr + GΔwr, G2 = A(ur0 ΔvrΔwr + vr0 ΔwrΔur + wr0 ΔurΔvr) + BΔurΔvr + CΔvrΔwr+ DΔwrΔur, and G3 = AΔurΔvrΔwr.

For the asymptotic decider test [], one uses bilinear interpolation to compute a value V = T00T11 - T01T10 - I(T00 + T11 - T01 - T10), for a given isovalue I, see Figure 2 (c). The test states that when V < 0, the isocontour is a rectangular hyperbola in the first and third quadrants, and when V > 0, then the isocontour is a rectangular hyperbola in the second and fourth quadrants.

DeVella′s necklace test states that if G[Fi] < 0 and Disc[Fi] > 0, then there exists a tunnel in the cell [] for values G[Fi] and Disc[Fi] obtained from the coefficients of the isosurface equation(Equation ( 7 )): G[Fi] = (AE - BD)(AF - BC)(AG-CD), Disc[Fi] = (AH)2 + (BG)2 + (CE)2 + (DF)2 - 2ABGH - 2ACEH - 2ADFH - 2BCEG - 2BDFG - 2CDEF + 4AEFG + 4BCDH.

## 9.  Acknowledgments

This work was supported by the National Science Foundation under contract ACI9624034 (CAREER Award), through the Large Scientific and Software Data Set Visualization (LSSDSV) program under contract ACI9982251, through the National Partnership for Advanced Computational Infrastructure (NPACI) and a large Information Technology Research (ITR) grant; the National Institutes of Health under contract P20 MH60975-06A2, funded by the National Institute of Mental Health and the National Science Foundation; and the U.S. Bureau of Reclamation. We thank the members of the Institute for Data Analysis and Visualization (IDAV) at the University of California, Davis.

## Bibliography

[Che95] Marching Cubes 33: Construction of Topologically Correct Isosurfaces Technical Report CERN CN 95-17,  CERN 1995.

[CHJ03] Iso-splatting: A Point-based Alternative to Isosurface Visualization Proceedings of the Eleventh Pacific Conference on Computer Graphics and Applications - Pacific Graphics 2003,  2003J. Rokne, W. Wang and R. Klein (eds.)pp. 325—334isbn 0-7695-2028-6.

[Gib98] Constrained Elastic Surface Nets: Generating Smooth Surfaces from Binary Segmented Data MICCAI '98: Proceedings of the First International Conference on Medical Image Computing and Computer-Assisted Intervention,  1998W. M. Wells, A. Colchester and S. Delp (eds.) Springer-Verlag London, UKpp. 888—898isbn 3-540-65136-5.

[GK03] Efficient Representation and Extraction of 2-Manifold Isosurfaces Using kd-Trees Proceedings of The Eleventh Pacific Conference on Computer Graphics and Applications - Pacific Graphics 2003,  J. Rokne, W. Wang and R. Klein (eds.) IEEE CS Press 2003pp. 364—376isbn 0-7695-2028-6.

[JLSW02] Dual contouring of hermite data Proceedings of the 29th annual conference on Computer graphics and interactive techniques - SIGGRAPH 2002,  2002Tom Appolloni (ed.) ACM Press pp. 339—346isbn 1-58113-521-1.

[KBSS01] Feature Sensitive Surface Extraction from Volume Data Proceedings of the 28th annual conference on Computer graphics and interactive techniques-SIGGRAPH 2001,  2001Eugene Fiume (ed.) ACM Press pp. 57—66isbn 1-58113-292-1.

[LB03] Improving the Robustness and Accuracy of the Marching Cubes Algorithm for Isosurfacing IEEE Transactions on Visualization and Computer Graphics 9 (2003)no. 116—29 IEEE Educational Activities Department Piscataway, NJ, USAissn 1077-2626.

[LC87] Marching cubes: A high resolution 3D surface construction algorithm Proceedings of the 14th annual conference on Computer graphics and interactive techniques - SIGGRAPH 1987,  1987M.C. Stone (ed.) ACM Press pp. 163—169isbn 0-89791-227-6.

[NH91] The asymptotic decider: resolving the ambiguity in marching cubes Proceedings of IEEE Conference on Visualization 1991,  IEEE Computer Society Press G. M. Nielson and L. J. Rosenblum (ed.)1991pp. 83—91.

[Nic93] A New Approach to solving the cubic: Cardan's solution revealed The Mathematical Gazette 77 (1993)354—359issn 0025-5572.

[Nie03] On Marching Cubes IEEE Transactions on Visualization and Computer Graphics 9 (2003)no. 3283—297issn 1077-2626.

[Nie04] Dual Marching Cubes VIS '04: Proceedings of the conference on Visualization '04,  2004Holly Rushmeier, Greg Turk and Jack Van Wijk (eds.) IEEE Computer Society Washington, DC, USApp. 489—496isbn 0-7803-8788-0.

[PFTV86] Numerical Recipes: The Art of Scientific Computing1st Edition Cambridge University Press 1986isbn 0-521-30811-9.

[PSL98] Interactive ray tracing for isosurface rendering VIS '98: Proceedings of the conference on Visualization '98,  1998David S. Ebert, Holly Rushmeier and Hans Hagen (eds.) IEEE Computer Society Press Los Alamitos, CA, USApp. 233—238isbn 0-8186-9176-X.

[SW04] Dual Marching Cubes: Primal Contouring of Dual Grids Pacific Conference on Computer Graphics and Applications,  2004D. Cohen-Or, H.-S. Ko, D. Terzopoulos and J. Warren (eds.)pp. 70—76isbn 0-7695-2234-3.