Boris Thibert

PROPERTIES OF GAUSS DIGITIZED SHAPES AND DIGITAL SURFACE INTEGRATION

This paper presents new topological and geometric properties of Gauss digitizations of Euclidean shapes, most of them holding in arbitrary dimension d. We focus on r-regular shapes sampled by Gauss digitization at gridstep h. The digitized boundary is shown to be close to the Euclidean boundary in the Hausdorff sense, the minimum distance

2 h being achieved by the projection map ξ induced by the Euclidean distance. Although it is known that Gauss digitized boundaries may not be manifold when d ≥ 3, we show that non-manifoldness may only occur in places where the normal vector is almost aligned with some digitization axis, and the limit angle decreases with h. We then have a closer look at the projection of the digitized boundary onto the continuous boundary by ξ. We show that the size of its non-injective part tends to zero with h. This leads us to study the classical digital surface integration scheme, which allocates a measure to each surface element that is proportional to the cosine of the angle between an estimated normal vector and the trivial surface element normal vector. We show that digital integration is convergent whenever the normal estimator is multigrid convergent, and we explicit the convergence speed. Since convergent estimators are now available in the literature, digital integration provides a convergent measure for digitized objects. KeywordsGauss digitization and geometric inference and digital integral and multigrid convergence and set with positive reach.

Introduction

Understanding what are the properties of real objects that can be extracted from their digital representation is a crucial task in knowledge representation and processing. In most applications, a real object or a scene is known only through some discrete finite representation, generally a digital image produced by some complex system, involving acquisition, sampling, quantization, and processing. This process is often called digitization or sampling and is realized by devices like CCD or CMOS cameras, document scanners, CT or MRI scanners. Since the digitization process aims to be as faithful as possible to the real data, it is very natural to look at topological and geometric properties that can be inferred from digital data for rather elementary digitization processes and classes of real Euclidean objects.

This paper focuses on several global and local topological and geometric properties that are preserved by Gauss digitization.

Definition 1 (Gauss digitization). Let h > 0 be a sampling grid step. The Gauss digitization of an Euclidean shape X ⊂ R d is defined as D h (X) := X ∩ (hZ) d (see Fig. 1).

It is thus one of the simplest conceivable digitization scheme. We study here more specifically the local connections between the boundary ∂X of the Euclidean shape and the boundary ∂ h X of its digitization (as an union of d -1-dimensional cubic faces, see below). It is clear that one cannot expect that many properties of real shapes be preserved by digitization for arbitrary digitization step h, just by some combinatorial argument. However one can expect that, as the grid step gets finer (h converges to 0), we can recover most of the properties of the real shape from its digitization. Indeed, the literature shows that topological properties may be preserved for fine enough digitization grids

X ∂X D h X h Q h D h X ∂ h X h Figure 1.
Illustration of the main definitions used troughout the paper: a shape X, its boundary ∂X, its Gauss digitization D h X, its h-cube embedding Q h D h X and its h-boundary ∂ h X.

for specific class of shapes, at least in dimension 2. For geometric properties, their "preservation" is rephrased in terms of accuracy of estimation. Thus, given some geometric estimator, the estimated quantity should tend towards the geometric quantity of the Euclidean shape as the digitization gridstep gets finer. The estimator is then said to be multigrid convergent with a speed depending on how the estimation error approaches zero. The objective of many works in the literature is to define geometric estimators and to prove their multigrid convergence. We review first previous works on topological and geometric properties inferred from digital data before describing our contributions in more details. Topological properties of digitizations. The inference of topological properties has been extensively studied especially in the 2D case, mainly with morphological tools. We may quote the seminal works of Pavlidis [START_REF] Pavlidis | Algorithms for graphics and image processing[END_REF] and Serra [START_REF] Serra | Image analysis and mathematical morphology[END_REF] who established the first homeomorphism theorems for sufficiently smooth shapes digitized on a square or hexagonal grid with Gauss digitization, provided the grid step is fine enough. A key ingredient for topology preservation independently discovered in their works was the R-regularity, later called par(R)-regularity. These results were extended to several other digitization schemes (square subset, intersection, v-area) by Latecki et al. [START_REF] Gross | Digitizations preserving topological and differential geometric properties[END_REF][START_REF] Latecki | Preserving topology by a digitization process[END_REF]. Along the same lines, a global digitization scheme called Hausdorff discretization was proposed in [START_REF] Ronse | Discretization in hausdorff space[END_REF][START_REF] Tajine | Topological properties of hausdorff discretization, and comparison to other discretization schemes[END_REF]. It was shown that connectivity is preserved by this scheme. Finally, Stelldinger and Köthe [START_REF] Stelldinger | Towards a general sampling theory for shape preservation[END_REF] achieves very general topology preservation theorems for arbitrary sampling grids, that applies not only to Gauss digitization but also to convolutions by a point spread function. It is worthy to note that their theorems are general enough to include most previously known homeomorphisms results [START_REF] Pavlidis | Algorithms for graphics and image processing[END_REF][START_REF] Serra | Image analysis and mathematical morphology[END_REF][START_REF] Gross | Digitizations preserving topological and differential geometric properties[END_REF][START_REF] Latecki | Preserving topology by a digitization process[END_REF].

Extending previous results to non R-regular shapes appears quite challenging. Giraldo et al. showed that finite polyhedra can be digitized such that the homotopy type is preserved [START_REF] Giraldo | Digitizations preserving shape[END_REF]. The more flexible R-stability property (a shape and its R-offset have same homotopy type) was proposed in [START_REF] Meine | A topological sampling theorem for robust boundary reconstruction and image segmentation[END_REF]. This approach allows topological stability even for plane partitions. Fewer works address the case of d-dimensional images, for d ≥ 3. One underlying reason is that topology preservation cannot be achieved in general already for d = 3. It is indeed easy to construct smooth sets, but with bad digitization at some arbitrary small step h. For instance, Stelldinger and Köthe ( [START_REF] Stelldinger | Towards a general sampling theory for shape preservation[END_REF], Theorem 3) exhibits a cylinder of radius R, the axis of which is aligned with the straight line z = 0 and x = y, and it contains the point (0, 0, -R), where 0 < h. The cylinder extremities are smoothed as spheres. Even for small h, its Gauss digitization induces a non-manifold digitized boundary. Worse, this issue arises for all classical digitization schemes. However, they show that objects keep identical homotopy tree through Gauss digitization ( [START_REF] Stelldinger | Towards a general sampling theory for shape preservation[END_REF], Theorem 1).

Several routes for solving the homeomorphism issue were proposed by Stelldinger et al. [START_REF] Stelldinger | Topological equivalence between a 3d object and the reconstruction of its digital image[END_REF]. A first idea is to refine the digitized object on a twice finer grid by majority interpolation, and this leads to a manifold digital surface close to the real object boundary ( Theorem 19). They also propose to recontruct from the digitized object an approximate surface based either on a union of ball, a modification of a marching-cubes algorithm, or a smoothing of the latter surface. Homeomorphism is achieved in all cases for fine enough grids. The projection map. Another path to handle topological or geometric inference problems is a more functional approach: the distance function to a shape and the associated projection map. It is a key tool since it encodes information on the shape and it is Hausdorff stable, whatever the dimension. The distance function of a compact set K is defined on R d by d K (x) := min{ x -y , y ∈ K}. The R-offset of K, denoted by K R , is the set whose points x satisfy d K (x) ≤ R. The medial axis MA(K) of K is the subset of R d whose points have at least two closest points on K.

Definition 2 (Projection map). The projection map of a compact set K is the map

ξ K : R d \ MA(K) → K that associates to any point x of R d \ MA(K) its unique closest point onto K.
The reach of K, denoted by reach(K), is the infimum of {d K (y), y ∈ MA(K)} [START_REF] Federer | Curvature measures[END_REF]. The projection map ξ K of a compact set K with positive reach is a useful tool because it allows to compare K with another shape lying in its neighborhood.

Note that the R-offsets allow to recover stable topological (and geometric) properties. If the shape K has positive reach and if a point cloud P is dense enough around K, then for some suitable values of R, the R-offsets of P are homotopy equivalent to K [START_REF] Chazal | Smooth manifold reconstruction from noisy and non-uniform approximation with guarantees[END_REF][START_REF] Niyogi | Finding the homology of submanifolds with high confidence from random samples[END_REF]. This result has been extended for digitizations close to K in the Hausdorff sense. They are shown to be homotopy equivalent to K, for suitable values of digitization step size [START_REF] Attali | Reconstructing shapes with guarantees by unions of convex sets[END_REF]. Global geometric properties of digitizations. Infering geometric properties of Euclidean object from their digitization has a long history. Until recently, most research efforts focused on global geometric properties. For instance, The area (in 2D) or volume (in 3D) may indeed be estimated just by counting the number of digital samples and this fact was known by Gauss and Dirichlet as reported for instance in [START_REF] Klette | Digital Geometry: Geometric Methods for Digital Picture Analysis[END_REF]. Further results show that volumes and also moments may be estimated by appropriate counting with even superlinear convergence for smooth enough classes of shapes [START_REF] Huxley | Exponential sums and lattice points[END_REF][START_REF] Klette | Multigrid convergence of calculated features in image analysis[END_REF].

It is harder to define length/perimeter estimators in 2D or area estimators in 3D with proven convergence. For length/perimeter, for specific classes of shapes, several approaches offer guarantees like segmentation into digital straight segments [START_REF] Klette | Multigrid convergence of calculated features in image analysis[END_REF], -sausage approach [START_REF] Klette | Digital Geometry: Geometric Methods for Digital Picture Analysis[END_REF], and minimum length polygon [START_REF] Sloboda | On piecewise linear approximation of planar Jordan curves[END_REF]. A more local approach based on tangent estimation and integration leads also to multigrid convergence with speed O(h 1 3 ) [START_REF] Lachaud | Fast, accurate and convergent tangent estimation on digital contours[END_REF][START_REF] Lachaud | Espaces non-euclidiens et analyse d'image : modèles déformables riemanniens et discrets, topologie et géométrie discrète[END_REF]. Few results exist for 3D area estimation. Most approaches try to assign weights to local configurations in order to minimize the maximal error [START_REF] Lindblad | Surface area estimation of digitized 3d objects using weighted local configurations[END_REF][START_REF] Ziegel | Estimation of surface area and surface area measure of three-dimensional sets from digitizations[END_REF], but such approaches cannot achieve multigrid convergence [START_REF] Kenmochi | Surface area estimation for digitized regular solids[END_REF]. Polyhedrization with digital planes for area estimation [START_REF] Klette | Digital planar segment based polyhedrization for surface area estimation[END_REF] is an interesting extension of 2D methods, but no theoretical guarantees have been established.

Finally three methods offer (some) theoretical guarantees. Area estimation by integration of normals, first proposed in [START_REF] Lenoir | Fast computation of the normal vector field of the surface of a 3-d discrete object[END_REF] and more formalized in [START_REF] Coeurjolly | Multigrid convergence and surface area estimation[END_REF], has the advantage of defining an elementary area measure, which in turn can provide the global area measure but may also be used for integration of other quantities. However their results rely on assumptions that are not satisfied by the Gauss digitization boundary. A second approach estimates the volume of an appropriate thickened version of the surface, and deduced the area [START_REF] Stelldinger | Topological equivalence between a 3d object and the reconstruction of its digital image[END_REF]. Their algorithm is not applicable as is on data since it requires to loop over finer and finer digitizations of the continuous object. Besides it is in fact very similar to Steiner tube formula dating from 1840. A third approach relies on Cauchy-Crofton integral formula and estimates area by statistical intersection of the volume with lines [START_REF] Liu | Surface area estimation of digitized 3d objects using quasi-monte carlo methods[END_REF]. It it important to note that all three methods do not provide an error bound. The speed of convergence of these estimators is thus unknown, even for specific classes of shapes.

Local geometric properties. It is often interesting to estimate more local geometric quantities like normal vector or tangent planes, curvatures or principal directions. Since accuracy is ambiguous at a given sampling, the definition of multigrid convergence is adapted to local geometric quantities (e.g., see [START_REF] Coeurjolly | Multigrid convergence of discrete geometric estimators[END_REF]). Several estimators are multigrid convergent: (i) digital straight segment recognition defines parameter-free convergent estimators of normal/tangent in 2D [START_REF] De Vieilleville | Maximal digital straight segments and convergence of discrete geometric estimators[END_REF][START_REF] Lachaud | Fast, accurate and convergent tangent estimation on digital contours[END_REF][START_REF] Lachaud | Espaces non-euclidiens et analyse d'image : modèles déformables riemanniens et discrets, topologie et géométrie discrète[END_REF], (ii) polynomial fitting induces convergent estimators of derivatives of any order in 2D [START_REF] Provot | Estimation of the derivatives of a digital function with a convergent bounded error[END_REF], (iii) binomial convolution leads also to convergent estimators of derivatives in 2D [START_REF] Esbelin | Convergence of binomial-based derivative estimation for c2-noisy discretized curves[END_REF][START_REF] Esbelin | Convergence of binomial-based derivative estimation for c 2 noisy discretized curves[END_REF]; (iv) the recently introduced integral invariants define multigrid convergent estimators of normals, 2D curvature, mean curvature [START_REF] Coeurjolly | Integral based curvature estimators in digital geometry[END_REF], and also 3D principal curvatures and principal directions [START_REF] Coeurjolly | Multigrid convergent principal curvature estimators in digital geometry[END_REF].

Note that the distance function to a shape and its projection map also encode information on the normals and on the curvatures. If K is a shape with positive µ-reach (a much less restrictive condition than positive reach), then the offset of point cloud approximating K provides estimation of the normals [START_REF] Chazal | Normal cone approximation and offset shape isotopy[END_REF] and of the curvature measures [START_REF] Chazal | Stability of curvature measures[END_REF][START_REF] Chazal | Boundary measures for geometric inference[END_REF] of K at a given scale. Voronoi covariance measure [START_REF] Mérigot | Voronoi-based curvature and feature estimation from point clouds[END_REF] may also be adapted to digital data to define multigrid convergent normal estimators in arbitrary dimension [START_REF] Cuel | Voronoi-based geometry estimator for 3d digital surfaces[END_REF]. Contributions. We establish both topological preservation and multigrid convergence results. After recalling useful notations and definitions in Section 2, Section 3 establishes elementary results on Gauss digitized sets. We connect in Lemma 1 two notions: the R-regularity of shapes known in digital geometry [START_REF] Gross | Digitizations preserving topological and differential geometric properties[END_REF][START_REF] Latecki | Preserving topology by a digitization process[END_REF][START_REF] Pavlidis | Algorithms for graphics and image processing[END_REF][START_REF] Serra | Image analysis and mathematical morphology[END_REF] and the reach of compact sets known in geometric measure theory [START_REF] Federer | Curvature measures[END_REF] and computational geometry. Such shapes have a good behaviour with respect to digitization. Then we establish that ∂X and ∂ h X are close to each other whatever the dimension (Theorem 1). This proximity is realized by the projection ξ of ∂ h X onto ∂X induced by the Euclidean distance.

We then address the homeomorphism problem between these two sets, which is caused by the possible non-manifoldness of the digitized boundary [START_REF] Stelldinger | Towards a general sampling theory for shape preservation[END_REF]. Although this problem is unavoidable starting from dimension 3, it is worth studying where non-manifoldness arises and if it is likely to arise often. With this information, it is then easier to take them into account, for instance to correct the digital dataset [START_REF] Siqueira | Topological repairing of 3d digital images[END_REF]. In Section 4, we show local sufficient conditions which guarantee that the digitized boundary is a manifold at this location (Theorem 2). They indicate that both sets ∂X and ∂ h X are "almost" homeomorphic, and that the area of non-homeomorphic places reduces generally toward 0 as the gridstep h gets finer and is reduced to 0. Furthermore, only places of ∂X with a normal very close to some axis direction may induce a non-manifold place in ∂ h X. This fact is illustrated on Fig. 2 as parts painted in dark grey on digitized boundary. Hence our approach is very different from the one of Stelldinger et al. [START_REF] Stelldinger | Topological equivalence between a 3d object and the reconstruction of its digital image[END_REF]. Instead of building a digitized surface different from the Gauss digitized boundary to get a homeomorphism, we characterize the rare places where the Gauss digitized boundary may not be a manifold.

Afterwards we establish in Section 5 several results related to the projection map between ∂X and ∂ h X. Even for smooth convex shapes, the projection map is not everywhere injective. However Theorem 3 shows that the size of the non-injective part on ∂X decreases linearly in h. Fig. 2 shows in light grey places where projection ξ might not be injective. Obviously, it includes zones in dark grey where the digitized boundary is not even a manifold.

Finally, using results from geometric measure theory, Section 6 shows the conditions under which digital integration on ∂ h X is multigrid convergent toward integration on ∂X, for an arbitrary integrable function from R d to R. Given some digital normal estimator, digital integration is defined as proposed in [START_REF] Lenoir | Fast computation of the normal vector field of the surface of a 3-d discrete object[END_REF][START_REF] Coeurjolly | Multigrid convergence and surface area estimation[END_REF] by summation over digital d -1-cells of the function value weighted by the inner product between trivial and estimated normal (see Definition 6). Theorem 4 demonstrates that digital integration is multigrid convergent toward usual integration as long as the normal estimator is multigrid convergent. The convergence speed is also fully explicited and is upper bounded on the one hand by the convergence speed of the normal estimator and on the other hand by the gridstep h. Since multigrid convergent normal estimators exist in arbitrary dimension [START_REF] Lachaud | Fast, accurate and convergent tangent estimation on digital contours[END_REF][START_REF] Coeurjolly | Multigrid convergent principal curvature estimators in digital geometry[END_REF][START_REF] Cuel | Voronoi-based geometry estimator for 3d digital surfaces[END_REF], our theorem proves that both local and global area estimation by digital integration is multigrid convergent, and it gives a well-defined measure on digitized boundary. ; their relative area is denoted by A nm . Zones in light grey (and dark grey) indicates the surface parts where projection ξ might not be an homeomorphism (Theorem 3); their relative area is denoted by A ξ . Clearly, both zones tends to area zero as the gridstep gets finer and finer, while parts where digitization might not be manifold are much smaller than parts where ξ might not be homeomorphic.

Preliminary notions and definitions

Given a compact shape X ⊂ R d , we wish to compare the topological boundary of X, denoted by ∂X, with the boundary of its Gauss digitization. As defined in the introduction, the Gauss digitization of X is a regular sampling of the characteristic function of X, with a parameterized sampling density h. Digitized sets are defined as subsets of (hZ) d . Since they have peculiar coordinates (multiple of h), points of such subsets will be called digital points. In order to define a digitized boundary, we have to see the digitized set as a union of cubes with edge length h. For some z ∈ (hZ) d , the closed d-dimensional axis-aligned cube of R d centered on z with edge length h is denoted by Q h z and called h-cube. The h-cube embedding of a digital set Z is naturally defined as

Q h Z := ∪ z∈Z Q h z .
Figure 3. Illustration of (left) the primal cubical grid F h in gray, (middle) the dual cubical grid Fh in cyan, (right) the Gauss digitization D h X in black and the h-boundary of X in red, seen as a subcomplex of F h . Note that some cells and their dual are shown between the two grids.

Definition 3 (h-boundary of X). The h-boundary of X, denoted by ∂ h X, is the topological boundary of the h-cube embedding of the Gauss digitization of X:

(1)

∂ h X := ∂ ∪ z∈D h X Q h z , where D h X is given in Definition 1.
The h-boundary of X is a d -1-dimensional staircase surface. This is the natural digital surface associated to ∂X at step h. We show in Theorem 1, Section 3.2, that ∂ h X gets closer to ∂X as h tends toward 0, whatever the dimension of the space. Primal cubical grid at step h. For several proofs, we need to consider the space as a cubical complex. Therefore, we subdivide the space into h-cubes. We need to distinguish cubes, faces, edges and so on. This is why we assign coordinates in ( h 2 Z) d to each cell of the space. To do so, we proceed in a standard manner by cartesian product. Let us associate to each t ∈ h 2 Z the set I h (t), such that for t ∈ hZ, I h (t) := [t -h 2 ; t + h 2 ], and I h (t) := {t} otherwise. Now, for arbitrary z ∈ ( h 2 Z) d , we set

I h (z) := I h (z 1 ) × . . . × I h (z d )
, where z i is the i-th coordinate of z.

Definition 4 (primal cubical grid). The set

F h := {I h (z)} z∈( h 2
Z) d tiles the Euclidean space R d with hypercubes and its faces. It is called the primal cubical grid at step h. Elements of F h are called cells.

The grid F h is a cell complex of dimension d, illustrated on Fig. 3, left. The partial order relation is defined as c 1 c 2 whenever c 1 ⊂ c 2 . The dimension of each cell c is the number of axes where the cell is not reduced to a point, and the set of k-dimensional cells is denoted by F k h . By construction, for any cell c of F h , there is exactly one z ∈ ( h 2 Z) d such that I h (z) is equal to c: the vector z forms the digital coordinates of the cell c. We also use the notation ċ to designate the centroid of the cell c. It is easily checked that they coincide, i.e. z = ċ.

By construction, the h-boundary of X is decomposable as a pure subcomplex of F h of dimension d -1 (see Fig. 3, right). Any one of its d -1-cell is bordered by two incident d-cells c, c in F d h . Their digital coordinates z and z are such that one is in D h X and the other not. Dual cubical grid at step h. It is obvious that we could have obtained a shifted cellular grid by inverting the role of multiples of h and half-multiples of h. Let us associate to each t ∈ h 2 Z the set Ĩh (t), such that for t ∈ hZ, Ĩh (t) := {t}, and Ĩh (t) := [t -h 2 ; t + h 2 ] otherwise. Similarly, we extend Ĩh to arbitrary z ∈ ( h 2 Z) d . Definition 5 (dual cubical grid). The set of sets Fh := { Ĩh (z)} z∈( h 2 Z) d is then called the dual cubical grid at step h, whose elements are called dual cells.

It is clearly a cell complex, with the same definitions of partial order and dimension. Digital coordinates and centroids are also defined similarly.

The sets F k h and Fd-k h have a natural duality isomorphism induced between cells and dual cells with identical coordinates. If we denote the duality operator on cells with the • operator, we clearly have for

c 1 , c 2 ∈ F h , c 1 c 2 ⇔ c2
c1 . The dual cubical grid and its duality with the primal cubical grid are illustrated on Fig. 3, left and middle. Sets with positive reach and properties of projection map. The projection map ξ is continuous on R d \ MA(K), and more precisely Proposition 1 (Theorem 4.8 of [START_REF] Federer | Curvature measures[END_REF]). Let K be a compact set with positive reach. Then for every p ∈ K and every α ∈ [0, 1[, the projection ξ K is 1 1-α Lipschitz in the ball centered on p with radius α • reach(K).

In the particular case where K = ∂X is the boundary of a compact domain of R d , we have the following equivalence: 1 , which means that it is of class C 1 and that the function that assigns the normal to ∂X to each point x on ∂X is Lipschitz.

Proposition 2 ([14]). Let X be a compact domain of R d . The reach of ∂X is strictly positive iff ∂X is a hyper surface of class C 1,
Remark that in Section 5 below, we will provide an explicit upper bound of the Lipschitz constant of the normal map (Lemma 5). Remark also that a manifold ∂X with strictly positive reach is thus of class C 2 almost everywhere. This is a consequence of Rademacher Theorem (3.1.6 in [START_REF] Federer | Geometric measure theory[END_REF]). In the following, we will denote by ξ = ξ ∂X the projection map on ∂X. R-regularity or par(R)-regularity. The R-regularity property was independently proposed by Pavlidis [START_REF] Pavlidis | Algorithms for graphics and image processing[END_REF] and Serra [START_REF] Serra | Image analysis and mathematical morphology[END_REF]. Gross and Latecki introduced the similar definition of par(R)-regularity in [START_REF] Gross | Digitizations preserving topological and differential geometric properties[END_REF], that is the shapes whose normal vectors do not intersect each other, when they are embedded as segments of length 2R. We prefer here to present the definition given in [START_REF] Latecki | Preserving topology by a digitization process[END_REF] with inside and outside osculating balls. A closed ball iob(x, R) of radius R is an inside osculating ball of radius R to ∂X at point

x ∈ ∂X if ∂X ∩ ∂iob(x, R) = {x} and iob(x, R) ⊆ X • ∪ {x}. A closed ball oob(x, R) of radius R is an outside osculating ball of radius R to ∂X at point x ∈ ∂X if ∂X ∩ ∂oob(x, R) = {x} and oob(x, R) ⊆ (R d \ X) ∪ {x}.
A set X is then par(R, +)-regular if there exists an outside osculating ball of radius R at each x ∈ ∂X. A set X is par(R, -)-regular if there exists an inside osculating ball of radius R at each x ∈ ∂X. The par(R)-regularity is the conjunction of these two properties. This definition implies the other definition.

First properties of the boundary of Gauss digitized sets

In this section, we show that the notion of reach, which is classical in geometric measure theory, and the notion of par(R)-regularity, which is known in digital geometry, are related (Lemma 1). We then show that the boundary of X is close to its h-boundary in the Hausdorff sense, and we give tight bounds on the distance (Theorem 1) for arbitrary dimensions. Hence, digitized surfaces tends to the original surface in the Hausdorff sense. Furthermore, the closest point is given by the projection map.

3.1. About R-regularity and positive reach. In the case where X is a d-dimensional object, the reach of ∂X and the R-regularity of X are related as follows.

Lemma 1. Let X be a d-dimensional compact domain of R d . Then reach(∂X) ≥ R ⇔ ∀R < R, X is par(R )-regular.
Proof. Suppose that the reach of ∂X is strictly less than R. We want to show that there exists R < R, such that X is not par(R )-regular. Since reach(∂X) < R, there exists a point x that has two closest points y 1 and y 2 on ∂X and such that d(x, ∂X) = R < R. For simplicity, we assume that x ∈ X. (If x is outside X, then the proof is similar.) Let R be such that R < R < R. We now proceed by contradiction: we assume that X is par(R )-regular and we are going to show that there does not

√ d 2 h-offset of ∂ h X (in gray).
exist any inside osculating ball to ∂X with center at x and radius R , contradicting the hypothesis that X is par(R )-regular. Note that the interior of the closed ball B x (R ) of center x and of radius R = x -y 1 = x -y 2 does not intersect ∂X, but B x (R ) intersects ∂X in at least the two points y 1 and y 2 . Then, the ball B x(R ) of center x = y 1 + R • (x -y 1 )/ x -y 1 and of radius R contains the point y 2 in its interior, thus cannot be an inside osculating ball to ∂X at y 1 .

Consider now any other ball B x(R ) of radius R whose center x does not belong to the straight line going through y 1 and x and such that y 1 ∈ ∂B x(R ). We want to prove that B x(R ) cannot be an inside osculating ball to ∂X at y 1 either. Since X is assumed to be par(R )-regular, there exists an outside osculating ball B x (R ) whose center x belongs to the straight line going through y 1 and x, as the two balls B x (R ) and B x (R ) are tangent at y 1 . But then the interior of the two balls B x(R ) and B x (R ) must intersect, which implies that B x(R ) does not belong entirely to the interior of X, since B x (R ) is an outside ball. So B x(R ) cannot be an inside osculating ball to ∂X at y 1 . This contradicts the fact that X is par(R )-regular.

Let us show the reverse. We suppose that the reach of ∂X is larger than R and are going to show that X is par(R )-regular for every R < R. Since reach(∂X) ≥ R, we know that ∂X is a (d -1)manifold of class C 1 . Let y ∈ ∂X. There exists a unit normal n y to ∂X at y. Furthermore, for any R < R, the point y + R • n y is at a distance R from ∂X. Hence the ball B y+R •ny (R ) only intersects ∂X at the point y. Similarly, the ball B y-R •ny (R ) also only intersects ∂X at the point y, which implies that X is par(R )-regular.

Remark 1. If X is a d-dimensional compact domain of R d whose boundary ∂X has a reach greater than R, then for R < R, any point x ∈ ∂X has an inside osculating ball of radius R and an outside osculating ball of radius R .

3.2.

Hausdorff distance between ∂X and its digital counterpart. We show in Theorem 1, below, that the boundary of X (in blue) and its digital counterpart ∂ h X (in red) are close in the Hausdorff sense, and this property is valid for arbitrary dimensions. This is illustrated on Fig. 4. Note that a 2D version of this theorem was given in [START_REF] Lachaud | Espaces non-euclidiens et analyse d'image : modèles déformables riemanniens et discrets, topologie et géométrie discrète[END_REF], Lemma B.9. For x ∈ ∂X, we denote by n(x, l) the segment of length 2l, centered on x and aligned with the normal vector to ∂X at x. Theorem 1. Let X be a compact domain of R d such that the reach of ∂X is greater than R. Then, for any digitization step 0 < h < 2R/ √ d, the Hausdorff distance between sets ∂X and ∂ h X is less than √ dh/2. More precisely:

∀x ∈ ∂X, ∃y ∈ ∂ h X, x -y ≤ √ d 2 h and y ∈ n(x, √ d 2 h), (2) ∀y ∈ ∂ h X, y -ξ(y) ≤ √ d 2 h. (3)
Remark that this bound is tight.

Proof. We first prove [START_REF] Chazal | Normal cone approximation and offset shape isotopy[END_REF]. Let x ∈ ∂X. Since ∂X has reach greater than R, there is an inside osculating ball of radius

√ d 2 h at x (from Remark 1 and h < 2R √ d ).
There is also an outside osculating ball of same radius at x. Let us denote by c i and c e their respective centers. Point c i (resp. c e ) belongs to at least one h-cube of center p i (resp. p e ), i.e. some Q h . Since p i is at a distance less than or equal to

√ d 2 h from c i (half-diameter of h-cube)
, point p i belongs to the inside osculating ball at x and is thus a point inside X • or is equal to x. Similarly point p e belongs to the complementary set of X or is equal to x. In the latter case, point c e is exactly in a corner of the h-cube Q h

x and we choose for p e another h-cube containing c e , hence p e = x and p e ∈ R d \ X.

The straight segment [c i c e ] is by definition the segment n(x,

√ d 2 h
). We show by contradiction that this segment intersects ∂ h X. Let D be the subset of h-cubes that intersect [c i c e ]. We already know that D contains at least two h-cubes, one of center p i that is in X, one of center p e that is outside X. By connectedness of segment [c i c e ], there is a covering sequence (P j ) j=0..l of h-cubes included in D so that: (i) P 0 has center p i , (ii) P l has center p e , (iii) ∀j, with 0 ≤ j < l and P j ∩ P j+1 = ∅. Since h-cubes are closed, it is easy to derive from (P i ) an enriched covering sequence (P j ) j=0..l of same extremeties such that any two consecutive h-cubes have a d -1-dimensional intersection. Since P 0 has center in X and P l has center outside X, there is an index k so that P k has center in X, and P k+1 has center outside X. By definition, P k ∩ P k+1 ⊂ ∂ h X. Now, [c i c e ] intersects both P k and P k+1 and, by convexity, their intersection. Let us denote by y this intersection. We have y

∈ P k ∩ P k+1 ⊂ ∂ h X. Since y ∈ [c i c e ] = n(x, √ d 2 h), y is at distance to x less than √ d 2 h
. We now prove [START_REF] Chazal | Stability of curvature measures[END_REF]. Let y ∈ ∂ h X. By the definition of h-boundary (cf. (1)), there must exist two h-cubes of center p 1 and p 2 such that p 1 ∈ X and p 2 ∈ X and they share a face (i.e. p 1 -p 2 1 = h). The closed straight segment [p 1 p 2 ] thus intersects ∂X at least once, say at x . By Pythagora's theorem, point x is at a distance less than √ d 2 h from y. Since this distance is smaller than the reach of ∂X, there is a unique point x onto ∂X that is closest to y. This implies that y -x ≤ y -x ≤ √ d 2 h. Furthermore, since ∂X is of class C 1 , the point y belongs to the line-segment normal to ∂X at x. Putting these two facts together gives y ∈ n(x,

√ d 2 h
). Clearly, this implies x = ξ(y) and (3).

Manifoldness of the boundary of Gauss digitized sets

In the whole section, the set X is a compact domain of R d , such that reach(∂X) is greater than some positive constant R. Hence, X is par(R )-regular for any 0 < R ≤ R (Lemma 1). Although Theorem 1 states that the h-boundary of X tends to the boundary of X in the Hausdorff sense, starting from d = 3 and as said in the introduction, the h-boundary of X may however not be a manifold. Focusing on d = 3, we thus exhibit local sufficient conditions which guarantee that the h-boundary is locally a 2-manifold (see Theorem 2 below). These conditions indicates that only places of ∂X with a normal very close to some axis direction may induce a non-manifold place in the h-boundary (dark grey zones in Fig. 2). Even better, if the shape is not flat at these places, these zones tend to area zero with finer digitization gridsteps.

Theorem 2 (Manifoldness sufficient condition). Let X be some compact domain of R 3 , with reach(∂X) greater than some positive constant R and h < 0.198R. Let y be a point of ∂ h X.

i) If y does not belong to some 1-cell of ∂ h X that intersect ∂X, then ∂ h X is homeomorphic to a 2-disk around y. ii) If y belongs to some 1-cell s of ∂ h X such that ∂X ∩ s contains a point P and if the angle α y between s and the normal to ∂X at P satisfies α y ≥ 1.260h/R, then ∂ h X is homeomorphic to a 2-disk around y.

The proof relies on the determination of necessary conditions for the presence of crossed configurations in the digitized set D h (X). A digital set without crossed configuration has the property to be well-composed [START_REF] Latecki | 3D well-composed pictures[END_REF]. And a well-composed set has a boundary that is a 2-manifold. The following subsections detail the steps of the proof of Theorem 2. 4.1. Terminology. Let 1 X be the indicator function of X. Hence, for any z ∈ (hZ) 3

, z ∈ D h (X) ⇔ 1 X (z) = 1. Any dual 3-cell ṽ of F3
h is a cube of side h and whose eight vertices (ṽ i ) 0≤i≤7 are points of (hZ) 3 , numbered according to the lexicographic ordering of their z, y, x coordinates. The 8configuration of X at ṽ is the 8-tuple 1 X (ṽ) := (1 X (ṽ 0 ), . . . , 1 X (ṽ 7 )). Let s be a dual 2-cell that is a face of ṽ. It is a square of side h whose four vertices s0 , s1 , s2 , s3 are points of (hZ) 3 , numbered counterclockwise when standing at the tip of the 1-cell s with maximal coordinate and looking at s. They form a subset of (ṽ i ) 0≤i≤7 . The 4-configuration of X at s is the 4-tuple 1 X (s) := (1 X (s 0 ), 1 X (s 1 ), 1 X (s 2 ), 1 X (s 3 )).

A crossed 8-configuration is any rotation or complementation of (1, 0, 0, 0, 0, 0, 0, 1) (there are 8 such configurations). A crossed 8-configuration at a dual 3-cell ṽ induces a non-manifold vertex in the hboundary of X, precisely at the primal 0-cell v. It corresponds locally to two cubes glued together only at one vertex. A crossed 4-configuration is either the 4-configuration (1, 0, 1, 0) or the 4-configuration (0, 1, 0, 1). It is obvious that a crossed 4-configuration at a dual 2-cell s induce a non-manifold edge in the h-boundary of X, precisely at the primal 1-cell s. It corresponds locally to two cubes glued together only along one edge. We recall (and adapt with our notations) Proposition 2.1 of [START_REF] Latecki | 3D well-composed pictures[END_REF].

Proposition 3 ([25]

). The h-boundary of X is a 2-dimensional manifold if and only if X has no crossed configurations in any dual 2-cell or 3-cell of Fh . (In this case, D h X is called a well-composed picture.) Non-manifoldness is thus determined by the presence of crossed configurations. We will thus exhibit sufficient conditions that prevent them to appear.

4.2.

Relations between crossed configurations and grid step. We study the presence of crossed configurations depending on whether the boundary ∂X intersects or not cells of the cubical grid sampled at step h. The first lemma is straightforward.

Lemma 2. If ∂X does not intersect a dual 2-cell s of F2

h , then the 4-configuration of X at s is not crossed.

Proof. Then s

⊂ R 3 \ ∂X = X • ∪ (R 3 \ X).
Since s is connected while the previous union is disjoint, we have two cases, either s ⊂ X • and the 4-configuration is (1, 1, 1, 1), or s ⊂ R 3 \ X and the 4configuration is (0, 0, 0, 0).

The second case tackled below is more involved. The idea is to look at how inner or outer osculating balls contains vertices of s or s. It appears that crossed 4-configurations cannot arise when h is small enough.

Lemma 3. Let h ≤ 0.198R. If ∂X intersects a dual 2-cell s of F2
h but does not intersect the corresponding primal 1-cell s, then the 4-configuration of X at s is not crossed.

Proof. This lemma is illustrated on Fig. 5. If all vertices of s are in X, or all vertices of s are outside X, then the 4-configuration of X at s is clearly not crossed and we are done. Hence at least one vertex of s, say s0 , is in X (but may be on ∂X) and at least one other vertex of s is outside X. We assume here that the primal 1-cell s (a segment of length h whose extremities are denoted by s 0 and s 1 ) lies outside X. Should the 1-cell s be completely inside X • , then we would reason on the vertex of s that lies outside X, and the reasonning would be symetrical. Without loss of generality, let s0 be this vertex in X, and let Q be the center of s0 s2 . The segment [s 0 Q] is a connected set that joins a point in X to a point in R 3 \ X (since Q ∈ s). Hence, there exists a point P ∈ [s 0 Q] ∩ ∂X. According to Remark 1, there is thus an inside osculating ball B in and an outside osculating ball B out of radius R at P . Let α be the angle between the normal n to ∂X at P and the segment s (oriented in the same direction). Let ñ be the projection of n onto the plane Π supporting s. The angle between ñ and the Left. The point s0 in black belongs to X • while the segment s in green lies in R 3 \ X. Point P in blue belongs to ∂X and its normal vector is also in blue. The angle between the normal at P and the axis of s is denoted by α. The inside osculating ball of radius R at P is drawn in orange. Middle. In this cross-section view (plane containing s0 and s), the projection of the normal at P makes an angle α 0 ≤ α with the the segment s. The trace of the inside osculating ball at P is drawn in orange. Right. In this transverse view (plane containing s), the trace of the outside osculating ball at P is also drawn in an orange dashed line. The angle β of the normal at P projected onto this plane is drawn in magenta. oriented segment s0 s2 is denoted by β. The angle α can be taken in [0, π 2 ], while β can be taken in [0, π] (negative β implies a reasonning on s3 instead of s1 ). We center the frame on Q, with x-axis aligned with [s 0 s2 ], z-axis aligned with segment s, y-axis aligned with [s 1 s3 ].

The idea of the proof is that, since the inside osculating ball B in at P does not touch s, the angle α may not be too small, which in turn prevents crossed configurations to occur. Indeed, this situation is depicted on Fig. 5. Setting the coordinates of P to be (-, 0, 0), the center C in of B in lies at (--R cos β sin α, -R sin β sin α, -R cos α). The vertex s 0 has coordinates (0, 0, -h 2 ). Since s 0 ∈ B in , we have

R 2 ≤ C in s 0 2 ⇔ hR cos α -2 R cos β sin α ≤ h 2 4 + 2 . ⇒ R(cos α -2 h cos β sin α) ≤ 3h 4 (4) (since 0 < < √ 2 2 h and h > 0) Since h > 0, cos β ≤ 1 and ≤ √ 2
2 h, we deduce that:

cos α - √ 2 sin α ≤ 3h 4R . (5)
Remark that the function g : [0, π 2 ] → R defined by g(α) := cos α -√ 2 sin α is decreasing and satisfies g(0) = 1. It follows that if h < cR, where c ≤ 4/3, then the angle α is greater than α , where g(α ) = 3c/4. The angle α may thus not be too small as the grid step gets finer.

To prove that the configuration at s is not crossed, it is sufficient to prove either that s1 ∈ B in or that s2 ∈ B out . Considering that s1 = (0, -

√ 2 2 h, 0) and s2 = ( √ 2
2 h, 0, 0), we derive in a similar manner the following relations:

s1 ∈ B in ⇔ C in s1 2 < R 2 ⇐ h R sin α < ( √ 2 sin β -2 h cos β) ⇐ h R sin α < √ 2 sin β - √ 2 cos β, (6)
since < h/ √ 2, and

s2 ∈ B out ⇔ C out s2 2 < R 2 ⇐ h R sin α < ( √ 2 2 cos β + h cos β) ⇐ h R sin α < √ 2 2 cos β, (7)
since 0 ≤ . It is thus enough for h/(R sin α) to be lower than the maximum of both bounds a(β

) := √ 2 sin β - √ 2 cos β and b(β) := √ 2
2 cos β of ( 6) and [START_REF] Coeurjolly | Integral based curvature estimators in digital geometry[END_REF]. One can show that there is a unique value β ∈ [0, π] such that a(β ) = b(β ). This value is given by β = tan -1 (3/2). It is easily seen that for every β ∈ [0, π], one has a(β

) = b(β ) ≤ max(a(β), b(β)). Since a(β ) = √ 26 13 , it follows h R sin α < √ 26 13 =⇒ s1 ∈ B in or s2 ∈ B out . (8)
We now wish to find the best constant c such that, when h < cR, either s1 ∈ B in or s2 ∈ B out , and thus the configuration is not crossed. We choose c such that 4 3 c = g(α ) and α is given by c sin α = √ 26 [START_REF] Esbelin | Convergence of binomial-based derivative estimation for c 2 noisy discretized curves[END_REF] . In that case, it follows from (5) that g(α) ≤ g(α ). Since g is decreasing, one has α ≥ α , and thus

c sin α ≤ c sin α = √ 26
13 . This implies by [START_REF] Coeurjolly | Multigrid convergent principal curvature estimators in digital geometry[END_REF] Numerical approximation gives c ≈ 0.198.

We turn ourselves to the last case, where we show that the direction of the normal to ∂X plays a role in the manifoldness of its digitized counterpart. Lemma 4. Assume ∂X intersects a primal 1-cell s of F 1 h at some point P . Let α be the angle between the normal n at P and the vector u aligned with direction s. Then the 4-configuration of X at s is not crossed whenever 1.260 h R < α. Proof. The idea is to measure the distance between vertices si (for i ∈ {0, 1, 2, 3}) and the center of the inside (resp. outside) osculating ball at P . Such osculating balls of radius R exist according to Remark 1. If this distance is smaller than R then we know that the value of 1 X (s i ) is 1 (resp. 0). Indices i are taken modulo 4. The distance of P to s is denoted by . Without loss of generality, the angle α is taken in [0, π 2 ]. Otherwise, a symmetric reasonning can be applied with the outside osculating ball. The frame denoted Π i is centered on P with x-axis directed as [s i si+2 ] and with z-axis directed as s, and oriented such that si has non positive z-coordinate. As in the proof of the previous lemma, let ñ be the projection of the outer normal at P onto the plane Π supporting s. The angle between ñ and the oriented segment [s i si+2 ] is denoted by β i .

Since α and β i represents the latitude and longitude of vector n, the center C in of the inside osculating ball has coordinates -R (sin α cos β i , sin α sin β i , cos α) in frame Π i . Furthermore, point si has coordinates (-√ 2 2 h, 0, -). Since the inside osculating ball is in X • , we deduce

1 X (s i ) = 1 ⇐ C in si 2 < R 2 ⇔ h 2 2 + 2 < √ 2Rh sin α cos β i + 2 R cos α ⇐ 3 4 h R < √ 2 sin α cos β i + 2 h cos α, (9) since ≤ h 2 . When angle β i ∈ [-π 4 , π 4 ], we have cos β i ≥ √ 2
2 . Inserting also ≥ 0 into (9) gives

1 X (s i ) = 1 ⇐ 3 4 h R < sin α ⇐ 1.179 h R < α, (10)
using α π/2 ≤ sin α and 3π 8 ≈ 1.1781. Clearly there is at least one

β j ∈ [-π 4 , π 4 ]
. Hence 1 X (s j ) = 1 for 1.179 h R < α. We prove either that the opposite vertex to sj on s is outside X, i.e. 1 X (s j+2 ) = 0, or that one of the neighboring vertex of sj is inside X, i.e. 1 X (s j+1 ) = 1 or 1 X (s j-1 ) = 1. We prove the case β j ∈ [0, π 4 ], hence we determine the bounds for which either 1 X (s j+2 ) = 0 or 1 X (s j+1 ) = 1. Negative values of β j are tackled similarly with 1 X (s j-1 ) = 1.

One easily checks that, in the frame Π j , sj+1 = (0, -

√ 2 2 h, -), sj+2 = ( √ 2
2 h, 0, -) and the center C out of the outside osculating ball has symmetric coordinates to C in , i.e. C out = -C in . With computations similar to [START_REF] Cuel | Voronoi-based geometry estimator for 3d digital surfaces[END_REF], we derive

1 X (s j+2 ) = 0 ⇐ C out sj+2 2 < R 2 ⇐ 3 4 h R < √ 2 sin α cos β j -2 h cos α. (11) 1 X (s j+1 ) = 1 ⇐ C in sj+1 2 < R 2 ⇐ 3 4 h R < √ 2 sin α sin β j + 2 h cos α. ( 12
)
It is sufficient to have either [START_REF] De Vieilleville | Maximal digital straight segments and convergence of discrete geometric estimators[END_REF] or [START_REF] Esbelin | Convergence of binomial-based derivative estimation for c2-noisy discretized curves[END_REF] to get a non crossed configuration. We look therefore at the maximum of both values. Denoting f (α, ν) := √ 2 sin α cos β j -2ν cos α and g(α, ν) := √ 2 sin α sin β j + 2ν cos α, we rewrite those equations as:

1 X (s j+2 ) = 0 or 1 X (s j+1 ) = 1 ⇐ 3 4 h R < max f α, h , g α, h ⇐ 3 4 h R < √ 2 2 f 2 α, h + g 2 α, h . ( 13
)
The last implication comes from the property that √ a 2 + b 2 / √ 2 ≤ max(a, b) holds not only for positive values a and b but in the more general case where -min(a, b) ≤ max(a, b). Here f may take negative values but, when negative, it is always smaller in absolute value than g. Simple calculations give:

f 2 (α, ν) + g 2 (α, ν) = 8 cos 2 αν 2 -2 √ 2 sin α cos α(cos β j -sin β j )ν + 2 sin 2 α ≥ 8 cos 2 αν 2 -2 √ 2 sin α cos αν + 2 sin 2 α, ( 14
) since β j ∈ [0, π 4 ]
. The last term is a degree 2 polynomial in ν that we denote h α (ν). It has discriminant -56 sin 2 α cos 2 α, which is non positive for arbitrary α ≥ 0. Hence, h α (ν) takes minimum value at

ν α = √ 2 8 tan α. Simple calculations lead to f 2 (α, ν) + g 2 (α, ν) > h α (ν α ) = 7 4 sin 2 α > 7 π 2 α 2 , ( 15
)
since sin α ≥ α/(π/2). Inserting inequality ( 15) into (13) implies:

1 X (s j+2 ) = 0 or 1 X (s j+1 ) = 1 ⇐ 3 4 h R < √ 2 2 √ 7 π α ⇐ 1.260 h R < α, ( 16
) since 3 √ 14π 28
≈ 1.2594. If both ( 10) and ( 16) hold, then the configuration at s is one of (1, 1, ?, ?), (1, ?, 0, ?) and circular permutations. Hence the configuration is not crossed when 1.260 h R < α.

To get an idea of the practical implication of previous Lemma, if one consider a shape with reach 1, then there might be a non-manifold zone on its digitization at gridstep 1 10 only at places where the normal makes an angle smaller than 7.5 • with one axis. For instance, this is less than 2.57% of the area on a sphere. We have now all the pieces to finish the proof of Theorem 2.

of Theorem 2. According to Proposition 3, the manifoldness of ∂ h X is determined by the absence of crossed configurations. Non manifoldness at a primal vertex v occurs only if the 8-configuration of X at ṽ is crossed. Theorem 13 of [START_REF] Stelldinger | Topological equivalence between a 3d object and the reconstruction of its digital image[END_REF] together with the equivalence of par-regularity and reach given by Lemma 1 show that h < 0.5R implies that the the 8-configuration is not crossed. Non manifoldness at a primal edge s occurs only if the 4-configuration of X at s is crossed. This case is fully studied in Lemma 2, Lemma 3 and Lemma 4. Non manifoldness at a primal 2-cell is impossible by construction. This concludes the proof.

Size of the non injective part

Here, the set X is a compact domain of R d , whose boundary ∂X has reach strictly greater than R. We assume that h ≤ R/ √ d, which implies by Theorem 1 that the Hausdorff distance between ∂X and ∂ h X is less than R/2. Therefore the projection map ξ on ∂X is well defined on ∂ h X. However, this map is not one-to-one in general.

The aim of this section is to show that the subset of ∂X for which ξ is not injective from ∂ h X, otherwise said the part of ∂X with multiplicity greater than one through projection, is small. We define the following set

mult(∂X) := {x ∈ ∂X, s.t. ∃y 1 , y 2 ∈ ∂ h X, y 1 = y 2 and ξ(y 1 ) = ξ(y 2 ) = x}. Theorem 3. If h ≤ R/ √ d, then one has Area (mult(∂X)) ≤ K 1 (h) Area (∂X) h,
where

K 1 (h) = 4 √ 3 d 2 R + O(h) ≤ 2 √ 3 d 2 4 d R .
Here and in the sequel, the constant appearing in O(h) only involves the dimension d and the reach R. Furthermore, the (d -1)-dimensional Hausdorff measure is denoted by Area and the d-dimensional Hausdorff measure is denoted by Vol.

Sketch of proof.

The assumption h ≤ R/ √ d implies by Theorem 1 that the Hausdorff distance between ∂ h X and ∂X is less than √ dh/2. In particular, one has for every y ∈ ∂ h X, y-ξ(y) ≤ √ dh/2. Furthermore, Theorem 1 also implies that the restriction of the projection map to ∂ h X is surjective. However, it may not be injective in general. We introduce the set mult(∂ h X) = ξ -1 (mult(∂X)). Clearly, the map ξ :

∂ h X \ mult(∂ h X) → ∂X \ mult(∂X)
is one-to-one. For any point x ∈ ∂X, we denote by n(x) the outward unit normal vector to ∂X at x and by n h (y) the outward unit normal vector to ∂ h X at y. Remark that n h (y) is defined almost everywhere for the (d -1) Hausdorff measure. If y belongs to the intersection of two or more (d -1)dual cells, then we can choose for n h (y) the outward unit normal to any of those cells. The outline of the proof is the following: i) We show that the scalar products between normals of ∂ h X and ∂X is always greater than -2 √ dh/R. ii) We show that mult(∂X) ⊂ ξ(P (h)), where

P (h) := {y ∈ ∂ h X, n(ξ(y)) • n h (y) ≤ 0}.
iii) We show that the jacobian of ξ at y is approximately |n(ξ(y)) • n h (y)|, hence the jacobian of its restriction to P (h) is in O(h). iv) We conclude that Area (mult(∂X)) is in O(h).

Angle relation between object boundary and its digitization.

Let X be a compact domain of R d , whose boundary ∂X has reach strictly greater than R. By Proposition 2, we know that ∂X is of class C 1,1 , meaning that the normal to ∂X is Lipschitz. We provide below an explicit upper bound of this Lipschitz constant. Lemma 5. For any x 1 , x 2 ∈ ∂X, one has

(17) n(x 1 ) -n(x 2 ) ≤ √ 3 R x 1 -x 2 .
Proof. For i = 1, 2 we denote by c i the center of the outside osculating ball of radius R to ∂X at the point x i , by c i the center of the inside osculating ball to ∂X at the point x i . Since the ball B c1 (R) is included in X and B c 2 (R) is included in the closure of R d \ X, their interior do not intersect and thus c 1 -c 2 ≥ 2R. From the fact that c i -

x i = R n(x i ), one has c 1 -c 2 = (c 1 -x 1 ) + (x 1 -x 2 ) + (x 2 -c 2 ) = R n(x 1 ) + (x 1 -x 2 ) + R n(x 2 ), which implies that c 1 -c 2 2 = 2R 2 + x 1 -x 2 2 + 2R 2 n(x 1 ) • n(x 2 ) + 2R(x 1 -x 2 ) • [n(x 1 ) + n(x 2 )] . The condition c 1 -c 2 2 ≥ 4R 2 thus implies that R 2 n(x 1 ) -n(x 2 ) 2 = 2R 2 [1 -n(x 1 ) • n(x 2 )] ≤ x 1 -x 2 2 + 2R (x 1 -x 2 ) • [n(x 1 ) + n(x 2 )] . (18) It remains to show that 2R (x 1 -x 2 ) • [n(x 1 ) + n(x 2 )] is bounded by 2 x 1 -x 2
2 , which will allow to conclude. Remark that the two points x 1 and x 1 + 2Rn(x 1 ) belong to the sphere ∂B c1 (R) and are diametrally opposed. Thus, since x 2 does not belong to the ball B c1 (R), one has

(x 2 -x 1 ) • (x 2 -(x 1 + 2R n(x 1 ))) ≥ 0 ⇔ (x 2 -x 1 ) • ((x 2 -x 1 ) -2R n(x 1 )) ≥ 0 ⇔ 2R (x 2 -x 1 ) • n(x 1 ) ≤ x 2 -x 1 2
Similarly, since x 2 does not belong to the ball B c 1 (R), the same inequality holds by replacing n(x 1 ) with -n(x 1 ) and thus

|2R (x 2 -x 1 ) • n(x 1 )| ≤ x 2 -x 1 2 .
Similarly,

x 1 does not belong to B c2 (R) ∪ B c 2 (R), which implies |2R (x 2 -x 1 ) • n(x 2 )| ≤ x 2 -x 1 2 .
Plugging these last two equations into (18) leads to

R 2 n(x 1 ) -n(x 2 ) 2 ≤ 3 x 2 -x 1 2 ,
which allows to conclude. Lemma 6. Let p ∈ X and q ∈ X, then there exists

x ∈ ∂X ∩ [pq] such that n(x) • - → pq ≥ 0.
Proof. First of all, X ∩ [pq] is not empty (it contains at least p) and is compact. In this compact set, we define x as the closest point to q. It is also clear that x ∈ ∂X. Assume that n(x) • -→ pq < 0, then the inside osculating ball at x of radius R intersect of (xq]. This is a contradiction since x was the closest point of X to q along this segment.

Lemma 7. For any y ∈ ∂ h X, the angle between the normal n h (y) of any (d-1)-cell of ∂ h X containing y and the normal of its projection x = ξ(y) onto ∂X satisfies:

n(x) • n h (y) ≥ - √ 3d R h.
Proof. Let x = ξ(y). If n(x)•n h (y) is positive, the result is obvious. We suppose now that n(x)•n h (y) < 0. Since y ∈ ∂ h X, it belongs to a primal 2-cell c, whose dual 1-cell c is a segment [pq], where p ∈ X and q ∈ X. Note that the normal n h (y) at y on ∂ h X points in the same direction as the vector -→ pq. Proof. Let x ∈ mult(∂X) and [ab] = n(x, √ dh/2) the segment centered in x of length √ dh and aligned with the normal n(x). We suppose that this segment touches several (d -1)-faces of ∂ h X and is not in the tangent plane of one of these faces (otherwise, the conclusions holds directly). To get the proof, it is sufficient to show that there is an orthonormal axis-aligned frame ( -→ e j ) j=1,...,d such that: (i) ∀j, with 1 ≤ j ≤ d, -→ ab • -→ e j ≥ 0, (ii) some intersected face of ∂ h X has a normal --→ e j2 . Let σ 1 , σ 2 be two d -1-faces of ∂ h X intersected by [ab]. We may consider the vector -→ ab to be in the first orthant of the space, with some choice of the reference frame ( -→ e j ) j=1,...,d . The segment [ab] crosses several cubes of F d h , from which one can extract a covering face-adjacent subsequence of cubes (c i ) i=1..m . Because -→ ab is in the first orthant, we have that ∀i, with 1 ≤ i < m, ∃k i ∈ {1, . . . , d},

----→ ċi ċi+1 ----→ ċi ċi+1 = + -→ e ki .
The faces σ 1 and σ 2 , being intersected by the segment, are the faces of some cubes c i1 and c i2 . Furthermore, the segment being not in their tangent planes, these faces are the intersection of consecutive cubes in the sequence (c i ), and we have σ 1 = c i1 ∩ c i1+1 and σ 2 = c i2 ∩ c i2+1 . We choose first i 1 < i 2 .

Two cases arise, either ċi1 ∈ X or not. In the first case, necessarily ċi1+1 ∈ X and the normal at σ 1 is then + --→ e ki 1 . Now since σ 2 ⊂ ∂ h X, either ċi2 or ċi2+1 belongs to X. Since i 1 < i 2 , there must be some i 3 , i 1 + 1 ≤ i 3 ≤ i 2 , with ċi3 ∈ X and ċi3+1 ∈ X. The face c i3 ∩ c i3+1 , which may be σ 2 , thus belongs to ∂ h X. Its normal vector is ---→ e ki 3 , which concludes this case. The other cases are solved identically.

= ξ |∂ h X is given by Jξ (y) = |n(ξ(y)) • n h (y)| K 2 (h)
where

K 2 (h) = 1 + O(h) ≤ 1 1 - √ d 2R h d-1 ≤ 2 d-1 .
Proof. First remark that if n(ξ(y)) • n h (y) = 0, then Jξ (y) = 0 and the result holds. If y ∈ ∂ h X is such that n(ξ(y)) • n h (y) = 0, then the map ξ is injective in a neighborhood of y. Furthermore, since ∂X is of class C 2 almost everywhere, we know that for almost every y ∈ ∂ h X such that n(ξ(y)) • n h (y) = 0, ∂X is of class C 2 at the point ξ(y). Let us take such a point y. It is known that ξ is differentiable at y and one has [32, Lemma 3, section 13. .

Hence, knowing that Jξ (y) = | det(Dξ (y))|, we get

|n(ξ(y)) • n h (y)| 1 1 + √ dh 2R d-1 ≤ Jξ (y) ≤ |n(ξ(y)) • n h (y)| 1 1 - √ dh 2R d-1
.

5.5.

Relating areas of continuous and digitized boundaries. We determine an explicit upper bound for the area of the digitized boundary ∂ h X with respect to the area of the continuous boundary ∂X. We denote by ∂X the -offset of ∂X (i.e., the Minkowski sum of ∂X with the ball of radius ), or equivalently

∂X := x ∈ R d , x -ξ(x) ≤ .
Lemma 10. Area (∂ h X) ≤ Area (∂X) K 3 (h), where

K 3 (h) = 4d 3 2 + O(h) ≤ 2 d+2 d 3 2 .
Proof. By Theorem 1, Equation (3), any point on ∂ h X is at distance lower than √ d 2 h from ∂X. Therefore, all faces of ∂ h X are included in the

√ d 2 h-offset of ∂X.
To get a set of cubes that contains all these faces, it suffices to take an offset twice bigger. Let us denote by F (h) the subset of the cellular grid F d h that lies in this offset ∂X √ dh , and by N (h) the number of (hyper)cubes of F (h). Every face of ∂ h X is some face of a cube of F (h). Hence, you may not have more faces in ∂ h X than they are faces of cubes of F (h). Since each cube has 2d faces, it follows that:

Area (∂ h X) ≤ 2d × h d-1 × N (h) From the fact that F (h) ⊂ ∂X √ dh , one has h d N (h) = Vol (F (h)) ≤ Vol ∂X √ dh ,
which implies with the previous equation that

Area (∂ h X) ≤ 2d × h d-1 × Vol ∂X √ dh h d ≤ 2d h Vol ∂X √ dh .
We put = √ dh. We are now going to bound the volume of ∂X . Weyl's tube formula expresses this volume as a polynomial in of degree d [START_REF] Weyl | On the volume of tubes[END_REF]. Since ∂X is of class C 2 almost everywhere, the coefficients are related to the principal curvatures but, here, every one of them can be upper bounded by 1/R. Hence, the volume is upper bounded as:

Vol (∂X ) ≤ 2Area (∂X) + d 1 1 R 2 + d 2 1 R 2 3 + . . . + d d 1 R d d+1 .
From this, we get that Vol (∂X ) ≤ Area (∂X) × 2( + O( 2)) and thus

Area (∂ h X) ≤ 2d h × 2 × √ dh + O(h 2 ) Area (∂X) ≤ 4d 3 2 + O(h) Area (∂X) .
One may also remark that since ≤ R, then we have an explicit upper bound Vol (∂X ) ≤ 2 d+1 Area (∂X) , which implies

Area (∂ h X) ≤ 2d h 2 d+1 Area (∂X) √ dh ≤ 2 d+2 d 3 2
Area (∂X) .

5.6. End of proof of Theorem 3. From Lemma 8, one has mult(∂X) ⊂ ξ(P (h)), where

P (h) := {y ∈ ∂ h X, n(ξ(y)) • n h (y) ≤ 0}.
Therefore Area (mult(∂X)) ≤ Area (ξ(P (h))). Let y ∈ P (h). By Lemma 7, one has

|n(ξ(y)) • n h (y)| ≤ √ 3d R h,
which implies by Lemma 9 that for almost every y ∈ P (h)

Jξ (y) ≤ √ 3d R h K 2 (h). Hence Area (mult(∂X)) ≤ √ 3d R h K 2 (h) Area (P (h)) .
Now, since P (h) ⊂ ∂ h X, one has by Lemma 10

Area (P (h)) ≤ Area (∂ h X) ≤ K 3 (h) Area (∂X) .
Putting this all together, one gets

Area (mult(∂X)) ≤ √ 3d R h K 2 (h) K 3 (h) Area (∂X) .
We conclude by letting

K 1 (h) = √ 3d R K 2 (h) K 3 (h).

Digital surface integration

In this section, we prove the convergence of a digital surface integral. Given a function f : R d → R, we let f ∞ := max x∈R d |f (x)| and denote Lipf := max x =y |f (x)-f (y)|/ x-y its Lipschitz constant, which can be infinite. We define the bounded-Lipschitz norm by f BL := f ∞ + Lip(f ). Given a normal estimator n defined on ∂ h X, we define the error of the normal estimation by nn est := sup y∈∂ h X n(ξ(y)) -n(y) .

We introduce the following digital surface integral. Definition 6. Let Z ⊂ (hZ) d be a digital set, with gridstep h > 0 between samples. Let f : R d → R be an integrable function and n be a digital normal estimator. We define the digital surface integral by DI

h (f, Z, n) := c∈F d-1 h ∩∂Q h Z h d-1 f ( ċ)|n( ċ) • n( ċ)|,
where ċ is the centroid of the (d -1)-cell c and n( ċ) is its trivial normal as a point on the h-boundary ∂ h X. The latter notation is valid only for cells of the primal cubical grid belonging to ∂ h X.

We prove the multigrid convergence of the digital surface integral toward the surface integral.

f (x)dx -DI h (f, D h (X), n) ≤ Area (∂X) f BL O(h) + O( n -n est ) .
Note that as before, the constant involved in the notation O(.) only depends on the dimension d and the reach R.

6.1. Multiplicity of the projection. We show in the section that the multiplicity of ξ is bounded almost everywhere for the (d -1)-Hausdorff measure. One introduces the subset C of ∂X as Let now x ∈ C. Then there exists y ∈ ξ -1 (x) such that the segment n(x, √ dh/2) is contained in a hyperplane P y orthogonal to n h (y). The number of intersections of n(x, √ dh/2) with the cells of F d-1 h that are not parallel to P y are bounded as previously by (d -1) √ d + 1 . For every y ∈ P y ∩ ξ -1 (x), one has n(ξ(y )) • n h (y ) = 0, hence the jacobian of ξ vanishes. Furthermore, in a neighborhood of x, C is included in ∂X ∩ P y which is a curve. Hence the area of C is equal to 0. 6.2. Proof of Theorem 4.

C := {ξ(y), s.t. y ∈ ∂ h X, n(ξ(y)) • n h (y) = 0}.
Step 1. We first show that

∂X f (x)dx = ∂X\mult(∂X) f (x)dx + K 1 (h)Area (∂X) f ∞ h. ( 19
)
We start by writing the integral of f as the sum of two other integrals:

∂X f (x)dx = ∂X\mult(∂X) f (x)dx + mult(∂X) f (x)dx.
According to Theorem 3 (Section 5), the second term is bounded by

mult(∂X) f (x)dx ≤ Area (mult(∂X)) f ∞ ≤ K 1 (h)Area (∂X) f ∞ h.
Step 2. The map ξ induces a bijection from ∂ h X \ mult(∂ h X) to ∂X \ mult(∂X). It is also a diffeomorphism since ∂ h X is within the reach of ∂X by Theorem 1. By the change of variable formula, one obtains: [START_REF] Klette | Digital Geometry: Geometric Methods for Digital Picture Analysis[END_REF] ∂X\mult(∂X)

f (x)dx = ∂ h X\mult(∂ h X)
f (ξ(y))Jξ(y)dy.

Step 3. We now want to show that

µ x f (x)dx ≤ Area (mult(∂X)) µ f ∞ ≤ K 1 (h)Area (∂X) µ f ∞ h.
Step 4. We now show that We then have (with Lemma 10) Step 5. We now show that We write the integral as a sum of integrals on each face of ∂ h X. ≤ Area (∂ h X) Lip(f )

√ dh + f ∞ n -n ∞ ≤ Area (∂X) K 3 (h) Lip(f ) √ dh + f ∞ n -n est .
End of proof. Putting together the equations ( 19), ( 20), ( 21), ( 22), ( 23) of Steps 1-5, one gets

∂X f (x)dx -DI h (f, D h (X), n) ≤ Area (∂X) (Lip(f ) + f ∞ )O(h) + f ∞ O( n -n est ) .
Experimental evaluation. We briefly evaluate numerically the digital surface integral formula for the purpose of area estimation of a 3D digital shape. Fig. 6 illustrates the area estimation error of digital surface integration for several digital normal estimators. Of course, the naive summation of the areas of each 2-cell leads to a non-convergent estimation that overestimates the true area by almost 45% (naive digital area). If the normal is estimated by averaging the trivial cell normals of cells at distance at most t (called trivial normal of radius t), then better area estimations are obtained (around 1% for t = 2). Still they are not convergent. If we use the exact ellipsoid normals (true normal) or convergent normal estimators like integral invariants (II, [START_REF] Coeurjolly | Integral based curvature estimators in digital geometry[END_REF][START_REF] Coeurjolly | Multigrid convergent principal curvature estimators in digital geometry[END_REF]) or Voronoi Covariance Measure (VCM, [START_REF] Cuel | Voronoi-based geometry estimator for 3d digital surfaces[END_REF]), then the digital surface integral appears convergent toward the true area. Even better, experimental convergence speed looks like O(h 2 ). Discussion. We have presented numerous properties of Gauss digitized sets in arbitrary dimension, with a special focus on the relations between the continuous boundary of the shape and the boundary of its digitization at some gridstep h. Although these sets are close in the Hausdorff sense through the projection map, they are not related by an homeomorphism starting from dimension 3. We have characterized precisely places where the digitized boundary is not a manifold in dimension 3. Their area is rapidly decreasing with the grid step (O(h 2 ) on non-flat parts). Furthermore, we have determined where the projection map is not a homeomorphism in arbitrary dimension, and it appears also that the problematic places on the shape boundary have an area that decreases toward zero in O(h). Thanks to this result, we have proven the validity of the digital surface integral as a multigrid convergent integral estimator, as long as the digital normal estimator is also multigrid convergent. Bounds have been made explicit and justify a posteriori previous papers using digital surface integration for area estimation [START_REF] Lenoir | Fast computation of the normal vector field of the surface of a 3-d discrete object[END_REF][START_REF] Coeurjolly | Multigrid convergence and surface area estimation[END_REF]. Experimental evaluations confirm this result. It remains to be understood why the convergence speed is better than expected. This observation seems related to the fact that places likely to induce a non homeomorphic projection are probably overestimated, and thus introduce a larger error on the integration. We are currently examining this issue. However, we cannot expect to achieve better than O(h 2 ) error since even true normals induce this error.

Figure 2 .

 2 Figure 2. Illustration of Theorem 2 and Theorem 3 on several Gauss digitizations of two polynomial surfaces (top row displays a Goursat's smooth cube and bottom row displays Goursat's smooth icosahedron). Zones in dark grey indicates the surface parts where the Gauss digitization might be non manifold (Theorem 2); their relative area is denoted by A nm . Zones in light grey (and dark grey) indicates the surface parts where projection ξ might not be an homeomorphism (Theorem 3); their relative area is denoted by A ξ . Clearly, both zones tends to area zero as the gridstep gets finer and finer, while parts where digitization might not be manifold are much smaller than parts where ξ might not be homeomorphic.

Figure 4 .

 4 Figure 4. Illustration of the fact that the boundary of X and the h-boundary of X are Hausdorff close, with distance no greater than √ d 2 h. On the left, ∂ h X lies in the

Figure 5 .

 5 Figure 5. Illustration of Lemma 3.Left. The point s0 in black belongs to X • while the segment s in green lies in R 3 \ X. Point P in blue belongs to ∂X and its normal vector is also in blue. The angle between the normal at P and the axis of s is denoted by α. The inside osculating ball of radius R at P is drawn in orange. Middle. In this cross-section view (plane containing s0 and s), the projection of the normal at P makes an angle α 0 ≤ α with the the segment s. The trace of the inside osculating ball at P is drawn in orange. Right. In this transverse view (plane containing s), the trace of the outside osculating ball at P is also drawn in an orange dashed line. The angle β of the normal at P projected onto this plane is drawn in magenta.

3 .Lemma 8 .

 38 Then we apply Lemma 6 for segment [pq], and we denote by x 2 the point of ∂X ∩[pq] such that n(x 2 )•n h (y) ≥ 0. By Theorem 1, equation (3), we have that x -y ≤ √ d 2 h. Since y ∈ c and x 2 ∈ [pq] = c, we also have y -x 2 ≤ √ d 2 h. We conclude by the triangle inequality that x -x 2 ≤ √ dh. Since h < R √ d , one has x -x 2 ≤ R, and one can apply Lemma 5 |n(x) • n h (y)| ≤ |(n(x) -n(x 2 )) • n h (y)| ≤ n(x) -n(x 2 ) Parameterization of mult(∂X). For every x ∈ mult(∂X), there exists y ∈ ∂ h X and a 2-cell c containing y, such that ξ(y) = x and n(x) • n h ( ċ) ≤ 0.

5. 4 .Lemma 9 .

 49 Jacobian of the projection. We consider here the restriction ξ := ξ |∂ h X of ξ to ∂ h X. Recall that the (d -1)-jacobian Jξ (y) of ξ at a point y measures the distortion of area induced by the map ξ near y, that is Jξ (y) := lim →0 Area (ξ (B(y, ))) Area (B(y, )) , where B(y, ) denotes the (d -1)-dimensional ball of radius centered at y on ∂ h X. For almost every y ∈ ∂ h X (for the (d -1)-Hausdorff measure), the (d -1)-jacobian of ξ

2 . 2 ] 1 ≤

 221 Dξ(y) = (Id ξ(y) -y -ξ(y) Dn(ξ(y))) -1 • π ξ(y) , where π ξ(y) is the orthogonal projection onto the plane tangent to ∂X at the point ξ(y), Id ξ(y) is the identity on the plane tangent to ∂X at the point ξ(y), and Dn is the differential of the normal map to ∂X. The same formula still holds if we replace ξ by its restriction ξ . The absolute value of the determinant of the restriction of π ξ(y) to the cell containing y is equal to |n(ξ(y))•n h (y)|. Furthermore, since the curvatures (that are the eigenvalues of Dn) are bounded by 1/R and y -ξ(y) | det((Id ξ(y) -y -ξ(y) Dn(ξ(y))) -1 )

Theorem 4 .

 4 Let X be a compact domain whose boundary has positive reach R. For h ≤ R √ d , the digital integral is multigrid convergent toward the integral over ∂X. More precisely, for any integrable function f : R d → R, one gets ∂X

Lemma 11 .

 11 One has the following properties • For every x ∈ ∂X \ C, the multiplicity µ x is less than µ:= d √ d + 1 . • For almost every point y ∈ ξ -1 (C) one has Jξ (y) = 0. • The area of C is equal to 0. Proof. Let x ∈ ∂X \ C and y ∈ ξ -1 (x). Then y belongs to the segment n(x, √ dh/2) centered in x, of length √ dh and aligned with the normal to ∂X at x. Since x / ∈ C, this segment is not contained in a plane orthogonal to n h (y). Since its length is less than √ dh, it cannot cross more than √ d + 1 cells of F d-1 h orthogonal to n h (y). The same bound holds for (d -1) other directions of the cells of F d-1 h . Hence µ x ≤ d √ d + 1 .

∂

  h X\mult(∂ h X) f (ξ(y))Jξ(y)dy = ∂ h X f (ξ(y))Jξ(y)dy + Area (∂X) µ f ∞ O(h).

( 21 )

 21 By Lemma 11 and the general coarea formula, one getsmult(∂ h X) f (ξ(y))Jξ(y)dy = mult(∂ h X)\ξ -1 (C) f (ξ(y))Jξ(y)dy = mult(∂X)\C

  ∂ h X f (ξ(y))Jξ(y)dy = ∂ h X f (ξ(y))|n(ξ(y)) • n h (y)|dy + f ∞ Area (∂X) O(h).(22) Lemma 9 implies that |Jξ (y) -|n(ξ(y)) • n h (y)| | = O(h).

  ∂ h X |f (ξ(y))| |Jξ (y) -|n(ξ(y)) • n h (y)| | dy ≤ f ∞ Area (∂ h X) O(h) ≤ f ∞ K 3 (h)Area (∂X) O(h).

f

  (ξ(y))|n(ξ(y)) • n h (y)|dy -DI h (f, D h (X), n) = Area (∂X) Lip(f )O(h) + f ∞ O( nn est ) .(23)

1 h

 1 ∂ h X f (ξ(y))|n(ξ(y)) • n h (y)|dy = c∈F d-∩∂ h X c f (ξ(y))|n(ξ(y)) • n h (y)|dy.For every square facec of ∂ h X (a d -1-cell of F d-1 h ), one has h d-1 f ( ċ)|n( ċ) • n( ċ)| = c f ( ċ)|n( ċ) • n( ċ)|dy.For every y ∈ c, one has n h (y) = n( ċ), and thenf (ξ(y))|n(ξ(y)) • n h (y)| -f ( ċ)|n( ċ) • n( ċ)| = f (ξ(y))|n(ξ(y)) • n( ċ)| -f ( ċ)|n(ξ(y)) • n( ċ)| + f ( ċ)|n(ξ(y)) • n( ċ)| -f ( ċ)|n( ċ) • n( ċ)| ≤ f (ξ(y)) -f ( ċ) n(ξ(y)) • n( ċ) + f ( ċ) |n(ξ(y)) • n( ċ)| -|n( ċ) • n( ċ)| ≤ f (ξ(y)) -f ( ċ) + f ( ċ) n(ξ(y)) -n( ċ) ≤ Lip(f ) ξ(y) -ċ + f ∞ nn est ≤ Lip(f ) √ dh + f ∞ nn est .Above, we use the relation that, for vectors a, b, u, ||a • u| -|b • u|| ≤ |(a -b) • u|. This relation comes from triangle inequalities. We deduce that (using alsoLemma 10) ∂ h X f (ξ(y))|n(ξ(y)) • n h (y)|dy -DI h (f, D h (X), n)

Figure 6 .

 6 Figure 6. Area estimation error of the digital surface integral (Definition 6) with several digital normal estimators. The shape of interest is 3D ellipsoid of half-axes 10, 10 and 5, for which the area has an analytical formula giving A ≈ 867.188270334505. The abscissa is the gridstep h at which the ellipsoid is sampled by Gauss digitization. For each normal estimator, the digital surface integral  is computed with f = 1, and the relative area estimation error | Â-A| A is displayed in logscale.

  that either s1 ∈ B in or s2 ∈ B out . A simple computation gives

		tan α =	3 √	52 26 + 52	√	2	.
	which implies that
	c =	4 2704 + (52 √ √ 26 2 + 3	√	26) 2

This work was partially supported by the ANR grants DigitalSnow ANR-11-BS02-009, KIDICO ANR-2010-BLAN-0205 and TopData ANR-13-BS01-0008. This is an author version of the paper accepted in JMIV in 2015.