Restricted Delaunay Triangulations and Normal Cycle

David Cohen-Steiner

Jean-Marie Morvan

ABSTRACT
We address the problem of curvature estimation from sampled smooth surfaces. Building upon the theory of normal cycles, we derive a denition of the curvature tensor for polyhedral surfaces. This denition consists in a very simple and new formula. When applied to a polyhedral approximation of a smooth surface, it yields an efcient and reliable curvature estimation algorithm. Moreover, we bound the difference between the estimated curvature and the one of the smooth surface in the case of restricted Delaunay triangulations.

Categories and Subject Descriptors

F.2.2 [Analysis of algorithms and problem complexity]: [Geometrical problems and computations, Computations on discrete structures]; G.2.3 [Discrete Mathematics]: Applications

General Terms
Algorithms, Theory

Keywords
curvature, geometric measure theory, mesh

Introduction
In many applications such as surface segmentation, anisotropic remeshing [21] or non-photorealistic rendering, a key step is Work partially supported by the IST Programme of the EU as a Shared-cost RTD (FET Open) Project under Contract No IST-2000-26473 (ECG Effective Computational Geometry for Curves and Surfaces). I.N.R.I.A., Sophia-Antipolis, France. E-mail: David.Cohen -Steiner@sophia.inria.fr I.N.R.I.A., Sophia-Antipolis, France. E-mail: JeanMarie.Morvan@sophia.inria.fr

to estimate the curvature of a smooth surface knowing only a polyhedral approximation of it. A lot of efforts have been devoted to this problem, leading to several curvature estimators, see [17] for a survey. Popular methods include quadric tting, where the estimated curvature is the one of the quadric that best ts the sample points locally. Other methods typically consider some denition of curvature that can be extended to the polyhedral setting. An example is [18], where the curvature is estimated using a discrete analog of the LaplaceBeltrami operator. The main shortcomings of these approaches are the lack of analysis of the quality of the obtained estimators, and also the lack of sound theoretical foundations. In this paper, we use the theory of normal cycles from differential geometry to dene curvature tensors for a general class of surfaces, including smooth and polyhedral ones. More precisely, we associate with each region a tensor which in the smooth case is the average of the curvature tensor over this region. The curvature tensor of a polyhedral approximation of a smooth surface then provides an estimator of the one of the smooth surface. Our main result is a bound on the difference between the estimated curvature and the actual one when the polyhedral approximation is chosen to be a Delaunay triangulation restricted to the surface. This case is important in practice, in particular when the triangulation is obtained by a Delaunay-based surface reconstruction algorithm. Under a local uniformity condition on the sampling, this bound implies that our estimator converges linearly with respect to the sampling density. To the best of our knowledge, only weaker results have been obtained in the past [19]. Our result can be viewed as a quantitative version of a theorem obtained by J.Fu [11] for gaussian and mean curvatures. This paper is organized in four sections. We rst introduce some notations and state the theorem (section 1). Then we present the theory of normal cycles (section 2) and how they can be used to deal with curvature tensors (section 3), which is a contribution of this paper. The proof of the theorem, based on geometric measure theory, is sketched in section 4.

1. STATEMENT OF THE THEOREM

In the sequel, we denote by M a surface in the three dimenPermission to make digital or hard copies of all or part of this work for 3 personal or classroom use is granted without fee provided that copies are sional oriented euclidean space . We assume for simplicity not made or distributed for prot or commercial advantage and that copies that M is the boundary of some compact set V 3 . bear this notice and the full citation on the rst page. To copy otherwise, to republish, to post on servers or to redistribute to lists, requires prior specic 1.1 Curvature measures permission and/or a fee. Let us rst recall some basic denitions and notations in SoCG03, June 8-10, 2003, San Diego, California, USA. the case where M is smooth. A good reference for these is Copyright 2003 ACM 1-58113-663-3/03/0006 ...$5.00.

312

[7]. The unit normal vector at a point p M pointing outward V will be refered to as n(p). Note that M is thereby oriented. Given a vector v in the tangent plane Tp M to M at p, the derivative of n(p) in the direction v is orthogonal to n(p) as n(q) has unit length for any q M . The derivative Dp n of n at p thus denes an endomorphism of Tp M , known as the Weingarten endomorphism. The Weingarten endomorphism can be shown to be symmetric ; the associated quadratic form is called the second fundamental form. Eigenvectors and eigenvalues of the Weingarten endomorphism are respectively called principal directions and principal curvatures. Both principal curvatures can be recovered from the trace and determinant of Dp n, also called mean and gaussian curvature at p. Our result does not involve curvatures at a single point, but rather curvature measures, which we dene here : D EFINITION 1. The gaussian curvature measure of M , G , V is the function that associates with every (Borel) set B 3 the quantity : G (B) = V G(p)dp
BM

HV (X, Y ) = 0 whenever X or Y is orthogonal to Tp M . The latter has already been considered by Taubin in [20]. We now introduce two matrix valued measures which are in some sense anisotropic versions of curvature measures : D EFINITION 3. The anisotropic curvature measures HV associate with every (Borel) set the 3 3 symmetric V and H bilinear form : HV (B) = HV (B) =
BM

HV (p)dp HV (p)dp

BM

Again, corresponding objects can be dened in the polyhedral case. If V is a polyhedron, we dene the discrete anisotropic curvature measures by : D EFINITION 4. (B) = V H (e)length(e B) e e
eE

where G(p) is the gaussian curvature of M at point p. Similarly, we dene the mean curvature measure H by : V H (B) = V H(p)dp
BM

HV (B) =

length(e B) [((e) sin(e))e+ e+ 2 eE +((e) + sin(e))e e ]

H(p) being the mean curvature of M at point p. Corresponding objects can be dened for triangulated surfaces. Assume now that V is a polyhedron with vertex set P and edge set E. D EFINITION 2. The discrete gaussian curvature measure of M , G , is the function that associates with every (Borel) V set B 3 the quantity : G (B) = V
pBP

where denotes a unit 3-vector with the same direction as e edge e, and e+ and e respectively denote the normalized sum and difference of unit normal vectors to triangles incident on e. If u and v are two vectors, u v is the bilinear form with matrix u.v t .

1.3 Theorem
We now go back to the case where M is smooth. In the sequel, P denotes a nite subset of M , and T is the Delaunay triangulation of P restricted to V , that is the union of Delaunay simplices the dual of which meet V . D EFINITION 5. P is said to be a -sample [2] of M if for all point p M , the ball B(p, lfs(p)) centered on p and with radius lfs(p) meets P. lfs, the local feature size [2], denotes the distance function to the medial axis of 3 \ M . T HEOREM 1. Let P an -sample of M with < 0.08. If B is the relative interior of a union of triangles of T , then : |G (B) G ((B))| K T V |H (B) H ((B))| K T V (B) H ((B))|| K T V ||H ||HT (B) HV ((B))|| K where for xed M K = O(
{tT, tB}

g(p)

(1)

where g(p) is the angle defect of M at point p, that is 2 minus the sum of angles between consecutive edges incident on p. Similarly, we dene the discrete mean curvature measure H V by : H (B) = V
eE

length(e B)(e)

(2)

|(e)| being the angle between the normals to the triangles of M incident on e. The sign of (e) is chosen to be positive if e is convex and negative if it is concave. In section 2 we will see where these formulas come from and why we use the same notation for continuous and discrete curvature measures.

1.2 Anisotropic curvature measures

In the case where M is smooth, the second fundamental form of M associates with each point p M a 2 2 symmetric bilinear form on Tp M , denoted by HV (p). The 2 2 symmetric bilinear form on Tp M having the same eigenvectors as HV (p) but with swapped eigenvalues will be denoted HV (p). These bilinear forms can be extended to 33 symmet ric bilinear forms HV (p) and HV (p) by setting HV (X, Y ) =

r(t)2 ) r(t))
{tT, tB, tB=}

+O(

r(t) being the circumradius of triangle t. In particular, when the sampling is locally uniform [3], K = O(area(B)+length(B))1
1 This result actually holds for any bounded aspect ratio triangulation projecting homeomorphically on M .

313

Here denotes the projection on M . Thanks to the assumption < 0.08, is indeed dened on T and is a homeomorphism from T to M as Amenta et al. showed [1].

B of M , then we have : V ol(V (B)) = area(B) + H (B) V 2 3 + G (B) V 2 3

2. NORMAL CYCLES AND CURVATURE MEASURES

Introduced by Wintgen and Z hle [15, 14], the theory of a normal cycles provides a unied way to dene curvature for both smooth and polyhedral surfaces2 . Here is a very crude overview of their approach. The rst observation is that the curvature measures of a smooth surface actually are byproducts of an object associated to the surface, called the normal cycle of the surface. More precisely, the curvature measures of a surface can be easily recovered from the normal cycle of the surface. Second, the denition of the normal cycle of a surface has a unique natural extension to the polyhedral case. Finally, the curvature measures of a polyhedron are dened to be the measures recovered from its normal cycle. Before explaining what a normal cycle is, we shortly review an early approach to curvature measures and the required background.

In the smooth case, the volume of V (B) is thus a polynomial in , and its coefcients are multiples of the curvature measures of B. H. Federer [9] actually showed that the volume of V (B) is always a polynomial in for < , even if the boundary of V is not smooth. The coefcients of this polynomial thus provide a way to generalize the denition of curvature measures as soon as is strictly positive. For instance, if V is a convex polyhedron, the obtained denitions agree with denition 2 3 . Unfortunately, this approach does not generalize to the case where V is a non convex polyhedron, for instance, as then equals 0 ; this is the reason why the theory of normal cycles was developed.

2.2 Background
The reader acquainted with differential calculus might want to skip this section. [8] provides a good introduction to the subject.

2.1 A rst approach

Historically, curvature measures were introduced by considering offsets of V , via the so-called tube formula. B

2.2.1 2-differential forms

Let S be a smooth manifold of dimension at least two embedded in some euclidean space k . If f is a vector eld on S, we denote by fx Tx S the vector associated with a point x S. 2-differential forms are, in a certain sense, 2dimensional analogs of vector elds : D EFINITION 6. A 2-differential form on S associates with every point x S a skew-symmetric bilinear form on Tx S, denoted by x . The following denition shows how a 2-differential form can be built from two vector elds : D EFINITION 7. The exterior product f g of two vector elds f and g on S, is the 2-differential form dened by : (f g)x (u, v) = (fx gx )(u, v) = for all x in S and (u, v) Tx S. Exterior products are special cases of 2-differential forms. However, they provide a good intuition of the general case : any 2-differential form can actually be written as a linear combination of exterior products of vector elds. It can be seen from the denition of an exterior product that if A is a linear transformation of the plane P spanned by u and v, then (f g)x (Au, Av) = det(A)(f g)x (u, v). In particular, (f g)x (u, v) = (f g)x (u , v ) for any two direct orthonormal frames (u, v) and (u , v ) of P . Note that this property extend to general 2-differential forms by linearity. Similarly, we have fx gx (u, v) = fx gx (u, v) for any couple of orthonormal frames (fx , gx ) and (fx , gx ) spanning the same oriented plane. Important examples of exterior products are
3 The case of a convex polyhedron is actually the rst considered historically, by Jacob Steiner [13].

2.2.1.1 Denition.

1111111111111111111111111111 0000000000000000000000000000 V (B)

1111111111111111111111111111 0000000000000000000000000000 1111111111111111111111111111 0000000000000000000000000000 1111111111111111111111111111 0000000000000000000000000000 V 1111111111111111111111111111 0000000000000000000000000000 1111111111111111111111111111 0000000000000000000000000000 1111111111111111111111111111 0000000000000000000000000000 M 1111111111111111111111111111 0000000000000000000000000000 1111111111111111111111111111 0000000000000000000000000000 1111111111111111111111111111 0000000000000000000000000000 1111111111111111111111111111 0000000000000000000000000000 1111111111111111111111111111 0000000000000000000000000000 1111111111111111111111111111 0000000000000000000000000000 1111111111111111111111111111 0000000000000000000000000000 1111111111111111111111111111 0000000000000000000000000000 1111111111111111111111111111 0000000000000000000000000000 1111111111111111111111111111 0000000000000000000000000000 1111111111111111111111111111 0000000000000000000000000000 1111111111111111111111111111 0000000000000000000000000000
Sk( \ V )
2

fx .u fx .v

gx .u gx .v

Figure 1: The tube formula

Let be the distance between M and the medial axis of the complement of V and / V = {p|p V d(p, V ) < } 3 that is the -offset of V minus V . The tube formula then reads : V ol(V ) = area(M ) + H (M ) V 2 3 + G (M ) V 2 3

for < . Moreover, this formula can be localized : if one only considers the part V (B) of V that projects on a subset
2 Normal cycles theory actually applies to a much more general class of objects, of any dimension and any codimension.

314

area forms. Area forms are a way to represent oriented surfaces as 2-differential forms. If T S is an oriented surface, then the area form of T is constructed as follows : for each point x T , pick a direct orthonormal frame of the tangent / plane Tx T , say (ux , vx ). For x T , set ux = vx = 0. The area form of T , denoted by aT , is the 2-differential form u v. Intuitively, area forms can be thought of as elds of surface elements : when applied to two vectors a and b in Tx S, aT x yields the signed area of the parallelogram spanned by the projections of a and b on Tx T .

The surface U that is setwise the same as T but with reverse orientation thus corresponds to the current T . Geometrically, integral 2-currents can be thought of as oriented surfaces with multiplicities. For instance, if T and T are two oriented surfaces such that orientations of T and T agree on T T , T +T can be represented as T T endowed with the same orientation as T and T , points in T T having a multiplicity equal to 2. If orientations of T and T do not agree, then summing T and T yields a cancellation on T T . T T T T +T T T

2.2.1.2 Integration.
2-differential forms can be integrated on oriented surfaces, in the same way vector elds can be integrated on oriented curves. To see how, let T be a oriented surface in S and, for each x T , let (ux , vx ) be a direct orthonormal frame of the tangent plane Tx T . The integral of a 2-differential form on T is dened to be :
T

For instance, one has T aT = area(T ), which is why area forms are called this way.

x (ux , vx )dx

Figure 2: Sum of integral currents.

2.2.1.3 Change of variable.

A change of variable is merely the data of a diffeomorphism : S S where S is the manifold where the new variables live. Using such a map, a 2-differential form on S can be transformed into a 2-differential form on S , by a process called pullback : D EFINITION 8. The pullback of by , denoted by is given by : x (u, v) = (x) (Dx (u), Dx (v)) for all x S and u, v Tx S . In a certain sense, pulling a 2-differential form back amounts to expressing it in terms of the new variables. The change of variable formula relates the integral of a 2-differential form with the one of its pullback. The result turns out to be particularly simple :
S

Now set S = 3 S 2 . S is obviously a subset of 3 3 . We will call the rst factor of the latter product the point space, Ep , and the second one the normal space, En . The reason for this is that an element of S can be thought of as a point in space together with a unit normal vector. If u is 3-vector, un will denote the vector (0, u) Ep En , and up the vector (u, 0) Ep En . Rigid motions of 3 can be naturally extended to S : if g is such a motion, one can set g (p, n) = (g(p), g (n)), where g is the rotation associated with g. We now dene two particular 2-differential forms on S : D EFINITION 9. Let (p, n) S and x, y 3 such that (x, y, n) is a direct orthonormal frame of 3 . We set :
H (p,n) G (p,n)

2.2.3 Invariant 2-forms

= xp y n + xn y p = xn y n

(S )

(3)

For example, if S = S = 2 and h an integrable function from S to , applying (3) to = haT yields = Jac()h aT : (3) thus generalizes the classical change of variable formula. For this formula to hold, need actually not be a diffeomorphism from S to S ; the only requirement is that should be a diffeomorphism from S to (S ).

One can actually check that these 2-forms do not depend on the choice of x and y. Moreover, they are invariant under rigid motions, that is satisfy g = for all rigid motion g. Geometric interpretations of these forms will be given in section 2.5.4. The dimension of the space of invariant forms is actually 4 [16].

2.3 Smooth case

The theory of normal cycles is inspired by the same ideas as the one presented in section 2.1, but transposed in a setting where they can be generalized : the theory of currents. Loosely speaking, normal cycles are a way to unfold offsets in a higher dimensional space : D EFINITION 10. The normal cycle N (V ) of V current associated with the set : ST V = {(p, n(p))|p M } Ep En endowed with the orientation induced by the one of M .
4 we will sometimes abuse the terminology and write the normal cycle of the oriented surface M instead. 4

2.2.2 Integral 2-currents

Integral 2-currents generalize oriented surfaces [12]. They can be formally dened as linear combinations of oriented surfaces with integral coefcients. In particular, any oriented surface T can be considered as an integral 2-current, which we will abusively also denote T . Integration of 2-differential forms is extended to integral 2-currents by linearity :
nT +pT

is the

=n

+p

315

M and ST V are obviously diffeomorphic via the map : i: M p ST V (p, n(p))

In particular, when V is convex and smooth, this denition agrees with the one given in the previous section. We now state a crucial property of the normal cycle, which we could have stated in the smooth case as well : the additivity. P ROPOSITION 3. Let V1 and V2 be two convex sets in 3 such that V1 V2 is convex. Then : N (V1 V2 ) + N (V1 V2 ) = N (V1 ) + N (V2 ) P ROOF. It is sufcient to show that the multiplicities of any point (p, n) in N (V1 V2 )+N (V1 V2 ) and N (V1 )+N (V2 ) agree. If p does not belong to V1 V2 , this is obvious. If p lies in V1 V2 , one concludes easily by noticing that N CV1 V2 (p) = N CV1 (p) N CV2 (p) and N CV1 V2 (p) = N CV1 (p) N CV2 (p).

The connection between normal cycles and curvature measures lies in the following lemma : L EMMA 2.
N(V ) G |i(BM ) H |i(BM )

= =

G (B) V H (B) V

N(V )

for all (Borel) set B 3 . Here |i(BM ) denotes the restriction of to i(B M ), that is the form which coincides with on i(B M ) and vanishes elsewhere. In words, curvature measures of a surface can be recovered by integrating specic differential forms on its normal cycle. P ROOF. By denition we have :
N(V ) G |i(BM ) =

V1

V1 V2

V2

G
i(BM )

The change of variable formula now states that :

i(BM )

N (V1 )

N (V2 )

G =

i G
(BM )

Figure 3: Additivity of the normal cycle In gure 3 normal cycles are graphically represented by their image under the map sending (p, n) Ep En to p + n.

To prove the rst claim, it is thus sufcient to show that : i G = GaM Let (u, v) be a direct orthonormal frame of Tx M , where x M . By denition, we have : Dx i(w) = wp + Dx n(w)n
G Expressing i(x) in the frame (un , v n , nn ), we get x

2.5 Polyhedral case

2.5.1 Denition
Once we know what the normal cycle of a convex is, there is at most one way of dening the normal cycle of a polyhedron while keeping the additivity property. Indeed, if one is given a triangulation of the polyhedron V into tetrahedra ti , i = 1..n, the normal cycle of V has to be :

(i G )x (u, v)

= =

un .(up + Dx n(u)n ) un .(vp + Dx n(v)n ) u.Dx n(u) u.Dx n(v)

vn .(up + Dx n(u)n ) vn .(vp + Dx n(v)n )

v.Dx n(u) = G(x) v.Dx n(v)

N (V )

=
n=1

(1)n+1
1i1 <..<in n

The proof of the second equality is similar. We omit it here as we will prove a stronger result in section 3.

N (n tij ) j=1

2.4 Convex case

When V is convex, a normal cycle can be dened even if M is not smooth. Indeed, in place of normal vectors, we can consider normal cones : D EFINITION 11. The normal cone N CV (p) of a point p V is the set of unit vectors v such that : q V 0 pq.v D EFINITION 12. The normal cycle of M is the current associated with the set {(p, n)|p V n N CV (p)} endowed with the orientation induced by the one of V .

by inclusion-exclusion. We will give in 2.5.3 a geometric description of the obtained current that does not depend on the chosen triangulation, so that N (V ) is well-dened in the polyhedral case.

2.5.2 Simplices
Let us now describe the normal cycle of the polyhedron V . The way it is dened suggest to look rst at the normal cycle of simplices. Remember that intuitively, these are unfolded versions of offsets of simplices. Just as their offsets, normal cycles of simplices can be decomposed into spherical parts, cylindric parts, and planar parts. The difference is that these parts now live in Ep En . We will say that a subset A of Ep En lies above a subset B 3 if the projection of A on the point space is included in B. Let us now describe in turn each type of part for a simplex S of varying dimension :

316

 spherical parts lie above vertices of S. They are subsets of {(p, n)| ||n|| = 1} where p is the considered vertex. If S is reduced to p, then the spherical part is a whole sphere. In case S is an edge, then each spherical part is a half sphere. When S is a triangle, they are spherical 2-gons, and if S is tetrahedron, they are the spherical triangles spanned by the normals of neighboring facets. Edges of these spherical polygons are dual to edges of S, and the external angle between two incident spherical edges equals the angle between corresponding dual edges. cylindric parts lie above edges of S. They are included in {(p, n)|p e, ||n|| = 1, n.e = 0}, e being the considered edge. If S is reduced to e, the cylindric part is a whole cylinder. If S is a triangle, it is a half-cylinder, and if S is tetrahedron, it is a portion of cylinder whose section is a circle arc joining the normals to incident facets. planar parts lie above facets. They have the form {(p, n)| p t} where t is the considered triangle and n a unit normal vector to t. If S is reduced to t, both possible orientation for n have to be taken into account, whereas if S is a tetrahedron, one should only consider the outward normal.

Figure 4: Normal cycle above a concave edge

L EMMA 4. V (p, h) = (St+ (p, h)) V where is the Euler characteristic and St+ (p, h) is the upper V star of p, that is the union of relative interiors of cells of V incident on p and lying in the half plane {x| 0} px.h P ROOF. One checks easily that both sides coincide when V is a simplex. The result then follows as both sides have the additivity property with respect to V . Note that this quantity, also called the index of p with respect to the direction h [6][5], is always smaller than 1 if p is regular in V .

2.5.3 General case

We can now go back to the case of a general polyhedron V . To begin with, for any point p lying in the interior of V , there is a triangulation of V such that p lies in the interior of a tetrahedron. Thus, there is nothing lying above the interior of V in N (V ). By a similar argument restricted to a face f of M = V , the part of N (V ) lying above f is the planar part f n, where n is the outward normal to f . The two remaining cases are slightly more involved.

2.5.4 Curvature measures for polyhedra

Curvature measures for a polyhedron V are dened by integration of corresponding invariant forms on N (V ), just like in the smooth case. Thanks to the structure of N (V ), it is sufcient to compute the integral of these forms on spherical, cylindric and planar parts : a tangent plane to a planar part is spanned by two vectors of E p . Applying G or H to a couple of two such vectors yields determinants with at least one zero column. Planar parts thus do not contribute to the curvature measures H and G , as could be expected. the tangent plane to a cylindric part at a point (p, n) is spanned by up and v n , where u is a vector parallel to the corresponding edge and v is orthogonal to u. For the same reason as above, G vanishes when applied to (up , v n ). u and v can be chosen so that (u, v, n) is H a direct orthonormal frame. Expressing (p,n) in the p n frame (u , v ), one obtains :
H (p,n) (up , v n ) =

2.5.3.1 Above edges.

If e is a convex edge of V , V can be triangulated in such a way that e is included in only one tetrahedron. Above e, N (V ) thus coincides with the normal cycle of this tetrahedron : we get a cylindric part delimited by the normals to faces incident on e, with multiplicity 1. If e is concave, one can nd a triangulation of V such that e is an edge of exactly two tetrahedra t and t . Above e, N (V ) is the sum of the cylindric parts of N (t) and N (t ) lying above e minus the cylindric part of N (t t ) lying above e. The above picture shows a cross section along a plane perpendicular to the concave edge e. As can be seen, we get again a cylindric part delimited by the normals to facets incident on e, but with a multiplicity equal to 1, that is with reverse orientation.

2.5.3.2 Above vertices.

Above a vertex p, the situation is more involved, as we obtain a linear combination of at least degree of p half-spheres, spherical 2-gons and spherical triangles. In magnitude, one can get arbitrarily large multiplicities if p is not supposed to be convex. A full description of the part of N (V ) lying above p can be given by computing the multiplicity V (p, h) (or (p, h) for short) of (p, h) in N (V ) for each unit vector h :

up .up up .v n

v n .up un .up + n n v n .v n u .v

v p .up =1 v p .v n

The integral of H over a subset of a cylindric part thus equals the area of this subset. a tangent plane to a spherical part is spanned by two vectors of E n . Thus, integrating H on a subset of such

317

a part yields 0. Integrating G yields the area of the subset, by a computation similar to the one given above. of a subset B is the The curvature measure sum of the areas of cylindric parts of N (V ) lying above B, weighted by their multiplicities. By the description of N (V ) given in 2.5.3.1, one obtains indeed the formula (2). G (B) is obtained by summing the areas of spherical parts V lying above B weighted by their multiplicities. Let us do the computation for parts lying above a vertex p B. V can be triangulated such that all tetrahedra incident on p share an edge pq. These tetrahedra can be numbered in a circular order around pq, say ti , i = 1..n. Let pi be the common vertex of ti , ti+1 and M , considering indices mod n, i be the angle pi1 ppi and i = pi pq. For each simplex S incident on p, the area of the spherical part SP (S) of N (S) lying above p is : H (B) V
3

3.1 Two more 2-differential forms

We now dene, for each couple of 3-vectors X and Y , two 2-differential forms on E p E n from which we will recover the second fundamental form. D EFINITION 13. Given a point (p, n) E p E n we set :
X,Y (p,n)

= (np p X p ) Y n = X p (nn n Y n )

(p,n) X,Y

where n and p respectively denote cross products in E n and in E p . Note that these two forms are bilinear in X and Y , but not symmetric. However, we will see that integrating them on normal cycles yields symmetric bilinear forms.

S = pq : 2 S = pi pq : 2( i ) as SP (S) is a spherical 2-gon with angle i at its vertices S = ti : 2 i i1 i by the formula giving the area of a spherical triangle as a function of its angles Let us now apply the inclusion-exclusion principle to nd the coefcient of each of the areas described above in the linear combination giving G (p). Areas of SP (ti ) appear once V each. Intersecting two tetrahedra ti1 and ti yields the triangle pi pq ; as these are obtained exactly once SP (pi pq) has coefcient 1. The remaining (n2 3n)/2 pairwise intersections all equal pq. For k 3, k-fold intersections also equal pq. Hence, the coefcient of SP (pq) is : n2 3n (1)k+1 + 2 k=3
n

3.2 Smooth case

L EMMA 5. If M is smooth, then : |i(BM ) X,Y
X,Y |i(BM )

N(V )

= HV (B)(X, Y ) = HV (B)(X, Y )

N(V )

n k

=1

P ROOF. As in the the proof of lemma 2, we perform a change of variable in the left-hand side. To compute i X,Y at a point p M , we consider the direct orthonormal frame (e1 , e2 , n) of 3 where e1 and e2 are principal directions and n = n(p). If the principal curvatures associated with e1 and e2 are respectively 1 and 2 , we have : X,Y i p (e1 , e2 ) = = (p,n) (ep + 1 en , ep + 2 en ) X,Y 1 1 2 2 X.e1 X.e2 det(n, Y, 1 e1 ) det(n, Y, 2 e2 )

Finally, we have :
G (p) V
n

=
i=1

(2 i i1 i )
i=1 n

2( i ) + 2

= = =

2
i=1

that is the classical denition of the angle defect at p. This computation thus agrees with the denition given by equation (1). Note that unlike the mean curvature measure, the gaussian curvature measure is independant on the orientation, in the sense that choosing the complement of V instead of V would yield the same measure.

X1 1 Y2 X2 2 Y1 2 X1 Y1 + 1 X2 Y2 HV (p)(X, Y )

where Xi and Yi , i = 1, 2, 3 are the components of X and Y in (e1 , e2 , n). The rst claim thus follows. The second one can be proved in a similar way.

3.3 Polyhedral case

The fact that integrals of the second fundamental form can be recovered from the normal cycle of a smooth surface by integration of 2-differential forms enables us to dene corresponding objects for polyhedral surfaces, as they also have a normal cycle. The next lemma justies the denition 4 : L EMMA 6. If V is a polyhedron, then : |i(BM ) X,Y
X,Y |i(BM )

3. THE SECOND FUNDAMENTAL FORM VIA THE NORMAL CYCLE

The concept of normal cycle was introduced to dene mean and gaussian curvature measures for a general class of objects. In this section, we show that it can actually provide a complete description of the curvature of an object. Not surprisingly, integrating invariant forms on the normal cycle yields integrals of invariants of the Weingarten endomorphism, namely its trace and determinant. The basic idea here is to integrate non invariant forms in order to obtain integrals of each coefcients in the Weingarten endomorphism matrix.

N(V )

= HV (B)(X, Y ) = HV (B)(X, Y )

N(V )

318

P ROOF. Clearly, X,Y vanishes on planar and spherical parts. Let e be an edge of M , or a segment included in such an edge, and CP be the cylindric part of N (V ) lying above e. The tangent plane to CP at (p, n) CP has a direct orthonormal frame of the form (up , v n ), where u is a unit vector parallel to e. We have : (p,n) (up , v n ) X,Y = = X.u 0 det(n, Y, 0) det(n, Y, v)

(X.u)(Y.u)

11111111111111111111111111 00000000000000000000000000 11111111111111111111111111 00000000000000000000000000 11111111111111111111111111 00000000000000000000000000 11111111111111111111111111 00000000000000000000000000 11111111111111111111111111 00000000000000000000000000 C 11111111111111111111111111 00000000000000000000000000 11111111111111111111111111 00000000000000000000000000 11111111111111111111111111 00000000000000000000000000 11111111111111111111111111 00000000000000000000000000 11111111111111111111111111 00000000000000000000000000 11111111111111111111111111 00000000000000000000000000 11111111111111111111111111 00000000000000000000000000 11111111111111111111111111 00000000000000000000000000 11111111111111111111111111 00000000000000000000000000 11111111111111111111111111 00000000000000000000000000 11111111111111111111111111 00000000000000000000000000 11111111111111111111111111 00000000000000000000000000 11111111111111111111111111 00000000000000000000000000 f 11111111111111111111111111 00000000000000000000000000 11111111111111111111111111 00000000000000000000000000 11111111111111111111111111 00000000000000000000000000 11111111111111111111111111 00000000000000000000000000 11111111111111111111111111 00000000000000000000000000 11111111111111111111111111 00000000000000000000000000 11111111111111111111111111 00000000000000000000000000 11111111111111111111111111 00000000000000000000000000 11111111111111111111111111 00000000000000000000000000 11111111111111111111111111 00000000000000000000000000 11111111111111111111111111 00000000000000000000000000 11111111111111111111111111 00000000000000000000000000 E 11111111111111111111111111 00000000000000000000000000 11111111111111111111111111 00000000000000000000000000 11111111111111111111111111 00000000000000000000000000 11111111111111111111111111 00000000000000000000000000 11111111111111111111111111 00000000000000000000000000

As a function of X and Y , (p,n) (up , v n ) is thus a symmetric X,Y bilinear form with 1 as unique non-zero eigenvalue and u (or e) as associated eigenvector. It thus equals /length(e)2 . e e Integration on CP yields :
CP

Figure 5: Proof of theorem 1

X,Y

(e) = e e length(e)

Then the mass M(T ) of T is dened to be the quantity :

n

M(T ) =
i=1

|i |area(Ti )

and the rst result follows. The derivation of the second one follows the same lines and is left to the reader. Both in the smooth and polyhedral case, anisotropic curvature measures generalize the mean curvature measure : indeed, the trace of HV (B)(X, Y ) or HV (B)(X, Y ) equals H (B). V

Intuitively, the mass of a current is its unsigned area, taking multiplicities into account. We now bound the mass of normal cycles in the case of restricted Delaunay triangulations. The following lemma is due to Amenta and Bern[2] : L EMMA 7. The angle between the normal to a triangle t in T and the normal to M at a vertex of t is O(r(t)) for xed M. In order to shorten notations we set : s(B) =
{tT, tB}

4. PROOF OF THE THEOREM

The idea behind the proof of theorem 1 is roughly as follows. Let E denote the part of N (T ) lying above B and D be the part of N (M ) lying above (B) (gure 5). Consider for simplicity that E and D are oriented surfaces, though it is not really accurate as they actually are currents. By lifting the projection to E p E n , one obtains a map f from E to D. Dene C to be the union of all line segments joining points of E with their image under f . C is a volume whose boundary is the union of E, D, and a surface A which is the union of all line segments joining points of E with their image under f . By applying Stokes theorem to C and a 2-differential form , one can express the difference between integrals of on E and D as the integral of on A plus an integral on C. In particular, when is a form associated with some curvature measure, this implies that the difference between the considered curvature measures of E and D is the sum of an integral on A plus an integral on C. In this particular case, the quantities to be integrated on E and D are bounded. Thus, to get a bound on the difference of curvature measures, it is sufcient to bound the volume of C and the area of A. To do so we rst bound the area of E and the length of E, which is the purpose of the next subsection. The desired bound then follows from the fact that f has bounded derivatives, and that the distance between points of E and their image under f is O(). Subsection 4.2 resumes the whole proof in a sketchy but more accurate way.

r(t)2 r(t)
{tT, tB, tB=}

sd(B) =

L EMMA 8. The mass of the part N (T )(BE n ) of N (T ) lying above B is O(s(B)). P ROOF. This mass can be decomposed in three terms : the mass lying above the interior of the triangles of T , M t , the mass lying above the interior of the edges of T , M e , and the mass lying above the vertices of T , M v . M t is merely the area of B, so it is O(s(B)). Let us now focus on M e . We have : Me =
e

|(e)| length(e) edge of B

4.1 Bounding the mass of normal cycles

We start with one more denition about currents : D EFINITION 14. Let T be an integral 2-current. Suppose T can be written as a linear combination of surface patches Ti whose relative interior are pairwise disjoint :
n

Let e be an edge of T and t, t be the triangles of T incident on e. The dihedral angle at e is O(r(t) + r(t )) by 7, as well as the length of e. Thus we also have M e = O(s(B)). The last quantity to consider is M v . Let u be a vertex of T , and ui , i = 1..n its neighbors in circular order. If ni is the unit normal to triangle uui ui+1 , then the mass lying above u is smaller than the sum of the areas of spherical triangles n(u)ni ni+1 . By lemma 7, the area of any such triangle is O((r(uui ui+1 ) + r(uui+1 ui+2 ))2 ). Summing on all u B, we get that M v = O(s(B)). We thus proved the announced claim. L EMMA 9. The mass of (N (T )(BE n )) is O(sd(B)).

T =
i=1

i Ti

319

P ROOF. This mass decomposes into the mass lying above triangles and the one lying above edges. The rst one is obviously O(sd(B)) and the second also by lemma 7. If the sampling is locally uniform in the sense of [3], then triangles of T are well-shaped. As a result, we have s(B) = O(area(B)) and sd(B) = O(length(B)).

4.2 Homotopying normal cycles

From now on, we will use the notations of [12] or [10]. Let f be dened by the following commutative diagram : SptN (V ) U E n
f

the polyhedral surface is the restricted Delaunay triangulation of a sampled smooth surface, we showed that these curvature tensors converge linearly towards the ones of the smooth surface, under a reasonable sampling condition. Choosing a small neighborhood of a given vertex as averaging region provides an estimation of the smooth curvature tensor at this vertex. In practice, this method is fast and provides very nice results. However, it still is not very clear how the neighborhood should be chosen. It would be interesting to try to choose it so as to minimize the theoretical error bound.

Acknowledgements
We are strongly indebted in Mariette Yvinec for reading the manuscript thoroughly.

p1

M U where U is the open set where is dened, p1 is the projection from E p E n to E p , and i(x) = (x, n(x)) for all x M . For simplicity we will denote the current N (T )(B E n ) by D and the current N (V )((B) E n ) by E, as in gure 5. Consider the afne homotopy h ([10] pp. 364) between f and the identity. As SptD is contained in p1 (U ), we can dene a 1 3-current C by : C = h ([0, 1] D) By the homotopy formula for currents, we have : C = f (D) D h ([0, 1] D) Now the image of the support of D by f is included in SptE and the one of SptD is included in SptE. By the constancy theorem ([10] 4.1.31), f (D) is a multiple of E. Moreover, comparing the multiplicities at points lying above images by of points in the interior of triangles contained in B shows that the proportionality constant is 1 : f (D) = E. As a consequence, the at norm of D E satises : F(D E) M(C) + M(h ([0, 1] D)) Also : M(C) M(h ([0, 1] D)) M(D)sup|f Id|sup||Df ||2 M(D)sup|f Id|sup||Df ||2

5. REFERENCES

by the inequality mentioned at page 364 in [10]. ||Df ||2 is less than ||D||2 ||Di||2 so sup||Df ||2 can be bounded by a function of the maximal curvature of M on (B). Let now x = (p, n) U E n . By denition f (x) = ((p), n((p))). It can easily be shown that |p (p)| = O(). Also by lemma 7 |n(p) n((p))| = O(). Plugging in the bounds on the masses obtained in the previous section, one gets that F(D E) K with the notations of the theorem. One concludes by noticing that the forms X,Y , X,Y are bounded and have bounded derivatives. In order to get the bound for the gaussian curvature measures, one needs rst to extend G to a bounded 2-form on E p E n having a bounded derivative, which can be done.

Conclusion
We presented a mathematically sound way to dene average curvature tensors for each region of a polyhedral surface, based upon the theory of normal cycles. Moreover, in the case where

[1] Nina Amenta, Sunghee Choi, Tamal Dey and Naveen Leekha. A simple algorithm for homeomorphic surface reconstruction. Proc. 16th Annu. ACM Sympos. Comput. Geom., pages 213-222, 2000. [2] Nina Amenta and Marshall Bern. Surface reconstruction by Voronoi ltering. Discrete Comput. Geom., 22(4):481-504, 1999. [3] S.Funke and E.A.Ramos Smooth-Surface Reconstruction in Near-Linear Time, to appear in SODA 2002. [4] H. Edelsbrunner and N. R. Shah. Triangulating topological spaces. Int. J. on Comp. Geom., 7:365-378, 1997. [5] Herbert Edelsbrunner, John Harer, and Afra Zomorodian. Hierarchical Morse complexes for piecewise linear 2-manifolds. In Proc. 17th Annu. ACM Sympos. Comput. Geom., pages 70-79, 2001. [6] Th. Banchoff, Critical points and curvature for embedded polyhedra, J. Diff. Geom 1 (1967) 245 256. [7] M. Berger, B. Gostiaux, G om trie diff rentielle : vari t s, e e e ee courbes et surfaces, Presses Universitaires de France. [8] H.Cartan, Cours de calcul diff rentiel, Hermann. e [9] H. Federer, Curvature measure theory, Trans. Amer. Math. Soc 93 (1959) 418 491. [10] H. Federer, Geometric Measure Theory, Springer-Verlag, New York, 1983. [11] J. Fu, Convergence of curvatures in secant approximations, J.Differential Geometry 37 (1993) 177 190. [12] F. Morgan, Geometric measure theory, Acad. Press, INC. 1987. [13] J. Steiner, Jber Preuss. Akad. Wiss. 114-118,(1840). In Gesammelte Werke, vol 2, New York, Chelsea 1971. [14] P. Wintgen, Normal cycle and integral curvature for polyhedra in Riemmannian manifolds, Differential Geometry (Gy. Soos and J. Szenthe, eds.), North-Holland, Amsterdam, 1982. [15] M. Z hle, Integral and current representations of Federers a curvature measures, Arch. Math. (Basel) 46, (1986), 557-567. [16] J.M. Morvan, On generalized curvatures, in preparation. [17] S. Petitjean, A survey of methods for recovering quadrics in triangle meshes, accepted. [18] M. Desbrun, M. Meyer, P. Schr der and A. Barr Discrete o differential-geometry operators in nD, preprint, the Caltech Multi-Res Modeling Group. [19] D. Meek and D. Walton, On surface normal and Gaussian curvature approximation given data sampled from a smooth surface, Computer-Aided Geometric Design 17, 521-543. [20] G. Taubin, Estimating the Tensor of Curvature of a Surface from a Polyhedral Approximation, Fifth International Conference on Computer Vision (ICCV95). [21] P. Alliez, D. Cohen-Steiner, M. Desbrun, O. Devillers and B. L vy, Anisotropic Polygonal Remeshing, to appear in e SIGGRAPH 2003.

320

Figure 6: Directions of minimal curvature estimated on a mesh of Michelangelos David. For each vertex, the averaging domain used for computations is the 2-ring of that vertex.

321