# Corrected Curvature Measures

Jacques-Olivier Lachaud, Pascal Romon, Boris Thibert

► To cite this version:

Jacques-Olivier Lachaud, Pascal Romon, Boris Thibert. Corrected Curvature Measures. Discrete and Computational Geometry, 2022, 68, pp.477-524. 10.1007/s00454-022-00399-4. hal-02193774

HAL Id: hal-02193774

https://hal.science/hal-02193774

Submitted on 28 Jul 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Corrected curvature measures*

J.-O. Lachaud†

P. Romon‡

B. Thibert§

July 25, 2019

## Abstract

This paper proposes a new mathematical and computational tool for infering 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 $R^3 \times S^2$. 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.

# 1 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 non-smooth 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 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.

<footnote>*This work has been partly funded by COMEDIC ANR-15-CE40-0006 research grant.</footnote>

<footnote>†Laboratoire de Mathématiques, Univ. Grenoble Alpes, Univ. Savoie Mont Blanc</footnote>

<footnote>‡Laboratoire d'Analyse et de Mathématiques Appliquées, Univ. Paris-Est</footnote>

<footnote>§Laboratoire Jean Kuntzman, Univ. Grenoble Alpes</footnote>
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 H. 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 [14]. The normal cycle of a general shape in $\mathbb{R}^n$ is a current in the Grassmann bundle $\mathbb{R}^n \times \text{Gr}(1, n)$ (where the $\text{Gr}(1, n) \simeq \mathbb{S}^{n-1}$ factor contains 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 [15, 9]. 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 extend 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 Section 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 $\mathbf{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 $\mathbf{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, \mathbf{u})$ where $S$ is a piecewise $C^1$-surface of $\mathbb{R}^3$ and $\mathbf{u}$ is a piecewise $C^1$ unit vector field and we call it the *corrected normal current* $N(S, \mathbf{u})$. Our corrected normal current $N(S, \mathbf{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, \mathbf{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 $\mathbf{u}$ is close to the normal vector field of $X$.

- We apply our results to the digital case and show 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 state-of-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 $\mathbb{R}^n \times \text{Gr}(d, n)$, where $\text{Gr}(d, n)$ denotes the Grassmannian manifold of unoriented $d$-planes in $\mathbb{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 [29, 28]. 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, 19], 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 Section 2, we formalize this new integral with currents in $\mathbb{R}^3 \times \mathbb{S}^2$, which we view as the set of couples (position, unit normal vector). The unit sphere $\mathbb{S}^2$ stands for the oriented Grassmannian $\text{Gr}(1, 3)$, the set of unit vectors in $\mathbb{R}^3$. Note that $\text{Gr}(1, 3)$ can be identified by duality to $\text{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, \mathbf{u})$ where $S$ is a surface and $\mathbf{u}$ a unit vector field along $S$ we associate a current $N(S, \mathbf{u})$ in $\mathbb{R}^3 \times \mathbb{S}^2$ (Definition 4). Intuitively we may view the current as a piecewise smooth oriented 2-surface in $\mathbb{R}^3 \times \mathbb{S}^2 \subset \mathbb{R}^6$. In particular, when the vector field $\mathbf{u}$ is the naive normal field of $S$ and $S = \partial V$ is the boundary of a domain of $\mathbb{R}^3$, then $N(S, \mathbf{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, \mathbf{u})$ gives us the flexibility we need.)

In Section 3, we define our *geometric measures*, or *corrected Lipschitz-Killing curvature forms* as the integral on $N(S, \mathbf{u})$ of the classical invariant forms on $\mathbb{R}^3 \times \mathbb{S}^2$ (Definition 5). 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). 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 $\mathbf{u}$.

We then provide in Section 4 a stability result for our corrected curvature measures (Theorem 1). More precisely, we show that the corrected curvatures of some $(S, \mathbf{u})$ approximate well the curvature measures of a surface $X = \partial V$ of
class $C^2$ provided that the surface $S$ is close to $X$ in the Hausdorff sense and that the vector field $\mathbf{u}$ approximates well the geometric normal $\mathbf{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, \mathbf{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, \mathbf{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, \mathbf{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 $\mathbf{u}$ is estimated by a multigrid convergent normal estimator (Theorem 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 Section 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 Section 7 and outline several research directions, especially the extension of our work to anisotropic curvature measures.

# 2 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*.

## 2.1 Corrected surface

We assume in the following that $S$ is a closed piecewise $C^1$ oriented surface of $\mathbb{R}^3$, namely $S = \cup_{i=1}^n 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 \neq j$, $S_{i,j} = S_i \cap 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 ori-
entation 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 $\mathbf{e}_{ij}$ along $S_{ij}$, which is tangent to $S_{i,j}$. If we are given a normal vector field $\mathbf{n}_i$ on $S_i$ respecting its orientation, then $\mathbf{n}_i \times \mathbf{e}_{ij}$ points toward the inside of $S_i$. Note that the adjacent face $S_j$ induces the opposite orientation on $S_{i,j}: \mathbf{e}_{ji} = -\mathbf{e}_{ij}$.

We say that $\mathbf{u}: S \setminus \{S_{i,j}, i \neq j\} \to \mathbb{S}^2$ is a corrected normal vector field on $S$ if: (i) its restriction $\mathbf{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 \cap S_j$, $\mathbf{u}_i$ and $\mathbf{u}_j$ are not antipodal; (iii) at any vertex $p = S_{i_1} \cap \dots \cap S_{i_k}$, the corrected normals $\mathbf{u}_{i_1}, \dots, \mathbf{u}_{i_k}$ lie in the same hemisphere. Note that, in this definition, we do not require $\mathbf{u}$ to be continuous or even defined on the edges.

**Definition 1 (Corrected surface).** We say that the couple $(S, \mathbf{u})$ is a corrected surface if $S \subset \mathbb{R}^3$ is a surface that satisfies the above assumptions and $\mathbf{u}$ is a corrected normal vector field on $S$.

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

*Example 2.* When $S$ is an oriented polyhedral surface, we may choose $\mathbf{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 $-\pi$ and $+\pi$ (resp. is a graph over the plane bounding the hemisphere).

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

## 2.2 Combinatorial structure

Starting from 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.<fnref n="1" />

The idea behind the construction is that each inner edge is blown into a strip (i.e. a combinatorial quadrilateral), while each inner vertex of degree $d$ is blown into a face with $d$ edges, see Figure 1. More precisely,

* the vertices of $V^*$ are the flags $\{(p, f) \in V \times F, p \in f\}$,

* the edges of $E^*$ are the flags $\{(e, f) \in E \times F, e \in f\}$,

<footnote n="1">¹See http://mathoverflow.net/questions/263452</footnote>
• the set $F^*$ is the set $\mathring{V} \cup \mathring{E} \cup F$, where $\mathring{V}$ (resp. $\mathring{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^* \in 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 \in V$, we denote it $S_p^*$ and let $S_1, \dots, 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^*, \dots, p_n^*$ where $p_i^* = (p, S_i)$.

• Whenever the face corresponds to an inner edge $S_{i,j} = S_i \cap S_j \in E$, and $S_{i,j}$ joins the vertices $p, q$, with the convention that $S_i$ induces the orientation $p \to q$ on $S_{i,j}$, and $S_j$ the opposite orientation. Then $S_{i,j}^*$ is the quadrilateral face joining $(p_i^*, p_j^*, q_j^*, q_i^*)$ where $p_k^* = (p, S_k)$ and $q_k^* = (q, S_k)$.

• Whenever the face $S_i^*$ corresponds to $S_i \in F$, and the vertices of $S_i$ (in order) are $p_1, \dots, p_n$, then $S_i^*$ is given by its vertices $p_1^*, \dots, 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 $\ell = d(p)$ edges that are common with faces $S_{i,j}^*$ and denoted by $s_1, \dots, s_\ell$; (i) we add a vertex $p^*$ at the interior of the face $S_p^*$; (ii) for every $1 \le i \le \ell$, we add the triangle with vertex $p^*$ and opposite edge $s_i$; (iii) we add the edges corresponding to these triangles, see Figure 1.

## 2.3 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 $\mathbb{R}^3 \times \mathbb{S}^2$. Its construction uses the corrected normal $\mathbf{u}$ and is done so that the corrected normal cone inherits the orientation of $S^*$.

**Definition 2.** Let $(S, \mathbf{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^*, \mathbf{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^*, \mathbf{u}) = NC(S_i, \mathbf{u})$ is the image of $S_i$ by the map $p \mapsto (p, \mathbf{u}(p))$;

• if $f^* = S_{i,j}^*$ corresponds to an edge $S_{i,j}$ of $S$, then

$$NC(S_{i,j}, \mathbf{u}) = \{(p, n), \text{ for } p \in S_{i,j} \text{ and } n \in \text{Arc}(\mathbf{u}_i(p), \mathbf{u}_j(p))\},$$

where $\text{Arc}(\mathbf{u}_i(p), \mathbf{u}_j(p))$ denotes the unique geodesic arc (on the sphere) between $\mathbf{u}_i(p)$ and $\mathbf{u}_j(p)$.
<figure><img src="image_2.png" /><figcaption>Figure 1: Expansion/cantellation at a degree three vertex, followed by subdivision of the face $S_p^*$ corresponding to an original vertex $p$. The vertex $p$ gives rise to a new face $S_p^*$ (further subdivided); edges $S_{i,j} = S_i \cap S_j$ give rise to new quadrilateral faces $S_{i,j}^*$; faces $S_i$ yield new faces $S_i^*$.</figcaption></figure>

• if $f^*$ is a triangle of $S_p^*$, then $\mathrm{NC}(f^*, \mathbf{u})$ denotes the spherical triangle with vertices $\mathbf{u}(p)$, $\mathbf{u}_i(p)$ and $\mathbf{u}_{i+1}(p)$ of area strictly less than $2\pi$, where $\mathbf{u}(p) \in \mathbb{S}^2$ is the normalized average of the $\mathbf{u}_i(p)$.

Remark that for each oriented face $f^*$ of $S^*$, the corrected normal cone $\mathrm{NC}(f^*, \mathbf{u})$ inherits from the orientation of $f^*$.

**Definition 3** (Corrected Normal Cone). The corrected normal cone $\mathrm{NC}(S, \mathbf{u})$ is the polygonal chain built as the sum of the $\mathrm{NC}(f^*, \mathbf{u})$, where $f^*$ ranges over the faces of $S^*$.

Remark that for every face $f^*$, $\mathrm{NC}(f^*, \mathbf{u})$ is a surface of class at least $C^1$ in $\mathbb{R}^3 \times \mathbb{S}^2$. The corrected normal cone $\mathrm{NC}(S, \mathbf{u})$ can be seen globally as the image by a continuous map of the combinatorial surface $S^*$ into $\mathbb{R}^3 \times \mathbb{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 1.** The corrected normal cone $NC(S_p^*, \mathbf{u})$ above a vertex $p$ of $S$ is made of $d(p)$ spherical triangles, each of them having $\mathbf{u}(p)$ as a vertex. Note that the algebraic sum of these oriented triangles do not depend on this arbitrary point $\mathbf{u}(p)$. The corrected normal cone $NC(S_p^*, \mathbf{u})$ can be seen as a set of point $n \in \{p\} \times \mathbb{S}^2$ with an integer multiplicity $\mu(p)$.

## 2.4 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 Section 4.3.

**Definition 4 (Corrected Normal Current).** Let $(S, \mathbf{u})$ be a corrected surface. The corrected normal current $N(S, \mathbf{u})$ of $(S, \mathbf{u})$ associates to every differential 2-form $\omega$ of $\mathbb{R}^3 \times \mathbb{S}^2$ the real number

$$N(S, \mathbf{u})(\omega) = \sum_{f^*} \int_{NC(f^*, \mathbf{u})} \omega$$

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.

**Proposition 1.** If $S$ has no boundary, then the corrected normal current $N(S, \mathbf{u})$ has no boundary.

*Proof.* By construction $NC(S, \mathbf{u})$ is the continuous image of a combinatorial surface $S^*$. This comes from the coincidence of incident faces along their boundaries. Now, if $S$ has no boundary, neither does $S^*$, nor $N(S, \mathbf{u})$. $\square$

**Remark 2.** Whenever $S = \partial V$ is a piecewise smooth or planar surface, which is the boundary of a domain in $\mathbb{R}^3$, and $\mathbf{u}$ is chosen as the unit normal on faces pointing toward the outside, then $N(S, \mathbf{u})$ coincides with the normal cycle $N(V)$ of $S$ [14, 36]. The corrected normal current thus generalizes the normal cycle.

**Remark 3.** For the sake of consistency and subsequent proofs, we have assumed that the vector field $\mathbf{u}$ has unit length. Note that it is sometimes useful to relax that requirement to nonzero vector fields, for example when using linear interpolation to 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 $\mathbf{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 $\mathbf{u}$ tends to the smooth normal $\mathbf{n}$, even if $\mathbf{u}$ is not of unit length. This case will be illustrated in a forthcoming article.
<figure><img src="image_3.png" /><figcaption>Figure 2: The corrected normal current of a planar curve, seen as a discrete curve in $\mathbb{R}^3$. The Grassmann bundle $\mathbb{R}^2 \times \text{Gr}(1, \mathbb{R}^2) = \mathbb{R}^2 \times \mathbb{S}^1$ is represented as $\mathbb{R}^3$ where the third coordinates is the angle $\theta$. 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\pi/2$ for some integer $k$); therefore the extra circular arcs at the vertices are vertical edges of length $\pm\pi/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 $\mathbb{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 into an edge.</figcaption></figure>

# 3 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 [14]). These invariant forms [34, 36] form a basis of the 2-forms in $\mathbb{R}^3 \times \mathbb{S}^2$ that are invariant by the action of the rigid motions of $\mathbb{R}^3$. They yield the area, the mean curvature and the Gaussian curvature.

## 3.1 Invariant forms, curvature measures, curvature estimators

Let us now recall the expression of the invariant differential 2-forms of $\mathbb{R}^3 \times \mathbb{S}^2$ (see [26] for more details). Let $(p, w)$ be any point in the *oriented* Grassmann bundle $\mathbb{R}^3 \times \mathbb{S}^2$. Clearly the tangent plane of $\mathbb{R}^3 \times \mathbb{S}^2$ at the point $(p, w)$ is a
space of dimension 5 spanned by the following vectors

$$ \varepsilon_1 = \begin{pmatrix} \mathbf{e}_1 \\ 0 \end{pmatrix}, \quad \varepsilon_2 = \begin{pmatrix} \mathbf{e}_2 \\ 0 \end{pmatrix}, \quad \varepsilon_3 = \begin{pmatrix} \mathbf{w} \\ 0 \end{pmatrix}, \quad \tilde{\varepsilon}_1 = \begin{pmatrix} 0 \\ \mathbf{e}_1 \end{pmatrix}, \quad \tilde{\varepsilon}_2 = \begin{pmatrix} 0 \\ \mathbf{e}_2 \end{pmatrix}, $$

where $\mathbf{e}_1$ and $\mathbf{e}_2$ are any vectors such that $(\mathbf{e}_1, \mathbf{e}_2, \mathbf{w})$ form a direct orthonormal frame of $\mathbb{R}^3$. The set of differential 2-forms of $\mathbb{R}^3 \times \mathbb{S}^2$ that are invariant under the action of rigid motions is a vector space of dimension four spanned by the following forms

$$ \omega_0 = \varepsilon_1^b \wedge \varepsilon_2^b, \quad \omega_1 = \varepsilon_1^b \wedge \tilde{\varepsilon}_2^b - \varepsilon_2^b \wedge \tilde{\varepsilon}_1^b, \quad \omega_2 = \tilde{\varepsilon}_1^b \wedge \tilde{\varepsilon}_2^b, \quad \omega_\Omega = \varepsilon_1^b \wedge \tilde{\varepsilon}_1^b + \varepsilon_2^b \wedge \tilde{\varepsilon}_2^b $$

where $\varepsilon^b(\xi)$ denotes the scalar product of $\varepsilon$ and $\xi$.

**Definition 5.** Let $k \in \{0, 1, 2, \Omega\}$. The $k^{th}$ Lipschitz-Killing corrected curvature measure of $S$ with respect to $\mathbf{u}$ associates to every Borel set $B$ the real number

$$ \mu_k(B) := \int_{NC(S, \mathbf{u}) \cap \pi^{-1}(B)} \omega_k, $$

where $\pi : \mathbb{R}^3 \times \mathbb{S}^2 \to \mathbb{R}^3$ denotes the projection on the position space.

**Remark 4.** When $S$ is smooth and $\mathbf{u} = \mathbf{n}$ the standard normal, $\mu_0$ is the area density dA, while $\mu_1$ and $\mu_2$ are respectively the mean curvature and Gaussian curvature densities, $-2H d\mathcal{H}^2$ and $K d\mathcal{H}^2$ (up to a scalar factor). This result is reproven below.

**Remark 5.** The measure $\mu_\Omega$ induced by the symplectic form $\omega_\Omega$ plays a different role. Instead of measuring metric data, it tests whether a surface in the $\mathbb{R}^3 \times \mathbb{S}^2$ is indeed the normal bundle of a surface in $\mathbb{R}^3$. Quite logically, it vanishes identically on smooth surfaces as soon as $\mathbf{u} = \mathbf{n}$. For an approximated surface, it measures how far the corrected normal $\mathbf{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:

**Definition 6.** For any Borel set $B$ such that $\mu_0(B) \neq 0$, the normalized corrected mean curvature $\hat{\mathbf{H}}^\mathbf{u}$ of $S$ with respect to $\mathbf{u}$ is defined as $\hat{\mathbf{H}}^\mathbf{u}(B) := \mu_1(B)/\mu_0(B)$, and the normalized corrected Gaussian curvature $\hat{\mathbf{G}}^\mathbf{u}$ of $S$ with respect to $\mathbf{u}$ is defined as $\hat{\mathbf{G}}^\mathbf{u}(B) := \mu_2(B)/\mu_0(B)$.

## 3.2 Calculation of the corrected curvature measures

Let $(S, \mathbf{u})$ be a corrected surface. We recall that $S = \cup_{i=1}^n 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 $\mathcal{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 $\mathbf{n}(p)$ the unit oriented normal vector, by $T_pS$ the plane tangent to $S$, by $\mathbf{e}_1(p)$ a vector of $T_pS$ given by $\mathbf{e}_1 = (\mathbf{n} \times \mathbf{u})/\|\mathbf{n} \times \mathbf{u}\|$ if $\mathbf{n}$ and $\mathbf{u}$ are not colinear, and given by one of the principal directions otherwise. We also introduce $\mathbf{e}_2 = \mathbf{u} \times \mathbf{e}_1$ and $\mathbf{e}'_2 = \mathbf{n} \times \mathbf{e}_1 \in T_pS$. Note that $(\mathbf{e}_1, \mathbf{e}_2, \mathbf{u})$ is a moving frame of $\mathbb{R}^3$ associated with $\mathbf{u}$, while $(\mathbf{e}_1, \mathbf{e}'_2, \mathbf{n})$ is a Darboux frame associated to the surface $S_i$.

At every point $p$ on an edge $S_{i,j}$, we denote by $\mathbf{e}(p) = \mathbf{e}_{ij}(p)$ the unit vector along the edge $S_{i,j}$ oriented<fnref n="2" /> as the boundary of $S_i$, by $\Psi(p) = \angle(\mathbf{u}_i(p), \mathbf{u}_j(p))$ the corrected dihedral angle between the faces $S_i$ and $S_j$, and $\mathbf{e}_1(p) = (\mathbf{u}_i \times \mathbf{u}_j)/\|\mathbf{u}_i \times \mathbf{u}_j\|$. (The last vector is only defined when $\mathbf{u}_i(p)$ and $\mathbf{u}_j(p)$ are not collinear; however, whenever that happens, $\Psi(p)$ vanishes and $NC(p, \mathbf{u})$ drops in dimension, hence does not contribute to the integrals. Without loss of generality, we will assume that $\mathbf{e}_1$ is always well-defined.)

In order to integrate differential forms over manifolds (more details are provided in Section 4.3), we will use the notion of Hausdorff measure. In the following, we denote by $\mathcal{H}^m$ the $m$-dimensional Hausdorff measure.

**Proposition 2.** *The corrected Lipschitz-Killing curvature measures of $S$ with respect to $\mathbf{u}$ associates to every Borel set $B$ the quantities*

$$
\begin{align*}
\mu_0(B) &= \int_{B \cap S} \langle \mathbf{u} | \mathbf{n} \rangle \, d\mathcal{H}^2, \\
\mu_1(B) &= \int_{B \cap S} \left( \langle \mathbf{du} \cdot \mathbf{e}'_2 | \mathbf{e}_2 \rangle + \langle \mathbf{u} | \mathbf{n} \rangle \langle \mathbf{du} \cdot \mathbf{e}_1 | \mathbf{e}_1 \rangle \right) \, d\mathcal{H}^2 \\
&\quad + \sum_{i \neq j} \int_{S_{i,j} \cap B} \Psi \langle \mathbf{e} | \mathbf{e}_1 \rangle \, d\mathcal{H}^1, \\
\mu_2(B) &= \int_{B \cap S} \left( \langle \mathbf{du} \cdot \mathbf{e}_1 | \mathbf{e}_1 \rangle \langle \mathbf{du} \cdot \mathbf{e}'_2 | \mathbf{e}_2 \rangle - \langle \mathbf{du} \cdot \mathbf{e}_1 | \mathbf{e}_2 \rangle \langle \mathbf{du} \cdot \mathbf{e}'_2 | \mathbf{e}_1 \rangle \right) \, d\mathcal{H}^2 \\
&\quad + \sum_{i \neq j} \int_{B \cap S_{i,j}} \langle (\mathbf{u}_i - \mathbf{u}_j) \times \mathbf{e}_1 | \mathrm{d}\mathbf{e}_1 \cdot \mathbf{e} \rangle \, d\mathcal{H}^1 \\
&\quad + \sum_{p \in B \cap \mathcal{E}(S)} \mathrm{Area}(NC(p, \mathbf{u})), \\
\mu_{\Omega}(B) &= \int_{B \cap S} \left( \langle \mathbf{du} \cdot \mathbf{e}'_2 | \mathbf{e}_1 \rangle - \langle \mathbf{u} | \mathbf{n} \rangle \langle \mathbf{du} \cdot \mathbf{e}_1 | \mathbf{e}_2 \rangle \right) \, d\mathcal{H}^2 + \int_{S_{i,j}} \langle \mathbf{e} | \mathbf{u}_i - \mathbf{u}_j \rangle \, d\mathcal{H}^1,
\end{align*}
$$

where $\mathrm{Area}(NC(p, \mathbf{u}))$ is the sum of the algebraic areas of each spherical triangle where the sign is induced by its orientation.

*Proof.* We recall that $\mathbf{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.

<footnote n="2">²so that $\mathbf{e}_{ji}(p) = -\mathbf{e}_{ij}(p)$.</footnote>
**Above the faces.** Remark that for any face $S_i$ of $S$ the map $\Gamma_u : p \in S_i \mapsto (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

$$ \mu_k(B \cap S_i) = \int_{NC(S_i, u) \cap \pi^{-1}(B)} \omega_k = \int_{\Gamma_u(B \cap S_i)} \omega_k = \int_{B \cap S_i} \Gamma_u^* \omega_k $$

Hence, the computation of the curvature measures amounts to computing the pull-back by $\Gamma_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_pS$ at $p$. Using the fact that $d\Gamma_u(p) = (\mathrm{Id}, du(p))$, one gets

$$ \begin{aligned} \Gamma_u^* \omega_0(e_1, e'_2) &= \omega_0(d\Gamma_u(p) \cdot e_1, d\Gamma_u(p) \cdot e'_2) \\ &= \varepsilon_1^b \wedge \tilde{\varepsilon}_2^b - \varepsilon_2^b \wedge \tilde{\varepsilon}_1^b \quad ((e_1, du \cdot e_1), (e'_2, du \cdot e'_2)) \\ &= \begin{vmatrix} \langle e_1 | e_1 \rangle & \langle e_1 | e_2 \rangle \\ \langle e'_2 | e_1 \rangle & \langle e'_2 | e_2 \rangle \end{vmatrix} \\ &= \langle e'_2 | e_2 \rangle \\ &= \langle u | n \rangle. \end{aligned} $$

Similarly, we evaluate $\Gamma_u^* \omega_1$ and $\Gamma_u^* \omega_2$ in the same basis.

$$ \begin{aligned} \Gamma_u^* \omega_1(e_1, e'_2) &= (\varepsilon_1^b \wedge \tilde{\varepsilon}_2^b - \varepsilon_2^b \wedge \tilde{\varepsilon}_1^b) ((e_1, du \cdot e_1), (e'_2, du \cdot e'_2)) \\ &= \begin{vmatrix} \langle e_1 | e_1 \rangle & \langle du \cdot e_1 | e_2 \rangle \\ \langle e'_2 | e_1 \rangle & \langle du \cdot e'_2 | e_2 \rangle \end{vmatrix} - \begin{vmatrix} \langle e_1 | e_2 \rangle & \langle du \cdot e_1 | e_1 \rangle \\ \langle e'_2 | e_2 \rangle & \langle du \cdot e'_2 | e_1 \rangle \end{vmatrix} \\ &= \begin{vmatrix} 1 & \langle du \cdot e_1 | e_2 \rangle \\ 0 & \langle du \cdot e'_2 | e_2 \rangle \end{vmatrix} - \begin{vmatrix} 0 & \langle du \cdot e_1 | e_1 \rangle \\ \langle u | n \rangle & \langle du \cdot e'_2 | e_1 \rangle \end{vmatrix} \\ &= \langle du \cdot e'_2 | e_2 \rangle + \langle u | n \rangle \langle du \cdot e_1 | e_1 \rangle \end{aligned} $$

$$ \Gamma_u^* \omega_2(e_1, e'_2) = \begin{vmatrix} \langle du \cdot e_1 | e_1 \rangle & \langle du \cdot e_1 | e_2 \rangle \\ \langle du \cdot e'_2 | e_1 \rangle & \langle du \cdot e'_2 | e_2 \rangle \end{vmatrix} \\ = \langle du \cdot e_1 | e_1 \rangle \langle du \cdot e'_2 | e_2 \rangle - \langle du \cdot e_1 | e_2 \rangle \langle du \cdot e'_2 | e_1 \rangle $$

Finally,

$$ \begin{aligned} \Gamma_u^* \Omega(e_1, e'_2) &= (\varepsilon_1^b \wedge \tilde{\varepsilon}_1^b + \varepsilon_2^b \wedge \tilde{\varepsilon}_2^b) ((e_1, du \cdot e_1), (e'_2, du \cdot e'_2)) \\ &= \begin{vmatrix} \langle e_1 | e_1 \rangle & \langle du \cdot e_1 | e_1 \rangle \\ \langle e'_2 | e_1 \rangle & \langle du \cdot e'_2 | e_1 \rangle \end{vmatrix} + \begin{vmatrix} \langle e_1 | e_2 \rangle & \langle du \cdot e_1 | e_2 \rangle \\ \langle e'_2 | e_2 \rangle & \langle du \cdot e'_2 | e_2 \rangle \end{vmatrix} \\ &= \begin{vmatrix} 1 & \langle du \cdot e_1 | e_1 \rangle \\ 0 & \langle du \cdot e'_2 | e_1 \rangle \end{vmatrix} + \begin{vmatrix} 0 & \langle du \cdot e_1 | e_2 \rangle \\ \langle u | n \rangle & \langle du \cdot e'_2 | e_2 \rangle \end{vmatrix} \\ &= \langle du \cdot e'_2 | e_1 \rangle - \langle u | n \rangle \langle du \cdot e_1 | e_2 \rangle \end{aligned} $$
**Above the edges.** Let $p: [0, L] \to S_{i,j}$ be the arc length parameterization of $S_{i,j}$ that satisfies $\dot{p}(t) = \mathbf{e}(p(t))$ (the dot denoting the derivative w.r.t. $t$). Then the corrected normal cone $NC(S, \mathbf{u}) \cap \pi^{-1}(S_{i,j})$ above the edge $S_{i,j}$ can be parameterized by

$$ \phi: \begin{array}{ccc} [0, L] \times [0, 1] & \to & NC(S, \mathbf{u}) \cap \pi^{-1}(S_{i,j}) \\ (s, t) & \mapsto & (p(t), \mathbf{v}(s, t)) \end{array} $$

with

$$ \mathbf{v}(s, t) = \cos(s\Psi(p(t)))\mathbf{u}_i(p(t)) + \sin(s\Psi(p(t))) \mathbf{e}_1(p(t)) \times \mathbf{u}_i(p(t)). $$

Note that the map $\phi$ 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 $(p, \mathbf{v})$ is $(\mathbf{e}_1, \mathbf{e}_2, \mathbf{v})$, where $\mathbf{e}_2 = \mathbf{v} \times \mathbf{e}_1$. Finally,

$$ \dot{\mathbf{v}} = \frac{\partial \mathbf{v}}{\partial t} = \cos(s\Psi) \, \mathrm{d}\mathbf{u}_i \cdot \mathbf{e} + \sin(s\Psi) \, \mathrm{d}(\mathbf{e}_1 \times \mathbf{u}_i) \cdot \mathbf{e} + s\lambda\mathbf{e}_2 $$

where $\lambda = \partial(\Psi \circ p)/\partial t$, while $\partial\mathbf{v}/\partial s = -\Psi\mathbf{e}_2$. We then get

$$
\begin{align*}
\mu_1(S_{i,j}) &= \int_{NC(S, \mathbf{u}) \cap \pi^{-1}(S_{i,j})} \varepsilon_1^b \wedge \tilde{\varepsilon}_2^b - \varepsilon_2^b \wedge \tilde{\varepsilon}_1^b \\
&= \int_0^L \int_0^1 (\varepsilon_1^b \wedge \tilde{\varepsilon}_2^b - \varepsilon_2^b \wedge \tilde{\varepsilon}_1^b) \left( \frac{\partial \phi}{\partial s}(s, t), \frac{\partial \phi}{\partial t}(s, t) \right) \mathrm{d}s \mathrm{d}t \\
&= \int_0^L \int_0^1 (\varepsilon_1^b \wedge \tilde{\varepsilon}_2^b - \varepsilon_2^b \wedge \tilde{\varepsilon}_1^b) ((0, -\Psi \mathbf{e}_2), (\mathbf{e}, \dot{\mathbf{v}})) \mathrm{d}s \mathrm{d}t \\
&= \int_0^L \int_0^1 \left( \begin{vmatrix} 0 & -\Psi \langle \mathbf{e}_2 | \mathbf{e}_2 \rangle \\ \langle \mathbf{e} | \mathbf{e}_1 \rangle & \langle \dot{\mathbf{v}} | \mathbf{e}_2 \rangle \end{vmatrix} - \begin{vmatrix} 0 & -\Psi \langle \mathbf{e}_2 | \mathbf{e}_1 \rangle \\ \langle \mathbf{e} | \mathbf{e}_2 \rangle & \langle \dot{\mathbf{v}} | \mathbf{e}_1 \rangle \end{vmatrix} \right) \mathrm{d}s \mathrm{d}t \\
&= \int_0^L \Psi(p(t)) \langle \mathbf{e}(p(t)) | \mathbf{e}_1(p(t)) \rangle \mathrm{d}t \\
&= \int_{S_{i,j}} \Psi(p) \langle \mathbf{e}(p) | \mathbf{e}_1(p) \rangle \mathrm{d}\mathcal{H}^1(p).
\end{align*}
$$

We may note that $\mu_0$ vanishes identically above the edges, since $\mathrm{Span}(\frac{\partial \phi}{\partial \dot{\psi}}, \frac{\partial \phi}{\partial 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 $\mu_1(B \cap S_{i,j})$, it amounts to
add the indicator function of B in the integrand. Similarly, one gets

$$
\begin{align*}
\mu_2(S_{i,j}) &= \int_{NC(S,u) \cap \pi^{-1}(S_{i,j})} \tilde{\boldsymbol{\varepsilon}}_1^b \wedge \tilde{\boldsymbol{\varepsilon}}_2^b \\
&= \int_0^L \int_0^1 \tilde{\boldsymbol{\varepsilon}}_1^b \wedge \tilde{\boldsymbol{\varepsilon}}_2^b ((0, -\Psi \mathbf{e}_2), (\mathbf{e}, \dot{\mathbf{v}})) \, ds dt \\
&= \int_0^L \int_0^1 \left| \begin{array}{cc} -\Psi \langle \mathbf{e}_2 | \mathbf{e}_1 \rangle & -\Psi \langle \mathbf{e}_2 | \mathbf{e}_2 \rangle \\ \langle \dot{\mathbf{v}} | \mathbf{e}_1 \rangle & \langle \dot{\mathbf{v}} | \mathbf{e}_2 \rangle \end{array} \right| \, ds dt \\
&= \int_0^L \int_0^1 \Psi \langle \dot{\mathbf{v}} | \mathbf{e}_1 \rangle \, ds dt \\
&= \int_0^L \int_0^1 \Psi \langle \cos(s\Psi) \mathbf{du}_i \cdot \mathbf{e} + (1 - \cos\Psi) \mathbf{d}(\mathbf{e}_1 \times \mathbf{u}_i) \cdot \mathbf{e} | \mathbf{e}_1 \rangle \, ds dt \\
&= \int_0^L \langle \sin\Psi \mathbf{du}_i \cdot \mathbf{e} + (1 - \cos\Psi) \mathbf{d}(\mathbf{e}_1 \times \mathbf{u}_i) \cdot \mathbf{e} | \mathbf{e}_1 \rangle \, dt \\
&= \int_0^L \left( \sin\Psi \langle \mathbf{du}_i \cdot \mathbf{e} | \mathbf{e}_1 \rangle + (1 - \cos\Psi) \langle \mathbf{d}(\mathbf{e}_1 \times \mathbf{u}_i) \cdot \mathbf{e} | \mathbf{e}_1 \rangle \right) \, dt
\end{align*}
$$

By differentiating the expression $\langle \mathbf{u}_i | \mathbf{e}_1 \rangle = 0$, we write $\langle \mathbf{du}_i \cdot \mathbf{e} | \mathbf{e}_1 \rangle = -\langle \mathbf{u}_i | \mathbf{d}\mathbf{e}_1 \cdot \mathbf{e} \rangle$. Similarly, one has $\langle \mathbf{d}(\mathbf{e}_1 \times \mathbf{u}_i) \cdot \mathbf{e} | \mathbf{e}_1 \rangle = -\langle \mathbf{e}_1 \times \mathbf{u}_i | \mathbf{d}\mathbf{e}_1 \cdot \mathbf{e} \rangle$. Using in addition that $\mathbf{u}_i = (\mathbf{e}_1 \times \mathbf{u}_i) \times \mathbf{e}_1$, we see that

$$
\begin{align*}
\mu_2(S_{i,j}) &= - \int_0^L \sin \Psi \langle \mathbf{u}_i | \mathbf{d}\mathbf{e}_1 \cdot \mathbf{e} \rangle + (1 - \cos \Psi) \langle \mathbf{e}_1 \times \mathbf{u}_i | \mathbf{d}\mathbf{e}_1 \cdot \mathbf{e} \rangle dt \\
&= - \int_0^L \langle \sin \Psi \mathbf{u}_i + (1 - \cos \Psi) (\mathbf{e}_1 \times \mathbf{u}_i) | \mathbf{d}\mathbf{e}_1 \cdot \mathbf{e} \rangle dt \\
&= - \int_0^L \langle (\sin \Psi \mathbf{e}_1 \times \mathbf{u}_i + (\cos \Psi - 1) \mathbf{u}_i) \times \mathbf{e}_1 | \mathbf{d}\mathbf{e}_1 \cdot \mathbf{e} \rangle dt \\
&= - \int_0^L \langle (\mathbf{u}_j - \mathbf{u}_i) \times \mathbf{e}_1 | \mathbf{d}\mathbf{e}_1 \cdot \mathbf{e} \rangle dt \\
&= - \int_{S_{i,j}} \langle (\mathbf{u}_j - \mathbf{u}_i) \times \mathbf{e}_1 | \mathbf{d}\mathbf{e}_1 \cdot \mathbf{e} \rangle d\mathcal{H}^1
\end{align*}
$$
Finally,

$$
\begin{align*}
\mu_{\Omega}(S_i, j) &= \int_{NC(S, \mathbf{u}) \cap \pi^{-1}(S_i, j)} \varepsilon_1^b \wedge \tilde{\varepsilon}_1^b + \varepsilon_2^b \wedge \tilde{\varepsilon}_2^b \\
&= \int_0^L \int_0^1 (\varepsilon_1^b \wedge \tilde{\varepsilon}_1^b + \varepsilon_2^b \wedge \tilde{\varepsilon}_2^b) ((0, -\Psi \mathbf{e}_2), (\mathbf{e}, \dot{\mathbf{v}})) \, ds dt \\
&= \int_0^L \int_0^1 \left( \begin{vmatrix} 0 & -\Psi \langle \mathbf{e}_2 | \mathbf{e}_1 \rangle \\ \langle \mathbf{e} | \mathbf{e}_1 \rangle & \langle \dot{\mathbf{v}} | \mathbf{e}_1 \rangle \end{vmatrix} + \begin{vmatrix} 0 & -\Psi \langle \mathbf{e}_2 | \mathbf{e}_2 \rangle \\ \langle \mathbf{e} | \mathbf{e}_2 \rangle & \langle \dot{\mathbf{v}} | \mathbf{e}_2 \rangle \end{vmatrix} \right) \, ds dt \\
&= \int_0^L \int_0^1 \langle \mathbf{e} | \Psi \mathbf{e}_2 \rangle \, ds dt = - \int_0^L \int_0^1 \langle \mathbf{e} | \frac{\partial \mathbf{v}}{\partial s} \rangle \, ds dt \\
&= \int_{S_{i,j}} \langle \mathbf{e} | \mathbf{u}_i - \mathbf{u}_j \rangle \, d\mathcal{H}^1
\end{align*}
$$

**Above the vertices.** When $p$ is a vertex, its contribution is restricted to $\mu_2$, since the position is constant in $NC(p)$. Since $\omega_2$ measures the area in the velocity component, the vertex adds a Dirac mass at $p$ to $\mu_2$, with coefficient equal to the (signed) area of $NC(p)$. $\square$

**Remark 6.** In the smooth case, there are no edges nor vertices, and taking $\mathbf{u} = \mathbf{n}$ (and hence $\mathbf{e}'_2 = \mathbf{e}_2$), the formula simplify. We recover

$$ \mu_0 = d\mathcal{H}^2, \mu_1 = \mathrm{Trace}(\mathrm{d}\mathbf{n}) \mathrm{d}\mathcal{H}^2 = -2H\mathrm{d}\mathcal{H}^2 \text{ and } \mu_2 = \det(\mathrm{d}\mathbf{n}) \mathrm{d}\mathcal{H}^2 = K\mathrm{d}\mathcal{H}^2 $$

while

$$ \mu_{\Omega} = (\langle \mathrm{d}\mathbf{n} \cdot \mathbf{e}_2 | \mathbf{e}_1 \rangle - \langle \mathrm{d}\mathbf{n} \cdot \mathbf{e}_1 | \mathbf{e}_2 \rangle) \mathrm{d}\mathcal{H}^2 = 0 $$

since $\mathrm{d}\mathbf{n}$ is symmetric on $T_pS$.

**Remark 7.** Similarly to [9], we can consider a vector-valued equivariant two-form on the Grassmann bundle defined for each fixed pair of vectors $\mathbf{X}, \mathbf{Y} \in \mathbb{R}^3$. This leads to an anisotropic curvature measure $\mu^{\mathbf{X},\mathbf{Y}}$ which converges to the second fundamental form by the same theorems in Section 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 (section 6).

**Proposition 3.** Let $S$ be a polyhedral surface and $\mathbf{u}$ a vector field constant on each planar face of $S$. Then the corrected Lipschitz-Killing curvature measures
<figure><img src="image_4.png" /><figcaption>Figure 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 $\frac{h}{m}$ and angle $\frac{\pi}{n}$. Any vertex is connected to the two nearby vertices above and the two nearby vertices below (except for the top and bottom levels).</figcaption></figure>

are

$$ \mu_0 = \sum_f \cos \alpha(f) d\mathcal{H}^2, $$

$$ \mu_1 = \sum_{S_{i,j}} \Psi \langle \mathbf{e} | \mathbf{e}_1 \rangle d\mathcal{H}^1 = \sum_{S_{i,j}} \frac{\Psi}{\sin \Psi} \det(\mathbf{e}, \mathbf{u}_i, \mathbf{u}_j) d\mathcal{H}_1, $$

$$ \mu_2 = \sum_p \text{area}(NC(p)) \delta_p, $$

$$ \mu_{\Omega} = \sum_{S_{i,j}} \langle \mathbf{e} | \mathbf{u}_i - \mathbf{u}_j \rangle d\mathcal{H}_1, $$

where $\alpha(f)$ is the angle between the corrected normal $\mathbf{u}$ and the naive normal $\mathbf{n}$, $\Psi$ is the oriented angle between the corrected normal $\mathbf{u}_i$ and $\mathbf{u}_j$ incident to the edge $S_{i,j}$ and $\mathbf{e}_1 = \frac{\mathbf{u}_i \times \mathbf{u}_j}{\|\mathbf{u}_i \times \mathbf{u}_j\|}$, while $\mathbf{e}$ is the oriented unit tangent to the edge.

## 3.3 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.

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 $\alpha$ between the corrected normal and the face's normal satisfies $\tan \alpha = \frac{2mr}{h} \sin^2 \frac{\pi}{2n}$, so that the corrected area of the lantern is

$$ \int_T \mu_0 = \cos \alpha \text{ area}(T) = \cos \alpha \times \frac{1}{2} \times 2r \sin \frac{\pi}{n} \times \frac{h}{m} \sqrt{1 + \tan^2 \alpha} = \frac{rh}{m} \sin \frac{\pi}{n}, $$

$$ \int_L \mu_0 = \sum_{i=1}^{m} \sum_{j=1}^{2n} \frac{rh}{m} \sin \frac{\pi}{n} \underset{n \to \infty}{\longrightarrow} 2\pi rh = \text{area}(C) \text{ independently of } m. $$

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 $\mathbf{u}_i, \mathbf{u}_j$ are horizontal (so that $\mathbf{e}_1 = \mathbf{u}_i \times \mathbf{u}_j$ is vertical) and at angular distance $\Psi = \pi/n$. For an edge $e$, the scalar product $\langle \mathbf{e} | \mathbf{e}_1 \rangle \ell(e)$ is then $h/m$, thus $\mu_1(e) = \frac{\pi h}{nm}$. Since there are exactly $2nm$ such edges, $\int_L \mu_1 = 2\pi h = 2 \int_C \frac{1}{2r} dA$. 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)$ convergent one for the area.

Now consider the same triangulation with a different corrected normal $\mathbf{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 $\pi$ on $C$: $\mathbf{u}(p) = r^{-1}\pi(p)$. On a face, $\mathbf{e}_1 = \mathbf{u} \times \mathbf{n}$ and $\mathbf{e}'_2 = \mathbf{n} \times \mathbf{e}_1$. Since $d\pi$ is the projection to the tangent plane to the cylinder, we have that $d\mathbf{u} \cdot \mathbf{e}_1 = r^{-1}\mathbf{e}_1$ and $d\mathbf{u} \cdot \mathbf{e}'_2 = r^{-1}\langle \mathbf{u} | \mathbf{n} \rangle \mathbf{e}_2$.

The field $\mathbf{u}$ being continuous on the edges and at the vertices, the curvature measures are carried only by the faces. Using Proposition 2, we get:

* the corrected area measure $\mu_0 = \langle \mathbf{u} | \mathbf{n} \rangle d\mathcal{H}^2$ is no other than the pullback of the area form on the cylinder by the horizontal projection, i.e. $\pi^*dA_C$; hence $\mu_0$ measures the area of the projected lantern $L$ onto the cylinder $C$, and globally gives exactly the area of the cylinder,

* for $\mu_1$, we have that

$$ \mu_1 = \langle d\mathbf{u} \cdot \mathbf{e}'_2 | \mathbf{e}_2 \rangle + \langle \mathbf{u} | \mathbf{n} \rangle \langle d\mathbf{u} \cdot \mathbf{e}_1 | \mathbf{e}_1 \rangle d\mathcal{H}^2 = \frac{2}{r} \langle \mathbf{u} | \mathbf{n} \rangle d\mathcal{H}^2 = \frac{2}{r} \mu_0; $$

hence we obtain exactly 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 $\mathbb{S}^2$; therefore $\mu_2$ vanishes identically as in the previous example.
# 4 Stability of the corrected curvature measures

We show in this section that the corrected curvatures of $(S, \mathbf{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 $\mathbf{u}$ approximates well the normal $\mathbf{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.

## 4.1 Background on sets with positive reach

The distance function $d_K$ of a compact set $K$ of $\mathbb{R}^d$ associates to any point $x$ of $\mathbb{R}^d$ its distance to $K$, namely $d_K(x) := \min_{y \in K} d(x, y)$, where $d$ is the Euclidean distance on $\mathbb{R}^d$. For a given real number $\varepsilon > 0$, we denote by $K^\varepsilon := \{x \in \mathbb{R}^d, d_K(x) \le \varepsilon\}$ the $\varepsilon$-offset of $K$. The Hausdorff distance $d_H(K, K')$ between two compact sets $K$ and $K'$ is the minimum number $\varepsilon$ such that $K \subset K'^\varepsilon$ and $K' \subset K^\varepsilon$.

The medial axis of $K$ is the set of points of $x \in \mathbb{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 $\text{reach}(K)$, is the infimum distance between $K$ and its medial axis. As a consequence, whenever the reach is positive and $\varepsilon < \text{reach}(K)$, the projection map

$$ \pi_K : K^\varepsilon \to K $$

is well-defined. It is well known that smooth compact submanifolds have positive reach [12, 26]. The following proposition will be useful for stating our main theorem.

**Proposition 4.** Let $X$ be a compact submanifold of class $C^2$ of $\mathbb{R}^d$. Then

$$ 0 < \text{reach}(X) \le 1/\rho_X, $$

where $\rho_X$ is the largest absolute value of the principal curvatures of $X$. Furthermore $\pi_X$ is differentiable on $X^\varepsilon$ for any $\varepsilon < \text{reach}(X)$ and

$$ \forall x \in X^\varepsilon \quad \|D\pi_X(x)\| \le \frac{1}{1 - \rho_X(\pi_X(x))\varepsilon}, $$

where $\rho_X(x)$ is the largest absolute value of the principal curvatures of $X$ at the point $x$.

## 4.2 Stability result

We provide in this section a stability result on the curvature measures, namely we show that the corrected curvature measures of $(S, \mathbf{u})$ approximate the curvature measures of $X$ provided that $S$ and $X$ are close in the Hausdorff sense and that $\mathbf{u}$ is close to the unit normal vector of $X$. In order to state the result,
we denote by $\mu_k^{S, u}$ the corrected curvature measures of $(S, u)$ and by $\mu_k^X$ the curvature measures of $X$, where $k \in \{0, 1, 2, \Omega\}$. 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$.

**Theorem 1.** Let $X$ be a compact surface of $\mathbb{R}^3$ of class $C^2$, of normal vector $n$, bounding a volume $V$, and $S = \cup_i S_i$ be a piecewise $C^{1,1}$ surface bounding a volume $W$, $u$ a corrected normal vector field on $S$. We assume the following two conditions:

* there exists an open set $U$ such that $\pi_X : U \cap S \to X$ is injective.

* $\varepsilon := d_H(S, X) < \text{reach}(X)$ is the position error.

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 = \cup_{i \in I} S_i$ of faces $S_i$ of $S$, one has

$$ \begin{aligned} |\mu_k^{S, u}(B) - \mu_k^X(\pi_X(B))| &\le 4 (\eta + \varepsilon) \left( \frac{\max(1, \rho_B)}{1 - \rho_B \varepsilon} \right)^2 \\ &\quad \left( (1 + L_u^2) (\mathcal{H}^2(B) + \mathcal{H}^1(\partial B)) + \frac{\delta_u^2 d_{\max}^3}{2\pi} N_v^S(B) \right. \\ &\qquad \left. + 2 \arcsin \left( \frac{\delta_u}{2} \right) \left( \sum_{S_{i,j} \subset B} \mathcal{H}^1(S_{i,j}) + d_{\max} N_v^S(\partial B) \right) \right), \end{aligned} $$

where $\delta_u = \sup_{i \sim j} \sup_{p \in S_{i,j}} \|u_i(p) - u_j(p)\|$ ($i \sim j$ whenever $S_i$ is adjacent to $S_j$), $\eta := \sup_{p \in S} \|u(p) - n(\pi_X(p))\|$ is the normal error, $\rho_B$ is the maximum absolute value of the principal curvatures of $X \cap \pi_X(B)$, $L_u$ is the maximum of the Lipschitz constants of $u$ over each face $S_i$, $d_{\max}$ is the maximum vertex degree, $N_v^S(B)$ and $N_v^S(\partial B)$ are the number of vertices of $S$ that respectively belong to $B$ and $\partial B$.

**Remark 8.** It may seem quite restrictive to evaluate the measures $\mu_k^{S, u}$ and $\mu_k^X$ 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 $\hat{H}^u(B)$ and $\hat{G}^u(B)$ of Definition 6. However, the faces can be subdivided at leasure 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 1 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_v^S(B)$ and $N_v^S(\partial 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.
**Corollary 1.** Under the same assumptions as in Theorem 1, if **u** is continuous over the whole surface S, one has:

$$|\mu_k^{S, \mathbf{u}}(B) - \mu_k^X(\pi_X(B))| \le 4 (\eta + \varepsilon) \left( \frac{\max(1, \rho_B)}{1 - \rho_B \varepsilon} \right)^2 (1 + L_{\mathbf{u}}^2) (\mathcal{H}^2(B) + \mathcal{H}^1(\partial B)).$$

**Corollary 2.** When S is a digital surface approximation with an estimated normal vector field **u** constant per face, we have the following simplified estimate

$$|\mu_k^{S, \mathbf{u}}(B) - \mu_k^X(\pi_X(B))| \le 4 (\eta + \varepsilon) \left( \frac{\max(1, \rho_B)}{1 - \rho_B \varepsilon} \right)^2 \left( \frac{2N_i + N_b}{4} h^2 + N_b h + 35 \delta_{\mathbf{u}}^2 \left( \frac{N_i}{2} + N_b \right) + 2 \arcsin \frac{\delta_{\mathbf{u}}}{2} (h(N_i + N_b) + 6N_b) \right),$$

where B is an union of surfels of edge size h, with N_i inner edges and N_b boundary edges.

*Proof.* In the digital case, a vertex has at most 6 neighbors so $d_{\max} \le 6$. (We then major $6^3/(2\pi)$ by 35.) Since **u** is constant per-face, we have $L_{\mathbf{u}} = 0$. Furthermore, the length of each edge is h and the area of each face is h$^2$. The formula follows. □

Estimation of curvatures on digital surfaces is studied in more details in Section 5. The remaining of this section is devoted to the proof of Theorem 1.

## 4.3 Notions of geometric measure theory

**Currents.** We denote by $D_m$ the set of $m$-differential forms on $\mathbb{R}^d$. This set can be endowed with the $C^\infty$-topology [31] and the set of $m$-currents on $\mathbb{R}^d$ is then by definition its topological dual. Therefore a current $T: D_m \to \mathbb{R}$ is a continuous linear map. The support of a current $T: D_m \to \mathbb{R}$ is the smallest closed set $K$ such that $\text{spt}(\omega) \cap K = \emptyset \Rightarrow T(\omega) = 0$. The boundary $\partial T$ of a $m$-current $T$ is the $(m-1)$-current defined by $\partial T(\omega) = T(d\omega)$, where d is the exterior derivative.

When $U$ and $V$ are open sets in Euclidean spaces, $f: U \to V$ is of class $C^\infty$, $T$ is a $m$-current on $U$ and $f_{|\text{spt}(T)}$ is proper, one define the $m$-current $f_{\#}T$ by $f_{\#}T(\omega) = T(f^*\omega)$ (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 $\mathbb{R}^d$ is $m$-rectifiable, then it is of class $C^1 \mathcal{H}^m$-almost everywhere and it is possible to integrate differential forms over it, and thus to define currents of the form

$$\omega \mapsto \int_S \omega = \int_S \omega(e_1, \dots, e_m) d\mathcal{H}^m(x),$$

where $(e_1, \dots, 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. A *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*<sub>*m*</sub>(ℝ<sup>*d*</sup>), one has

$$
T(\omega) = \int_S \mu \omega,
$$

where $\mu(x) \in \mathbb{Z}$ is the multiplicity and satisfies $\int_S \mu(x) d\mathcal{H}^m(x) < \infty$. When $\mu(x) = k$ is constant, we denote $T = k \cdot 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 $\partial 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

$$
M(T) = \sup\{T(\omega), \omega \in D_m(\mathbb{R}^d) \text{ and } \forall x \in \mathbb{R}^d \|\omega(x)\|^* \le 1\},
$$

where $\|\omega(x)\|^* = \sup\{\omega(x)(\xi_1, \dots, \xi_m), \xi_i \in \mathbb{R}^d \text{ and } \|\xi_i\| \le 1\}$. The flat norm $\mathcal{F}(T)$ of $T$ is given by

$$
\mathcal{F}(T) = \inf \left\{ M(A) + M(B), T = A + \partial B, \begin{array}{l} A \text{ is } m\text{-rectifiable} \\ B \text{ is } (m+1)\text{-rectifiable} \end{array} \right\}.
$$

**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 2 (Constancy Theorem).** Let $X$ be an $m$-dimensional oriented submanifold of $\mathbb{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 $\lambda$ such that

$$
T(\omega) = \lambda \int_X \omega.
$$

**4.4 Proof of Theorem 1**

The proof is based on the *homotopy formula for currents* and is in the same
spirit than the proof of [9]. We recall that the projection map $\pi_X : X^\varepsilon \to X$ is
well defined and differentiable since $\varepsilon < \text{reach}(K)$. We first build the Lipschitz
map

$$
f: \begin{array}{rcl} X^{\epsilon} \times \mathbb{R}^3 & \to & \operatorname{spt}(N(X)) \\ (p, u) & \mapsto & (\pi_X(p), \mathbf{n}_X(\pi_X(p))). \end{array}
$$

**Lemma 1.**

$f_{\#} N(S, \mathbf{u}) = N(X).$
*Proof.* By definition, $f_{\sharp}N(S, \mathbf{u})$ is a 2-current supported in the set spt($N(X)$). Since $N(S, \mathbf{u})$ has no boundary, $f_{\sharp}N(S, \mathbf{u})$ also has no boundary. Furthermore, spt($N(X)$) is of class $C^1$, so by the Constancy Theorem (Theorem 2), one has $f_{\sharp}N(S, \mathbf{u}) = \lambda \cdot N(X)$, where $\lambda$ is an integer. We know by assumption that there exists an open set $U$ such that $\pi_X: U \cap S \to X$ is injective. This implies that the restriction of $f$ to spt($N(S, \mathbf{u})$) is one-to-one onto its image, so the constant $\lambda = 1$. $\square$

We denote by $D = N(S, \mathbf{u})_{\perp}(B \times \mathbb{R}^3)$ the restriction of the current $N(S, \mathbf{u})$ to $B \times \mathbb{R}^3$. By Lemma 1, one has

$$f_{\sharp}D = f_{\sharp}(N(S, \mathbf{u})_{\perp}(B \times \mathbb{R}^3)) = N(X)_{\perp}(\pi_X(B) \times \mathbb{R}^3).$$

Let $h$ be the affine homotopy between $f$ and the identity:

$$h: \begin{array}{rcl} [0,1] \times X^{\varepsilon} \times \mathbb{R}^3 & \to & \mathbb{R}^3 \times \mathbb{R}^3 \\ (t,x) & \mapsto & (1-t)(x) + tf(x). \end{array}$$

In the remaining of the proof, in order to shorten equations, we put $x = (p, \mathbf{u})$. Using the *affine homotopy formula* for locally Lipschitz map (see [23, page 187] or [13, 4.1.14]), one has

$$N(S, \mathbf{u})_{\perp}(B \times \mathbb{R}^3) - N(X)_{\perp}(\pi_X(B) \times \mathbb{R}^3) = D - f_{\sharp}D = \partial E - F,$$

where $E = h_{\sharp}(D \times [0, 1])$ and $F = h_{\sharp}(\partial D \times [0, 1])$ are respectively 3-current and 2-current that satisfy

$$M(E) \leq \sup_{x \in \text{spt}(D)} \|x - f(x)\| \cdot \sup_{x \in \text{spt}(D)} (1, \|Df(x)\|^2) \cdot M(D),$$

$$M(F) \leq \sup_{x \in \text{spt}(D)} \|x - f(x)\| \cdot \sup_{x \in \text{spt}(D)} (1, \|Df(x)\|) \cdot M(\partial D),$$

where the norm of the linear map $Df(x)$ is $\|Df(x)\| = \sup_{\|h\|=1} Df(x) \cdot h$. By definition of the flat norm one has

$$\begin{align} \mathcal{F}(N(S, \mathbf{u})_{\perp}(B \times \mathbb{R}^3) - N(X)_{\perp}(\pi_X(B) \times \mathbb{R}^3)) & \nonumber \\ & \leq M(E) + M(F) \nonumber \\ & \leq \sup_{x \in \text{spt}(D)} \|x - f(x)\| \cdot \sup_{x \in \text{spt}(D)} (1, \|Df(x)\|^2) \cdot (M(D) + M(\partial D)). \tag{1} \end{align}$$

The following lemmas bound the terms involved in the right hand side term of the previous equation.

**Lemma 2.** For every $x \in \text{spt}(D)$, one has

$$\|x - f(x)\| \leq \sqrt{\eta^2 + \varepsilon^2} \leq \eta + \varepsilon \quad \text{and} \quad \max(1, \|Df(x)\|^2) \leq \left( \frac{\max(1, \rho_B)}{1 - \rho_B \varepsilon} \right)^2.$$
*Proof.* Clearly, by definition of the Euclidean norm, one has

$$ \|x - f(x)\| \le \sqrt{\eta^2 + \varepsilon^2} \le \eta + \varepsilon. $$

Using the fact that $\|D\pi_X(x)\| \le 1/(1 - \rho_X(\pi_X(x))\varepsilon)$ for every $x$ in $X^\varepsilon$ and that $\|D\mathbf{n}_X(p)\| \le \rho_X(p)$, one has by composition:

$$ \|Df(x)\| \le \max \left( \frac{1}{1 - \rho_B \varepsilon}, \frac{\rho_B}{1 - \rho_B \varepsilon} \right) = \frac{\max(1, \rho_B)}{1 - \rho_B \varepsilon}. $$

The conclusion follows from the fact that the upper bound is greater than 1. $\square$

**Lemma 3.**

$$ M(D) \le (1 + L_{\mathbf{u}}^2) \mathcal{H}^2(B) + 2 \arcsin(\frac{\delta_{\mathbf{u}}}{2}) \sum_{S_{i,j} \subset B} \mathcal{H}^1(S_{i,j}) + \frac{\delta_{\mathbf{u}}^2 d_{\max}^3}{2\pi} N_v^S(B). $$

*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_i^0, \mathbf{u})$ above the interior $S_i^0$ of the face $S_i$ is parameterized by $\Gamma_{\mathbf{u}} : p \in S_i^0 \mapsto (p, \mathbf{u}(p))$. Since $\mathbf{u}$ is $L_{\mathbf{u}}$-Lipschitz on each face $S_i$, the map $\Gamma_{\mathbf{u}}$ is $\sqrt{1 + L_{\mathbf{u}}^2}$-Lipschitz. The mass above $S_i^0$ is given by the General Area-Coarea Formula [25, 3.13]

$$ \int_{NC(S_i^0, \mathbf{u})} d\mathcal{H}^2 = \int_{S_i^0} J_2(\Gamma_{\mathbf{u}}) d\mathcal{H}^2 \le \int_{S_i^0} \sqrt{1 + L_{\mathbf{u}}^2} d\mathcal{H}^2 = (1 + L_{\mathbf{u}}^2) \mathcal{H}^2(S_i) $$

where $J_2(\Gamma_{\mathbf{u}})$ denotes the 2-Jacobian. Summing over all the faces $S_i^0$ of $S$, one has

$$ M^f \le (1 + L_{\mathbf{u}}^2) \mathcal{H}^2(S). $$

**Above the edges of $S$.** We recall that $\pi : \text{spt}(N(s, \mathbf{u})) \to \mathbb{R}^3$, defined by $\pi(p, n) = p$, denotes the restriction to the support of $N(S, \mathbf{u})$ of the projection onto the position component. Note that the support of $N(S, \mathbf{u})$ above an edge $S_{i,j}$ is given by $\text{spt}(N(S, \mathbf{u}) \perp (S_{i,j} \times \mathbb{R}^3)) = \pi^{-1}(S_{i,j})$. The Coarea formula thus gives:

$$ \begin{aligned} M(N(S, \mathbf{u}) \perp (S_{i,j} \times \mathbb{R}^3)) &= \int_{\pi^{-1}(S_{i,j})} d\mathcal{H}^2 \\ &= \int_{S_{i,j}} \left( \int_{\pi^{-1}(p)} \frac{1}{J_1(\pi)} d\mathcal{H}^1 \right) d\mathcal{H}^1(p) \\ &= \int_{S_{i,j}} \mathcal{H}^1(\pi^{-1}(p)) d\mathcal{H}^1(p) \\ &= \int_{S_{i,j}} 2 \arcsin \left( \frac{\|u_i(p) - u_j(p)\|}{2} \right) d\mathcal{H}^1(p) \\ &\le 2 \arcsin(\frac{\delta_{\mathbf{u}}}{2}) \mathcal{H}^1(S_{i,j}). \end{aligned} $$
The third line follows from the fact that the 1-Jacobian of $\pi$ satisfies $J_1(\pi) = 1$ and that $\pi^{-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

$$M^e \le 2 \arcsin\left(\frac{\delta_{\mathbf{u}}}{2}\right) \sum_{S_{i,j} \subset B} \mathcal{H}^1(S_{i,j}).$$

Above the vertices of $S$. Let $p \in B$ be a vertex of $S$. The support of the corrected cycle $N(S, \mathbf{u})_{\perp}(\{p\} \times \mathbb{R}^3)$ above $p$ obviously lies in the set $\{p\} \times \mathbb{S}^2$. By construction, the multiplicity $\mu(\mathbf{n})$ of a point $(p, \mathbf{n})$ lying in this support is an integer that satisfies $|\mu(\mathbf{n})| \le \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 $\|\mathbf{u}_i(p) - \mathbf{u}_j(p)\| \le \delta_{\mathbf{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, \mathbf{u})_{\perp}(\{p\} \times \mathbb{R}^3)$ is at most $L = \deg(p)\delta_{\mathbf{u}}$. By the isoperimetric inequality, this implies that the area of this normal component is at most $L^2/(2\pi) \le \deg(p)^2 \delta_{\mathbf{u}}^2/(2\pi)$. Therefore, one has:

$$M(N(S, \mathbf{u})_{\perp}(\{p\} \times \mathbb{R}^3)) = \int_{\mathbb{S}^2} |\mu(n)| d\mathcal{H}^2(n) \le \frac{\delta_{\mathbf{u}}^2 \deg(p)^2}{2\pi} \deg(p).$$

Summing over all the vertices of $S \cap \partial B$, one has

$$M^v \le \frac{\delta_{\mathbf{u}}^2}{2\pi} \sum_{\text{vertices } p} \deg(p)^3 \le \frac{\delta_{\mathbf{u}}^2 d_{\max}^3}{2\pi} N_v^S(B).$$

□

**Lemma 4.**

$$M(\partial D) \le \sqrt{1 + L_{\mathbf{u}}^2} \mathcal{H}^1(\partial B) + 2d_{\max} \arcsin\left(\frac{\delta_{\mathbf{u}}}{2}\right) N_v^S(\partial B).$$

*Proof.* Since $B = \cup_{i \in I} S_i$ is a connected union of faces $S_i$, its boundary is a union of edges $S_{i,j}$ where $i \in I$ and $j \notin I$. The boundary $\partial D$ of the restricted current $D = N(S, \mathbf{u})_{\perp}(B \times \mathbb{R}^3)$ is supported on curves of two different kinds.

Above the relative interior of each $S_{i,j} \subset \partial B$, the support $D$ is parameterized by $\Gamma_{\mathbf{u}_j} : S_{i,j}^0 \to \mathbb{R}^3 \times \mathbb{R}^3$. One then has

$$\begin{align*}
M((\partial D)_{\perp}(S_{i,j}^0 \times \mathbb{R}^3)) &= \int_{\Gamma_{\mathbf{u}_j}(S_{i,j}^0)} d\mathcal{H}^1 \\
&= \int_{S_{i,j}^0} J_1(\Gamma_{\mathbf{u}_j}) \, d\mathcal{H}^1 \\
&\le \int_{S_{i,j}^0} \sqrt{1 + L_{\mathbf{u}}^2} \, d\mathcal{H}^1 \\
&\le \sqrt{1 + L_{\mathbf{u}}^2} \mathcal{H}^1(S_{i,j}).
\end{align*}$$
Above a vertex $p$ of $S \cap \partial 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 $\mathbf{u}_{j_1}(p)$ and $\mathbf{u}_{j_2}(p)$, and is of length less than $2 \arcsin(\|\mathbf{u}_{j_1}(p) - \mathbf{u}_{j_2}(p)\|/2) \le 2 \arcsin(\delta_{\mathbf{u}}/2)$. Therefore

$$M((\partial D)_{\perp}(\{p\} \times \mathbb{R}^3)) \le 2 \deg(p) \arcsin(\frac{\delta_{\mathbf{u}}}{2}).$$

The result follows by summing over all the vertices and edges of $\partial B$. $\square$

End of proof of Theorem 1. Plugging the results of Lemmas 2, 3 and 4 with Equation (1), one gets:

$$\begin{aligned} & \mathcal{F}(N(S, \mathbf{u})_{\perp}(B \times \mathbb{R}^3) - N(X)_{\perp}(\pi_X(B) \times \mathbb{R}^3)) \\ & \le (\eta + \varepsilon) \left( \frac{\max(1, \rho_X)}{1 - \rho_X \varepsilon} \right)^2 \\ & \quad \left( (1 + \mathrm{L}_{\mathbf{u}}^2) (\mathcal{H}^2(B) + \mathcal{H}^1(\partial B)) + \frac{\delta_{\mathbf{u}}^2 \mathrm{d}_{\mathrm{max}}^3}{2\pi} N_v^S(B) \right. \\ & \quad \left. + 2 \arcsin(\frac{\delta_{\mathbf{u}}}{2}) \left( \sum_{S_{i,j} \subset B} \mathcal{H}^1(S_{i,j}) + \mathrm{d}_{\mathrm{max}} N_v^S(\partial B) \right) \right), \end{aligned}$$

For every form $\omega \in \{\omega_0, \omega_1, \omega_2, \Omega\}$, one has $\|\omega\| \le 2$ and $\|d\omega\| \le 4$. Let $i \in \{0, 1, 2\}$. One has

$$\begin{aligned} & |\mu_k^S(B) - \mu_k^X(\pi_X(B))| \\ & = |N(S, \mathbf{u})_{\perp}(B \times \mathbb{R}^3)\omega_k - N(X)_{\perp}(\pi_X(B) \times \mathbb{R}^3)\omega_k| \\ & \le \sup(\|\omega_k\|, \|d\omega_k\|) \mathcal{F}(N(S, \mathbf{u})_{\perp}(B \times \mathbb{R}^3) - N(X)_{\perp}(\pi_X(B) \times \mathbb{R}^3)) \\ & \le 4 \mathcal{F}(N(S, \mathbf{u})_{\perp}(B \times \mathbb{R}^3) - N(X)_{\perp}(\pi_X(B) \times \mathbb{R}^3)), \end{aligned}$$

which implies the result.

# 5 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 $\mathbb{R}^3$ whose boundary $X := \partial V$ is of class $C^3$. We assume that the reach of $X$ is greater than $\rho > 0$. Let $h > 0$ be the sampling gridstep of the regular grid $h\mathbb{Z}^3$. The Gauss digitization of $V$ is defined as $G_h(V) := V \cap h\mathbb{Z}^3$. Digital sets are subsets of $h\mathbb{Z}^3$. Let us denote $Q_z^h$ the axes-aligned cube of edge length $h$ centered on a point $z \in h\mathbb{Z}^3$. Digital sets can thus be seen as a union of such cubes in $\mathbb{R}^3$. We can now defined the digitized surface $X_h$ of $X := \partial V$ at step $h$ as the topological boundary of the Gauss digitization of $V$, seen as a union of cubes:

$$X_h := \partial(\bigcup_{z \in G_h(V)} Q_z^h).$$

We call *surfel* any square of edge length $h$ that is the face of some $Q_z^h$, with $z \in 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) an open set $U$ should exists where the projection $\pi_X : U \cap S \to 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 Theorem 2 of [22]:

**Theorem 3.** Let $h < 0.198\rho$, letting $y \in 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* of $X_h$ and $s \cap X = \emptyset$, or (iii) $y$ belongs to a *linel* of $X_h$ and $\exists p \in s \cap X$ but the angle $\alpha_y$ between $s$ and the normal to $X$ at $p$ satisfies $\alpha_y \ge 1.260h/\rho$.

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), Theorem 3 of [22] 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 $\pi_X : U \cap X_h \to 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_p^R$ the ball centered at $p$ and of radius $R$.
**Lemma 5.** Let $p$ be an arbitrary point of $X_h$. Let $B$ be the surfels of $X_h$ that lie in $X_h \cap \mathbb{B}_p^R$. Then the numbers $N_i$ of inner linels (the ones bordering two surfels of $B$) and $N_b$ of boundary linels of $B$ (the ones bordering only one surfel of $B$) follow these bounds for radius $R = Kh^\alpha$:

$$N_i = \mathcal{O}(h^{4\alpha-3} + h^{2\alpha-2}), \quad N_b = \mathcal{O}(h^{\alpha-1} + h^{7\alpha-3}).$$

The proof of Lemma 5 relies on the following intermediary lemma.

**Lemma 6.** Under the same assumptions as in Lemma 5, one has

$$X_h \cap \partial \mathbb{B}_p^R \subset (X \cap \partial \mathbb{B}_p^R)^{\mathcal{O}(h+h^{3\alpha})}.$$

of Lemma 6. Let $x \in X_h \cap \partial \mathbb{B}_p^R$. We denote by $\varepsilon$ the Hausdorff distance between $X$ and its Gauss digitization $X_h$. It is known that $\varepsilon = \mathcal{O}(h)$. We put $p' = \pi(p)$ and $x' = \pi(x)$. We denote by $\mathcal{C}$ the geodesic starting at the point $p'$ and passing through $x'$. We first assume that $x' \in \mathbb{B}_p^R$. In that case, the geodesic curve $\mathcal{C}$ is extended until it reaches a point $\tilde{x} \in \partial \mathbb{B}_p^R$. Note that such a curve always exist if $R$ is small enough. In the following of the proof, we denote by $\mathcal{C}_{a,b}$ the curve between any two points $a$ and $b$ and by $\ell(\mathcal{C}_{a,b})$ its length. Since the length of a curve is greater than its chord, we have

$$\ell(\mathcal{C}_{p',x'}) \geq \|p' - x'\| \geq \|p - x\| - \|x - x'\| - \|p' - p\| \geq R - 2\varepsilon.$$

If $R$ is small enough, then the line segment $[p', \tilde{x}]$ belongs to the offset $X^r$ of $X$, where $r = \text{reach}(X)$ is the reach of $X$. By using the Lipschitz property of the projection map $\pi$ onto $X$ (see Proposition 4), one gets

$$\ell(\mathcal{C}_{p', \tilde{x}}) \leq \ell(\pi([p', \tilde{x}])) \leq \frac{1}{1 - \tilde{\varepsilon}/r} \|p' - \tilde{x}\|,$$

where $\tilde{\varepsilon} = \max_{w \in [p', \tilde{x}]} \|w - \pi(w)\| = \mathcal{O}(R^2)$. Since $\|p' - \tilde{x}\| \leq \|p - \tilde{x}\| + \|p' - p\| \leq R + \varepsilon$, one gets

$$\ell(\mathcal{C}_{x', \tilde{x}}) = \ell(\mathcal{C}_{p', \tilde{x}}) - \ell(\mathcal{C}_{p', x'}) \leq (1 + \mathcal{O}(R^2))(R + \varepsilon) - R + 2\varepsilon = \mathcal{O}(R^3 + R^2\varepsilon + \varepsilon)$$

Since the curve is longer than its chord, $R = \mathcal{O}(h^\alpha)$, $\varepsilon = \mathcal{O}(h)$, we have

$$\|x' - \tilde{x}\| = \mathcal{O}(h^{3\alpha} + h^{2\alpha+1} + h) = \mathcal{O}(h^{3\alpha} + h).$$

If $x' \notin \mathbb{B}_p^R$, the same result holds with a similar proof. In that case the geodesic curve $\mathcal{C}$ does not have to be extended and $\tilde{x} = \mathcal{C} \cap \mathbb{B}_p^R$. We deduce that $\|x - \tilde{x}\| = \mathcal{O}(h^{3\alpha} + h)$, which ends the proof. $\square$

**Proof of Lemma 5.** Let us first bound $N_i$. Since $X$ is of class $C^2$, the intersection of $X$ with a ball of radius $R = Kh^\alpha$ is contained in a box of dimensions $(2R) \times (2R) \times \mathcal{O}(R^2)$. It is known that the digitized surface $X_h$ is at an Hausdorff distance less than $\frac{\sqrt{3}}{2}h$ from $X$ ([22], Theorem 1). Hence the set $B$ is
included in a set of dimensions $(2R + 2h) \times (2R + 2h) \times (O(R^2) + 2h)$. Since $\alpha \in (0, 1)$, we have $h$ in $O(h^\alpha)$. Furthermore, since $R = Kh^\alpha$, the set $B$ is thus included in a domain of volume $O(h^\alpha h^\alpha h^{\min(2\alpha, 1)})$. This implies that the number of voxels intersecting $B$ is less than $O(h^{\min(4\alpha, 2\alpha+1)-3})$. This also implies that $N_i = O(h^{\min(4\alpha, 2\alpha+1)-3})$.

Let us now bound $N_b$. Every point $x \in \partial B$ belongs to a surface that intersects $X_h \cap \partial \mathbb{B}_p^R$, so is at a distance less than $\sqrt{2h}/2$ from $X \cap \partial \mathbb{B}_p^R$. By Lemma 6, one has

$$ \partial B \subset (X \cap \partial \mathbb{B}_p^R)^{O(h+h^{3\alpha})}. $$

Since the length of $X \cap \partial \mathbb{B}_p^R$ is in $O(R)$, the volume of $(X \cap \partial \mathbb{B}_p^R)^{O(h+h^{3\alpha})}$ is $O(R(h+h^{3\alpha})^2)$. The number of boundary edges $N_b$ is of the same order than the number of cubes of size $h$ intersecting this volume so

$$ N_b = \frac{O(R(h+h^{3\alpha})^2)}{h^3} = O(h^{\alpha-1} + h^{7\alpha-3}). $$

□

We may now state our convergence result for the normalized corrected curvatures onto digitized surfaces.

**Theorem 4.** Let $V$ be a compact domain of $\mathbb{R}^3$ whose boundary $X := \partial V$ is of class $C^{3, \tau}$. Let $X_h$ be the boundary of the Gauss digitization of $V$ with step $h$. Suppose that the normal estimator satisfies $\delta_u = O(h^\beta)$ with $\beta \le 1$ and $\eta = O(h^{\frac{2}{3}})$. Let $p \in X_h$ and $B$ be the set of surfels of $X_h$ contained in the ball centered at $p$ and of radius $Kh^\alpha$ (for arbitrary $K > 0$ and $\alpha \in (0, 1)$). Then

$$ |\hat{H}_u(B) - H(\pi_X(p))| = O(h^{\gamma'}) \quad (2) $$

$$ |\hat{G}_u(B) - G(\pi_X(p))| = O(h^{\gamma'}), \quad (3) $$

where $\gamma' = \min((1 + \tau)\alpha, 2\beta - 4/3, \beta - \alpha - 1/3, 2\alpha + 2\beta - 7/3, 5\alpha + \beta - 7/3)$.

**Remark 9.** Note that for any $\beta \in ]2/3, 1[$, there exists $\alpha \in ]7/6 - \beta, \beta - 1/3[$ such that $\gamma' > 0$, which implies that we have the convergence of the pointwise mean and Gaussian curvature measures.

*Proof.* This proof relies on Corollary 2. We know that the Hausdorff $\varepsilon$ between $X$ and the boundary of its Gauss discretization $X_h$ is no greater than $\frac{\sqrt{3}}{2}h$, so $\varepsilon = O(h)$. By plugging $\varepsilon = O(h)$, $\delta_u = O(h^\beta)$, $\eta = O(h^{\frac{2}{3}})$ with $\beta \le 1$ in Corollary 2, one gets

$$ \Delta := \mu_0^{S, u}(B) - \mu_0^X(\pi_X(B)) = O\left(h^{\frac{2}{3}}\left(N_i h^{2\beta} + N_b h^\beta\right)\right). $$

Furthermore, since we are in the hypotheses of Lemma 5, we have bounds for $N_i = O(h^{4\alpha-3} + h^{2\alpha-2})$ and $N_b = O(h^{\alpha-1} + h^{7\alpha-3})$. This leads to

$$ \Delta = h^{2\alpha} O\left(h^{2\alpha+2\beta-\frac{7}{3}} + h^{2\beta-\frac{4}{3}} + h^{\beta-\alpha-\frac{1}{3}} + h^{5\alpha+\beta-\frac{7}{3}}\right) = h^{2\alpha} O(h^{\gamma}), $$
where

$$ \gamma = \min \left( 2\alpha + 2\beta - \frac{7}{3}, 2\beta - \frac{4}{3}, \beta - \alpha - \frac{1}{3}, 5\alpha + \beta - \frac{7}{3} \right) $$

Using the fact that there exists a constant C such that $\mu_0^X(\pi_X(B)) \ge Ch^{2\alpha}$, one gets

$$ \mu_0^{S, \mathbf{u}}(B) = \mu_0^X(\pi_X(B))(1 + \mathcal{O}(h^\gamma)). $$

Similarly, for $k = 1, 2$

$$ \mu_k^{S, \mathbf{u}}(B) = \mu_k^X(\pi_X(B)) + \mathcal{O}(h^{2\alpha+\gamma}). $$

Finally, we can relate our normalized mean curvature with the ratio of the mean curvature measure and the area measure:

$$ \hat{H}_{\mathbf{u}}(B) = \frac{\mu_1^{S, \mathbf{u}}(B)}{\mu_0^{S, \mathbf{u}}(B)} = \frac{\mu_1^X(\pi_X(B)) + \mathcal{O}(h^{2\alpha+\gamma})}{\mu_0^X(\pi_X(B))(1 + \mathcal{O}(h^\gamma))} = \frac{\mu_1^X(\pi_X(B))}{\mu_0^X(\pi_X(B))} + \mathcal{O}(h^\gamma). \quad (4) $$

It remains to relate this ratio of measures to the mean curvature at $p' := \pi_X(p)$. Recall that the surface is of class $C^{3,\tau}$, $\tau \le 1$. From the proof of the preceding Lemma, we know that there are two real numbers $R_1$ and $R_2$ with $D_{R_1}(p') \subset \pi_X(B) \subset D_{R_2}(p')$, where $D_r(p')$ is the geodesic disk centered at $p'$ of radius $r$, with $R_1, R_2 = \mathcal{O}(h^\alpha)$ and $R_2 - R_1 = \mathcal{O}(h + h^{3\alpha})$. Then we can write in geodesic polar coordinates

$$ \int_{D_{R_1}(p')} H dA = \int_0^{R_1} \int_0^{2\pi} H(r, \theta)r(1 + \mathcal{O}(r)) dr d\theta. $$

Using the $C^{1,\tau}$ regularity of $H$ we have, as $r \to 0$,

$$ H(r, \theta) = H(p') + ar \cos \theta + br \sin \theta + \mathcal{O}(r^{1+\tau}) $$

for some constant coefficients $a, b$. Therefore,

$$ \begin{align*} \int_{D_{R_1}(p')} H dA &= H(p') \operatorname{area}(D_{R_1}(p')) \\ &\quad + \int_0^{R_1} \int_0^{2\pi} (a \cos \theta + b \sin \theta + \mathcal{O}(r^\tau))(1 + \mathcal{O}(r)) r^2 dr d\theta \\ &= H(p') \operatorname{area}(D_{R_1}(p')) + \int_0^{R_1} \int_0^{2\pi} (a \cos \theta + b \sin \theta + \mathcal{O}(r^\tau)) r^2 dr d\theta \\ &\quad \text{because } \tau \le 1 \\ &= H(p') \operatorname{area}(D_{R_1}(p')) + \int_0^{R_1} \int_0^{2\pi} \mathcal{O}(r^{2+\tau}) dr d\theta \\ &= H(p') \operatorname{area}(D_{R_1}(p')) + \mathcal{O}(R_1^{3+\tau}) \\ &= H(p') \operatorname{area}(D_{R_1}(p')) + \mathcal{O}(h^{(3+\tau)\alpha}) \end{align*} $$
Similarly the difference between the integrals over the two disks satisfies

$$ \left| \int_{R_1}^{R_2} \int_{0}^{2\pi} (H(r, \theta) - H(p')) r(1 + \mathcal{O}(r)) dr d\theta \right| = \mathcal{O}(R_2^3 - R_1^3) = \mathcal{O}(h^{2\alpha+1} + h^{5\alpha}) $$

so that

$$ \begin{aligned} & \mu_1^X(\pi_X(B)) - H(p')\mu_0^X(\pi_X(B)) \\ &= \int_{D_{R_1}(p')} H dA - H(p') \text{area}(D_{R_1}(p')) + \mathcal{O}(h^{2\alpha+1} + h^{5\alpha}) \\ &= \mathcal{O}(h^{(3+\tau)\alpha} + h^{2\alpha+1} + h^{5\alpha}) \\ &= \mathcal{O}(h^{(3+\tau)\alpha} + h^{2\alpha+1}) \end{aligned} $$

and

$$ \begin{aligned} \mu_1^X(\pi_X(B)) &= H(p')\mu_0^X(\pi_X(B)) + \mathcal{O}(h^{(3+\tau)\alpha} + h^{2\alpha+1}) \\ &= H(p')\mu_0^X(\pi_X(B))(1 + \mathcal{O}(h^{(1+\tau)\alpha} + h)). \end{aligned} $$

Finally,

$$ \frac{\mu_1^X(\pi_X(B))}{\mu_0^X(\pi_X(B))} = H(p') + \mathcal{O}(h^{\min(1,(1+\tau)\alpha)}) \quad (5) $$

Combining (4) and (5) yields

$$ \begin{aligned} \hat{H}_{\mathbf{u}}(B) &= H(p') + \mathcal{O}(h^{\min(1,(1+\tau)\alpha,2\beta-4/3,\beta-\alpha-1/3,2\alpha+2\beta-7/3,5\alpha+\beta-7/3)}) \\ &= H(p') + \mathcal{O}(h^{\min((1+\tau)\alpha,2\beta-4/3,\beta-\alpha-1/3,2\alpha+2\beta-7/3,5\alpha+\beta-7/3)}). \end{aligned} $$

The same holds for the Gaussian curvature. □

<figure><img src="image_5.png" /></figure>

**Remark 10.** The above bound $\eta = \mathcal{O}(h^{2/3})$ has been established in [19] for the normal vector estimator based on digital integral invariants. Note that $\beta = 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 $\beta = 1$). However, we have run simulations and both Voronoi Covariance measure [24, 11] and digital Integral Invariant [19] appear to be Lipschitz normal estimator (see Figure 4).

**Remark 11.** If we choose $\mathbf{u}$ as the normals estimated by digital integral invariant normal estimator and we assume $\beta = 1$, then the previous theorem implies that, when the set $B$ of surfels is taken in a ball of center $p$ and radius $Kh^{1/3}$, then the mean corrected curvature $\hat{H}_{\mathbf{u}}(B)$ (resp. the Gaussian corrected curvature $\hat{G}_{\mathbf{u}}(B)$) tends to the mean curvature $H(\pi_X(p))$ (resp. to the Gaussian curvature $G(\pi_X(p))$) with a speed $\mathcal{O}(h^{1/3})$. In the terminology of [18], $\hat{H}_{\mathbf{u}}$ (resp. $\hat{G}_{\mathbf{u}}$) is a multigrid convergent mean (resp. Gaussian) curvature estimator.
<figure><img src="image_6.png" /><figcaption>Figure 4: Experimental evaluation of the Lipschitz property of VCM and II normal vector estimators as a function of the gridstep $h$. Asymptotic behavior is observed when reading the graph from right to left. It is observed that $\delta_u \le O(h)$ for both estimators.</figcaption></figure>

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 $Kh^{1/3}$ with a speed at least $O(h^{1/3})$, we will run our experiments with a much smaller radius $Kh^{2/3}$ (see Figures 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.

# 6 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 $\hat{H}$ and the normalized Gaussian curvature $\hat{G}$ on digital surfaces (see Definition 6). Indeed, their canonical normal vectors can
$\ell_\infty$- and $\ell_2$ errors of estimators $\hat{H}^u$ according to $\alpha$ parameter (for $u$=VCM)

<figure><img src="image_7.png" /><figcaption>Figure 5: Experimental evaluation of the asymptotic error of mean curvature along a digitized polynomial shape (“goursat”, see next section) as a function of the gridstep $h$. Convergence behavior is observed when reading the graph from right to left, where smaller values of $h$ are located. We check two different exponents $\alpha$ of the computation radius $Kh^\alpha$: $\alpha = \frac{1}{3}$ and $\alpha = \frac{1}{2}$. Thick lines represent $\ell_\infty$-error, thin lines represent $\ell_2$-error. Theoretically $\alpha = \frac{1}{3}$ should be the best choice and implies a $\ell_\infty$-error of $O(h^{\frac{1}{3}})$. In practice, the error is smaller and seems below $O(h^{\frac{1}{2}})$. Furthermore, even better results are achieved when choosing $\alpha = \frac{1}{2}$.</figcaption></figure>

only take six different values. Last, we compare our new approach with the state-of-the-art method of [7, 19], which is based on digital integral invariants. We also show that the normal cycle approach [34, 36, 8, 9] is neither accurate nor convergent for digital surfaces.

## 6.1 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.
<figure><img src="image_8.png" /><figcaption>Figure 6: Experimental evaluation of the asymptotic error of Gaussian curvature along a digitized polynomial shape ("goursat", see next section) as a function of the gridstep $h$. Convergence behavior is observed when reading the graph from right to left, where smaller values of $h$ are located. We check two different exponents $\alpha$ of the computation radius $Kh^\alpha$: $\alpha = \frac{1}{3}$ and $\alpha = \frac{1}{2}$. Thick lines represent $\ell_\infty$-error, thin lines represent $\ell_2$-error. Theoretically $\alpha = \frac{1}{3}$ should be the best choice and implies a $\ell_\infty$-error of $O(h^{\frac{1}{3}})$. In practice, the error is smaller and seems below $O(h^{\frac{1}{2}})$. Furthermore, even better results are achieved when choosing $\alpha = \frac{1}{2}$.</figcaption></figure>

**Ground-truth surfaces.** We shall test shapes whose boundary is at least twice differentiable. As a representative example, we choose the “Goursat” implicit polynomial shape $X := \{(x, y, z) \in \mathbb{R}^3, P(x, y, z) \ge 0\}$, with $P(x, y, z) := 4 - 0.015(x^4 + y^4 + z^4) + x^2 + y^2 + z^2$. Its minimal, mean and maximal mean curvatures are respectively approximately $-0.1070, 0.0956, 0.3448$. Its minimal, mean and maximal Gaussian curvatures are respectively approximately $-0.0337, 0.0080, 0.1189$. Its reach is greater than 2.9.

**Input digitized surfaces.** We digitize a shape $X$ using Gauss digitization $G_h(X) := X \cap h\mathbb{Z}^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
<figure><img src="image_9.png" /><figcaption>Figure 7: Illustration of “Goursat” polynomial surface, from left to right: smooth surface, its corresponding digitized surfaces at $h = 0.5$, $h = 0.25$, $h = 0.125$.</figcaption></figure>

forms a digital surface (e.g. see [22]). We denote it by $\partial_h X$. As illustrated on Figure 7, digitized surfaces $\partial_h X$ tends toward the smooth surface $\partial X$ in Hausdorff distance. In fact they are a Hausdorff approximation of $\partial X$ at distance less than $\frac{\sqrt{3}}{2}h$ [22]. 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 $h$.

**Input corrected normal vector field u.** We must estimate the field $\mathbf{u}$ solely from the input digitized surface $\partial_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 $\mathbf{u}$:

* Trivial Normal estimator (TN): this estimator just replicates $\mathbf{n}$ (i.e. $\mathbf{u} = \mathbf{n}$). We use it in experiments since our measures become then equivalent to the normal cycle [34, 36, 8, 9].

* Digital Integral Invariant normal estimator (II): this estimator provides a convergent normal vector field $\mathbf{u}$ for a certain family of parameters [19]. It depends on a radius of integration parameter $r := kh^\alpha$. As shown by our experiments, the values $\alpha = 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 $\mathbf{u}$ for a certain family of parameters [11, 10]. It depends on a radius of integration parameter $r := kh^\alpha$, which we set exactly as parameter $r$ of II, and on a distance of computation $R := Kh^\alpha$, where $K = 10$ gives good results.

**Estimated curvatures.** We will estimate the accuracy and stability of the following curvatures estimators:

* Normalized corrected mean curvature ($\hat{H}^\mathbf{u}$) and Gaussian curvature ($\hat{G}^\mathbf{u}$): since they are ratios of measures, we must choose a Borel set on which
measures are computed. For any point $p \in \partial_h X$, we simply compute $\hat{H}^u(p) := \hat{H}(B_\rho(p)) = \mu_1(B_\rho(p))/\mu_0(B_\rho(p))$ and $\hat{G}^u(p) := \hat{G}(B_\rho(p)) = \mu_2(B_\rho(p))/\mu_0(B_\rho(p))$, where $\rho := mh^\beta$. As shown by our experiments, the values $\beta = \frac{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^{\frac{1}{2}}$ is used).

* Digital integral invariant mean and Gaussian curvature estimators ($\hat{H}^{II}$) and ($\hat{G}^{II}$): both are parameterized by the radius $r'$ of integration. Experiments show that for $r' = k'h^{\alpha'}$, we must set $\alpha' = \frac{1}{3}$ (and no greater value) to get convergence and set $k' = 6$ to minimize estimation errors.

* Normal Cycle mean and Gaussian curvatures ($\hat{H}^{NC}$) and ($\hat{G}^{NC}$): they are defined similarly with $\hat{H}^u$ and $\hat{G}^u$ except that we use the Trivial Normal estimator to compute $u$, i.e. $\hat{H}^{NC} := \hat{H}^n$ and $\hat{G}^{NC} := \hat{G}^n$.

**Measuring $\ell_2$- and $\ell_\infty$-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 \in \partial X$ that is closest to $p$. Note that $q = \pi(p)$ since $p$ is in the reach of $\partial X$. For instance, letting $\sigma$ be the set of centroids of the surfel of $\partial_h X$, we define the errors between $\hat{H}$ and $H$ as:

$$ \ell_2(\hat{H}^u) := \sqrt{\frac{1}{|\sigma|} \sum_{p \in \sigma} (\hat{H}^u(p) - H(\pi(p)))^2} \\ \ell_\infty(\hat{H}^u) := \max_{p \in \sigma} |\hat{H}^u(p) - H(\pi(p))|. $$

Note that when $\ell_\infty$-error tends to zero as $h$ tends to zero implies the classical multigrid convergence. It also implies that $\ell_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 $\delta$ is the discrete distance of the voxel to its nearest boundary point (minimum distance is 1), then this voxel has a probability of $p^\delta$ to be flipped inside/out.

## 6.2 Evaluation of mean curvature estimator $\hat{H}^u$

**Asymptotic behavior of $\hat{H}^u$ and normal cycle $\hat{H}^{NC}$**. First we measure both the $\ell_\infty$- and $\ell_2$-errors of our proposed mean curvature estimator $\hat{H}^u$. We further test several normal estimators for $u$: II, VCM, and TN. Results are displayed on Fig. 8. Graphs show an experimental convergence speed of $O(h^{\frac{2}{3}})$, II and VCM normals. It shows also that the normal cycle method is not convergent, since $\hat{H}^u$ with TN corresponds to the normal cycle mean curvature $\hat{H}^{NC}$.
<figure><img src="image_10.png" /><figcaption>Figure 8: (Left) Asymptotic errors of mean curvature estimations by $\hat{H}^u$ along the “goursat” shape, when changing the digital normal vector estimator: II and VCM are convergent normal estimators (see text), TN is the canonical normal vector, and it corresponds to the Normal Cycle mean curvature estimator $(\hat{H}^{NC})$. (Right) Respective asymptotic errors of mean curvature estimations by $\hat{H}^u$ and $\hat{H}^{II}$ along the “goursat” shape. In both figures, $\ell_\infty$-errors are drawn with thick lines, $\ell_2$-errors are drawn with thin lines.</figcaption></figure>

**Comparative evaluation of $\hat{H}^u$ and $\hat{H}^{II}$**. We now compare the $\ell_\infty$- and $\ell_2$-errors of $\hat{H}^u$ with the ones of $\hat{H}^{II}$, which is the state-of-the-art method for digital surface curvature estimation. We use the II normal estimator for $\hat{H}^u$ with $r = 3h^{\frac{1}{2}}$. Results are shown on Fig. 8. First it confirms that $\hat{H}^{II}$ is not convergent if its parameter $r$ follows some $\Theta(h^{\frac{1}{2}})$. Secondly we do obtain a convergence speed of $\Theta(h^{\frac{1}{3}})$ for bigger radii $r = \Theta(h^{\frac{1}{3}})$. Last our new estimator $\hat{H}^u$ has a much faster convergence, approximately $\Theta(h^{\frac{2}{3}})$, despite the fact that both the integration radius $r$ for $u$ and the integration radius $\rho$ for $\hat{H}^u$ are much smaller, i.e. $3h^{\frac{1}{2}}$.

We further illustrate the differences of the two estimators $\hat{H}^u$ and $\hat{H}^{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 $\hat{H}^u$, we measure the $\ell_\infty$-error when changing $k$ in the radius $r := kh^{\frac{1}{2}}$ and changing $m$ in the radius $\rho := mh^{\frac{1}{2}}$. For $\hat{H}^{II}$ we simply change $k$ in the radius $r := kh^{\frac{1}{3}}$. See Figure 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) \in [2.5, 4] \times [2.5, 6]$.
<figure><img src="image_11.png" /><figcaption>Figure 9: Illustration of the accuracy of mean curvature estimators $\hat{H}^u$ and $\hat{H}^{II}$ at different resolutions: 1st row is ground-truth $H$, 2nd row is our approach $\hat{H}^u$, 3rd row is its local estimation error $|\hat{H}^u - H|$, 4th row is integral invariant approach $\hat{H}^{II}$, 5th row is its local estimation error $|\hat{H}^{II} - H|$. Curvatures are displayed with colors from dark blue (-0.3) to white to red (0.3), except a black band at [-0.01; 0.01]. Errors are displayed with colors from white (0) to red (0.025) to black (0.05).</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 10: Stability of $\hat{H}^u$ and $\hat{H}^{II}$ with respect to parameters and gridstep. Estimator $\hat{H}^u$ is parameterized by the integration radius $r := kh^{1/2}$ of its normal estimator $\mathbf{u}=\mathrm{II}$ and by the measure radius $\rho := mh^{1/2}$. We plot $\ell_2$- and $\ell_\infty$-errors of $\hat{H}^u - H$ as a function of $k$ and $m$. Estimator $\hat{H}^{II}$ is parameterized by the integration radius $r' := k'h^{1/3}$. We plot $\ell_2$- and $\ell_\infty$-errors of $\hat{H}^{II} - H$ as a function of $k'$ (displayed as 2d plot to make easier the comparison).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 11: Comparison of mean curvature estimators $\hat{H}^u$, $\hat{H}^{NC}$ and $\hat{H}^{II}$ on classical digital shapes at two different resolutions. We choose parameters according to section 6.1, setting $h = 1$ for “dragon 256³” and naturally $h = 0.5$ for “dragon 512³”. Curvatures are displayed with colors from dark blue (-0.3) to white to red (0.3), except a black band at [-0.01; 0.01].</figcaption></figure>

**Evaluation of $\hat{H}^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 on Figure 11 and Figure 12. It is clear that both $\hat{H}^{NC}$ and $\hat{H}^{II}$ oscillate around the correct solution (see black random or *moiré* patterns on both figures, whichever the resolution). On the contrary, $\hat{H}^u$ is stable in zero=curvature regions (like on octaflower) while accurately delineating the small scales of the chinese dragon.

**Robustness to noise of $\hat{H}^u$.** Last we have checked the robustness to noise of our mean curvature estimator $\hat{H}^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
<figure><img src="image_14.png" /><figcaption>Figure 12: Comparison of mean curvature estimators $\hat{H}^u$, $\hat{H}^{NC}$ and $\hat{H}^{II}$ on classical digital shapes at two different resolutions. We choose parameters according to section 6.1, setting $h = 0.5$ for “octaflower 256³” and naturally $h = 0.25$ for “octaflower 512³”. Curvatures are displayed with colors from dark blue (-0.2) to white to red (0.2), except a black band at [-0.005; 0.005].</figcaption></figure>
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].

6.3 Evaluation of Gaussian curvature estimator $\hat{G}^u$

Asymptotic behavior of $\hat{G}^u$ and normal cycle $\hat{G}^{NC}$. Then we measure both the $\ell_\infty$- and $\ell_2$-errors of our proposed Gaussian curvature estimator $\hat{G}^u$. We further test several normal estimators for $u$: II, VCM, and TN. Results are displayed on Fig. 14. Graphs show again an experimental convergence speed of $O(h^{3/2})$ for II and VCM normals. It shows also that the normal cycle method is not convergent, since $\hat{G}^u$ with TN corresponds to the normal cycle Gaussian curvature $\hat{G}^{NC}$.

Comparative evaluation of $\hat{G}^u$ and $\hat{G}^{II}$. We compare the $\ell_\infty$- and $\ell_2$-errors of $\hat{G}^u$ with the ones of $\hat{G}^{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 $\hat{G}^u$ with respect to $\hat{G}^{II}$: much faster convergence speed with small integration radii.

We further illustrate the differences of the two estimators $\hat{G}^u$ and $\hat{G}^{II}$ by displaying the estimated Gauss curvatures and the localization of errors on several digitizations of “goursat” on Fig. 15. 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.

**Evaluation of $\hat{G}^u$ on digital shapes** We also check Gaussian curvature estimators $\hat{G}^u$, $\hat{G}^{NC}$ and $\hat{G}^{II}$ on the same classical shapes as in the previous subsection. Results are displayed on Figure 16 and Figure 17; keeping the same parameters as above. We also observe the oscillations and *moiré* patterns in $\hat{H}^{II}$ estimations. Furthermore the Normal Cycle estimator $\hat{G}^{NC}$ is clearly incorrect and gives only extremal results. This is because $\hat{G}^{NC}$ takes into account solely six possible normals.

**Robustness to noise of $\hat{G}^u$.** Last we have checked the robustness to noise of our Gaussian curvature estimator $\hat{G}^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].
<figure><img src="image_15.png" /><figcaption>Figure 13: Robustness to noise of our mean curvature estimator $\hat{H}^u$. Noise with parameter $p$ (see text) was added to the digital shape “octaflower 256³”. The gridstep was set to $h = 0.5$ for parameterizing the estimators. Curvatures are displayed with colors from dark blue (-0.2) to white to red (0.2), except a black band at [-0.005; 0.005].</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 14: (Left) Asymptotic errors of Gaussian curvature estimations by $\hat{G}^u$ along the “goursat” shape, when changing the digital normal vector estimator: II and VCM are convergent normal estimators (see text), TN is the canonical normal vector, and it corresponds to the Normal Cycle Gaussian curvature estimator ($\hat{G}^{NC}$). (Right) Respective asymptotic errors of Gaussian curvature estimations by $\hat{G}^u$ and $\hat{G}^{II}$ along the “goursat” shape. In both figures, $\ell_\infty$-errors are drawn with thick lines, $\ell_2$-errors are drawn with thin lines.</figcaption></figure>

## 6.4 Practical computation times of our curvature estimators.

We note first that the computation times of the measures $\mu_0, \mu_1, \mu_2, \mu_\Omega$ and $\mu_{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 $\rho$ (2% when $N \approx 1e^3$, 0.01% when $N \approx 4e^6$). Secondly computation times for Gaussian curvature estimator $\hat{G}^u$ or principal directions is almost the same as the one of $\hat{H}^u$ and $\hat{G}^{II}$, and are thus not displayed.

We have plotted the measured computation times of our mean curvature estimator $\hat{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, 19]. We have not compared with the execution times of Normal Cycle method, since $\hat{H}^{NC}$ is not accurate. Note that using VCM normal estimator instead of II normal estimator to compute $\mathbf{u}$ is faster for large digital shapes, but requires more memory.

To sum up, our method has an experimental complexity in $\Theta(n^{3/2})$ compared to Integral Invariant method that has an experimental complexity in $\Theta(n^2)$. Our code could be further optimized in the integration step of the measures $\mu_k$ by reusing the results of a neighboring point, as it is already done in Integral
<figure><img src="image_17.png" /><figcaption>Figure 15: Illustration of the accuracy of Gaussian curvature estimators $\hat{G}^u$ and $\hat{G}^{II}$ at different resolutions: 1st row is ground-truth $G$, 2nd row is our approach $\hat{G}^u$, 3rd row is its local estimation error $|\hat{G}^u - G|$, 4th row is integral invariant approach $\hat{G}^{II}$, 5th row is its local estimation error $|\hat{G}^{II} - G|$. Curvatures are displayed with colors from dark blue (-0.1) to white to red (0.1), except a black band at [-0.002; 0.002]. Errors are displayed with colors from white (0) to red (0.0125) to black (0.025).</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure 16: Comparison of Gaussian curvature estimators $\hat{G}^u$, $\hat{G}^{NC}$ and $\hat{G}^{II}$ on classical digital shapes at two different resolutions. We choose parameters according to section 6.1, setting $h = 1$ for “dragon 256³” and naturally $h = 0.5$ for “dragon 512³”. Curvatures are displayed with colors from dark blue (-0.05) to white to red (0.05), except a black band at [-0.0005; 0.0005].</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure 17: Comparison of Gaussian curvature estimators $\hat{G}^u$, $\hat{G}^{NC}$ and $\hat{G}^{II}$ on classical digital shapes at two different resolutions. We choose parameters according to section 6.1, setting $h = 0.5$ for “octaflower 256³" and naturally $h = 0.25$ for “octaflower 512³". Curvatures are displayed with colors from dark blue (-0.05) to white to red (0.05), except a black band at [-0.0005; 0.0005].</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure 18: Robustness to noise of our Gaussian curvature estimator $\hat{G}^u$. Noise with parameter $p$ (see text) was added to the digital shape “octaflower 256³”. The gridstep was set to $h = 0.5$ for parameterizing the estimators. Curvatures are displayed with colors from dark blue (-0.05) to white to red (0.05), except a black band at [-0.0005; 0.0005].</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure 19: Computation times of curvature estimators $\hat{H}^u$ and $\hat{H}^{II}$ as a function of the number of surfels of the input digitized boundary (left) or as a function of the $l_2$-error (right).</figcaption></figure>

Invariant method. Last, for a given accuracy, our method is even much faster (5000 times faster than II to achieve 0.0065 $l_2$-accuracy).

# 7 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 Section 4. These works will be the focus of a forthcoming paper. Preliminary results are displayed on Figure 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.

# References

[1] William K Allard, *On the first variation of a varifold*, Annals of mathematics (1972), 417–491.
<figure><img src="image_22.png" /><figcaption>Figure 20: Computation of principal directions (represented as sticks) and principal curvatures (same colormap as in "goursat" figures) with a corrected anisotropic curvature measure. Top row: first (left) and second (right) principal directions with **u** = VCM. Second row: expected first (left) and second (right) principal directions.</figcaption></figure>
[2] Frederick J Almgren, *Plateau's problem: an invitation to varifold geometry*, vol. 13, American Mathematical Soc., 1966.

[3] Blanche Buet, Gian Paolo Leonardi, and Simon Masnou, *A varifold approach to surface approximation*, Arch. Ration. Mech. Anal. **226** (2017), no. 2, 639–694 (English).

[4] F. Cazals and M. Pouget, *Estimating differential quantities using polynomial fitting of osculating jets*, Computer Aided Geometric Design **22** (2005), no. 2, 121–146.

[5] Frédéric Chazal, David Cohen-Steiner, André Lieutier, and Boris Thibert, *Stability of curvature measures*, Computer Graphics Forum **28** (2009), no. 5, 1485–1496.

[6] D. Coeurjolly, M. Foare, P. Gueth, and J.-O. Lachaud, *Piecewise smooth reconstruction of normal vector field on digital data*, Comput. Graph. Forum **35** (2016), no. 7, 157–167, Proc. of Pacific Graphics 2016.

[7] D. Coeurjolly, J.-O. Lachaud, and J. Levallois, *Multigrid convergent principal curvature estimators in digital geometry*, Computer Vision and Image Understanding **129** (2014), 27 – 41, Special section: Advances in Discrete Geometry for Computer Imagery.

[8] David Cohen-Steiner and Jean-Marie Morvan, *Restricted delaunay triangulations and normal cycle*, Proceedings of the nineteenth annual symposium on Computational geometry, ACM, 2003, pp. 312–321.

[9] —, *Second fundamental measure of geometric sets and local approximation of curvatures*, Journal of Differential Geometry **74** (2006), no. 3, 363–394.

[10] L. Cuel, J.-O. Lachaud, Q. Mérigot, and B. Thibert, *Robust geometry estimation using the generalized voronoi covariance measure*, SIAM Journal on Imaging Sciences **8** (2015), no. 2, 1293–1314.

[11] L. Cuel, J.-O. Lachaud, and B. Thibert, *Voronoi-based geometry estimator for 3d digital surfaces*, Proc. Int. Conf. on Discrete Geometry for Computer Imagery (DGCI’2014), Sienna, Italy (Elena Barcucci, Andrea Frosini, and Simone Rinaldi, eds.), Lecture Notes in Computer Science, vol. 8668, Springer International Publishing, 2014, pp. 134–149 (English).

[12] Herbert Federer, *Curvature measures*, Transactions of the American Mathematical Society **93** (1959), no. 3, 418–491.

[13] —, *Geometric measure theory*, Springer, 2014.

[14] Joseph HG Fu, *Curvature measures of subanalytic sets*, American Journal of Mathematics (1994), 819–880.
[15] Joseph HG Fu et al., *Convergence of curvatures in secant approximations*, Journal of Differential Geometry **37** (1993), no. 1, 177–190.

[16] Rocio Gonzalez-Diaz, Maria-Jose Jimenez, and Belen Medrano, *3d well-composed polyhedral complexes*, Discrete Applied Mathematics **183** (2015), 59–77.

[17] Klaus Hildebrandt, Konrad Polthier, and Max Wardetzky, *On the convergence of metric and geometric properties of polyhedral surfaces*, Geometriae Dedicata **123** (2006), no. 1, 89–112.

[18] Reinhard Klette and Azriel Rosenfeld, *Digital geometry: Geometric methods for digital picture analysis*, Elsevier, 2004.

[19] J.-O. Lachaud, D. Coeurjolly, and J. Levallois, *Robust and convergent curvature and normal estimators with digital integral invariants*, Modern Approaches to Discrete Curvature (L. Najman and P. Romon, eds.), Lecture Notes in Mathematics, vol. 2184, Springer International Publishing, Cham, 2017, pp. 293–348.

[20] J.-O. Lachaud, X. Provençal, and T. Roussillon, *An output-sensitive algorithm to compute the normal vector of a digital plane*, Theor. Comput. Sci. **624** (2016), 73–88.

[21] __________, *Two plane-probing algorithms for the computation of the normal vector to a digital plane*, Journal of Mathematical Imaging and Vision **59** (2017), no. 1, 23–39.

[22] J.-O. Lachaud and B. Thibert, *Properties of gauss digitized shapes and digital surface integration*, Journal of Mathematical Imaging and Vision **54** (2016), no. 2, 162–180 (English).

[23] Fang-Hua Lin and Xiaoping Yang, *Geometric measure theory: an introduction*, Science Press, 2002.

[24] Q. Mérigot, M. Ovsjanikov, and L. Guibas, *Voronoi-based curvature and feature estimation from point clouds*, Visualization and Computer Graphics, IEEE Transactions on **17** (2011), no. 6, 743–756.

[25] Frank Morgan, *Geometric measure theory: a beginner's guide*, Academic press, 2016.

[26] Jean-Marie Morvan, *Generalized curvatures*, Springer, 2008.

[27] Jean-Marie Morvan and Boris Thibert, *Approximation of the normal vector field and the area of a smooth surface*, Discrete & Computational Geometry **32** (2004), no. 3, 383–400.

[28] H. Pottmann, J. Wallner, Q. Huang, and Y. Yang, *Integral invariants for robust geometry processing*, Computer Aided Geometric Design **26** (2009), no. 1, 37–60.
[29] H. Pottmann, J. Wallner, Y. Yang, Y. Lai, and S. Hu, *Principal curvatures from the integral invariant viewpoint*, Computer Aided Geometric Design **24** (2007), no. 8-9, 428-442.

[30] Marcelo Siqueira, Longin Jan Latecki, Nicholas Tustison, Jean Gallier, and James Gee, *Topological repairing of 3d digital images*, Journal of Mathematical Imaging and Vision **30** (2008), no. 3, 249-274.

[31] Michael Spivak, *Calculus on manifolds: a modern approach to classical theorems of advanced calculus*, CRC Press, 2018.

[32] Peer Stelldinger, Longin Jan Latecki, and Marcelo Siqueira, *Topological equivalence between a 3d object and the reconstruction of its digital image*, IEEE transactions on pattern analysis and machine intelligence **29** (2007), no. 1, 126-140.

[33] T. Surazhsky, E. Magid, O. Soldea, G. Elber, and E. Rivlin, *A comparison of gaussian and mean curvatures estimation methods on triangular meshes*, Robotics and Automation, 2003. Proceedings. ICRA '03. IEEE International Conference on, vol. 1, 2003, pp. 1021-1026.

[34] Peter Wintgen, *Normal cycle and integral curvature for polyhedra in riemannian manifolds*, Differential geometry **21** (1982).

[35] G. Xu, *Convergence analysis of a discretization scheme for gaussian curvature over triangular surfaces*, Computer Aided Geometric Design **23** (2006), no. 2, 193-207.

[36] Martina Zähle, *Curvatures and currents for unions of sets with positive reach*, Geometriae Dedicata **23** (1987), no. 2, 155-171.

J.-O. Lachaud

Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS

Laboratoire de Mathématiques jacques-olivier.lachaud@univ-smb.fr

P. Romon

Université Paris-Est

Laboratoire d'Analyse et de Mathématiques Appliquées

Université Paris-Est Marne-la-Vallée, Université Paris-Est Créteil, CNRS

pascal.romon@upem.fr

B. Thibert

Université Grenoble Alpes

Laboratoire Jean Kuntzman

boris.thibert@imag.fr