Corrected Curvature Measures

This paper proposes a new mathematical and computational tool for inferring the geometry of shapes known only through approximations such as triangulated or digital surfaces. The main idea is to decouple the position of the shape boundary from its normal vector field. To do so, we extend a classical tool of geometric measure theory, the normal cycle, so that it takes as input not only a surface but also a normal vector field. We formalize it as a current in the oriented Grassmann bundle R3×S2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbb {R}^3 \times \mathbb {S}^2$$\end{document}. By choosing adequate differential forms, we define geometric measures like area, mean and Gaussian curvatures. We then show the stability of these measures when both position and normal input data are approximations of the underlying continuous shape. As a byproduct, our tool is able to correctly estimate curvatures over polyhedral approximations of shapes with explicit bounds, even when their natural normal are not correct, as long as an external convergent normal vector field is provided. Finally, the accuracy, convergence and stability under noise perturbation is evaluated experimentally onto digital surfaces.


Introduction
We address in this paper the problem of defining a notion of curvature on non-smooth sets, with the goal of estimating the curvatures of a smooth surface through a nonsmooth approximation such as a triangulation or a digitization. In order to tackle this problem, there are two main issues: (i) the first one is to define consistent notions of curvatures; (ii) the second is to show that these newly defined curvatures are stable, namely that one can can estimate the curvatures of a smooth object from an approximation. This problem has various applications in computer science, in particular in geometry processing, computer graphics or digital imaging where discrete surfaces are ubiquitous.
There is a long history on the generalization of the curvatures in non-smooth geometry. One important work on this topic is the seminal paper of Federer [12], who first defined curvature measures for sets with positive reach. This generalizes the notion of Gaussian and mean curvatures for convex and smooth objects, but it unfortunately does not apply to triangulations. Using the notion of normal cycle introduced by Wintgen [34], this notion was then extended to a wider class of objects including triangulations, digitized objects and subanalytic sets [15]. The normal cycle of a general shape in R n is a current in the Grassmann bundle R n × Gr(1, n) (where the Gr(1, n) S n−1 factor stands for the normal cone) that encodes the geometry of the shape and allows to define curvature measures.
The stability of curvature measures has been well investigated the last 25 years. It is known that the curvature measures of a smooth object can be approximated by the ones of a triangulation, provided that the points and the normals of the triangulation are close to the ones of the smooth object [9,14]. There exist different kinds of stability or convergence results for curvatures or curvature measures, including anisotropic curvature measures, which were introduced to estimate principal curvatures and directions [9,26]. We may also quote the work of [5] which extends the normal cycle to arbitrary cloud of points by using offset surfaces, also with stability results. However these approaches do not provide a sound definition of curvatures when the naive normals do not converge toward smooth normals [5,17]. This is the case for instance for the famous counterexample of the Schwarz lantern (see Sect. 3.3) that converges to a cylinder in the Hausdorff sense, but whose normals diverge. It also fails utterly in the case of digital approximations of surfaces, since the naive normal vectors take only six different possible values, parallel to the axes.
The key idea of this article is to replace the normal vector field of a surface S by another vector field u which we assume to be geometrically more meaningful. For instance if S is a digitization of a smooth surface X , one may take for u a local average of the naive normals of X . The contribution of the paper is therefore the following: -We extend the notion of normal cycle of a surface to a couple (S, u) where S is a piecewise C 1 -surface of R 3 and u is a piecewise C 1 unit vector field and we call it the corrected normal current N(S, u). Our corrected normal current N(S, u) allows to define new curvature measures that we call corrected curvature measures. -We derive stability results for our corrected curvature measures with explicit bounds. In particular, we show that the corrected curvature measures of (S, u) approximate well Federer's curvature measures of a smooth surface X , provided that S is close to X in the Hausdorff sense and that u is close to the normal vector field of X . -We apply our results to the digital case and show that that it allows to define convergent pointwise estimators for the mean and Gaussian curvature with explicit convergence rates. We show that our estimators outperform also in practice stateof-the-art methods like digital integral invariants [7].

Alternative Approach with Varifolds
Another notable mathematical tool for representing generic shapes is the varifold, which has been introduced to solve shape optimization problems like Plateau's problem [2]. A d-varifold is a Radon measure on R n × Gr(d, n), where Gr(d, n) denotes the Grassmannian manifold of unoriented d-planes in R n . Its first variation is indeed related to the mean curvature vector field [1]. Varifolds were recently proposed for surface approximation and mean curvature estimation in [3]. Both approaches, the one with (corrected) normal cycles and the one with varifolds, are similar in that they rely on the Grassmann bundle. They differ however since N(S, u) is an (oriented) integral current which possesses an additional combinatorial structure, which can be integrated against ambient invariants forms, yielding all the geometric curvatures, whereas the varifold defined in [3] is by definition unoriented and yields a vector valued mean curvature. Our framework is therefore less universal (for instance it does not encompass point clouds, which are a strong focus of the varifold approach). However, our explicit yet flexible construction allows us to use the homotopy lemma and compute actual convergence rates, which are further refined in a few test cases and which result in extremely efficient approximation rates. Furthermore, our results concern both Gaussian and mean curvature and can be extended to deal with principal curvatures and directions.

Related Works in Geometry Processing
Since accurate geometry estimations is an important step in many geometry processing tasks like sharp features extraction, surface approximation, shape segmentation, shape matching and identification, mesh denoising, the problem of having stable yet accurate normal and curvature estimations on discrete surfaces has been widely studied. We mention here some representative approaches.
The first methods for estimating curvatures on polyhedral meshes (generally triangulated surfaces) were based on local formula using a simple neighborhood of vertices (e.g. see [33] for a survey). However it was quickly shown that even the classical angle defect method for Gaussian curvature estimation does not converge in most of the cases [35]. Polynomial fitting was also popular to estimate curvatures, but, again, even the popular osculating jets of [4] are convergent under hypotheses that are not met in practice (like vertices lying on the ideal continuous surface).
For triangulated meshes, it was noted that convergence of normals was essential for getting convergence of area or other quantities [27]. Other authors then show the equivalence of convergence of normal fields, metric tensors, area, and Laplace-Beltrami operator [17]. This induces the weak convergence of mean curvatures, but only if the naive normals of the mesh are convergent.
Recognizing the weaknesses of local fitting approaches, integral methods (like in geometric measure theory) were explored also in this field. Integral invariants were studied as a way to estimate curvatures onto a mesh, by a proper eigendecomposition of a local covariance matrix [28,29]. Unfortunately, this method is sensitive to errors in position. The Voronoi covariance measure is another way to compute geometric information from arbitrary compacts [24], with stability results. It indeed carries information related to curvature, but it is unclear if it induces pointwise convergence of curvature estimates.
Our work is significantly different from all the previous approaches. Our definition of curvatures is valid for surfaces that are union of cells homeomorphic to a disk, equipped with a vector field that only needs to be C 1 per cell. We provide stability results for curvature measures with respect to a smooth surface. And given a convergent corrected normal vector field, we exhibit pointwise convergence for mean and Gaussian curvatures onto digital surfaces. Note that there exists convergent normal vectors onto digital surfaces [11,22], so our mathematical framework is effective and gives convergent mean and Gaussian curvatures. Last but not least, our experiments show that they are not only convergent but outperform the state-of-the-art.

Detailed Outline of the Paper
In Sect. 2, we formalize this new integral with currents in R 3 × S 2 , which we view as the set of couples (position, unit normal vector). The unit sphere S 2 stands for the oriented Grassmannian Gr (1,3), the set of unit vectors in R 3 . Note that Gr(1, 3) can be identified by duality to Gr (2,3), the set of oriented planes. The advantage in using currents is that they encompass both discrete and continuous objects. To each couple (S, u) where S is a surface and u a unit vector field along S we associate a current N(S, u) in R 3 ×S 2 (Definition 2.8). Intuitively we may view the current as a piecewise smooth oriented 2-surface in R 3 × S 2 ⊂ R 6 . In particular, when the vector field u is the naive normal field of S and S = ∂ V is the boundary of a domain of R 3 , then N(S, u) corresponds to the normal cycle of V [26]. Then it is known that the integral of invariant forms over the normal cycle yields the area, mean and Gaussian curvature integrals over S. (However, the local convergence is only true when the normal cone approximates the smooth normal, whereas using N(S, u) gives us the flexibility we need.) In Sect. 3, we define our geometric measures, or corrected Lipschitz-Killing curvature forms as the integral on N(S, u) of the classical invariant forms on R 3 × S 2 (Definition 3.1). We further show that these measures have explicit and simple expressions in the specific but crucial case of polyhedral surfaces and a corrected normal vector field that is constant per face (Proposition 3.8). They reduce to finite sums of quantities attached to each discrete cell of S, i.e., face for area measure, edge for mean curvature measure, and vertex for Gaussian curvature measure. We illustrate with the case of the Schwarz lantern, known to be problematic even for the area, and we show how our corrected measures give very accurate or even exact curvature measures for several natural choices of the corrected normal field u.
We then provide in Sect. 4 a stability result for our corrected curvature measures (Theorem 4.2). More precisely, we show that the corrected curvatures of some (S, u) approximate well the curvature measures of a surface X = ∂ V of class C 2 provided that the surface S is close to X in the Hausdorff sense and that the vector field u approximates well the geometric normal n of X . This result relies on several notions of Geometric Measure Theory: it uses the notion of reach [12] introduced by Federer that allows the projection of the support of N(S, u) onto the one of the normal cycle N(V ); the affine homotopy formula [13] for currents is a key tool that allows to bound the flat norm between the normal cycle N(V ) and the corrected normal cycle N(S, u); the result also makes use of the Constancy Theorem [25] that allows to solve multiplicity issues that appear when projecting the support of N(S, u) onto the one of N(V ).
As a byproduct, in the difficult case of digital surface approximations to smooth surface, we show that our normalized definitions of mean and Gaussian curvatures converge pointwise to the mean and Gaussian curvatures of the continuous surface S, provided u is estimated by a multigrid convergent normal estimator (Theorem 5.4). The theoretical convergence speed is O(h 1/3 ) if h is the sampling step of the digitized surface, and equals the best known bound for curvature estimation O(h 1/3 ) for digital integral invariant method [7].
Following this stability and convergence results, we confront theory with practice in Sect. 6 and present an experimental evaluation of our new curvature estimator (defined from curvature measures) in the case of digital surfaces, which is a good testbed for our framework since naive normals of S never converge toward the normals of the continuous surface X (a problem known as metrication in the case of the area measure). We show that our estimators outperform also in practice state-of-the-art methods like digital integral invariants [7] and convergence speeds in O(h 2/3 ) are reached in practice.
We conclude in Sect. 7 and outline several research directions, especially the extension of our work to anisotropic curvature measures.

Normal Cycle Corrected by a Vector Field u
We introduce in this section the notion of corrected normal current for piecewise C 1 surfaces endowed with a vector field in the three dimensional space, called the corrected normal.

Corrected Surface
We assume in the following that S is a closed piecewise C 1 oriented surface of R 3 , namely S = n i=1 S i where each S i is a compact surface homeomorphic to a disc. Furthermore, we assume that we have the following combinatorial structure: (i) each S i is called a face of S; (ii) for every i = j, S i, j = S i ∩ S j is either empty or a point or a non-degenerated connected C 1 curve; when not empty or reduced to a point, it is called an edge of S; (iii) the intersection of three (or more) faces of S is either empty or equal to a point that is called a vertex in the latter case.
Note that the orientation of S is equivalent to having an abstract orientation of the combinatorial structure. Indeed, each face S i induces naturally an orientation of its bounding curves S i, j , where j ranges over the faces S j adjacent to S i . This orientation turns S i, j into a directed edge, and can be represented in space by the unit vector field e i j along S i j , which is tangent to S i, j . If we are given a normal vector field n i on S i respecting its orientation, then n i × e i j points toward the inside of S i . Note that the adjacent face S j induces the opposite orientation on S i, j : e ji = −e i j .
We say that u : S \ {S i, j : i = j} → S 2 is a corrected normal vector field on S if: (i) its restriction u i to the relative interior of the face S i is of class C 1 , and extends continuously to S i ; (ii) along any curve S i, j = S i ∩ S j , u i and u j are not antipodal; (iii) at any vertex p = S i 1 ∩ . . . ∩ S i k , the corrected normals u i 1 , . . . , u i k lie in the same hemisphere. Note that, in this definition, we do not require u to be continuous or even defined on the edges. Definition 2.1 (corrected surface) We say that the couple (S, u) is a corrected surface if S ⊂ R 3 is a surface that satisfies the above assumptions and u is a corrected normal vector field on S.

Example 2.2
When S consists of a single C 1 face, we can choose the vector u to be the usual normal vector n. This generalizes easily to a piecewise C 1 surface. The non-degeneracy across vertices or edges prohibits cusps or cuspidal edges.

Example 2.3
When S is an oriented polyhedral surface, we may choose u to be the (unit) normal vector on each face. Here again, non-degeneracy across edges (resp. vertices) means that the dihedral angles are between −π and +π (resp. is a graph over the plane bounding the hemisphere).

Example 2.4
We can also consider any C 1 unit vector field u along S, even when S is a polyhedral surface, e.g. S can be a digital surface and u can be a bilinear interpolation at the normals given at the vertices (in which case u usually ceases to have unit length throughout, hence should be normalized).

Combinatorial Structure
Starting with the combinatorial structure of the surface S, i.e., the set of its vertices V , edges E and oriented faces F, together with their incidence relations, we define a new (abstract) combinatorial surface S * over which the corrected normal current will be constructed. A similar construction exists for convex polyhedra, and is known under various names: expansion, cantellation, or complete truncation. 1 The idea behind the construction is that each inner edge is blown up into a strip (i.e., a combinatorial quadrilateral), while each inner vertex of degree d is blown up into a face with d edges, see Fig. 1. More precisely, -the set F * is the setV ∪E ∪ F, whereV (resp.E) denotes the set of inner vertices (resp. edges).
The incidence relations are completely determined by describing the faces as an ordered list of vertices. Since a face f * ∈ F * is either a vertex p, an edge e of a face f of S, we consider separately all three cases.
-Whenever the face corresponds to an inner vertex p ∈ V , we denote it S * p and let S 1 , . . . , S n be an ordered list of the faces incident to p (the order is induced from the orientation and is unique, up to circular permutation). Then S * p is described by its vertices p * 1 , . . . , p * n where p * i = ( p, S i ). -Whenever the face corresponds to an inner edge S i, j = S i ∩ S j ∈ E, and S i, j joins the vertices p, q, with the convention that S i induces the orientation p → q on S i, j , and S j the opposite orientation. Then S * i, j is the quadrilateral face joining -Whenever the face S * i corresponds to S i ∈ F, and the vertices of S i (in order) are p 1 , . . . , p n , then S * i is given by its vertices p * 1 , . . . , p * n where p * k = ( p k , S i ). We slightly modify this definition by dividing all the faces S * p into triangles. More precisely, let S * p be a face of S * corresponding to a vertex p of S. The boundary of this face is composed of = d( p) edges that are common with faces S * i, j and denoted by s 1 , . . . , s ,: (i) we add a vertex p * at the interior of the face S * p ; (ii) for every 1 ≤ i ≤ , we add the triangle with vertex p * and opposite edge s i ; (iii) we add the edges corresponding to these triangles, see Fig. 1.

Corrected Normal Cone
We need to define the notion of corrected normal cone which is the image by a continuous map of the combinatorial surface S * into R 3 × S 2 . Its construction uses the corrected normal u and is done so that the corrected normal cone inherits the orientation of S * . Definition 2.5 Let (S, u) be a corrected surface and S * be the combinatorial surface of S whose orientation is inherited from S. The corrected normal cone NC( f * , u) of a face f * of S * is defined by: -if f * = S * i corresponds to a face S i in S, then NC( f * , u) = NC(S i , u) is the image of S i by the map p → ( p, u( p)); where Arc(u i ( p), u j ( p)) denotes the unique geodesic arc (on the sphere) between u i ( p) and u j ( p); -if f * is a triangle of S * p , then NC( f * , u) denotes the spherical triangle with vertices u( p), u i ( p), and u i+1 ( p) of area strictly less than 2π , where u( p) ∈ S 2 is the normalized average of the u i ( p).
Remark that for each oriented face f * of S * , the corrected normal cone NC( f * , u) inherits from the orientation of f * . Remark that for every face f * , NC( f * , u) is a surface of class at least C 1 in R 3 × S 2 . The corrected normal cone NC(S, u) can be seen globally as the image by a continuous map of the combinatorial surface S * into R 3 × S 2 . However it may not be an embedding nor even an immersion over edges and vertices; indeed, the image may have multiplicity above the vertices, and singularities above the edges.

Remark 2.7
The corrected normal cone NC(S * p , u) above a vertex p of S is made of d( p) spherical triangles, each of them having u( p) as a vertex. Note that the algebraic sum of these oriented triangles do not depend on this arbitrary point u( p). The corrected normal cone NC(S * p , u) can be seen as a set of point n ∈ {p} × S 2 with an integer multiplicity μ( p).

Corrected Normal Current
We may now define our variant of normal cycle, defined as a current with support given by the corrected normal cone. For the reader unfamiliar with the notion of currents, we recall briefly the main notions in Sect. 4.3.

Definition 2.8 (corrected normal current)
Let (S, u) be a corrected surface. The corrected normal current N(S, u) of (S, u) associates to every differential 2-form ω of R 3 × S 2 the real number where f * ranges over all the faces of S * .
The following proposition is an obvious consequence of the construction but is the heart of the notion introduced in this paper. Thanks to this property, the corrected normal current is globally coherent. In particular, this property will be central in the proof of stability results.
θ Fig. 2 The corrected normal current of a planar curve, seen as a discrete curve in R 3 . The Grassmann bundle R 2 × Gr(1, R 2 ) = R 2 × S 1 is represented as R 3 where the third coordinates is the angle θ . The smooth planar curve is a circular arc lifting to a piece of a helix. Its approximation by a digital curve can be lifted as the normal cycle (in purple), in which case the normal on each edge follows the axes (hence lies at height kπ/2 for some integer k); therefore the extra circular arcs at the vertices are vertical edges of length ±π/2. On the contrary, corrected normals on the digital curve are closer to the smooth normals; as a result, the corresponding lift (in blue) is also closer to the helix in R 3 , and the vertical edges are shorter.
Obviously, the corrected current is closer to the smooth lift. Note that the combinatorial structure of the current is extremely simple: each vertex is blown up into an edge construct a smooth vector field from few samples. Indeed unit length is not required to define the corrected normal cone over a face of S i . However over an edge S i, j or a vertex, it asks for different interpolation formulas, which we will not delve into in this article. Nevertheless, if u is continuous, edge and vertex contributions have zero Lebesgue measure, and the integral formulas below can be defined. Moreover they will converge to the same limit as u tends to the smooth normal n, even if u is not of unit length. This case will be illustrated in a forthcoming article.

Corrected Curvature Measures
Our goal is to obtain geometric information on the surface S, namely its area and curvatures, which are independent of the position of S in space. It is classically known in the smooth case that the area and curvature measures can be computed by integrating the invariant forms on the normal cycle; it also extends to discrete surfaces (see [15]). These invariant forms [34,36] form a basis of the 2-forms in R 3 × S 2 that are invariant by the action of the rigid motions of R 3 . They yield the area, the mean curvature and the Gaussian curvature.

Invariant Forms, Curvature Measures, Curvature Estimators
Let us now recall the expression of the invariant differential 2-forms of R 3 ×S 2 (see [26] for more details). Let ( p, w) be any point in the oriented Grassmann bundle R 3 × S 2 .
Clearly the tangent plane of R 3 × S 2 at the point ( p, w) is a space of dimension 5 spanned by the following vectors: where e 1 and e 2 are any vectors such that (e 1 , e 2 , w) form a direct orthonormal frame of R 3 . The set of differential 2-forms of R 3 × S 2 that are invariant under the action of rigid motions is a vector space of dimension four spanned by the following forms: where ε (ξ ) denotes the scalar product of ε and ξ .
The kth Lipschitz-Killing corrected curvature measure of S with respect to u associates to every Borel set B the real number where π : R 3 × S 2 → R 3 denotes the projection on the position space.

Remark 3.2 When
S is smooth and u = n the standard normal, μ 0 is the area density d A, while μ 1 and μ 2 are respectively the mean curvature and Gaussian curvature densities, −2H dH 2 and K dH 2 (up to a scalar factor). This result is reproven below.

Remark 3.3
The measure μ induced by the symplectic form ω plays a different role. Instead of measuring metric data, it tests whether a surface in the R 3 × S 2 is indeed the normal bundle of a surface in R 3 . Quite logically, it vanishes identically on smooth surfaces as soon as u = n. For an approximated surface, it measures how far the corrected normal u is from being normal to the faces. While pointwise never true, its averaged value yields a measure of the quality of the choice of corrected normal.
In order to recover classical geometric invariants, we further define:

Calculation of the Corrected Curvature Measures
Let (S, u) be a corrected surface. We recall that S = n i=1 S i is a union of faces of class C 1 , that the edges S i, j are of class C 1 and the set E(S) of edges of S is finite. The surface S is oriented. In order to state the proposition we need to introduce the following notations.

Notations
At every point p of the interior of a face S i , we denote by n( p) the unit oriented normal vector, by T p S the plane tangent to S, by e 1 ( p) a vector of T p S given by e 1 = (n × u)/ n × u if n and u are not collinear, and given by one of the principal directions otherwise. We also introduce e 2 = u × e 1 and e 2 = n × e 1 ∈ T p S. Note that (e 1 , e 2 , u) is a moving frame of R 3 associated with u, while (e 1 , e 2 , n) is a Darboux frame associated to the surface S i .
At every point p on an edge S i, j , we denote by e( p) = e i j ( p) the unit vector along the edge S i, j oriented 2 as the boundary of S i , by Ψ (p) = ∠(u i ( p), u j ( p)) the corrected dihedral angle between the faces S i and S j , and e 1 ( p) = (u i ×u j )/ u i ×u j . (The last vector is only defined when u i ( p) and u j ( p) are not collinear; however, whenever that happens, Ψ (p) vanishes and NC( p, u) drops in dimension, hence does not contribute to the integrals. Without loss of generality, we will assume that e 1 is always well defined.) In order to integrate differential forms over manifolds (more details are provided in Sect. 4.3), we will use the notion of Hausdorff measure. In the following, we denote by H m the m-dimensional Hausdorff measure.

Proposition 3.5 The corrected Lipschitz-Killing curvature measures of S with respect to u associates to every Borel set B the quantities
where Area (NC( p, u)) is the sum of the algebraic areas of each spherical triangle where the sign is induced by its orientation.
Proof We recall that u is of class C 1 . The measures have absolutely continuous and atomic parts; by additivity of the integral, we may consider separately the contributions of faces, edges and vertices.
Above the Faces. Remark that for any face S i of S the map Γ u : p ∈ S i → ( p, u( p)) is a diffeomorphism between S i and the corrected normal cone NC(S i , u). The change of variable formula implies that Hence, the computation of the curvature measures amounts to computing the pullback by Γ u of the corresponding curvature forms. In order to do that, we consider the orthonormal frame (e 1 , e 2 ) of the tangent plane T p S at p. Using the fact that dΓ u ( p) = (Id, du( p)), one gets Similarly, we evaluate Γ * u ω 1 and Γ * u ω 2 in the same basis: Finally, = e 1 | e 1 du · e 1 | e 1 e 2 | e 1 du · e 2 | e 1 + e 1 | e 2 du · e 1 | e 2 e 2 | e 2 du · e 2 | e 2 Above the Edges. Let p : [0, L] → S i, j be the arc length parameterization of S i, j that satisfiesṗ(t) = e( p(t)) (the dot denoting the derivative w.r.t. t). Then the corrected normal cone NC(S, u) ∩ π −1 (S i, j ) above the edge S i, j can be parameterized by Note that the map φ preserves the orientation. In order to shorten the equations below, we often (but not always) omit to specify that the quantities are considered at the point p or p(t). Note also that the reference frame at We then get We may note that μ 0 vanishes identically above the edges, since Span (∂φ/∂ψ, ∂φ/∂t) projects to a line on the position component. Note that this calculation can be easily localized over a Borel set B: if we calculate μ 1 (B ∩ S i, j ), it amounts to add the indicator function of B in the integrand. Similarly, one gets By differentiating the expression Similarly, one has d(e 1 × u i ) · e | e 1 = − e 1 × u i | de 1 · e and thus We note that the angle ψ between the two vectors u i and u j belongs to Therefore u i , u j , and e 1 × u i are all orthogonal to e 1 , so they belong to the same plane and we have This leads to Finally, v | e 2 ds dt Above the Vertices. When p is a vertex, its contribution is restricted to μ 2 , since the position is constant in NC( p). Since ω 2 measures the area in the velocity component, the vertex adds a Dirac mass at p to μ 2 , with coefficient equal to the (signed) area of NC( p).

Remark 3.6
In the smooth case, there are no edges nor vertices, and taking u = n (and hence e 2 = e 2 ), the formulas simplify. We recover since dn is symmetric on T p S.

Remark 3.7
Similarly to [9], we can consider a vector-valued equivariant two-form on the Grassmann bundle defined for each fixed pair of vectors X, Y ∈ R 3 . This leads to an anisotropic curvature measure μ X,Y which converges to the second fundamental form by the same theorems in Sect. 4. Study of this measure allows to compute and approximate principal curvatures and principal directions. This will be the topic of a forthcoming article.
In the case of polyhedral surfaces with a normal constant per face, we have the following simplifications, which we will use in our experiments (Sect. 6).
where α( f ) is the angle between the corrected normal u and the naive normal n, Ψ is the oriented angle between the corrected normal u i and u j incident to the edge S i, j and e 1 = (u i × u j )/ u i × u j , while e is the oriented unit tangent to the edge.

The Schwarz Lantern
To illustrate the remarkably fast rate of convergence of our approach, let us apply it to the Schwarz lantern L, a C 0 approximation by triangles of a cylinder C of radius r and height h, known for its failure to converge in curvature and even area (see Fig. 3.3). The failure occurs because the normals do not converge to the cylinder's (see [26] for the definition of the Schwarz lantern and an analysis of the problem). We will use two different choices of corrected normals on the Schwarz lantern. Fig. 3 The Schwarz lantern. The vertical cylinder of height h and radius r is cut horizontally along m + 1 evenly spaced circles; each circle is replaced by a regular n-gon (here m = 6 and n = 7). However, consecutive polygons are not parallel but rather obtained one from another by a screw motion of vertical translation h/m and angle π/n. Any vertex is connected to the two nearby vertices above and the two nearby vertices below (except for the top and bottom levels) We first consider the corrected normals to be constant on each triangle T and equal to the cylinder's at the top (or bottom) of the triangle. Then the angle α between the corrected normal and the face's normal satisfies tan α = (2mr/h) sin 2 (π/(2n)), so that the corrected area of the lantern is Edges are either horizontal, in which case the corrected normal is identical on the faces above and below, or slanted, in which case the two normals u i , u j are horizontal (so that e 1 = u i × u j is vertical) and at angular distance Ψ = π/n. For an edge e, the scalar product e | e 1 (e) is then h/m, thus μ 1 (e) = π h/(nm). Since there are exactly 2nm such edges, Finally, at a vertex p, the corrected normals are all horizontal, so that NC( p) has measure zero, as expected. We have thus an exact result for the mean and Gaussian corrected curvature measures and a O(1/n 2 ) convergent one for the area. Now consider the same triangulation with a different corrected normal u, defined as follows: at every vertex, it is the same (horizontal) normal as the one to the cylinder; along edges and faces, it is defined by linear interpolation followed by normalization to one. It is easy to see that these normals coincide with the normal to the cylinder of the horizontal projection π on C: u( p) = r −1 π( p). On a face, e 1 = u × n and e 2 = n × e 1 . Since dπ is the projection to the tangent plane to the cylinder, we have that du · e 1 = r −1 e 1 and du · e 2 = r −1 u | n e 2 .
The field u being continuous on the edges and at the vertices, the curvature measures are carried only by the faces. Using Proposition 3.5, we get: -the corrected area measure μ 0 = u | n dH 2 is no other than the pullback of the area form on the cylinder by the horizontal projection, i.e., π * d A C ; hence μ 0 measures the area of the projected lantern L onto the cylinder C, and globally gives exactly the area of the cylinder; -for μ 1 , we have that hence we obtain exactly twice the mean curvature density of the cylinder; -the corrected Gaussian curvature measure is also a pullback; but we may directly notice that the corrected normal is horizontal, hence stays on the equator of S 2 ; therefore μ 2 vanishes identically as in the previous example.

Stability of the Corrected Curvature Measures
We show in this section that the corrected curvatures of (S, u) approximate well the curvature measures of a surface X of class C 2 provided that the surface S is close to X in the Hausdorff sense and that the vector field u approximates well the normal n X of X . In order to be able to compare the two surfaces S and X , one first need to recall the notion of reach introduced by Federer [12]. Remark that this notion was originally introduced in order to define the notion of curvature measures.

Background on Sets with Positive Reach
The distance function d K of a compact set K of R d associates to any point x of R d its distance to K , namely d K (x) := min y∈K d(x, y), where d is the Euclidean distance on R d . For a given real number ε > 0, we denote by K ε := {x ∈ R d : d K (x) ≤ ε} the ε-offset of K . The Hausdorff distance d H (K , K ) between two compact sets K and K is the minimum number ε such that K ⊂ K ε and K ⊂ K ε . The medial axis of K is the set of points of x ∈ R d such that the distance d(x, K ) is realized by at least two points y and y in K . The reach of K , denoted by reach(K ), is the infimum distance between K and its medial axis. As a consequence, whenever the reach is positive and ε < reach(K ), the projection map is well defined. It is well know that smooth compact submanifolds have positive reach [12,26]. The following proposition will be useful for stating our main theorem. Proposition 4.1 Let X be a compact submanifold of class C 2 of R d . Then where ρ X is the largest absolute value of the principal curvatures of X . Furthermore π X is differentiable on X ε for any ε < reach(X ) and where ρ X (x) is the largest absolute value of the principal curvatures of X at the point x.

Stability Result
We provide in this section a stability result on the curvature measures, namely we show that the corrected curvature measures of (S, u) approximate the curvature measures of X provided that S and X are close in the Hausdorff sense and that u is close to the unit normal vector of X . In order to state the result, we denote by μ S,u k the corrected curvature measures of (S, u) and by μ X k the curvature measures of X , where k ∈ {0, 1, 2, }. Note that the curvature measures of X coincide with the corrected curvature measures of (X , n) where n is the geometric oriented unit normal of X . Then the corrected curvature measures of (S, u) are close to the curvature measures of X in the following way: for any connected union B = i∈I S i of faces S i of S, one has  and μ X k on union of faces instead of generic Borel sets. In particular, these faces may be large, which contradicts the idea behind the normalized corrected curvatureĤ u (B) and G u (B) of Definition 3.4. However, the faces can be subdivided at measure in order to give neighborhoods as small as necessary (extra edges and vertices will carry no curvature). Working on faces ensures a minimum of regularity in our theorem, without loss of generality. Theorem 4.2 allows us to estimate the rate of convergence of curvature measures for a sequence of approximations (S k , u k ) together with their corrected normals vector fields. However, some terms may prove difficult, especially N S v (B) and N S v (∂ B) which will tend to infinity very fast. We give below two corollaries of the previous stability result. The first one is specific to the case where the corrected normal field u is continuous. The second case concerns digital surface approximations, where u is constant per face.  H 1 (∂ B)).

Corollary 4.5 When S is a digital surface approximation with an estimated normal
vector field u constant per face, we have the following simplified estimate: where B is an union of N s surfels of edge size h, with N b boundary edges.
Proof In the digital case, a vertex has at most six neighbors so d max ≤ 6. (We then major 6 3 /(2π) by 35.) Since u is constant per-face, we have L u = 0. Furthermore, the length of each edge is h, the area of each face is h 2 , N S v (B) ≤ 4N s and the number of edges is less than 4N s . The formula follows.
Estimation of curvatures on digital surfaces is studied in more details in Sect. 5. The remaining of this section is devoted to the proof of Theorem 4.2.

Currents
We denote by D m the set of m-differential forms on R d . This set can be endowed with the C ∞ -topology [31] and the set of m-currents on R d is then by definition its topological dual. Therefore a current T : D m → R is a continuous linear map. The support of a current T : D m → R is the smallest closed set K such that spt(ω) ∩ K = ∅ ⇒ T (ω) = 0. The boundary ∂ T of a m-current T is the (m − 1)-current defined by ∂ T (ω) = T (dω), where d is the exterior derivative.
When U and V are open sets in Euclidean spaces, f : U → V is of class C ∞ , T is a m-current on U and f |spt(T ) is proper, one define the m-current f T by f T (ω) = T ( f * ω) (see [13, 4.1.9]). This notion is extended when f is a locally Lipschitz map (see [13, 4.1.14]).
Note that if a subset S of R d is m-rectifiable, then it is of class C 1 H m -almost everywhere and it is possible to integrate differential forms over it, and thus to define currents of the form where (e 1 , . . . , e m ) is any direct orthonormal basis of the tangent space at x.

Integral Currents
There exist several categories of currents, but we are mainly mentioning the one we use in this paper. An m-current T is said to be rectifiable if its support S = spt(T ) is m-rectifiable, compact and oriented and if for every ω ∈ D m (R d ), one has where μ(x) ∈ Z is the multiplicity and satisfies S μ(x) dH m (x) < ∞. When μ(x) = k is constant, we denote T = k · S. Note that a current can be rectifiable without having its boundary rectifiable. A current T is said to be integral if both T and ∂ T are rectifiable.

Flat Norm and Mass
We can also define semi-norms over the set of currents. The mass M(T ) of a m-current T is given by

Constancy Theorem
A key tool in the proof is the Constancy Theorem (see [13, 4.1.14] or [25, 3.13]). This theorem is important since it implies that the multiplicity of a current supported on a C 1 submanifold is constant. We state it for integral currents in the case where there is no boundary, even though it is true in a much more general setting.
Theorem 4.6 (Constancy Theorem) Let X be an m-dimensional oriented submanifold of R d of class C 1 with no boundary. If T is an integral current supported in X with no boundary, then there exists an integer λ such that

Proof of Theorem 4.2
The proof is based on the homotopy formula for currents and is in the same spirit as the proof of [9]. We recall that the projection map π X : X ε → X is well defined and differentiable since ε < reach(K ). We first build the Lipschitz map

Lemma 4.7 f N(S, u) = N(X ).
Proof By definition, f N(S, u) is a 2-current supported in the set spt (N(X )). Since N(S, u) has no boundary, f N(S, u) also has no boundary. Furthermore, spt (N(X )) is of class C 1 , so by the Constancy Theorem (Theorem 4.6), one has f N(S, u) = λ · N(X ), where λ is an integer. We know by assumption that there exists an open set U such that π X : U ∩ S → X is injective. This implies that the restriction of f to spt (N(S, u)) is one-to-one onto its image, so the constant λ = 1.
We denote by D = N(S, u) (B ×R 3 ) the restriction of the current N(S, u) to B ×R 3 . By Lemma 4.7, one has Let h be the affine homotopy between f and the identity: In the remainder of the proof, in order to shorten equations, we put x = ( p, u). Using the affine homotopy formula for locally Lipschitz map (see [23, p. 187] or [13, 4.1.14]), one has are respectively 3-current and 2-current that satisfy where the norm of the linear map By definition of the flat norm one has ).
(1) The following lemmas bound the terms involved in the right hand side term of the previous equation.

Lemma 4.8
For every x ∈ spt(D), one has Proof Clearly, by definition of the Euclidean norm, one has We can write f = g • π X • p 1 , where p 1 : X ε × R 3 → X ε is the projection and g : X → spt(N(X )) is given by g( p) = ( p, n X ( p)). Since for any p ∈ X , one has Using the fact that Dπ X (x) ≤ 1/(1 − ρ X (π X (x))ε) for every x in X ε , one has by composition: The conclusion follows from the fact that the upper bound is greater than 1.
Proof We can decompose the mass of D into three terms: the mass M f above the relative interior of the faces S i of S, the mass M e above the relative interior of the edges of S i, j , and the mass M v above the vertices of S.
Above the Faces of S. We recall that the corrected normal cone NC(S 0 i , u) above the interior S 0 i of the face S i is parameterized by Γ u : p ∈ S 0 i → ( p, u( p)). Since u is L u -Lipschitz on each face S i , the map Γ u is 1 + L 2 u -Lipschitz. The mass above S 0 i is given by the General Area-Coarea Formula [25, 3.13] NC(S 0 i ,u) where J 2 (Γ u ) denotes the 2-Jacobian. Summing over all the faces S 0 i of S, one has Above the Edges of S. We recall that π : spt(N(s, u)) → R 3 , defined by π( p, n) = p, denotes the restriction to the support of N(S, u) of the projection onto the position component. Note that the support of N(S, u) above an edge S i, j is given by spt (N(S, u) . The Coarea Formula thus gives: The third line follows from the fact that the 1-Jacobian of π satisfies J 1 (π ) = 1 and that π −1 ( p) is an arc of circle of length 2 arcsin ( u i ( p) − u j ( p) /2). Summing over all the edges S i, j one has Above the Vertices of S. Let p ∈ B be a vertex of S. The support of the corrected cycle N(S, u) ({ p} × R 3 ) above p obviously lies in the set { p} × S 2 . By construction, the multiplicity μ(n) of a point ( p, n) lying in this support is an integer that satisfies |μ(n)| ≤ deg( p), where deg( p) is the degree of p, namely the number of edges S i, j incident to p. By definition, one has u i ( p)−u j ( p) ≤ δ u for every pair of adjacent faces i and j containing p, which implies that the length of the boundary of the normal component of the support of N(S, u) ({ p} × R 3 ) is at most L = δ u deg( p). By the isoperimetric inequality on the sphere, since the support of the normal component above the point p lies in half a sphere, its area is at most L 2 /(2π) ≤ δ 2 u deg( p) 2 /(2π). Therefore, one has Summing over all the vertices of S ∩ ∂ B, one has

Lemma 4.10
Proof Since B = i∈I S i is a connected union of faces S i , its boundary is a union of edges S i, j where i ∈ I and j / ∈ I . The boundary ∂ D of the restricted current D = N(S, u) (B × R 3 ) is supported on curves of two different kinds.
Above the relative interior of each S i, j ⊂ ∂ B, the support D is parameterized by Above a vertex p of S ∩ ∂ B, the support of the D is the union of at most deg( p) arcs of circles. Each such arc of circle is a geodesic between two vectors u j 1 ( p) and u j 2 ( p), and is of length less than 2 arcsin ( u j 1 ( p) − u j 2 ( p) /2) ≤ 2 arcsin(δ u /2). Therefore The result follows by summing over all the vertices and edges of ∂ B.

Convergence of Curvatures on Digital Surfaces
Digital surfaces come from the sampling of Euclidean shapes in a regular grid. They appear naturally when analyzing 3D images (coming from tomography or MRI). But for estimating curvatures, the digital surfaces are the most challenging among all discrete surfaces. As boundaries of union of cubes, their geometric normals can take only six possible directions. We show in this section how our normalized corrected curvatures provide accurate pointwise approximations of the curvatures of the original Euclidean shape. More precisely, we show their multigrid convergence and give an explicit speed of convergence.
We begin by recalling useful notions of digital geometry and establishing a few lemmas linking the local geometry of the digitized surface and the local geometry of the continuous surface. We will then use the stability result of the previous section to establish multigrid convergence results.
In this section, let V be a compact domain of R 3 whose boundary X := ∂ V is of class C 3 . We assume that the reach of X is greater than ρ > 0. Let h > 0 be the sampling gridstep of the regular grid hZ 3 . The Gauss digitization of V is defined as G h (V ) := V ∩ hZ 3 . Digital sets are subsets of hZ 3 . Let us denote Q h z the axes-aligned cube of edge length h centered on a point z ∈ hZ 3 . Digital sets can thus be seen as a union of such cubes in R 3 . We can now defined the digitized surface X h of X := ∂ V at step h as the topological boundary of the Gauss digitization of V , seen as a union of cubes: We call surfel any square of edge length h that is the face of some Q h z , with z ∈ G h (V ) and that is included in X h . We call linel any edge of a surfel.
In order to apply the stability result presented in the previous section, a few requirements are necessary: (i) the object S := X h should be a surface (i.e., a 2-manifold), and (ii) there exists an open set U of the form U = π −1 X (V ) such that U ∩ S = ∅ and π X : U ∩ S → X is injective.
Concerning point (i), unfortunately, X h is not a 2-manifold in general, even for smooth and convex shapes V [32]. However X h is almost a 2-manifold since the places where it is not a 2-manifold tends to zero quickly as h tends to zero. We recall [19,Thm. 2]: Theorem 5.1 Let h < 0.198ρ, letting y ∈ X h , then the digital surface X h is locally homeomorphic to a 2-disk around y if either (i) y does not belong to a linel of X h , or (ii) y belongs to a linel s of X h and s ∩ X = ∅, or (iii) y belongs to a linel s of X h and there exists p ∈ s ∩ X but the angle α y between s and the normal to X at p satisfies α y ≥ 1.260h/ρ. So X h may not be a manifold only when the normal of X is very close to one of the axes. It is possible to fix locally the manifoldness of digital surfaces by making it well-composed, either by subsampling the grid and doing majority interpolation [32], or repairing the digital surface by adding voxels [30] or by splitting vertices and edges [16]). All these transformations affect a very small part of the digital surface according to the previous theorem. Therefore from now on we shall assume that X h is a 2-manifold.
As for point (ii), [19,Thm. 3] gives a positive answer to the existence of an open set U where the projection is injective. Indeed the area on X , where the projection π X : U ∩ X h → X is not injective is proportional to Area(S)h and tends to zero.
The following results relate combinatorial properties of X h to geometric properties. In all the sections, we denote by B R p the ball centered at p and of radius R.

Lemma 5.2 Let p be an arbitrary point of X h . Let B be the surfels of X h that lie in
Then the numbers N s of surfels and N b of boundary linels of B (the ones bordering only one surfel of B) follow these bounds for radius R = K h α : The proof of Lemma 5.2 relies on the following intermediary lemma.

Proof of Lemma 5.3
Let x ∈ X h ∩∂B R p . We denote by ε the Hausdorff distance between X and its Gauss digitization X h . It is known that ε = O(h). We put p = π( p) and x = π(x). We denote by C the geodesic starting at the point p and passing through x . We first assume that x ∈ B R p . In that case, the geodesic curve C is extended until it reaches a pointx ∈ ∂B R p . Note that such a curve always exist if R is small enough. In the following of the proof, we denote by C a,b the shortest path on X between any two points a and b and by (C a,b ) its length. Since the length of a curve is greater than its chord, we have If R is small enough, then the line segment [ p ,x] belongs to the offset X r of X , where r = reach(X ) is the reach of X . By using the Lipschitz property of the projection map π onto X (see Proposition 4.1), one gets Since the curve is longer than its chord, If x / ∈ B R p , the same result holds with a similar proof. In that case the geodesic curve C does not have to be extended andx = C ∩ B R p . We deduce that x −x = O(h 3α + h), which ends the proof.

Proof of Lemma 5.2 Let us first bound
. It is known that the digitized surface X h is at a Hausdorff distance less than h √ 3/2 from X [19, Thm. 1]. Hence the set B is included in a set of dimensions (2R + 2h) × (2R + 2h) × (O(R 2 ) + 2h). Since α ∈ (0, 1), we have h in O(h α ). Furthermore, since R = K h α , the set B is thus included in a domain of volume O(h α h α h min(2α,1) ). This implies that the number of voxels intersecting B is less than O(h min(4α,2α+1)−3 ). This also implies that N s = O(h min(4α,2α+1)−3 ).
Let us now bound N b . Every point x ∈ ∂ B belongs to a surfel that intersects X h ∩ ∂B R p , so is at a distance less than h √ 2/2 from X ∩ ∂B R p . By Lemma 5.3, one has . The number of boundary edges N b is of the same order than the number of cubes of size h intersecting this volume so We may now state our convergence result for the normalized corrected curvatures onto digitized surfaces. > 0 and α ∈ (0, 1)). Then
Proof This proof relies on Corollary 4.5. We know that the Hausdorff ε between X and the boundary of its Gauss discretization X h is no greater than h Furthermore, since we are in the hypotheses of Lemma 5.2, we have bounds for Using the fact that there exists a constant C such that μ X 0 (π X (B)) ≥ Ch 2α , one gets Similarly, for k = 1, 2, Finally, we can relate our normalized mean curvature with the ratio of the mean curvature measure and the area measure: It remains to relate this ratio of measures to the mean curvature at p := π X ( p). Recall that the surface is of class C 2,1 . From the proof of the preceding lemma, we know that there are two real numbers R 1 and R 2 with Then we can write in geodesic polar coordinates Since H is K -Lipschitz on D R 1 ( p ) for some K > 0 one has Similarly, the difference between the integrals over the two geodesic disks satisfies Finally, Combining (2) and (3) The same holds for the Gaussian curvature.

Remark 5.6
The above bound η = O(h 2/3 ) has been established in [22] for the normal vector estimator based on digital integral invariants. Note that β = 1 is the optimal convergence rate, since it is the one obtained by taking the ground truth normal, i.e., taking the normal at the projection on X . There is yet no formal proof that any digital normal vector estimator is Lipschitz (which implies β = 1). However, we have run simulations and both Voronoi Covariance measure [11,24] and digital Integral Invariant [22] appear to be Lipschitz normal estimator (see Fig. 4).

Remark 5.7
If we choose u as the normals estimated by digital integral invariant normal estimator and we assume β = 1, then the previous theorem implies that, when the set B of surfels is taken in a ball of center p and radius K h 1/3 , then the mean corrected curvatureĤ u (B) (resp. the Gaussian corrected curvatureĜ u (B)) tends to the mean curvature H (π X ( p)) (resp. to the Gaussian curvature G(π X ( p))) with a speed O(h 1/3 ). In the terminology of [18],Ĥ u (resp.Ĝ u ) is a multigrid convergent mean (resp. Gaussian) curvature estimator.
In the following section, we perform a comparative evaluation of our curvature estimators. Although we have checked that they are indeed convergent for a radius K h 1/3 with a speed at least O(h 1/3 ), we will run our experiments with a much smaller radius K h 1/2 (see Figs. 5 and 6). Indeed it provides not only faster estimators but also a slightly better error bound in practice, closer to O(h 2/3 ). The reason is still investigated.

Experimental Evaluation on Digital Surfaces
We present here a series of experiments which demonstrates that the corrected measures provide accurate and stable curvature information. We do not evaluate experimentally the accuracy of measures themselves (which are already established through our main theorem), but the much more difficult problem of pointwise mean and curvature inference. As a stressful testbed that maximizes the difficulty of estimating curvatures, we evaluate the accuracy of the normalized corrected mean curvatureĤ and the normalized Gaussian curvatureĜ on digital surfaces (see Definition 3.4). Indeed, their canonical normal vectors can only take six different values. Last, we compare our new approach with the state-of-the-art method of [7,22], which is based on digital integral invariants. We also show that the normal cycle approach [8,9,34,36] is neither accurate nor convergent for digital surfaces.

Methodology of Evaluation on Digitizations of Polynomial Surfaces
We evaluate the accuracy of our geometric estimators on the digitization of implicitly defined polynomial shapes X , in order to have ground-truth curvatures. Let us detail our methodology for evaluating our new curvature estimators.

Input Digitized Surfaces
We digitize a shape X using Gauss digitization G h (X ) := X ∩ hZ 3 at several gridsteps h. If we see the discrete set G h (X ) as a union of axis-aligned cubes of edge-length h centered on those points, its topological boundary is then a union of axis-aligned squares of edge-length h that forms a digital surface (e.g. see [19]). We denote it by ∂ h X . As illustrated in Fig. 7, digitized surfaces ∂ h X tends toward the smooth surface ∂ X in Hausdorff distance. In fact they are a Hausdorff approximation of ∂ X at distance less than h √ 3/2 [19]. However their natural normal vectors n do not tend toward the normal vectors of the smooth surface, since they can take only six different values whatever is h.

Input Corrected Normal Vector Field u
We must estimate the field u solely from the input digitized surface ∂ h X . We will use several normal vector estimators in the experiments, in order to show the importance of having a convergent estimator but also to show that our method gives stable results for any convergent u: -Trivial normal estimator (TN): this estimator just replicates n (i.e., u = n). We use it in experiments since our measures become then equivalent to the normal cycle [8,9,34,36]. -Digital Integral Invariant normal estimator (II): this estimator provides a convergent normal vector field u for a certain family of parameters [22]. It depends on a radius of integration parameter r := kh α . As shown by our experiments, the values α = 0.5 and r = 3 provide both good and stable results. -Voronoi Covariance Measure normal estimator (VCM): this estimator provides a convergent normal vector field u for a certain family of parameters [10,11]. It depends on a radius of integration parameter r := kh α , which we set exactly as parameter r of II, and on a distance of computation R := K h α , where K = 10 gives good results.

Estimated Curvatures
We will estimate the accuracy and stability of the following curvatures estimators: -Normalized corrected mean curvature (Ĥ u ) and Gaussian curvature (Ĝ u ): since they are ratios of measures, we must choose a Borel set on which measures are computed. For any point p ∈ ∂ h X , we simply computeĤ u ( p) : where ρ := mh β . As shown by our experiments, the values β = 1/2 and m = 3 provide both good and stable results. We use the corrected normal field u given either by II or VCM normal estimators (in both cases, r = 3h 1/2 is used). -Digital integral invariant mean and Gaussian curvature estimators (Ĥ II ) and (Ĝ II ): both are parameterized by the radius r of integration. Experiments show that for r = k h α , we must set α = 1/3 (and no greater value) to get convergence and set k = 6 to minimize estimation errors.
-Normal Cycle mean and Gaussian curvatures (Ĥ NC ) and (Ĝ NC ): they are defined similarly withĤ u andĜ u except that we use the Trivial Normal estimator to compute u, i.e.,Ĥ NC :=Ĥ n andĜ NC :=Ĝ n .

Measuring 2 -and ∞ -Errors
Curvatures are estimated at the centroid p of each surfel element (the squares that form the digitized shape boundary), and are compared to the curvatures of the point q ∈ ∂ X that is closest to p. Note that q = π( p) since p is in the reach of ∂ X . For instance, letting σ be the set of centroids of the surfel of ∂ h X , we define the errors betweenĤ and H as: Note that when ∞ -error tends to zero as h tends to zero implies the classical multigrid convergence. It also implies that 2 -error tends to zero.

Robustness to Noise
In practical application, input data are rarely perfect digitizations and may be corrupted by noise. We have used a noise model parameterized by a probability p, which perturbates the input voxels according to their distance to the exact digitized boundary. More precisely, if δ is the discrete distance of the voxel to its nearest boundary point (minimum distance is 1), then this voxel has a probability of p δ to be flipped inside/out.

Asymptotic Behavior ofĤ u and Normal CycleĤ NC
First we measure both the ∞ -and 2 -errors of our proposed mean curvature estimator H u . We further test several normal estimators for u: II, VCM, and TN. Results are displayed in Fig. 8. Graphs show an experimental convergence speed of O(h 2/3 ), II and VCM normals. It shows also that the normal cycle method is not convergent, sincê H u with TN corresponds to the normal cycle mean curvatureĤ NC .

Comparative Evaluation ofĤ u andĤ II
We now compare the ∞ -and 2 -errors ofĤ u with the ones ofĤ II , which is the state-of-the-art method for digital surface curvature estimation. We use the II normal estimator forĤ u with r = 3h 1/2 . Results are shown in Fig. 8. First it confirms that H II is not convergent if its parameter r follows some (h 1/2 ). Secondly we do obtain (Right) Respective asymptotic errors of mean curvature estimations byĤ u andĤ II along the "Goursat" shape. In both figures, ∞ -errors are drawn with thick lines, 2 -errors are drawn with thin lines a convergence speed of (h 1/3 ) for bigger radii r = (h 1/3 ). Last our new estimator H u has a much faster convergence, approximately (h 2/3 ), despite the fact that both the integration radius r for u and the integration radius ρ forĤ u are much smaller, i.e., 3h 1/2 .
We further illustrate the differences of the two estimatorsĤ u andĤ II by displaying the estimated mean curvatures and the localization of errors on several digitizations of "Goursat" on Fig. 9. It is clear that errors are mostly localized on places of extremal curvatures, but our estimator converges much faster everywhere visibly and does not oscillate around the correct value.
Last, we measure the stability of both estimators with respect to their parameters. For H u , we measure the ∞ -error when changing k in the radius r := kh 1/2 and changing m in the radius ρ := mh 1/2 . ForĤ II we simply change k in the radius r := kh 1/3 . See Fig. 10. First, the results show that the exponents chosen for the gridstep are consistent (best results are achieved for the same constant at a finer scale). It confirms that integral invariant methods require a radius of integration that is much larger asymptotically. Secondly it shows that our method is more stable with respect to parameter settings. We have almost the same errors in the range (k, m) ∈ [2.5, 4] × [2.5, 6].

Evaluation ofĤ u on Digital Shapes
We also run our method on digitizations of classical shapes ("dragon" and "octaflower"), trying several parameters (here we checked several initial gridsteps) such that Integral Invariant and Normal Cycle methods give the best possible results. Outputs are displayed in Figs. 11 and 12. It is clear that bothĤ NC andĤ II oscillate around the correct solution (see black random or moiré patterns on both figures, whichever the resolution). On the contrary,Ĥ u is stable in zero-curvature regions (like on the octaflower) while accurately delineating the small scales of the Chinese dragon.  Fig. 10 Stability ofĤ u andĤ II with respect to parameters and gridstep. EstimatorĤ u is parameterized by the integration radius r := kh 1/2 of its normal estimator u = II and by the measure radius ρ := mh 1/2 . We plot 2 -and ∞ -errors ofĤ u − H as a function of k and m. EstimatorĤ II is parameterized by the integration radius r := k h 1/3 . We plot 2 -and ∞ -errors ofĤ II − H as a function of k (displayed as 2d plot to make easier the comparison)

Robustness to Noise ofĤ u
Last we have checked the robustness to noise of our mean curvature estimatorĤ u . Experiments show that our method is mostly sensitive to the quality of the input corrected normal vector field u. On the one hand VCM normal estimator is relatively accurate but more unstable than II normal estimator. On the other hand, II normal estimator is inaccurate at places with sharp features (as one can see on the outer parts of the screw of "octaflower"). We should certainly in this case use smarter normal estimators, like the AT normal vector method of [6], which is able to compute piecewise smooth normal vector field, or a normal estimator which takes into account the digital nature of the input shape, like the plane-probe algorithms of [20,21].

Comparative Evaluation ofĜ u andĜ II
We compare the ∞ -and 2 -errors ofĜ u with the ones ofĜ II , which is also the state-of-the-art method for digital surface curvature estimation. We pick the same parameterization as in the previous paragraph for all estimators. We observe the same behaviour ofĜ u with respect toĜ II : much faster convergence speed with small integration radii. We further illustrate the differences of the two estimatorsĜ u andĜ II by displaying the estimated Gauss curvatures and the localization of errors on several digitizations of "Goursat" in Fig. 15 curvatures, but our estimator converges much faster everywhere visibly and does not oscillate around the correct value.

Evaluation ofĜ u on Digital Shapes
We also check Gaussian curvature estimatorsĜ u ,Ĝ NC , andĜ II on the same classical shapes as in the previous subsection. Results are displayed in Figs. 16 and 17; keeping the same parameters as above. We also observe the oscillations and moiré patterns inĤ II estimations. Furthermore the Normal Cycle estimatorĜ NC is clearly incorrect and gives only extremal results. This is becauseĜ NC takes into account solely six possible normals.

Robustness to Noise ofĜ u
Last we have checked the robustness to noise of our Gaussian curvature estimatorĜ u . Experiments show that our method is mostly sensitive to the quality of the input corrected normal vector field u. On the one hand VCM normal estimator is relatively accurate but more unstable than II normal estimator. On the other hand, II normal estimator may be inaccurate at places with sharp features (as one can see on the outer parts of the screw of "octaflower"). As already said in the previous section, we should certainly in this case use smarter normal estimators [6,20,21]. Error |G −Ĝ u | II methodĜ II [7] Error |G −Ĝ II | Fig. 15 Illustration of the accuracy of Gaussian curvature estimatorsĜ u andĜ II at different resolutions: 1st row is ground-truth G, 2nd row is our approachĜ u , 3rd row is its local estimation error |Ĝ u − G|, 4th row is integral invariant approachĜ II , 5th row is its local estimation error |Ĝ II − G|. Curvatures are displayed with colors from dark blue (−0.1) to white to red (0.

Practical Computation Times of Our Curvature Estimators
We note first that the computation times of the measures μ 0 , μ 1 , μ 2 , μ , and μ X,Y per cell are linear with the number of cells and are extremely fast. Their running times are negligible with respect to the estimation of corrected normal vector field or their integration in a ball a radius ρ (2% when N ≈ 1e 3 , 0.01% when N ≈ 4e 6 ). Secondly, computation times for Gaussian curvature estimatorĜ u or principal directions is almost the same as the one ofĤ u andĜ II , and are thus not displayed.
We have plotted the measured computation times of our mean curvature estimator H u , first as a function of the number of surfels of the digitized boundary, and after as a function of the accuracy (see Fig. 19). Running times for all estimators were measured with a mono-threaded CPU implementation on an average server (Intel Xeon Gold 6128 processor, 3.4 GHz, cache memory 19.25 Mb, each thread is evaluated at 6785.92 bogomips). In both cases, our approach clearly outperforms the digital Integral Invariant approach of [7,22]. We have not compared with the execution times of Normal Cycle method, sinceĤ NC is not accurate. Note that using VCM normal estimator instead of II normal estimator to compute u is faster for large digital shapes, but requires more memory.
To sum up, our method has an experimental complexity in (n 3/2 ) compared to Integral Invariant method that has an experimental complexity in (n 2 ). Our code could be further optimized in the integration step of the measures μ k by reusing the results of a neighboring point, as it is already done in Integral Invariant method. Last, for a given accuracy, our method is even much faster (5000 times faster than II to achieve 0.0065 2 -accuracy).

Conclusion
We have proposed a sound mathematical framework for defining area and curvature measures over rather general surfaces. We have also shown that these measures are stable with an error proportional to the sum of the position error and the normal error. This framework induce sound definitions of curvatures on polyhedral surfaces, which are easy and fast to compute. We have evaluated extensively the numerical accuracy of our approach on digital surfaces, which are polyhedral surfaces with bad naive normals. It shows that our method is also effective and accurate in practice. We are currently adapting the anisotropic curvature measures introduced in [8,9] (see also textbook [26]) to our framework. These new measures give estimates of principal curvatures and principal directions, which are also stable and convergent by the same principles stated in Sect. 4. These works will be the focus of a forthcoming paper. Preliminary results are displayed in Fig. 20. We are also currently investigating smooth corrected normal vector field obtained by interpolation, which may reveal to be even more accurate than piecewise constant corrected normal vector field, both theoretically and practically.