Shrinkwrap: An efficient adaptive algorithm for triangulating an iso-surface.
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 deformations to this triangulation to transform it into the required surface. These deformations leave the topology invariant, so the final iso-surface should be homeomorphic with a sphere. The 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. A quantitative analysis of the accuracy of the algorithm is given along with an empirical comparison with earlier algorithms.
Introduction
For a scalar function V defined on IR 3 , an iso-surface 1 is the collection of points r where V (r) takes a given
value,
. Iso-surfaces 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]).
In order to visualise iso-surfaces, ray tracing is an often used technique ([Nishimura 85]) which produces high-quality
images. In the case of all but trivial cases, however, it is an arduous task to compute points on the
iso-surface and the computational burden of ray tracing often imposes restraints on its application.
An alternative to ray tracing is the conversion of the iso-surface into polygon meshes which are rendered after-
wards; an additionaladvantage 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]).
Also called implicit surface
Most currently existing techniques for the polygonisationof iso-surfaces are based on data structures that allow
spatial indexing: either a voxel-based structure ([Bloomentha 88]) or the hash-table structure of ([Wyvill 86])
may be used.
Some inherent disadvantages of these data structures exist:
First, the data structure comprises a partitioning of the space rather a tesselation of the surfaces to be polyg-
onalised. Especially in the case of animation (e.g. in the computer animation "The great train rubbery",
[Wyvill 88]), this is likely to cause geometric artifacts that are fixed with respect to space, thus moving in
an incoherent way over every moving surface.
Second, there is an apparent mismatch between the number of triangles that is generated by these algorithms
and the complexity of the surface that is approximated: even relatively smooth and flat segments of an iso-surface
usually result in large amounts of facets. Bloomental ([Bloomentha 88]) uses an adaptive version of
the spatial indexing data structures, an octree in order to reduce the amount of polygons produced in tesselating
an iso-surface. This indeed reduces the amount of polygons generated, but full advantage of large cells can
only be taken if the flat regions of the surface happen to fall entirely within the appropriate octants. The same
observation applies to the approach to adaptive tesselation using the geometry of a tetrahedron rather than a
cubical structure.
Third, in earlier approaches using space partitioning with packed cubes, a value is sampled at each cube vertex
and a check is made to see if the vertex is inside or outside the surface ([Lorensen 87]). A table determines
the polygons which are used to replace each cube. Unfortunately this method has various ambiguous cases
where there are alternative configurations for certain in-out combinations. Ning and Bloomenthal ([Ning 93])
explore this problem and point out tht each cube may be subdivided in 5 or 6 tetrahedra which may be replaced
unambiguously by triangles. However, such an algorithm produces many unnecessary triangles. The
algorithm presented in this paper does not suffer from such ambiguities.
In this paper, a different approach to the tesselation of a class of iso-surfaces is taken. This class is defined in
2.1. The proposed approach is adaptive to the local behavior of the surface rather than being imposed by an
octtree with a priori defined cutting planes; this causes the tesselation to move along with the surface in the
case of (smooth) animations. Finally, as a side effect of the algorithm we propose, a coordinate system can
easily be introduced on the surface which can be useful e.g. for applying surface texturing.
In section 2 the intuition behind the algorithm is sketched, and the algorithm proper is presented. It makes use
of the notion of "acceptable edge", and this will be discussed in a more quantitative manner in relation to the
accuracy of the algorithm in section 3. The results are summarised in section 4.
An algorithm to arrive at an adaptive triangulation for an iso-surface
First we define the class of iso-surfaces for which our algorithm should produce triangulations. Next the al-
gorithm, 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 in 3-D space such that a given function V
. We consider the
class of functions V where V
. The summation over i indicates that the function is composed
of a number of components with relative strength (weight) ae i . Components can be thought to be generated by
different types of geometric primitives: points, line segments or convex polygons.
If a component is generated by a point, then R i is the location of this point, and the set ae i
designates
a sphere with radius ae i =V 0
round R i .
The element can also be a line segment, say ab. In that case R i depends on is the projection of r onto ab
and the set ae i
designates a cylinder with hemi-spherical caps with radius ae i =V 0
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 R i . In these cases the designated surface is an offset surface of the
original polygon with cylindrically rounded edges. Another way to envisage the designated surface is as the
Minkowsky sum of a sphere with radius ae i =V 0
and the geometric object (point, line segment or convex poly-
gon).
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 skeletal element ([Bloomentha 90]).
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.
2.2 Intuitive introduction of the shrinkwrap algorithm
The basic idea that underlies the algorithm is that the iso-surface is to be sampled with sufficient density in
order to capture all shape detail; further, that these samples are connected by edges to form a triangular mesh.
In voxel-based methods, such edges result from intersecting the surface with voxel boundaries, and hence
their directions all lay in one of three orthogonal planes (the XY, YZ or ZX planes of some world coordinate
system). These directions have no apparent relations with the iso-surface, and therefore many unnecessarily
short edges, and consequently, small triangles, result. In the algorithm proposed here, the edges should adjust
themselves to the shape of the surface. This means that e.g. in a cylindrical part of the surface, there should be
relatively long edges directed more or less along the axis of the cylinder and relatively short edges in directions
perpendicular to the axis. This allows for an adaptive triangulation. In order to let edges gradually adjust
themselves to the shape of the surface, we develop an iterative approach where the surface develops in several
steps from a sphere to the final shape. Indeed, the original shape (a sphere) has uniform curvature, and hence
adaptivity plays no role there; the triangulation of a sphere in order to achieve an initial estimate for the mesh
topology is simple. There are two natural ways to have an iso-surface develop itself from a sphere: One method
operates by first collapsing the entire skeleton into a point and gradually expanding back ("inflating") to the
desired skeleton geometry. The second method ("shrinking"), which allows a slightly simpler mathematical
analysis, and which will be chosen therefore for our algorithm, is based on the following observation.
Consider color plate 1. It depicts a 2-dimensional cross section through a distribution of some point skeleton
colors indicate the function value V (r) in a point. The iso-contour we are interested in is at the
boundary between the yellow and the magenta regions. We observe that iso-surfaces
with have a shape which is less involved, whereas iso-surfaces with are more involved. An
extreme case of the first example is which produces a sphere with an infinitely large radius.
So an algorithm could start by setting V 0
to a value close to 0, and providing a triangulation for the resulting
sphere. This triangulation should consist of approximately equilateral triangles, since the curvature is the same
anywhere.
Next the value of V 0
should be 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 should move towards
the new surface. This is discussed in section 2.3. Also, in order to meet with accuracy requirements, to be
discussed in 3.3, it might be necessary to split edges and triangles. But we note that edges and triangles are
only split when necessary, so there should be fewer unnecessarily small triangles at the end of the porcess than
with voxel-based methods.
Of course, this gradually shrinkingof the triangulation will only work as long as the topological structure of the
surface stays equivalent (homeomorphic) to the sphere that we start with. The issue of topologicalchanges during
the shrink process has been studied in ([Bottino 95]) and will be published elsewhere ([AB96]). However
to make this paper self-contained, we give a coarse outline of the proposed technique to deal with topological
changes (ruptures and holes) in Appendix B (topological changes).
2.3 getting the vertices onto the surface
Using the stepwise approach, and assuming the difference in iso-values V 0
fromone step to the next sufficiently
small 2 , we will use a Newton-Raphson method to displace vertices. So we make use of a first order 3 Taylor
expansion to compute a first estimate for the new location of a vertex when increasing
with an amount \DeltaV
to
\DeltaV , and, when necessary, iterate. Assume r is on the
. Next we look for
a new location, r
Taylor expansion round
or
\DeltaV  (ffi \Delta rV (r)):
Of course, this does not tell us in which direction the step ffi should be taken. A reasonable choice seems to be
to set
which gives
\DeltaV
So the new location is
\DeltaV rV (r)
(1)
2 and assuming that no topological changes occur between this step and the next
3 See section 4.2 for a proposal for a more sophisticated approach.
2.4 The shrinkwrap algorithm
We are now able to write down the total shrinkwrap algorithm. Vertices are defined as tuples (r; E; V; d) 2
the r-attribute is the location; E is the gradient rV (r); V is the value of the function
(r), and d is the displacement vector that should apply to this vertex in order to get it to the surface with next
higher iso-value.
Edges contain two references e1 and e2 to the vertices in the extremes, and two references to the two adjacent
triangles. Furthermore, an edge has a boolean n to indicate if it is non-acceptable (the acceptability of an edge
is discussed in section 3; for now it suffices to observe that edges should be in some way close to the underlying
surface in order to make quantitative statements about the accuracy of the triangulation; this can be obtained
by splitting edges that are non-acceptable).
The difference \DeltaV is an entire fraction of V 0
, say
the value will be
used. When using this convention, the interpretation of the value ae i is simply "the radius of the offset surface
if component i was the only one skeleton element".
The issue of the number of steps N steps in relation to the robustness of the algorithm is discussed in section
3.4.
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 V 0
1=N steps . All vertices v have pointing radially outwards, and v:d is
proportional to v:E=(v:E \Delta v:E) (as derived in 2.3).
The global structure of the algorithm (in pseudo-Pascal) reads as follows:
begin
while V0 < 1.0 do
begin
{ all vertices are on the surface;
for every vertex v we
have v.V=V(v.r)
v.E=grad V(v.r)
{.and all edges are acceptable}
{ all vertices are on the surface;
all edges are acceptable,
and V0=1.0 }
The statement S1 reads:
for all vertices v do
repeat
v.r:=v.r+v.d;
v.E:=grad V(v.r);
until |v.d|< Epsilon;
The statement S2 reads:
unlabel all triangles;
for all edges c do
begin
c.n:="c is non-acceptable";
{ Use the analysis of
section 3.2 and 3.3 to decide
acceptability.}
while edges are unacceptable do
begin
for all edges c with c.n=true do
begin
create a new vertex w for edge
c by splitting the edge;
move w to the surface similar as
in S1;
label the two triangles
adjacent to e; create two
new edges for c's
fragments, c1 and c2;
c1.n:="c1 is non-acceptable";
c2.n:="c2 is non-acceptable";
remove edge c;
{ all initially unacceptable edges
have been split, but new unacceptable
edges may have been introduced }
for all labeled triangles t do
begin
split triangle t
create 2, 3 or 4 new triangles;
unlabel the new triangles;
create 1, 2 or 3 new edges, ci;
for the new edge(s) ci, do ci.n:="ci
is non-acceptable";
remove triangle t;
{ all triangles bounded by initially
unacceptable edges
have been split,
but new unacceptable edges may
have occurred }
{ end of the while loop; all
edges are accpetable}
In the next two sections we discuss how to split edges and how to split triangles.
2.5 How to split edges
We have not defined yet what an acceptable edge is. For now it suffices to state that an acceptable edge should
be short enough so that the surface cannot bend away too much between the two extremes of the edge. Con-
versely, an unacceptable edge is an edge which is too long. So we see that the remedy to an unacceptable
edge is to split it, and to make sure that the new midpoint is again on the iso-surface. A naive way to do so is
depicted in the left column of figure 1 (figs. 1a-1d). The original edge is
in fig. 1a; M 1
)=2.
The array of dotted curves indicate the direction of the gradient of the function V (r) in the neighbourhood of
the iso-surface; the thick curve represents the iso-surface proper. If we move the point M 1
in accordance with
the local gradient, we arrive at M 0, as indicated by the dash-dotted arrow. Note that although M 0will be close
to the surface, since we only use a linear approximation for it will not lie on the surface in general.
In order to get it closer to the surface, we have to iterate as in statement S1 above. Now M 0will be a new
vertex. The edge
is very likely acceptable, but it may be much shorter than needed. On the other
hand, probably still unacceptable. As shown in fig. 1b, we therefore have to repeat the process on
edge
is the M 0from the previous phase), which yields M 0. As a result, we end up with a series of
unnecessary short edges, as depicted in fig. 1d. The main cause for this unfortunate behaviour is that we use
information about the geometry of the function V (r) and its gradients in the points M 1
, ., which
are possibly far from the surface. Evaluating the gradients in these points may yield misleading information
on the geometry of the iso-surface, causing a slow convergence and many unnecessary short edges.
A more efficient splitting strategy therefore should make use of reliable information only, that is information
in points that are already on the surface. In fig. 1e, the same configuration is shown as in fg. 1a. The dashed
thick curve is a curve which passes through the extremes A and B and is perpendicular to the normal vectors
nA and nB , respectively. Moreover, it is the smoothest curve with these properties; Appendix A contains a
derivation of an analytical expression for this curve. This curve serves to approximate a curve that lies in the
iso-surface
through A and B, i.e. the thick solid curve in the figure. Based on the dashed thick
curve, we propose point M , i.e. its parametric midpoint as a next point to evaluate the function's gradient.
Point M in fig. 1e is likely to be closer to the iso-surface than M 1
in fig. 1a, so the gradient computed in M
is likely to be more adequate to get acceptable edges than the gradient in M 1
. In this case,
might be acceptable (the edge fig. fig. 1f) might need one more
subdivision as depicted in fig. 1f.
2.6 How to split triangles
In case one or more edges are unacceptable, they have to be split. Fig. 2 shows a splitting scheme which
illustrates how a triangle can be subdivided into smaller triangles. In the top row one of the edges is subdivided;
the bottom row shows the case of three subdivided edges. In the case of two subdivided edges, two possible
schemes exist; in case we choose the first alternative; otherwise we choose the
second one.
3 Robustness and accuracy
In order to make quantitative statements about the behaviour of the shrinkwrap algorithm, we have to address
several issues:
ffl a. what assumptions have to be made about the iso-surface V
ffl b. based on the assumptions of (a), how can we assure that the shrinkwrap process captures all large-scale
structure of the developing iso-surface, in other words how can we prevent a situation like figure
3 happening (this regards the robustness of the algorithm);
ffl c. assuming we can guarantee robustness with respect to (b), how can we enforce bounds on the difference
between the triangulation and the underlying surface (this regards the accuracy of the algorithm);
ffl d. what is the minimal value for N steps ;
These items will be discussed in the subsequent sub-sections.
3.1 Requirements of the iso-surface
The shrinkwrap algorithm constructs a discrete model of a continuous object by means of sampling. This
means that the sampling density should be sufficiently high in order to capture all shape detail of the surface.
Since in shrinkwrap the samples are vertices of a piecewise planar mesh, all curved regions of the surface
require samples to be sufficiently close together. If the local curvature radius of the surface is known to be
nowhere smaller than a given value fi, then we can define a sample density (that is, the maximal radius of
a sphere round any sample such that no other sample lies within that sphere) that is guaranteed sufficient to
capture all surface detail.
So we assume that for every value of V 0;k
that is to be used as an iso-value in the k-th step of shrinkwrap, a
value of fi k is available such that the radius of curvature of the surface V
is nowhere less than fi k .
The value of the local curvature in a point r, V
can be computed using standard calculus, by fitting
a quadratic polygonomial V quad with suitable coefficient matrix
A which approximates V (r) in the neighbourhood of r, and deducing the curvature of this quadric using the
coefficients in A. However, for an arbitrary iso-surface V
it is in general not possible to compute the
globally smallest curvature radius analytically. Exhaustively sampling the space in the vicinity of the surface
to estimate this curvature would be computationally prohibitive. Instead, the values of fi k should be provided
beforehand by the user. The requirement that the values fi k be known beforehand seems to be impractically
restrictive. However, three observations apply.
ffl First, this requirement is not unique for shrinkwrap. The correctness of all purely point-sampling-based
methods for visualising iso-surfaces, such as ray tracing, uniform voxel-space algorithms and occtree-
or tetrahedron-based algorithms, relies on the assumption that the curvature radius of the surface is
bounded by a given constant. More sophisticated sampling techniques, such as interval artihmetic and
affine arithmetic may provide useful here, but these techniques apply in the context of shrinkwrap as
well.
ffl For a function V (r) as defined in 2.1 that consists of one single component, the fi k;i for that component
i can be straightforwardly given. Indeed, since the iso-surface is then an offset surface 4 with
radius ae i =V 0;k
defined on the skeleton element, fi
. Superposition of several elements
(all with positive ae i ) usually causes the minimal curvature radius to become larger, and in these cases
\Delta) is appropriate. Care should be taken, however for configurations such
as depicted in figure 4 where the value for fi k for the higher values V 0;k
should be taken smaller. Finally,
notice that the user should only be concerned about the value fi k for the value of V for the final
step. For an earlier step k 0 , we can set fi k
ffl In case the curvature radius of the surface should be locally smaller than the value for fi k provided by
the user, then this only deteriorates the quality of the approximation in the vicinity of this region with
high curvature: the rest of the surface (with sufficient large curvature radius) is not affected.
4 in other words: since it is the Minkowski sum of the skeletal element and a sphere with radius ae i =V 0;k
it follows that for convex
skeletal elements, the curvature radius is nowhere smaller than the radius of this sphere.
3.2 Capturing large-scale structure
Assume that the iso-surface has nowhere a smaller curvature radius than fi k . Consider a triangle ABC of
adjacent samples A, B, and C that are all on the surface. Let the gradients in these points be nA , nB , nC ;
they are normalised such that jn A 1. Consider the points a
and figure 5). These points will be called the adjoint points of A, B, and C , respectively.
Because the surface has nowhere a curvature radius less than fi k , it stays outside ("above" in figure 5) the three
spheres with radius fi k centered around a, b and c. Now if in the triangle abc the circle sectors with radius fi k
and centres a, b and c cover triangle abc (that is, no point of the triangle abc falls outside all three circle sectors),
the iso-surface cannot penetrate the triangle abc. Indeed: in order to penetrate the triangle abc there should be
a region in abc which is outside the spheres with radius fi k and centres a, b, and c as is illustrated for a 2-D
case in figure 5b. Now since the iso-surface cannot penetrate triangle abc, it can not move too far from triangle
ABC, because ABC and abc cannot be too far from each other. Hence the portion of the iso-surface that is
approximated by ABC is bounded from below by triangle abc. More precisely: since no point of ABC is
further away than fi k from triangle abc, the entire triangle ABC can be nowhere further away than fi k from
the iso-surface. In other words: there cannot be a portion of the surface that "escapes" between the samples
C as in figure 3.
Now a sufficient (although in general too restrictive) condition for abc to be covered by three circles with radius
round its corners a, b, and c is that
3 and similar for bc and ca, as follows from the geometry
of an equilateral triangle (see figure 5c).
A similar construction should be made for a triangle a 0 b 0 c 0 which is on the other side of ABC so that a
etc., to make sure the iso-surface cannot "escape" in the upward direction. So the portion of the
iso-surface approximated by ABC is bounded between abc and a 0 b 0 c 0 .
Thus in order to guarantee that the large-scale structure of the surface is captured, an edge AB should be
sufficiently 5 short that
3: (2)
3: (3)
This is the first condition (robustness) for acceptability for an edge. This condition should hold throughout
the shrinkwrap process for every iteration k. Notice that we used no assumptions about the iso-surface except
that its curvature is bounded by fi k , so these conditions are valid for all types of skeletal elements.
3.3
In this section we derive a condition for the length of an edge AB, such that no point in the interior of a triangle
ABC is further away than a predefined distance fflfi k , ffl !! 1, from the iso-surface. Similar as in the previous
section, we construct a plane that bounds the iso-surface from below (and one to bound the surface from
above). This time the plane has to be at distance fflfi k from ABC. This means that we have to introduce points
a", b" and c" with A is the angle between the gradient nA and the plane
through ABC. Similar for b" and c". The triangle a"b"c" is also intersected by the three spheres, centered
round a, b, and c with radius fi k , and we have to derive a condition such that the intersected sectors cover the
entire area of a"b"c" in order to enforce the iso-surface to be close to ABC . Therefore we consider the plane
ff through A and a, perpendicular to the plane through ABC . See figure 6. The point labeled m indicates
the intersection of ff with a"b"c" and the sphere round a, and it can be seen that a circular sector with radius
5 observe that this condition always can be satisfied: if we take A and B closer together, as is the case when splitting edges, nA and
nB have to get closer together, too.
round a" in a"b"c" falls entirely within the intersecting sector between a"b"c" and the sphere round
a with radius fi k . The distance can be computed from the cosine rule in triangle a"ma. We
have that
and
so
This gives for x:
sin
sin
or
sin
sin
By substituting for should find that so the +-sign option in 7 applies.
We note that the radius of the intersecting sectors is bounded by a number which is dependent of ffi A , and
this angle can of course be different from ffi B and ffi C . To facilitate our further analysis, we will assume that
derive a safe condition stating that the iso-surface does not penetrate a"b"c".
Notice that replacing ffi A , by their maximum means that we should actually compute the ffi M for
every triangle before we can apply the acceptability test.
Then again we have, from the geometry of the equilateral triangle, that ja" \Gamma b"j ! x
3 and similar for
b"c" and "c"a" where x follows from 7. A similar argumentation holds for the plane a 000 b 000 c 000 where a
So in order to guarantee that the surface is nowhere further away from the triangular mesh than fflfi k , an edge
AB should be sufficiently 6 short that
Now have arrived at the second condition for acceptability for an edge (accuracy). This condition should only
hold towards the end of the shrinkwrap process when V 0;k
approaches the final value of 1. A useful strategy
therefore is to decrease ffl gradually from 1 to the desired final value during iteration. This causes the first
iterations to take place with relatively few triangles.
6 observe that this condition also always can be satisfied for the same reasons as in section 3.2
Notice that conditions 9 and 10 are less critical than 2 and 3: failure of 9 or 10 does not jeopardize the global
structure of the triangulation. Therefore, instead of the global value for fi k , we may dare to use a local approximation
of fi k , which may be different for each triangle ABC , derived from the angles between nA , nB and
nC and the plane through A, B, and C. In this way, again triangles in relatively flat areas of the surface are
allowed to be larger which results in an even better adaptation. The experiments in the results-section were
made using this optimisation.
3.4 The number of steps
One parameter of the shrinkwrap algorithm has received no attention yet: the choice for N steps . On the one
hand, N steps should be as small as possible to give a small number of evaluations and hence an efficient al-
gorithm. On the other hand, a large value of \DeltaV may impede the Newton Raphson iteration in statement S1
from converging. Therefore a possible implementation could be to apply relatively large steps of \DeltaV and as
soon as a convergence problem occurs, go back one iteration and reduce \DeltaV .
However in practice it works better to take a value of \DeltaV slightly smaller than the maximally allowable value
for Newton Raphson convergence: even if convergence occurs, it may happen that the distribution of small
and large triangles is not very well adapted to the curvature of the surface if the vertices have to move a large
distance between each iteration and the next one (this is e.g. the case with the penguin in table 1 for 3 and
4 iterations). Moreover, due to edge splitting underway, it is obvious that during the first N steps \Gamma 3 or so
iterations significantly less vertices have to be updated than during the last few iterations. This means that the
amount of processing is very much concentrated in the last 3 or 4 iterations, whatever the number of iterations
is. So if N is a minimal value for N steps such that convergence marginally occurs, then N
much safer value which only costs a small amount of additional computing effort (this can be seen from table
1). This issue is discussed further in section 4.
4 Results, future research
The algorithm described has been implemented. In color plate 2, its qualitative behaviour is depicted: while
increasing
in 9 steps from 0.1 to 1.0, we see a penguin emerge from what starts of as a feature-less large
spherical shape. Note how gradually the details become visible, first the most protruded ones (for this reason
we equipped the penguin with an oversized bill), later on the more subtle ones. Color plate 3 shows the final
object together with the triangle mesh. Here, the adaptiveness of the algorithm is clearly visible, e.g. the bill
consists of mostly very slender triangles, whereas in the spherical part of the head we find more obtuse ones.
Also, in the concave regions, having a relatively high curvature, the triangle mesh has a much higher resolution
than elsewhere.
Some of the quantitative features are presented in the graph in figure 7. Here, we see the increase of the number
of triangles while iterating the penguin. It follows that the majority of the triangles is only created towards the
final (few) iterations. As a result, the efficiency of the algorithm is not as bad as one might expect from an
iterative approach: the work of the first few iterations is very small compared to the work in the final phase
of the process. For instance, when we apply 7 steps, the total number of function evaluations called is 7640,
whereas the final version contains 2152 triangles; that is 3.55 function evaluation per triangle. If we apply
smaller steps for V 0
, the result is of course worse; e.g. with 10 iterations we have 9389 function evaluations
and 2118 triangles in the final object, so 4.43 overhead.
So, as can be seen in table 1, the efficiency of the method is not extremely sensitive to the number of iterations
we use.
4.1 Comparison with earlier algorithms
In this version of the algorithm we do not attempt to reproduce the precise shapes as produced by other al-
gorithms; to provide a full comparison, the faithfulness to the surface of the polygonal approximation would
have to be measured. This could be done by taking samples for each triangle and computing the error from
the true surface. The sum of such an error could be averaged over all the triangles. Polygonisation algorithms
for iso-surfaces can be divided into two classes, adaptive and non-adaptive. For a given number of triangles
an adaptive algorithm should have a smaller error than for a non-adaptive algorithm.
In order to find the surface, polygonisation algorithms such as the one described in this paper, must evaluate a
function V (r) for a number of points in space. An exhaustive comparitive analysis is beyond the scope of this
paper, however it is likely (and certainly true for the algorithms we have tested) that the performance of such
algorithms depends on the number of such function evaluations (FE's), averaged over the number of triangles.
per triangle - FEPT).
Shrinkwrap was compared with an earlier algorithm as indicated below:
method 1: Shrinkwrap, as described in this document.
ffl method 2: Soft Objects, based on ([Wyvill 86]).
The voxel-based algorithm that we consider for comparison performs some processing to find a point of the
surface by intersecting the surface with one of the voxel boundaries. The point of intersection is then cal-
culated. Since this can be done by a variety of techniques (e.g. regula falsi, binary search) this part of the
algorithm should be the same for each, thus it is not counted in the total number of FEPTs. Another area
where the algorithms differ is the method chosen for computing the normal. In the algorithm of this paper the
normal is computed as a by-product of the search technique. The technique from ([Wyvill 86]) takes components
computed from V (X is 0:01 the length of a voxel side. In fact algorithm 2
without computing normals has an FEPT count about 20 times less than with the normal calculation.
Three typical modelling situations were chosen (a single sphere, where there is no adaptivity; 2 spheres that
blend together: and three blended line segments) to be tested along with a more complicated model, the pen-
guin. Although the primitives produce similar models for both methods, the penguin surface depends on the
style of blending chosen as well as the precise position for the final iso-surface. Thus figures for the penguin
can only be thought of as a guide.
The data was used by both algorithms and the following results computed and shown in Table 1. The numbers
are FEPTS (Function Evaluations Per Triangle).
When the penguin is triangulated with 3 or 4 iterations, we find that the bill is not completely triangulated due
to convergence failures with the Newton Raphson process. With 5 or more iterations, the triangulation is OK.
Notice that the number of FEPTS increases sublinearly with the number of iterations.
4.2 Future research
Experiments show that the version of the shrinkwrap algorithm as discussed in this paper seems to be a promising
approach to efficient adaptive triangulation of a relevant class of iso-surfaces. It is therefore interesting to
look at some obvious improvements of the algorithm.
ffl the skeleton elements need not necessarily be limited to points, line segments and convex polygons. In
any geometric primitive that allows a differentiable distance function is appropriate.
Algorithm Sphere 2 Spheres 3 lines penguin
Shrinkwrap 8.03 8.56 3.22 2.03 (N steps = 3; incomplete triangulation)
Shrinkwrap 2.46 (N steps = 4; incomplete triangulation)
Shrinkwrap 2.91 (N steps = 5; triangulation OK)
Shrinkwrap 3.25 (N steps = 6; triangulation OK)
Shrinkwrap 3.55 (N steps = 7; triangulation OK)
Shrinkwrap 3.83 (N steps = 8; triangulation OK)
Shrinkwrap 4.14 (N steps = 9; triangulation OK)
Shrinkwrap 4.43 (N steps = 10; triangulation OK)
Shrinkwrap 5.78 (N steps = 15; triangulation OK)
Objects 12.56 12.58 12.61 12.59

Table

1: Comparison of two polygonisation algorithms in terms of FEPTS
ffl the method to move vertices onto the surface for every next value of V 0;k
uses a Newton-Raphson al-
gorithm. This algorithm computes the distance a vertex should move along a straight line (namely the
gradient direction). If the distance between the iso-surface for k and k + 1 is too large (that is, if \DeltaV
is too large), this straight line can be a bad approximation of the curved field line that the vertex should
ideally follow. Therefore, the steps \DeltaV cannot be too large, and hence N steps cannot be too small. The
algorithm therefore could be made much more efficient if the algorithm to move the vertices were not
based on a first order Taylor expansion, but on a second order expansion. In that case, the curvature of
the ideal trajectory of the moving vertex could be approximated better, and a vertex could move over a
longer distance in one (curved) step. We are currently investigating this option.
ffl the criterion for acceptability of an edge is too strict. Indeed, irrespective of the local curvature, an edge
(and hence, a triangle) is split as soon as one of the accuracy conditions is violated. It is, however, only
essential that the surface is checked in samples that are not further away than the indicated distances;
it is not necessary that the triangles are split if samples in the interior of such a triangle would indicate
that the iso-surface is sufficiently close to the plane of the triangle. So a further reduction of the number
of triangles in flat regions could be obtained by a slightly more sophisticated splitting criterion.
ffl along the same line, we can even afford to leave out sampling the surface with a sample density in accordance
with the criteria 2 and 3 if a test could indicate that the entire region of an iso-surface within a
triangle is sufficiently close to the plane of that triangle. Such a test should employ interval arithmetic
or affine arithmetic on the bounding box of the triangle where a coordinate transformation should be
applied such that the plane of the triangle becomes parallel to one of the coordinate planes. Then a very
thin bounding box results, and the range of values of V (r) on this box could be established with very
few evaluations only.
ffl as indicated in Appendix B, the restriction of iso-surfaces that are homeomorphic to sphere can be easily
lifted.
5 Conclusion
We have presented a new algorithm for polygonizing an iso-surface, defined by a set of skeletal elements.
Whereas earlier approaches depend on a grid structure, the approach presented here is adaptive in the sense
that the tiling is at arbitrary locations on the surface. The algorithm produces triangles of different sizes and
shapes according to the local curvature of the surface offering an approximation with a known accuracy. The
advantages of this technique are summarized as follows:
ffl Adaptive mesh
due to motion in fixed grid
ffl Vertex Normals are calculated essentially for 'free'
ffl Faster (fewer FEPTs) than previous adaptive algorithms.
A disadvantage of this algorithm, is that it can only cope with a single closed surface without holes. However,
since the submission of the first version of this paper, we have devised an extension of this approach to solve
this problem; this is outlined in Appendix B.
6

Acknowledgements

We would like to thank Jules Bloomenthal and Geoff Wyvill for their continuing interest and enthusiasm for
this research. Also, the other members of the U. of Calgary graphics lab. (the GraphicsJungle) who have
helped with the ongoing research on iso-surfaces. This research was partly sponsored by the Natural research
and Engineering Council of Canada. One of the authors (CWAMvO) thanks the Department Board of Mathematics
and Computing Science of EUT for giving him the opportunity of his sabbatical leave.

Appendix

A: a smooth interpolation curve
The method for splitting edges as explained in 2.5 is based on finding a curve that should lay close to the iso-surface
and that passes through two given points, say p
and p 3
(the A and B from section 2.5). Since it should
be close to the iso-surface, it should also be perpendicular to the normals in these two points, say n
and n 3
We approximate this curve by a cubic Bezier curve, and we demand it to be optimally smooth. Call this curve
, and p 3
. Thus the extremes of the edge are p 0
and p 3
; the
normal vectors in these points are n 0
and n 3
. The unknown variables are p 1
and
. The boundary conditions
of perpendicularity to the normal vectors are expressed by
and
We solve this problem by minimising the average curvature of the curve subject to 11 and 12 by introducing
the Lagrange multipliers  0
and  3
. Then the minimisation problem is expressed as "minimise \Phi(p
where
(n 3
subject to 11 and 12.
Let us express the curve p(t) into the cubic Bernstein polynomials (they provide a complete base for cubic
polynomials, and cubic polynomialsare known to minimise the average curvature as expressed in the integral):
This yields for the integral (apart from a multiplicative factor)
Z 1j
We proceed by demanding @\Phi=@p 1
and @\Phi=@p 2
to vanish; this gives
To establish the values of the multipliers, (15) is substituted back, yielding the following set of linear equations
to be solved:
Observe that, in the case of normalised normal vectors, the coefficient matrix of this set takes the form
where OE is the angle between the two normal vectors. This matrix is always non-singular.
Finally, the desired midpoint M is found by taking
3:

Appendix

B: topological changes
The shrinkwrap algorithm as described in this paper assumes that the surface for V
to a sphere. It is not a priori possible to deduce from a given skeleton if this will be the case; moreover, the
applicability of shrinkwrap would be significantly improved if it could deal with arbitrary topological configu-
rations. Since the submission of the first version of this paper, an extension to shrinkwrap has been developed
which allows such generalisation. It will be submitted for publication elsewhere; however we give a coarse
outline of the proposed method here to make this paper self contained. The modifications to shrinkwrap are
as follows:
ffl a topological change occurs in a point (a saddle point) where rV in principle, for every
vertex r a search could be done to a nearby point r saddle with rV (r saddle using also 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 saddle point is near a given vertex;
ffl once a saddle point is found, its function value is computed:
ffl this value V saddle is inserted in the list of iso-values that is used for the shrink wrap iterations. Most
saddle will be the next higher iso-value adjacent to the current one;
ffl when the vertices are displaced to arrive on the surface with iso-value V saddle , make sure that a vertex
lands in the saddle point;
ffl in this saddle point, distinguish between a rupture and a hole. This distinction can be made on the basis
of a local approximation of the function near r saddle as a quadratic polynomial in x, y, and z. The
coefficients of this polynomial can be inspected in order to make the required distinction;
ffl disconnect the edges to the vertex r saddle and re-connect them in order to form the right topological
configuration;
ffl proceed with the iteration either to the next V 0;k
or the next V saddle for an optional next saddle point.



--R

Kees van Overveld
A Generalization of Algebraic Surface Drawing.
Polygonisation of Implicit Surfaces.
Interactive Techniques for Implicit Modeling.
How To Shrinkwrap a Saddle Point.
Marching Cubes: A High Resolution 3D Surface Construction Algorithm.
Discrete smooth interpolation in geometric modelling.
An evaluation of implicit surface tilers.

Boolean operations in solid modelling boundary evaluation and merging algorithms.


A coherant Projection Approach for Direct Volume
Data Structure for Soft Objects.
The Great Train Rubbery.
--TR

--CTR
Charlie C. L. Wang, Direct extraction of surface meshes from implicitly represented heterogeneous volumes, Computer-Aided Design, v.39 n.1, p.35-50, January, 2007
implicit surfaces;iso-surfaces;polygonization