GRAPP 2006
Topologically Accurate Dual Isosurfacing Using Ray Intersection
urn:nbn:de:0009611700
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 lookup table for all configurations. We use thelookup table to optimize the rayintersection 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.
Keywords: Isosurface, Dual Contouring, Ray Intersection, Trilinear Interpolation
Subjects: Computer Graphics, Threedimensional Computer Graphics, Ray Tracing
Given a scalar field discretized on a threedimensional 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 twomanifold isosurface for each cell independently, in a piecewise fashion [ LC87 ]. 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 lookup 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 lookup 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 [ Che95 ] extended the casetable of 14 “signconfigurations” to obtain 33 “topological configurations” to deal with complex topologies like tunnels. This configuration set was further reduced to 31 distinct topological configurations in [ LB03 ], which are analytically discussed in [ Nie03 ]. We also use the asymptotic decider test [ NH91 ] and DeVella′s necklace test [ Nie03 ] 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 undersampling 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 [ Gib98 ], the “Extended MC” method [ KBSS01 ], and dual contouring using Hermite interpolation [ JLSW02 ] 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. [ GK03 ] extended dual contouring [ JLSW02 ] 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 [ SW04 ] and then implement MC on the new grid. Nielson [ Nie04 ] improved the MCsurface in a twostep procedure of generating a “MCPatch” surface and a “MCdual” surface, which is dual to the MCpatch surface. MCpatch surface is the MCsurface after eliminating edges of triangles within the cell.
Figure 1. Better triangles in isosurface obtained using dual contouring (left) than MC method (right) for the fuel dataset (from http://www.volvis.org). The dual contouring result is obtained from the isosurface points computed using our method.
We present a dual isosurfacing method based on ray intersection to determine isosurface points [ CHJ03 ] [ PSL98 ]. 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 rayisosurface 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 lookup table based on all 31 topological configurations of the cells [ LB03 ]. Lopes et. al. [ LB03 ] 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 rayintersection method. We formulate certain observed properties of the cell diagonals used as rays, and enlist the rules used to optimize the rayintersections. 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.
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 [x_{min},x_{max}] x [y_{min}, y_{max}] x [z_{min},z_{max}] ∈ R^{3} , i. e., .
With (u_{i},v_{j},w_{k}) 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
where U_{i} = uu_{i}, V_{j} = v  v_{j}, and W_{k} = w  w_{k}.
A hexahedral cell is represented by vertices with
minimum coordinates (u_{0},v_{0},w_{0}) = (0,0,0) and maximum
coordinates (u_{1},v_{1},w_{1}) = (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 rayintersection 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
G_{3}t^{3} + G_{2}t^{2} + G_{1}t + G_{0} = 0. (4)
The computation of its coefficients G_{0},G_{1},G_{2} and G_{3} is explained in Appendix A.
We solve Equation ( 4 ) using the analytical Cardan′s method [ Nic93 ] in case the equation has points of inflexion or the numerical NewtonRaphson method [ PFTV86 ] 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 D_{3} for Equation ( 4 ) is evaluated based on the point of inflexion (x_{N},y_{N}) [ Nic93 ]. We need to compute
If G_{3} = 0, then the cubic equation reduces to a quadratic one, and the quadratic discriminant D_{2} determines the nature of the roots.
.
Table 1. Roots of governing cubic equation of trilinear function depending on its degree.
Degree 
Discriminant 
Nature of Roots 
Cubic 
D_{3} > 0 
1 real, 


2 imaginary 
(G_{3} ≠ 0) 
y_{N} = D_{3} = 0 
3 equal real 

y_{N} ≠ D_{3} = 0 
1 distinct real, 


2 equal real 

D_{3} < 0 
3 distinct real 
Quadratic 
D_{2} < 0 
2 distinct real 
(G_{3} = 0, 
D_{2} = 0 
2 equal real 
G_{2} ≠ 0) 
D_{2} < 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.
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 samesign ended rays (+/+ or /) or differentsign ended rays (+/ or /+). In the former case, there are either zero or two solutions to the cubic equation for rayisosurface 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 samesign endpoints as “(+/+) or (/) edges,” and differentsign 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 samesign ended diagonals and four (+/) edges is called an ambiguous face. An ambiguous face separated with respect to a specific sign [ Che95 ] [ NH91 ] 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 [ NH91 ]. A face separated with respect to positive sign will be denoted as “(+) ambiguousface,” and that with respect to negative sign will be denoted as “() ambiguousface.” Equations for the asymptotic decider are given in Appendix A.
Figure 2. (a) Sign of ray/edge (b) Sign of cell (c) Ambiguous face with the notations used for vertices and their scalar values in bilinear interpolation. Solid contours occur when the face is separated with respect to positive sign, and dotted contours occur when with respect to negative sign (d) Notation for cell vertices for trilinear interpolation. Black and white vertices are positive and negative, respectively.
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 casetable [ Che95 ] [ Nie03 ] [ LB03 ] 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 [ NH91 ], and the sign of the diagonals′ vertices. However, the condition of “onepointperdisjointcomponent” 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 signconfigurations (numbered 0 to 13) of the MCcase 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 [ Che95 ] [ LB03 ], of the form “x,” “x.y” or “x.y.z,” where “x” is the sign configuration, “y” the topological configuration, and subconfiguration “z”. Nielson [ Nie03 ] showed analytically that there can be just three “levels of characterization” for the cell configurations. We use DeVella′s necklace test [ Nie03 ] (which leads to a positive result exclusively in the presence of tunnels) in a few cases to distinguish its different topological subconfigurations. 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 signconfiguration 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 subconfigurations [ Nie03 ]. 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.
Figure 3. Topological configurations with no ambiguous faces: trilinear model in yelloworange, isopoints in black, positive nodes in blue, rays in cyan.
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 samesign ended diagonal intersects the isosurfaces, and the vertices on the diagonal may be internally connected [ Nie03 ] to form a tunnel. Thus, there are two subconfigurations 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.
Figure 4. Topological configurations with one ambiguous face: trilinear model in yelloworange, isopoints in black, positive nodes in blue, rays in cyan.
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 samesign 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.
Figure 5. Topological configurations with two ambiguous faces: trilinear model in yelloworange, isopoints in black, positive nodes in blue, rays in cyan.
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 subconfigurations, 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.
Figure 6. Topological configurations with three ambiguous faces: trilinear model in yelloworange, isopoints in black, positive nodes in blue, rays in cyan.
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 samesigned 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 samesigned 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.
Figure 7. Topological configurations with six ambiguous faces: trilinear model in yelloworange, isopoints in black, positive nodes in blue, rays in cyan.
Figure 8. (a) Associations of isosurface point (in red) with vertex A (in blue), which is the endpoint of the intersecting diagonal closer to the isosurface (curved lines), and the edges “e” containing A; (b) Configurations where isosurface point is associated with vertex A on the intersecting diagonal, and additional associations between isosurface point and vertices, B, are made. Blue cell vertices are associated with the red isosurface point, the bold dotted line is the intersecting diagonal, and the curved boundary is the isosurface.
Our algorithm involves two major steps: (1) isosurface point computation and (2) triangulating the isosurface points.
We generate a lookup 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 lookup 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 endpoint 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 nonoverlapping triangles. This strategy is similar to the connectivity used in earlier dual isosurfacing works [ Gib98 ] [ JLSW02 ], 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 [ LB03 ].
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.
Figure 9. Trilinear surface (left), MC result (middle), and Dual Isosurfacing result (right) for 2x2x2 synthetic dataset with a tunnel and no ambiguous faces.
Figure 10. Trilinear surface (left), MC result (middle), and Dual Isosurfacing result (right) for 2x2x2 synthetic dataset with a tunnel and one ambiguous face.
Figure 11. Trilinear surface (left), MC result (middle), and Dual Isosurfacing result (right) for 2x2x2 synthetic dataset with a tunnel and two ambiguous faces.
Figure 12. Trilinear surface (left), MC result (middle), and Dual Isosurfacing result (right) for 2x2x2 synthetic dataset with a tunnel and three ambiguous faces.
Figure 13. Trilinear surface (left), MC result (middle), and Dual Isosurfacing result (right) for 2x2x2 synthetic dataset with a tunnel and six ambiguous faces.
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 stateoftheart 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 
Figure 14. Dual contouring result (left), MC result (middle), and comparisons of both surfaces by superposition for regular datasets of size 64x64x64 of Neghip (top) and Buckyball (second), and of size 256x256x256 (rendered in 1:4) of Bonsai (third) and Skull (bottom). (Datasets available at http://www.volvis.org.)
The polygonization step still needs to be improved. Using efficient data structures like kdtrees [ GK03 ], 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 [ Nie03 ], and employ the designated diagonals for intersection.
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.
The parametric equation for the isosurface (Equation (
2 )) in Section 2 can be rewritten as:
F_{i}(u,v,w,I) = A_{uvw} + B_{uv} + C_{vw} + D_{wv} + E_{u} + F_{v} + G_{w} + H, (7)
where I is the given isovalue and T_{ijk} = f(u_{i},v_{j},w_{k})
denote the vertices of a cell, see Figure 2 (d). The values
of coefficients are :
A = T_{000} + T_{010}  T_{110} + T_{100} + T_{001}  T_{011} + T_{111}  T_{101},
B = T_{000}  T_{010}  T_{100} +T_{110},
C = T_{000}  T_{010}  T_{100} + T_{110},
D = T_{000}  T_{100}  T_{001} + T_{101},
E = T_{000} + T_{100},
F = T_{000} + T_{010},
G = T_{000} + T_{001}, and
H = T_{000}  I.
Let (u_{r0
},v_{r0
},w_{r0
}) and (u_{r1
},v_{r1
},w_{r1
}) be the
endpoints of the ray (in parametric form). The
intersection point is given, in parametric form, by
(u^{*},v^{*},w^{*}) = (u_{r0
} + tΔu_{r},v_{r0
}+tΔv_{r},w_{r0
}+tΔw_{r}), where
Δu_{r} = u_{r1
}  u_{r0
}
, Δv_{r} = v_{r1
}  v_{r0
}
, Δw_{r} = w_{r1
}  w_{r0
}
, and
parameter t, 0 ≤ t ≤ 1, computed using the ray equation
(Equations ( 3 ) and ( 4 )). The coefficients G_{0}
, G_{1}
, G_{2}
,
and G_{3}
, of the cubic equation (Equation ( 4 )) are
obtained by using Equation ( 7 ), F_{i}(u^{*},v^{*},w^{*},I) = 0.
G_{0} = F_{i}(u_{r0
},v_{r0
},w_{r0
},I),
G_{1} = A(u_{r0
}v_{r0
}Δw_{r} + v_{r0
}w_{r0
}Δu_{r} + w_{r0
}u_{r0
}Δv_{r})
+B(v_{r0
}Δu_{r} + u_{r0
}Δv_{r}) + C(w_{r0
}Δr_{r} + v_{r0
}Δw_{r}) + D(u_{r0
}Δw_{r} + w_{r0
}Δu_{r}) + EΔu_{r} + FΔv_{r} + GΔw_{r},
G_{2} =
A(u_{r0
}Δv_{r}Δw_{r} +
v_{r0
}Δw_{r}Δu_{r} +
w_{r0
}Δu_{r}Δv_{r}) +
BΔu_{r}Δv_{r} +
CΔv_{r}Δw_{r}+
DΔw_{r}Δu_{r}, and
G_{3} = AΔu_{r}Δv_{r}Δw_{r}.
For the asymptotic decider test [ NH91 ], one uses bilinear interpolation to compute a value V = T_{00}T_{11}  T_{01}T_{10}  I(T_{00} + T_{11}  T_{01}  T_{10}), 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[F_{i}] < 0 and Disc[F_{i}] > 0, then there exists a tunnel in the cell [ Nie03 ] for values G[F_{i}] and Disc[F_{i}] obtained from the coefficients of the isosurface equation(Equation ( 7 )): G[F_{i}] = (AE  BD)(AF  BC)(AGCD), Disc[F_{i}] = (AH)^{2} + (BG)^{2} + (CE)^{2} + (DF)^{2}  2ABGH  2ACEH  2ADFH  2BCEG  2BDFG  2CDEF + 4AEFG + 4BCDH.
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 MH6097506A2, 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.
[Che95] Marching Cubes 33: Construction of Topologically Correct Isosurfaces, Technical Report CERN CN 9517, CERN, 1995.
[CHJ03] Isosplatting: A Pointbased Alternative to Isosurface Visualization, Proceedings of the Eleventh Pacific Conference on Computer Graphics and Applications  Pacific Graphics 2003, 2003, J. Rokne, W. Wang and R. Klein (eds.), pp. 325—334, isbn 0769520286.
[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 ComputerAssisted Intervention, 1998, W. M. Wells, A. Colchester and S. Delp (eds.), SpringerVerlag, London, UK, pp. 888—898, isbn 3540651365.
[GK03] Efficient Representation and Extraction of 2Manifold Isosurfaces Using kdTrees, 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, 2003, pp. 364—376, isbn 0769520286.
[JLSW02] Dual contouring of hermite data, Proceedings of the 29th annual conference on Computer graphics and interactive techniques  SIGGRAPH 2002, 2002, Tom Appolloni (ed.), ACM Press, pp. 339—346, isbn 1581135211.
[KBSS01] Feature Sensitive Surface Extraction from Volume Data, Proceedings of the 28th annual conference on Computer graphics and interactive techniquesSIGGRAPH 2001, 2001, Eugene Fiume (ed.), ACM Press, pp. 57—66, isbn 1581132921.
[LB03] Improving the Robustness and Accuracy of the Marching Cubes Algorithm for Isosurfacing, IEEE Transactions on Visualization and Computer Graphics (2003), no. 1, 16—29, IEEE Educational Activities Department, Piscataway, NJ, USA, issn 10772626.
[LC87] Marching cubes: A high resolution 3D surface construction algorithm, Proceedings of the 14th annual conference on Computer graphics and interactive techniques  SIGGRAPH 1987, 1987, M.C. Stone (ed.), ACM Press, pp. 163—169, isbn 0897912276.
[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.), 1991, pp. 83—91.
[Nic93] A New Approach to solving the cubic: Cardan's solution revealed, The Mathematical Gazette (1993), 354—359, issn 00255572.
[Nie03] On Marching Cubes, IEEE Transactions on Visualization and Computer Graphics (2003), no. 3, 283—297, issn 10772626.
[Nie04] Dual Marching Cubes, VIS '04: Proceedings of the conference on Visualization '04, 2004, Holly Rushmeier, Greg Turk and Jack Van Wijk (eds.), IEEE Computer Society, Washington, DC, USA, pp. 489—496, isbn 0780387880.
[PFTV86] Numerical Recipes: The Art of Scientific Computing, 1st Edition, Cambridge University Press, 1986, isbn 0521308119.
[PSL98] Interactive ray tracing for isosurface rendering, VIS '98: Proceedings of the conference on Visualization '98, 1998, David S. Ebert, Holly Rushmeier and Hans Hagen (eds.), IEEE Computer Society Press, Los Alamitos, CA, USA, pp. 233—238, isbn 081869176X.
[SW04] Dual Marching Cubes: Primal Contouring of Dual Grids, Pacific Conference on Computer Graphics and Applications, 2004, D. CohenOr, H.S. Ko, D. Terzopoulos and J. Warren (eds.), pp. 70—76, isbn 0769522343.
Fulltext ¶
 Volltext als PDF ( Size 10.5 MB )
License ¶
Any party may pass on this Work by electronic means and make it available for download under the terms and conditions of the Digital Peer Publishing License. The text of the license may be accessed and retrieved at http://www.dipp.nrw.de/lizenzen/dppl/dppl/DPPL_v2_en_062004.html.
Recommended citation ¶
Jaya SreevalsanNair, Lars Linsen, and Bernd Hamann, Topologically Accurate Dual Isosurfacing Using Ray Intersection. JVRB  Journal of Virtual Reality and Broadcasting, 4(2007), no. 4. (urn:nbn:de:0009611700)
Please provide the exact URL and date of your last visit when citing this article.