Technical Report 2009-002

Polygonization of implicit surfaces using Delaunay
triangulation

Pierre-Alain Fayolle
August 19, 2009

Computer Graphics Laboratory
The University of Aizu
Tsuruga, Ikki-Machi, Aizu-Wakamatsu City
Fukushima, 965-8580 Japan

Technical Report 2009-002
Title:
Polygonization of implicit surfaces using Delaunay triangulation

Authors:
Pierre-Alain Fayolle
Key Words and Phrases:
Polygonization, implicit surfaces, mesh generation, mesh optimization

Abstract:
We present in this report two algorithms for polygonizing implicit surfaces based on Delaunay
triangulation. To improve the quality of the obtained triangle mesh, we propose a mesh
optimization algorithm for: retrieving sharp features, regularizing mesh triangles, andminimizing the approximation error. We show how to combine an adaptive selection of vertices,with a Delaunay triangulation and mesh optimization to produce good quality triangle meshes.

Report Date:
8/19/2009

Written Language:
EnglishAny Other Identifying Information of this Report:

Distribution Statement:
First Issue: 10 copies

Supplementary Notes:

Computer Graphics Laboratory
The University of Aizu
Aizu-Wakamatsu
Fukushima 965-8580
JapanPolygonization of implicit surfaces using Delaunay
triangulation
Pierre-Alain Fayolle
07/15/2009

Contents
1 Introduction
1.1 Related works . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Overview and main contributions . . . . . . . . . . . . . . . .

2
3
5

2 The first polygonization algorithm
2.1 Projection on the implicit surface . . . . . . . . . . . . . . . .
2.1.1 Projection for general implicit surface . . . . . . . . .
2.1.2 First order approximation . . . . . . . . . . . . . . . .
2.2 Adding ghost points . . . . . . . . . . . . . . . . . . . . . . .
2.3 Removing external tetrahedrons . . . . . . . . . . . . . . . . .
2.4 Mesh optimization . . . . . . . . . . . . . . . . . . . . . . . .
2.4.1 Optimization of the connectivity . . . . . . . . . . . .
2.4.2 Sharp features extraction . . . . . . . . . . . . . . . .
2.4.3 Mesh regularization . . . . . . . . . . . . . . . . . . .
2.4.4 Projection . . . . . . . . . . . . . . . . . . . . . . . . .
2.5 Thin triangles removal . . . . . . . . . . . . . . . . . . . . . .
2.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5
6
6
7
8
8
9
10
10
11
12
12
12

3 Adaptive polygonization of implicit surfaces
13
3.1 Previous works . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.2 Main algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.2.1 Adding tetrahedron circumcenters . . . . . . . . . . . 15
3.2.2 Adding triangle centroids . . . . . . . . . . . . . . . . 15
3.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
4 Results and discussions

16

5 Conclusion

20
Abstract1

We present in this report two algorithms for polygonizing implicit
surfaces based on Delaunay triangulation. To improve the quality ofthe obtained triangle mesh, we propose a mesh optimization algorithm
for: retrieving sharp features, regularizing mesh triangles, and minimizing the approximation error. We show how to combine an adaptive
selection of vertices, with a Delaunay triangulation and mesh optimization to produce good quality triangle meshes.

1

Introduction

In this report, we present two algorithms for the polygonization of implicit
surfaces based on the Delaunay triangulation algorithm. Fig. 1 is illustratingsome steps of the second algorithm.

Figure 1: Some steps of the adaptive algorithm presented in section 3.
Implicit surfaces [6], or F-Reps [3] represent surfaces as the set of points
in space taking a given function value (isovalue), for example the surface
defined by the isovalue 0 of f in R 3 is: S = {p ∈ R3 : f (p) = 0}. Taking
the inequality instead of the equality defines more generally a solid: V ={p ∈ R3 : f (p) >= 0}.
This representation is interesting because of its ability to describe complicated shapes and handle easily topological changes during simulation [14].

2Visualization is typically done by ray-tracing [11], i.e. shooting rays on
a two dimensional grid and checking for their intersection with the object,
or by polygonization [18, 4, 8]. Polygonization consists in generating anapproximation of the surface f = 0 by a set of triangles. It is usually doneby sampling the function on a regular grid and inspecting each cell to see
if it intersects with the surface, which is done by looking at the sign of thevalues at each corner of the cell.
In this paper, we propose a new approach for the polygonization of implicit surfaces: random points are projected on the surface f = 0 and atetrahedralization of the this point-set is generated by the Delaunay triangulation algorithm. The sign of the function is used to discard outside
tetrahedrons. Triangles on the surface are extracted and their shape is improved by using information from the function and the current triangular
mesh.

1.1

Related works

The Marching Cube algorithm [18] works by inspecting each cells of a regular grid and check for a potential intersection of the surface with the cell.
An intersection is detected when the values at each corner of the current
cell have different signs. A lookup table is used to retrieve the intersection
topology based on the configuration of the signs in the cell. It is a popular algorithm for extracting a triangulated surface from a regular grid ofvalues but suffers from several drawbacks such as: topological ambiguities,
exhaustive investigation of cells, smoothing of the surface’s sharp features,among others. For an in-depth review of the Marching Cube algorithm,
its improvements and extensions, the reader is referred to the survey from
Newman and Yi [17].
Pasko et al [4] proposed a similar algorithm, which also exhaustively
inspects all cells in the grid, but is free of topological ambiguity in contrast to
the original Marching Cube algorithm. The ambiguity in a cell is removed bynoticing that the intersection of the implicit surface with the cell corresponds
to the branches of an hyperbola.
An algorithm for retrieving the sharp features of the surface was proposed by Kobbelt et al. [12]. Their algorithm extends the Marching Cube
algorithm by using the local distance field information and its gradient tocompute and insert into the mesh additional sample points lying on the surface’s feature. Ohtake et al [19, 20] used a set of post-processing steps toimprove the triangular mesh generated by the Marching Cube algorithm.
Their work attempt to handle sharp features correctly as well as to regularize the mesh by relocating the mesh vertices based on the current meshinformation and the implicit surface (function value and gradient).
In general, these algorithms tend to produce an excessive number oftriangles, especially in the regions of low curvature, because of their use of3

a uniform grid. Moreover, the shapes of the triangles are not controlled andcan degenerate even when a post-processing optimization is used.
A different approach is the Dual Contouring algorithm introduced by Ju
et al [15], which is designed for extracting surfaces with adaptive resolution(in contrast to using a uniform grid) while reproducing sharp features (basedon normal information). However, the original Dual Contouring algorithm
is not guaranteeing intersection-free surfaces, and it was later extended by
Ju and Udeshi [16] to resolve this issue. In these methods, the size of thetriangle mesh is controlled by the octree’s depth. There is no mechanism inplace to control the shape of triangles.
A different family of algorithm for the polygonization of implicit surfaces
rely on surface tracking, i.e. these methods start from a seed triangle on
the surface and grow by iteratively adding new triangles to approximate thesurface. An example of surface tracking method is the Marching Triangles
algorithm [2].
Starting from a seed triangle on the surface, the algorithm iteratively
creates new triangles covering the surface by: adding a vertex in the planeof an existing triangle, projecting it to the implicit surface and deciding if
the newly formed triangle should be added. The later decision is based on
the Delaunay surface constraint proposed by Boissonnat [9].
The original Marching Triangles algorithm is not crack-free as noticed
by Akkouche and Galin. They proposed to extend it by: improving thetriangle creation step, especially the projection of the additional point on thesurface, and to tessellate the cracks. Additionally, an incremental algorithm
is proposed to integrate the improved marching triangles algorithm in aninteractive implicit surface editing environment.
However, none of these algorithms seem to be able to handle properly
sharp features of the surface.
Another type of algorithm was proposed by Desbrun et al [13]. Seedpoints on the bounding volume of an implicit surface migrate to the surface
to generate a mesh approximating the surface. The algorithm is designed
to work with implicit surfaces defined by skeletal elements only. The mainstrength of the algorithm is to be fast and to allow interactive update of
the mesh when adding or removing primitives. This is due to the fact thatthe algorithm operates on elements locally and not globally on the implicit
surface. The difficulty of the algorithm is to attach the triangles patches
together. The algorithm may fail to tessellate correctly regions where many
primitives blend. As mentioned earlier, it is also restricted to a single classof implicit surfaces: skeletal elements.

4

De Figueiredo et al. [7] proposed an algorithm inspired by physical systems to solve the problem of polygonization of implicit surfaces. First, aparticle system is constrained on the implicit surface by solving a partial
differential equation. Particles tend to be attracted by area of high curvatures. The sampling is then improved and made uniform by introducing
repelling forces between each particles and solving until equilibrium. An
improved algorithm for sampling particles on an implicit surface was presented by Witkin and Heckbert in [5]. Finally, a triangular approximationof the surface is obtained by keeping from the Delaunay triangulation of thepoint-set the triangles that approximate the tangent plane at each vertex.
If the initial set of particles is not properly scattered, this algorithm may
miss disjoint pieces and holes of the implicit surfaces. Scattering the points
by solving the partial differential equation is time consuming, as is the regularization by introducing repelling forces. Finally, sharp features are nothandled properly by the algorithm.

1.2

Overview and main contributions

In this paper, we present first a simple algorithm using Delaunay triangulation for the polygonization of an implicit surface f given a vertex budget
n. The algorithm consists in a few simple steps: 1) generate n random
points and project them on the boundary surface f = 0; 2) generate a set of
tetrahedrons using the Delaunay triangulation algorithm; 3) discard outside
tetrahedrons and generate triangles on the surface; 4) optimize the mesh(connectivity, surface approximation, triangles shape and sharp features).
An adaptive version of the algorithm is then presented. Starting from asmall sample of random points, it iteratively adds points in area of highcurvature while optimizing the generated triangle mesh.
The presented algorithms are relatively simple to implement, give a direct control on the number of vertices on the mesh, handle sharp featuresand generate relatively well shaped triangles. In contrast, to marching cube
algorithms, the function is sampled only on the vertices and not on every
cells’ corners (including in empty cells away from the surface). Additionally,
the adaptive method provides a better distribution of the triangles’ size byadding smaller triangles only in the high curvature areas.

2

The first polygonization algorithm

In this section, we describe the first version of the polygonization algorithm.
The outline of the algorithm is given first. The details for each steps are
given in the following subsections.
Polygonize (Polygonization of an implicit surface). Given a function f :R3 → R, a number of points n and a rectangular box bbox bounding the5space {(x, y, z) : f (x, y, z) ≥ 0}, compute a polygonal approximation of thesurface {(x, y, z) : f (x, y, z) = 0}.
1. [Generate points.] Generate n points p in the bounded space[bbox(1, 1), bbox(2, 1)] × [bbox(1, 2), bbox(2, 2)] × [bbox(1, 3), bbox(2, 3)].
Points can be randomly generated or lying on a regular grid.
2. [Projection.] Project the points p on the surface f = 0.
3. [Add ghost points.] Generate a small number m << n of ghost points
g lying on regular grids, on each faces of a slightly enlarged box containing all the mesh vertices. Add the ghost points g to the original
point-set p: p ← p + g.
4. [Delaunay triangulation.] Compute the Delaunay triangulation t of
the point-set p.
5. [Discard tetrahedrons.] Remove the tetrahedrons from t outside of thevolume {(x, y, z) : f (x, y, z) ≥ 0}.
6. [Generate triangles.] Find the set of triangles tri on the surface. A
triangle is lying on the surface if it is not shared by tetrahedrons.
7. [Mesh optimization.] Optimize the triangular mesh (tri, p) by repeating a set of simple steps that will optimize the connectivity, extract
sharp features, regularize the mesh and maintain vertices close to thesurface.
8. [Remove.] Remove triangles with an angle below some user defined
threshold.
The algorithm (steps 1 to 6) is illustrated in two dimensions by Fig. 2.
Steps 2, 3, 5, 7 and 8 are detailed in the following subsections.

2.1

Projection on the implicit surface

Given a point p and a function f , we are searching for the projection p sof the point p on the implicit surface f = 0. If f is the signed distancefunction, then the projection: p ← p − f (p)∇f (p) is exact. If f is not thesigned distance function, it is possible to project the point on the surface by
a numerical algorithm or to approximate its projection.
2.1.1

Projection for general implicit surface∆p = p − ps is parallel to ∇f (ps ) and at ps , f (ps ) = 0, gives the following
system of equations:f (p + ∆p) = 0∆p + t∇f (p + ∆p) = 0
6(1)

Figure 2: Illustration of the algorithm Polygonize, steps 1 to 6 in two
dimensions.
By eliminating t, we can get a system in ∆p only:f (ps )R(ps ) =  fy (xs − x) − fx (ys − y)  = 0
fz (xs − x) − fx (zs − z)(2)where fx , fy and fz are the partial derivatives of f with respect to x, y andz at ps .
Eq. (2) is solved for the unknown vector column p s with the damped
Newton method using as initial guess the current point p: let J be theJacobian matrix of the system R, pk+1 = pk − αJ −1 (pk )R(pk ) is iterated
until R(pk ) < .
2.1.2

First order approximation

The previous projection can be time consuming and we may want to avoid
it. It is possible to obtain an approximation of the projection. Starting witha first order Taylor expansion: f (ps ) = f (p + ∆p) ≈ f (p) + ∆p.∇f (p),and taking the gradient gives: ∇f (p + ∆p) ≈ ∇f (p). By combining with∆p + t∇f (p + ∆p) = 0, we obtain: ∆p + t∇f (p) = 0. f (p s ) = 0 gives:f (p) − tk∇f (p)k2 = 0
and:
∆p =

∇f (p)k∇f (p)k2

7

(3)

(4)

The first order projection of p to the implicit surface f = 0 is thus given
by:
p ←p−

∇f (p)k∇f (p)k2(5)

Practically, one iteration is not sufficient to reach the surface as illustrated
by Fig. 3. We can iterate the projection (5) until f (p) < .
2
1.5
1
0.5
0
−0.5
−1
−1.5
−2
−2

−1.5

−1

−0.5

0

0.5

1

1.5

2

Figure 3: Error committed when projecting a point on an implicit surface
with a first order projection. The original point is represented by the circle.
The projected point is represented by the cross. The implicit curve is:x2 + y 2 + xy − 0.5x2 y 2 − 0.25.

2.2

Adding ghost points

The step (3) of the Polygonize is adding a small number of points on aregular grid. The rectangular grid is obtained by slightly enlarging in eachdirections the rectangular box enclosing all mesh vertices (see Fig. 2). We
found that adding these extra points increase the quality of the result.

2.3

Removing external tetrahedrons

The Delaunay triangulation of a point-set gives a triangulation of the convex
hull of the point-set. Tetrahedrons outside of the domain need to be peeledoff (see Fig. 4 for an illustration in two dimensions). Distinguishing between
inside and outside tetrahedrons can be done by inspecting the sign of thefunction f at the centroid c of each tetrahedron t: if f (c) < 0, then thetetrahedron is outside, otherwise it is inside.
Practically, we may want to keep a slightly thicker solid. One possible
approach is to discard tetrahedron for which f (c) < −. Another approach
8

5

5

4.5

4.5

4

4

3.5

3.5

3

3

2.5

2.5

2

2

1.5

1.5

1

0

2

4

1

6

0

2

4

6

Figure 4: Delaunay triangulation of a point-set covers the convex hull. Outside tetrahedrons (triangles in two dimensions) need to be removed.

is to compute the circumcenter c and the circumradius r of the tetrahedront, and keep t if it is outside and if d(c)/r < α, where α is a user definedconstant (smaller than 1) and d is the distance from c to the surface. If f is
not the signed distance function, then d(c) can be approximated by kc−c s k,
where cs is the projection of c to f = 0 and is computed as explained in(2.1).

2.4

Mesh optimization

After the step (6) of the algorithm Polygonize, a triangular mesh approximating the surface f = 0 is obtained. However, the triangles are not necessarily well shaped and (in case the surface has sharp features) some sharp
features may not have been extracted correctly.
Optimization of the mesh is done using an iterative scheme which resembles what was proposed by Ohtake et al [20]. This iterative scheme consistsof the following steps:
Mesh optimization . Optimize the triangular mesh (tri, p) by repeating a set of simple steps that will optimize the connectivity, extract sharp
features, regularize the mesh and maintain vertices close to the surface:• [Optimize connectivity.] if needed: recompute the Delaunay triangulation of p, remove the external tetrahedrons and identify the triangles
on the surface: tri.
• [Sharp features extraction.] Move the triangle vertices to align thetriangle normal with ∇f , the gradient of f .
9• [Regularization.] Regularize the shape of the triangles.
• [Projection.] Project the points p on the surface f = 0 as in step 2 ofPolygonize.
The first step modifies the connectivity of the triangle mesh while theother steps are vector fields Fi that modifies the mesh vertices: pnew ←Fi (pold ).
2.4.1

Optimization of the connectivity

The mesh connectivity is optimized by computing the Delaunay triangulation of the point-set p1 . This is not needed for the first iteration of the meshsmoothing (since the mesh just got created), nor is it needed if the relative
movement of each vertices is small (below a user defined threshold). Similar
as in the triangle mesh creation, external tetrahedrons need to be removed(step 5) and triangles on the surface need to be extracted (step 6). This canbe summarized by the following algorithm.
Optimize connectivity . Given a point-set p, a triangular mesh tri, andthe current iteration index iter. Let p 0 = p.
1. [Exit ?] if kp − piter−1 k < ttol exit.
2. [Delaunay triangulation.] Compute the Delaunay triangulation t of
the point-set p.
3. [Discard tetrahedrons.] Remove the tetrahedrons from t outside of thevolume {(x, y, z) : f (x, y, z) >= 0}.
4. [Generate triangles.] Find the set of triangles tri on the surface. A
triangle is lying on the surface if it is not shared by tetrahedrons.
2.4.2

Sharp features extraction

For extracting the sharp features, we follow the method of Ohtake et al(section 2.2 in [20]): mesh vertices are moved in order to align the triangle
normals with the gradient of f .
Let N (p, f ) be the transformation applied to a mesh vertex p, then themodified vertex pnew is set to: pnew ← pold + N (pold , f ).
∇f (x,y,z)
Let n(T ) be the normal to the triangle T , m(x, y, z) = |∇f(x,y,z)| be theunit gradient to f , and A(T ) the area of the triangle T . The transformation
N is defined by:
N (p, f ) = P
1

1X

1ring(p) A(T ) 1ring(p)A(T )v(T )the point-set contains the mesh vertices and the ghost points.

10

(6)

−
−
where v(T ) = (→pc.m(c)).m(c) is the projection of →pc on the m(c) direction,c is the centroid of T and 1ring(p) is a function returning the set of triangles
T with p as a vertex.
As noted by Ohtake et al in [20], this step is actually doing more than
only restoring sharp features: it also improves the position of the meshvertices, as illustrated by the figure 5 in [20].
2.4.3

Mesh regularization

The modified position pnew of a mesh vertex pold is obtained by:1

pnew ← P

X

1ring(pold ) A(T ) 1ring(pA(T )c(T )

(7)

old )

where 1ring(p) is a function returning the set of triangles T with p as avertex, A(T ) is the area of the triangle T and c(T ) is the centroid of thetriangle T .
Intuitively, if we consider a set of planar triangles around a vertex p,then applying Eq. (7) to p should move it to balance the area of eachtriangles. Especially, if the triangle vertices (other than p) lay on a circle,
then p should be moved to the center of the circle as illustrated in Fig. 5.

2.5

2.5

2

2

1.5

1.5

1

1

0.5

0.5

0

0

−0.5

−0.5

−1

−1

−1.5

−1.5

−2

−2

−2.5

−2.5
−1

0

1

−1

0

1

Figure 5: Illustration of the mesh regularization procedure Eq. 7. p is moved
to the center of the circle in order to balance the area of each triangle.
Updating mesh vertices according to Eq. (7) destroys sharp edges. It ispreferable to identify sharp features and mark them as fixed, before invoking the mesh regularization procedure. Identifying sharp features is done
following a method proposed by Kobbelt et al [12]. A mesh vertex p istagged as belonging to a sharp feature if: min i,j (ni .nj ) < α, where ni is the11gradient of f at the centroid ci of the i-th triangle in the one ring of p, andα a user controlled parameter.
2.4.4

Projection

Eq. (7) tends to slightly move the mesh vertices away from the surface. The
mesh vertices are projected back to the surface using the methods describedin section 2.1. Even if the function f is not a signed distance function,
convergence of the methods 2.1.1 or 2.1.2 is practically fast, since the vertices
are initially close to the surface.

2.5

Thin triangles removal

At the end of the step (7) of the algorithm Polygonize, it is possible thatthe triangular mesh contains some thin triangles with very small angles.
This will most likely happen when the surface contains sharp features, because mesh vertices on sharp features are marked and fixed during the meshregularization.
To resolve this issue, an additional step is added to iterate through thetriangles and remove those with angle(s) below some user defined threshold.
Removal is done by collapsing the edge opposite to the small angle andmerging its incident vertices (see Fig. 6 top). If the triangle is flat with twoangles below the threshold, we collapse the smaller of the two edges.

Figure 6: Removing thin triangles is done by collapsing edges opposite tosmall angles within a triangle. The top row illustrates the case where thetriangle has one small angle and the bottom row illustrates the case where
the triangle has two small angles.

2.6

Discussion

If the point-set is well sampled, then the polygonal approximation will begood (in term of topology and distance). The work of Amenta and Bern[1] characterizes good surface samplings. If the result of Polygonize is notsatisfactory, it is always possible as a first and rough solution to increase
the number of points n.
In term of complexity, steps (1) to (6) (included) have a worst casecomplexity of O(n2 ), because of the Delaunay triangulation. One iterationof the mesh optimization algorithm, step (7), has a worst case complexity12O(n2 ) (sharp features extraction and regularization are both O(n)). If the
optimization is iterated N times then the total worst case complexity isO(N n2 ).
Algorithm Polygonize can be improved: the initial sampling of the npoints is random, and sharp features (if any) tend to attract points during the projection (see Figure). Even if the mesh optimization (especiallystep 2 and 3) tries to improve the vertices locations, there may still be too
many triangles in area of low curvature and not enough in area of high curvature. It is possible to resolve that issue by: collapsing edges in area oflow curvature and using an adaptive subdivision in area of high curvature.
Another possibility explored in the next section and inspired by the work of
Boissonnat and Oudot [10] is to use an adaptive algorithm.

3

Adaptive polygonization of implicit surfaces

In this section an adaptive version of the algorithm Polygonize is described.
This algorithm extends the previous one by adding points in area insufficiently sampled, while maintaining an optimized mesh.

3.1

Previous works

Boissonnat and Oudot gave in [10] the theoretical background and an algorithm for the creation of a good surface sampling. Their algorithm proceedby iteratively adding the centers of the bad surface Delaunay balls to thecurrent point-set, until all Delaunay ball become good. A surface Delaunayball is defined in their work as any ball circumscribing a Delaunay facet fand centered at the intersection of the surface with the Voronoi edge dual
to f . Whether a surface Delaunay ball is good or not is decided based on
the ratio of the circumradius with a 1-lipschitz function. In their analysis,
the distance to the medial axis is used, which is rather difficult to computepractically.
In the following algorithm, the idea of adding points to the point-set is
kept but the criteria for adding points is simplified. In addition, at each iteration, the current mesh is optimized (following the step (7) of Polygonize
for a small number of iterations).

3.2

Main algorithm

Adaptive polygonization . Given a function f : R 3 → R, a number of points n and a rectangular box bbox bounding the space {(x, y, z) :f (x, y, z) ≥ 0}, compute a polygonal approximation of the surface {(x, y, z) :f (x, y, z) = 0}.

13

1. [Generate points.] Generate n points p in the bounded space[bbox(1, 1), bbox(2, 1)] × [bbox(1, 2), bbox(2, 2)] × [bbox(1, 3), bbox(2, 3)].
Points can be randomly generated or lying on a regular grid.
2. [Projection.] Project the points p on the surface f = 0.
3. [Add ghost points.] Generate a small number m << n of ghost points
g lying on regular grids, on each faces of a slightly enlarged box containing all the mesh vertices. Add the ghost points g to the original
point-set p: p ← p + g.
4. [Loop.] Repeat the following steps:(a) [Delaunay triangulation.] Compute the Delaunay triangulation tof the point-set p.
(b) [Discard tetrahedrons.] Remove the tetrahedrons from t outsideof the volume {(x, y, z) : f (x, y, z) ≥ 0}.
(c) [Add tetrahedron circumcenters.] Select tetrahedron circumcenters c such that d(c)/r < β, where d is the unsigned distance tothe surface, r is the circumradius and β a user defined constant.
Project c on the surface f = 0, and add to the list of additional
points.
(d) [Generate triangles.] Find the set of triangles tri on the surface.
A triangle is lying on the surface if it is not shared by tetrahedrons.
(e) [Mesh optimization.] Optimize the triangular mesh (tri, p) by
repeating (for a small number of iterations) a set of simple steps
that will: extract sharp features, regularize the mesh and maintain vertices close to the surface.
(f) [Add triangle centroids.] Select triangle centroids c such thatd(c) < h0 and r < Rmin where r is the triangle circumradius,h0 and Rmin are user defined constant. Project c on the surface
f = 0 and add to the list of additional points.
(g) [Exit ?] if the number of mesh vertices is greater than n or there
is no more points to add then exit the loop
5. [Mesh optimization.] Optimize the triangular mesh (tri, p) by repeating a set of simple steps that will optimize the connectivity, extract
sharp features, regularize the mesh and maintain vertices close to thesurface.
6. [Remove.] Remove triangles with an angle below some user defined
threshold.

14

2
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0
−0.5

0

0.5

1

1.5

2

Figure 7: Because of an insufficient sampling during the first iteration of thealgorithm, the surface may cross the inside of a tetrahedron.

Most of the components of this algorithm have already been discussed in
section 2. Step (4.e) and (5) differ as no connectivity optimization is done in(4.e). There is no need to call the connectivity optimization in (4.e) because
a Delaunay triangulation is already performed at each iteration of the loop(4) (in step (4.a)).
3.2.1

Adding tetrahedron circumcenters

Let c be the circumcenter of a tetrahedron t. c is selected if the ratio
between the distance from c to the surface and the circumradius r is lessthan a threshold. The goal of this rule is to try to capture undersampledareas, which appear at the first iterations of the algorithm, when the numberof points is small. Fig. 7 shows the type of configuration that we wantto capture. Selected circumcenters are then projected on the surface andmarked as additional points.
3.2.2

Adding triangle centroids

Let c be the circumcenter of a triangle tri. ct is selected if its distanceto the surface is greater than a given threshold and the circumradius r is
greater than a given threshold. Selecting triangle circumcenters in the area
where the approximation error is high, allow to properly adapt the surface
approximation. Limiting the size of the circumradius, allow to avoid too
small triangles in area of high curvature. Selected triangle circumcenters
are then projected on the surface and marked as additional points.

15

3.3

Discussion

The algorithm Adaptive polygonization improves the algorithm Polygonize by adding new vertices where the geometric approximation is insufficient instead of starting with a randomly distributed set of points. Thismeans that Adaptive polygonization should allow fewer triangles than
Polygonize in area of low curvature and more triangles in area of high curvature. Given a number of vertices n, it does a better usage of the vertices.
In term of complexity, the loop (4) for has a worst case complexity ofO(n3 ) (in the case where one point is added at each iteration of the loop (4)).
Step (5) has a worst case complexity of O(N n 2 ), where N is the number of
iterations in (5) (see section 2.6).

4

Results and discussions

It is possible to polygonize rather complicated shapes as illustrated by Fig. 8.
This shape contains multiple disjoint components (the sphere in the middle,
the eight spherical parts around the core sphere and the outside torii), andsharp features (see the right image in Fig. 8). The polygonal model in Fig.
8 was obtained by using the adaptive algorithm (section 3.2). 10 steps of theloop (4) were needed to reach the final sampling. 50 iterations of the mesh
optimization were used to optimize the mesh quality. The mesh is made of
29,998 vertices and 59,980 faces. The average radius ratio 2 is 1.13 with astandard deviation of 0.17. The maximum ratio is 4.2.

Figure 8: Polygonization of a complex model (’core’) with several disjoint
components, and sharp edges. See how sharp features are properly retrieved(right). The polygonal model is made of 29,998 vertices.
2

the radius ratio is defined as the ratio of the radius of the circumscribing circle to theradius of the inscribed circle.

16

Both the algorithm presented in section 2, and the adaptive algorithm
presented in section 3.2, produces good results, with good mesh quality(distributions of radius ratio and minimum angle are presented) as illustrated Fig. 9 (for the algorithm section 2) and Fig. 10 (for the adaptivealgorithm section 3.2). The resulting polygonal models shown in Fig. 9 and
10 are made of approximately 7000 vertices.

Figure 9: Result of the polygonization of a torus with the algorithm presented in section 2. Distributions of radius ratio and minimum angle illustrate the mesh quality.

Figure 10: Result of the polygonization of a torus with the adaptive algorithm presented in section 3. Distributions of radius ratio and minimum
angle illustrate the mesh quality.
In general, the adaptive algorithm does a better treatment of positioning
the vertices as illustrated in Fig. 11. The model on the left has been createdwith the algorithm from section 2, while the model on the right has beencreated with the adaptive algorithm. The top row illustrates a view of the
top of the object, which is flat. See how the adaptive algorithm adapts
the vertex distribution to the zone of high curvature. On the other hand,17algorithm Polygonize use lots of vertices on the flat zone, while area ofhigher curvature are more sparse (see Fig. 11 bottom left).

Figure 11: Illustration of the vertex distribution with the algorithm Adaptive polygonization (top and bottom rows, right column) against the vertex distribution with the algorithm Polygonize (top and bottom rows, left
column) for the model ’sand’. See how vertices are concentrated on the areaof high curvature with the adaptive algorithm.
Table 1 shows some qualitative results of the polygonization of several
models using both the adaptive algorithm from section 3.2 and the algorithm
Polygonize from section 2. For both algorithms, 50 iterations of the mesh
optimization loop were performed.
Effect of fixing vertices on sharp features during the mesh regularization Before optimizing the position of each vertices according to eq.
7, we advised to identify vertices on sharp features and mark them, suchthat they remain fixed. Fig. 12 shows a cube polygonized with (left) andwithout (right) fixing the vertices on sharp features.
Fixing vertices on sharp features has a negative effect on the triangles
shape, as it forbids the triangles around such vertices to be regularized. A
practical approach is to allow the full mesh optimization, without fixing
vertices, for some iterations, then allow to fix vertices on sharp featuresduring the remaining iterations. See the results in fig. 13 where the left
picture corresponds to 50 mesh optimization iterations with fixed vertices on18Model

#v

#tris
1996

max / averageradius ratio
1.65 / 1.07

min / averageminimum angle
30 / 47.6

spheresphereadaptive
torus
torus
adaptivecylindercylinderadaptivecore
core
adaptivesand
sand
adaptive

1000
500
6998

996
13996

1.7 / 1.08
2.2 / 1.09

30 / 47.4
25 / 46.4

7219
3270

14438
6538

2.5 / 1.07
6.04 / 1.18

21 / 48.7
15.11 / 44.5

4467
29396

8934
58780

5.02 / 1.15
6.2 / 1.13

15.06 / 45.3
15.08 / 44.7

21769
27755

43526
55540

5.5 / 1.15
7.45 / 1.13

15.1 / 44.7
15 / 45.5

28103

56346

7.1 / 1.14

15 / 44.8

Table 1: Max/average radius ration and min/average minimum angle obtained on a series of models.

sharp features, while the right picture corresponds to 15 iterations without
fixing vertices and 35 iterations with fixing the vertices, as they are reaching
sharp features.
Limitations The limitations of the first algorithm (section 2) have been
already mentioned in section 2.6 and addressed with the second algorithm(section 3).
Both algorithms rely on an implementation of the Delaunay triangulation
algorithm. In our experiments, we used Matlab R , and its implementation of
the Delaunay triangulation. The Delaunay triangulation is the bottleneckof the algorithm in terms of speed. A fast and robust implementation ofDelaunay triangulation is required.
Both algorithms are not very fast. It is possible to speed up the adaptative algorithm by omitting the Delaunay triangulation during the mesh
optimization (step 5). Practically, we did not find that removing the Delaunay triangulation from the mesh optimization, in step 5, had any noticeable
effect on the final result. Figure 14 shows the minimum angle and radiusratio distributions for the model core with (left) and without (right) the
Delaunay triangulation in step 5. On the other hand, for the algorithm
polygonize (section 2), removal of the Delaunay triangulation in the mesh
optimization results in a noticeable lower quality for the final mesh.

19

Figure 12: The impact of fixing vertices on sharp features is illustrated on
the left; by allowing vertices to move following Eq. 7 sharp features are notretrieved properly (right) despite the sharp feature extraction step of section
2.4.2
The second bottleneck of the algorithms is the number of function (implicit surface) evaluations. If the gradient of the function is unknown, it
has to be computed numerically by finite differences. In total, there are
4 function calls at each projection (there is one projection at each triangle
mesh creation, and one projection at each mesh optimization). The gradientof the function is also needed in the mesh optimization for retrieving sharp
features.
Both algorithms tend to produce triangles with one (or two) very small
angle(s) close to sharp edges and corners. This is caused by the fact thatvertices on sharp features are fixed during the mesh optimization. As discussed above, it is possible to first optimize the mesh without fixing vertices
on sharp features, and then apply again the mesh optimization but with
fixed vertices on sharp features. We also proposed to remove skinny triangles with small angles (see section 2.5). For the threshold angle, we used 15
degrees for model with sharp features (like for example: core, sand), and 20
degrees for others (like for example: sphere, torus).

5

Conclusion

We have proposed two algorithms for the polygonization of implicit surfaces
based on Delaunay triangulation. We have illustrated the performance of
the algorithms on several examples, including complex implicit surfaces withsharp features.

20

Figure 13: Early fixing of vertices on sharp features forbids good triangle
shapes close to the sharp features. Left model was obtained with 50 iterations of the mesh optimization with vertices fixed on sharp features. Rightmodel was obtained with 15 iterations of mesh optimization without locking
vertices, and 35 iterations with locked vertices.

References[1] Amenta A. and Bern M. Surface reconstruction by voronoi filtering.
In ACM Symposium on Computational Geometry proceedings, pages
39–48, 1998.
[2] Hilton A., Stoddart A., Illingworth J., and Windeatt T. Marching triangles : range image fusion for complex object modelling. In International
Conference on Image Processing proceedings, 1996.
[3] Pasko A., Adzhiev V., Sourin A., and Savchenko V. Function representation in geometric modeling: concepts, implementation and applications. The Visual Computer, 8(11):429–446, 1995.
[4] Pasko A., Pilyugin V., and Pokrovskiy V. Geometric modeling in theanalysis of trivariate functions. Computers and Graphics, 12:457–465,
1988.
[5] Witkin A. and Heckbert P. Using particles to sample and control implicit surfaces. In ACM Siggraph 1994 proceedings, pages 260–268, 1994.
[6] Jules Bloomenthal, editor. Introduction to Implicit Surfaces. MorganKaufmann, 1997.
[7] De Figueiredo L. H., De Miranda Gomes J., Terzopoulos D., and Velho
L. Physically-based methods for polygonization of implicit surfaces. In
Graphics Interface 92 proceedings, pages 250–257.
[8] Bloomenthal J. Polygonization of implicit surfaces. Computer Aided
Geometric Design, 5(4):341–355, 1988.
21

[9] Boissonnat J.-D. Geometric structures for three-dimensional shape representation. ACM Transactions on Graphics, 3(4), 1984.
[10] Boissonnat J.-D. and Oudot S. An effective condition for sampling
surfaces with guarantee. Technical Report RR-5064, INRIA, December
2003.
[11] Hart J.C. Ray tracing implicit surfaces. Technical Report EECS-93014, WSU, 1993.
[12] Kobbelt L., Botsch M., Schwanecke U., and Seidel H.-P. Feature sensitive surface extraction from volume data. In SIGGRAPH’01 Proceedings, pages 57–66.
[13] Desbrun M., Tsingos N., and Gascuel M.-P. Adaptive sampling of
implicit surfaces for interactive modeling and animation. Computer
Graphics Forum, 15(5):319–325, 1996.
[14] Osher S. and Fedkiw R. Level Set Methods and Dynamic Implicit Surfaces. Springer, 2002.
[15] Ju T., Losasso F., Schaefer S., and Warren J. Dual contouring on
hermite data. In ACM Siggraph 2002 proceedings.
[16] Ju T. and Udeshi T. Intersection-free contouring on an octree grid. In
Pacific Graphics 2006 proceedings.
[17] Newman T. and Yi H. A survey of the marching cubes algorithm.
Computers and Graphics, 30:854–879, 2006.
[18] Lorensen W. and Cline H. Marching cubes: A high resolution 3d surface
construction algorithm. Computer Graphics, 21(4), 1987.
[19] Ohtake Y. and Belyaev A. Dual/primal mesh optimization for polygonized implicit surfaces. In ACM Solid Modeling proceedings, pages
171–178.
[20] Ohtake Y., Belyaev A., and Pasko A. Dynamic mesh optimization forpolygonized implicit surfaces with sharp features. The Visual Computer, 19(2-3):115–126, 2003.

22

Figure 14: Minimum angle and ratio distributions for the model core obtained with the adaptative algorithm with (left) and without (right) Delaunay triangulation in the mesh optimization step 5. Similar distributions
indicate that the effect of the Delaunay triangulation in the mesh optimization may not be so important.

23