How to Shrinkwrap through a Critical Point: an Algorithm for the Adaptive Triangulation of Iso-Surfaces with Arbitrary Topology Andrea Bottino, Wim Nuij and Kees van Overveld Department of Mathematics and Computing Science Eindhoven University of Technology September 20, 1996
Abstract An algorithm is presented which generates a triangular mesh to approximate an iso-surface. It starts with a triangulation of a sphere and next applies a series of adjustments to this triangulation to transform it into the required surface. It is based on the shrinkwrap method ([Overveld 93]), but unlike the original shrinkwrap algorithm, it allows for arbitrary topological genus of the resulting iso-surface. As the original shrinkwrap algorithm, the new algorithm is adaptive in the sense that the lengths of the sides of the triangles in the mesh vary with the local curvature of the underlying surface.
keywords: iso-surfaces, implicit surfaces, polygonization, bifurcations
1 Introduction
()
For a scalar function V defined on IR3 , an iso-surface1 is the collection of points r where V r takes a given value, V r V0 . They play a significant role in computer visualisation. Various researchers have used skeletal iso-surfaces for modelling: ([Blinn 82], [Bloomentha 90], [Wyvill 86], [Nishimura 85]) and several of the surfaces that are of interest in pure mathematics, physics, or chemistry are iso-surfaces ([Wilhelms 91]). It is often necessary to convert the iso-surface into a polygon mesh which is rendered afterwards; the advantage of this approach is the availability of a full 3-D approximation of the iso-surface which allows for fast viewing from arbitrary directions, as well as the application of post-processing techniques such as mesh reduction ([Turk 92]), relaxation ([Mallet 92]) or free-form deformation ([Sederberg 86]). Also, the resulting polygon mesh is a closed manifold, so it may be used in B-rep-based CSG operations ([Requicha 83]). In this paper we focus on a recent technique to obtain such meshes, namely the shrinkwrap algorithm ([Overveld 93]). Shrinkwrap is adaptive to the local shape of the surface. This means that relatively flat regions are approximated with large triangles whereas regions with a smaller curvature radius are tessellated with small triangles. As a result, for a given geometric tolerance, the total number of required triangles is lower than could be obtained from e.g. a space partitioning method with a uniform voxel
( )=
Address: P.O.Box 513, 5600 MB, Eindhoven, The Netherlands. Email:
[email protected] 1 Also
called implicit surface
1
structure. Also, when shrinkwrap is used to construct a tessellation for an iso-surface, this tessellation will move along with the surface in the case of (smooth) animations, unlike tessellations that are constructed with uniform voxel structures. In order to make this paper self contained, in section 2 a short review of shrinkwrap is presented. The original shrinkwrap algorithm cannot deal with surfaces that are not topologically equivalent to a sphere. Therefore a technique is discussed in section 3 which is an extension to shrinkwrap. The results are summarised in section 4.
2
The shrinkwrap algorithm
First we define the class of iso-surfaces for which shrinkwrap should produce triangulations. Next an outline of the algorithm, together with an intuitive motivation are given, and some aspects are discussed in further detail.
2.1
The definition of the iso-surface
( )=
An iso-surface is the collection of points r 2 IR3 such that a given function V r V0 . We consider the i class of functions V where V r . The summation over i indicates that the function is comi jr,Ri j posed of a number of components with relative strength (weight) i . Components can be thought to be generated by different types of geometric primitives: points, line segments or convex polygons. In view of the obvious analogy with electrical charges we shall often refer to the function V as the ’potential’ . i If a component is generated by a point, then Ri is the location of this point, and the set jr,R V0 ij designates a sphere with radius i=V0 round Ri . The element can also be a line segment, say ab. In that case Ri depends on r; Ri is the projection of r onto ab and the set jr,iRij V0 designates a cylinder with hemi-spherical caps with radius i =V0 and ab as the axis. Similarly, elements can be triangles or other convex polygons in which case a projection of r onto that polygon has to take place in order to obtain Ri. In these cases the designated surface is an offset surface of the original polygon with cylindrically rounded edges. The collection of points, line segments and convex polygons that define V r is called the skeleton; each geometric object in the skeleton is called a skeleton element. The addition of the several components results in an iso-surface which is a smoothly blended union of the several iso-surfaces associated with the individual skeleton elements.
( )=P
=
=
()
2.2
Outline of the shrinkwrap algorithm
Shrinkwrap starts by setting V0 to a value close to 0, and providing a triangulation for the resulting isosurface. As can be easily seen, this is approximately a sphere with very large radius. This triangulation consists of approximately equilateral triangles. Next the value of V0 is increased, in a number of steps, and with each step the surface shrinks a bit towards its final shape and size. With every shrink step, the vertices of the triangulation are moved towards the new surface. This is done as follows: We make use of a first order Taylor expansion to compute a first estimate for the new location of a vertex when increasing V0 to V0 V , and, when necessary, iterate. Assume r is on the V V0 surface: V r V0. Next we look for a new location, r rV r , such that
( )=
+
+
()
V (r + rV (r)) = V0 + V:
2
=
= 0 gives: V (r + rV (r)) = V (r) + (rV (r) rV (r)) + O(jj2) = V0 + V;
Taylor expansion around
or
V = (rV (r) rV (r)):
This gives
=
V (rV (r) rV (r)) :
Also, to meet with accuracy and robustness requirements, discussed in [Overveld 93], it may be necessary to split edges. An edge is split if the angle between the gradients in the endpoints is larger than a certain threshold. If an edge is split the adjacent triangles should be split as well. It should be noted that a small angle between gradients will not assure a good planar approximation, but an error would not remain unnoticed: in the iteration unacceptable edges would be created that cannot be divided into acceptable ones. Initially, the set of vertices consists of the vertices in the initial object (a triangulated sphere with more or less equilateral triangles); this object is assumed to be sufficiently large to be outside the entire iso-surface even for iso-value V0 =Nsteps . The global structure of the shrinkwrap algorithm (in pseudo-Pascal) reads as follows:
=1
begin V0:=0; while V0 < 1.0 do begin V0:=V0+DeltaV; move all vertices to the surface; repeat check the acceptability of the edges; split unacceptable edges; split the adjacent triangles; until there are no more unacceptable edges end; end;
3 Extending shrinkwrap for topological changes
()
The shrinkwrap algorithm as described above assumes that the gradient rV x is not zero for all x such that V x . It is not a priori possible to deduce from a given skeleton if this will be the case; moreover, the applicability of shrinkwrap significantly improves if it can deal with arbitrary topological configurations. In order to devise the required modifications to shrinkwrap, we start by observing that a topological change or bifurcation occurs in a critical point, i.e. a point rc where rV rc . Such a bifurcation can be the occurence of a hole (as during the formation of a torus) or the occurrence of a rupture (e.g. when the iso-surface consists of several disjoint components). The proposed modification therefore consists of the following steps:
() 1
( )=0
( )=0
In principle, for every vertex r a search could be done to a nearby point rc with rV rc using a Newton Raphson iteration. Of course, it is useless to actually perform such an iteration for every vertex, and a condition is derived which can be used for a fast check to test if a critical point is near a given vertex r; 3
once a critical point is found, its function value is computed: Vc
= V (rc);
this value Vc , augmented by a small amount, say 0.001, is inserted in the list of iso-values that is used for the shrinkwrap iterations. By adding this small amount to the potential we are sure that the corresponding iso-surface has changed topologically while the triangulation is still sufficiently close to allow a correction. in this critical point, distinguish between a rupture and a hole. This distinction can be made on the basis of a local approximation of the function near rc as a quadratic polynomial in x, y , and z . The coefficients of this polynomial can be inspected in order to make the required distinction;
disconnect the edges to the vertex rc and re-connect them in order to form the right topological configuration;
proceed with the iteration to the next value of the potential V .
The above issues will be explained in more detail in the next sections.
3.1
Finding critical points
In principle, for every shrinkwrap iteration, we could start a Newton-Raphson iteration at every vertex position to see if a solution for rV r can be found in the neighbourhood. Of course, most of these tests will have a negative result, and a large amount of computational effort would be wasted. Therefore, we first propose a simple condition which states that no critical point is near and next we discuss how the Newton-Raphson iteration is executed for those vertices that don’t satisfy the condition.
( )=0
3.1.1 A sufficient condition for the absence of a critical point No topological changes occur as long as the shrinking iso-surface passes no critical point, that is a point where rV , so it is rather important to know whether there are such critical points. Moreover, we only sample the surface in a finite number of points, but we would like to be sure that the triangulated mesh based upon these points is a good approximation of the iso-surface. The two estimates we shall derive next will be helpfull. Suppose we are given a function V on IR3 , and the points a, b, and c on the iso-surface V V0. Let r be a point in the convex hull of a, b, and c and define the function ga by ga t V a t r , a , and similarly gb and gc . Using integration by parts we have:
=0
()= ( + (
V (r) = V (a) + rV (r) (r , a) , Using r
Z1 0
=
))
tga00(t)dt
= a + b + c, where + + = 1, addition results in V (r) = V (a) + V (b) + V (c) ,
Z1 0
(ga00 + gb00 + gc00)tdt:
The quantity jga00j is bounded by M jr , aj2, where M is a upperbound for the second derivative of V . (For the upperbound M we can take the largest eigenvalue of the Jacobian of rV ). A rather involved computation shows that
1 3
jr , aj2 + jr , bj2 + jr , cj2 L2 ; where L is the diameter of the triangle abc (this inequality is sharp), so we have 4
jV (r) , V0j 61 M L2:
Next we consider the rate of change of V in the direction normal to the triangle to prove that V attains the value V0 in the vicinity of the triangle abc, and that the gradient differs from zero. Let the normal vector to the plane abc be denoted by n. Then
rV (r + tn) n rV (a) n , M jr , aj , M tjnj; and consequently
1 6
1 2
V (r + tn) , V0 , M L2 + trV (a) n , tM L , M t2 :
Substituting t
= L=p3, we can simplify this last estimate to p
3) , V0 pL3 [rV (a) n , 2M L]: p Hence if rVp(a) n > 2M L then V (r + Ln= 3) , V0 is positive. In a similar way it is proved that V (r , Ln= 3) , V0 is negative, so somewhere inbetween V (r + tn) is equal to V0 . Obviously, in this case we also have rV (r + tn) n > 0, so there is a smooth surface V (r) = V0 near the triangle abc and no critical point in the neighbourhood. We conclude that a critical point can only be expected where the gradient in the endpoints of the edges is small. For the condition rV (a) n > 2M L to be V (r + Ln=
of practical use, we have to assume a value for M . Since M cannot be computed in closed form for arbitrary surfaces, a user–provided value should be used instead. An obvious application is a potential consisting of a sum of functions of =r-type. In this case bounds for the second derivatives of the contributing functions sum up to a (conservative) estimate for M.
1
3.1.2 Iterating towards a critical point In the course of the shrinkwrap algorithm we maintain a set S , consisting of all vertices v of the polygonal mesh with jrV v j < . This set can be made arbitrarely large by choosing a greater value of , and its elements can be taken as starting points for a Newton-Raphson method to find the critical points. This method gives a very efficient means of converging to a root, if there is a sufficiently good initial guess. It can also spectacularly fail to converge, indicating (though not proving) that no root exists nearby. The set S as defined before can be considered as a good initial guess. Newton-Raphson method tries to reach the roots of an equation by adjusting the argument, iterating the following formula:
()
xn+1
= xn , J ,1 rV (xn)
(1)
(see [Mizrami 82] and [Press 86a]) where J is the Jacobian of rV and x0 2 S . If this process converges a critical value Vc is found. This value is added to the list of function values for which a skeleton is computed.
5
3.2
Establishing the topological nature of the bifurcation
Using the eigenvectors of the Jacobian in the critical point, we can introduce a coordinate system wich describes clearly the behaviour of the iso-surface near the critical point. In these coordinates we can write:
V
2
2
2
= ( xa2 + yb2 , zc2 ) + Vc
(2)
If is negative, the iso-surface will change from a one sheet hyperboloid for V < Vc into a two sheet hyperboloid for V > Vc – a ’rupture’ occurs, while for positive the change is the other way around – a ’hole’ is formed. See fig 1,2, and 3. The plane z , which we will call the cutting plane, separates the splitting surfaces generated by the critical point ([Press 86b]). The definition of this plane is essential for the process of fixing the topology (see x3.3.2). In general, the Jacobian J of rV has to be computed by numerical differentiation, which is not a very stable process, but the sign of the determinant of J , and hence the sign of will be reliable.
=0
A
B
Figure 1: Iso-potential surface for V0 < Vc
A
B
Figure 2: Iso-potential surface for V0
= Vc
When reaching a critical value of the potential the iso-surface possesses a bifurcation point, and increasing the potential still a little bit more will result in a smooth surface of a different topological nature. 6
A
B
Figure 3: Iso-potential surface for V0 > Vc As a consequence some edges in the neighbourhood of the critical point ’lose their ground’: before there was a curve on the surface near the edge, but now these edges connect different components of the surface, or – in the case of a hole – edges connect two points on the same component but the curve on the surface connecting the same points has to make a large detour. In order to establish the nature of the bifurcation we have to find these edges. Let us call those edges red edges, to make the talking easy. In the neighbourhood of a critical point the potential can be approximated by a quadratic function 2. In [Nuij 95] we explain how we can determine whether two points lay on different parts of a two sheet hyperboloid or the edge connecting them crosses the interior of a one sheet hyperboloid. In the current implementation we compute the critical value Vc , and this makes the test simpler. Let ab be an edge of the triangulation, so V a V b > Vc . As V is quadratic, its restriction to the line through a and b has an extremal value in the point a b = . Hence in the case of a rupture we have: the edge ab connects b two components if and only if V a+ 2 < Vc . For the other case the test is more complicated, due to the fact that the distinction between an edge which is close to the surface and an edge which crosses a hole b cannot be made uniquely. See figure 4 Therefore we decide that an edge ab is red iff V a+ 2 < Vc and the gradients in the endpoints make a large angle, say more then = . Clearly, only edges with endpoints in the set S need to be considered.
( )= ( ) ( +)2 ( )
2
edge close to surface
?
( )
edge crossing the hole
Figure 4: Crossing the hole After determining the set of red edges surrounding the critical point in the structure, we can build a set composed by red triangles. A red triangle is a triangle of the mesh in which at least one of its edges is red. Removing the red triangles cuts holes in the triangulated surface; our task is to fix these holes so that the triangulation matches the new topological structure of the iso-surface. First we have to determine the borders that enclose the red triangles on the surface. To do that we introduce a neighbour relation between the vertices of the triangulation: two vertices are neighbours iff they are connected by a non-red edge of a red triangle. The vertices that have exactly two neighbours form two chains on opposite sides of the cutting plane. In the sequel we will denote these chains by T1 and T2. Due to numerical errors 7
not all vertices have either two or zero neighbours. How to correct the coloring in that case is described in [Bottino 95].
Revs Red edges
Figure 5: The kind of singularity is reflected in the structure of the red edges. When a rupture occurs, the two borders of the rupture are connected together only by red edges, and all red edges intersect the cutting plane. When the singularity is a hole, all red edges connect two vertices of the same border, and no red edge intersects with the cutting plane. For a sketch of the two situations we refer to fig 5.
8
3.3
Adjusting the topological structure of the mesh in a critical point
Now we know whether we have to make a hole in the surface or to tear it apart. There are two different ways to proceed from here. One thing in common for the two methods is that before creating a new topology, we have to destroy the old one. That is, we have to remove red edges, red triangles, and all the vertices that become isolated because of the removal of all the surrounding triangles. 3.3.1 Fixing a rupture For each border we have to connect together vertices belonging to the same subset. A straightforward way to do this is: 1. calculate the center of the border (Ci
= Pv 2T j
i
vj n
where n
= jTij )
2. shift Ci to the iso-potential surface 3. connect all the vertices of the border with it We will obtain n new triangles and 1 new vertex. This method produces a low number of new triangles but sometimes produces very stretched ones. An alternative method which gives more regular triangles is discussed in [Bottino 95]. 3.3.2 Fixing a hole Now the situation is a little more difficult. We have to connect vertices belonging to T1 only with vertices of T2. Here the connecting edges should not cross the hole. Therefore we try to make the angle between the crossing edge and the separating plane as large as possible. The main idea is to connect vertices following the same orientation in the neighbourhood chains T1 and T2. Since the two chains can have different number of elements, we link a vertex with the nearest vertex on the other chain.
4
Results
A quantitative analysis of the overload introduced by the new features is not straightforward, since no comparison can be made on the same instances with the earlier Shrinkwrap. So one must be satisfied with computing the mean increase in terms of IFEPTs (Implicit Functions Evaluated Per Triangle). The results are not yet comparable because the number of IFEPTs depends on the number of iterations done; the number of iterations is increased by two for each critical value. The overload introduced is mostly due to the application of the Newton-Raphson method, used to find the critical points. In Table 1 can be seen which percentage of the number of function evaluations is spent in the Newton-Raphson method. The old Shrinkwrap gives an IF EP T ' : ; now we have IF EP T ' : (+40.67%) when Newton computation is off and IF EP T ' : (+108.42%) when it is on. The examples used in table 1 contain usually only one critical point except points that consists of four equal point sources on the vertices of a square and that creates one hole in the center of the square and four ruptures on the edges of the square. Two observations can be made. First, the instances in the first three lines all generate holes; their relatively little overhead is due to the fact that a hole usually gives rise to a much lower number of red edges than a rupture. Second, the number of computations of the Newton-Raphson method penalizes
8 31 17 32
9
11 69 4
Instance 4lines 1torus 4points 2triangles 2points 1line-1point
Triangles 746 458 2366 640 1204 606
Newt.on 7668 5230 43618 11875 26152 14267
IFEPTs 10.28 11.42 18.44 18.55 21.72 23.54
Newt.off 6582 4475 32035 8107 16679 6969
IFEPTs 8.82 9.77 13.54 12.67 13.85 11.50
Difference -14.20 % -14.44 % -26.57 % -32.80 % -36.22 % -48.85 %
Table 1: considerably the performances. This number depends mainly on the dimension of the set S ; by reducing this set, the overload introduced is reduced. The default threshold for S is maybe expensive, but has the advantage of being able to deal with any input instance, irrespective of its complexity.
5
Conclusion
We present an improvement of the Shrinkwrap algorithm. Now a large variety of final iso-potential surface can be polygonized even if we have to change the topology of the starting surface (the original sphere). All new features respect the basic philosophy of the primary Shrinkwrap, and they can be made ”transparent” to the user, since the list of the iso-values used to reach the final surface is automatically changed in order to take care of the critical points introduced by the instances.
References [Blinn 82]
James Blinn. A Generalization of Algebraic Surface Drawing. ACM Transactions on Graphics, 1:235, 1982.
[Bloomentha 90] Jules Bloomenthal and Brian Wyvill. Interactive Techniques for Implicit Modeling. Computer Graphics, 24(2):109–116, 1990. [Bottino 95]
Andrea Bottino. How To Shrinkwrap a Saddle Point. Master’s thesis, Eindhoven University of Technology, Eindhoven, the Netherlands, October 1995.
[Mallet 92]
J.L. Mallet. Discrete smooth interpolation in geometric modelling. Computer Aided Design, 24(4):178–191, April 1992.
[Mizrami 82]
Abe Mizrami and M. Sullivan. Calculus and analytical geometry. Wadsworth Publishing Company, 1982.
[Nishimura 85]
H. Nishimura, A. Hirai, T. Kawai, T. Kawata, I .Shirakawa, and K. Omura. Object Modelling by Distribution Function and a Method of Image Generation. Journal of papers given at the Electronics Communication Conference ’85, J68-D(4), 1985. In Japanese.
[Nuij 95]
W. Nuij. Bridging a gap with an old theorem, in: E.h.l. aarts, h.m.m. ten eikelder, c. hemerik, m. rem (eds), simplex sigillum veri: Een liber amicorum voor prof. dr. f.e.j. kruseman aretz. Technische Universiteit Eindhoven, ISBN 90-386-0197-2, pages 261–264, 1995. 10
[Overveld 93]
C.W.A.M. van Overveld and Brian Wyvill. Shrinkwrap: an adaptative algorithm for polygonizing an implicit surface. The University of Calgary, Department of computer science, Research Report No. 93/514/19, March 1993.
[Press 86a]
William H. Press, Brian P. Flannery, Saul A. Teukolsky, and William T. Vetterling. Numerical Recipes. Cambridge University Press, Cambridge, New York, Victoria, 1986. page 270–273
[Press 86b]
William H. Press, Brian P. Flannery, Saul A. Teukolsky, and William T. Vetterling. Numerical Recipes. Cambridge University Press, Cambridge, New York, Victoria, 1986. page 335–380
[Requicha 83]
A.A.G. Requicha and H.B. Voelcker. Boolean operations in solid modelling boundary evaluation and merging algorithms. Proceedings of the IEEE, 3(7):30–44, October 1983.
[Sederberg 86]
T.W. Sederberg and S.R. Parry. Free-form deformation of solid geometric models. Computer Graphics (Proc. SIGGRAPH 86), 20(4):151–160, Aug 1986.
[Turk 92]
Greg Turk. Re-tiling polygonal surfaces. Computer Graphics (Proc. SIGGRAPH 92), 26(2):55–64, July 1992.
[Wilhelms 91]
Jane Wilhelms and Allan van Gelder. A coherant Projection Approach for Direct Volume Rendering. Computer Graphics (Proc. SIGGRAPH 91), 25(4):275–284, August 1991.
[Wyvill 86]
Geoff Wyvill, Craig McPheeters, and Brian Wyvill. Data Structure for Soft Objects. The Visual Computer, 2(4):227–234, February 1986.
11
