Thomas Caissard
email: thomas.caissard@liris.cnrs.fr

David Coeurjolly
email: david.coeurjolly@liris.cnrs.fr

Jacques-Olivier Lachaud
email: jacques-olivier.lachaud@univ-savoie.fr

Tristan Roussillon
email: tristan.roussillon@liris.cnrs.fr

Laplace-Beltrami Operator on Digital Surfaces

Keywords: Laplace-Beltrami Operator, Digital Surface, Discrete Geometry, Differential Geometry

published or not. The documents may come

Introduction

Computer graphics, and particularly the field of geometry processing, revolves around studying discrete embedded surfaces (in many cases 2D surfaces in 3D). The Laplace-Beltrami operator (the Laplacian on a manifold) is a fundamental tool in geometry as it holds many properties of the surface. Eigenfunctions of the operator form a natural basis for square integrable functions on the manifold, in the same manner as Fourrier harmonics for functions on a circle. It is used for example as a basis for Functional Maps [START_REF] Ovsjanikov | Computing and processing correspondences with functional maps[END_REF] or mesh compression [START_REF] Lévy | Spectral Mesh Processing[END_REF]. Other applications are for example surface fairing, mesh smoothing, remeshing or feature extractions (see [START_REF] Lévy | Spectral Mesh Processing[END_REF]). The operator is also related to diffusion and the heat equation on a surface and connected to a large field of classical mathematics linking geometry of manifold to properties of the heat flow (see for example [START_REF] Rosenberg | The Laplacian on a Riemannian Manifold[END_REF]).

Many characterizations of discrete surfaces exist such as triangular, quadrangular meshes (or more generally simplicial complexes), points clouds, etc. Our model of surface comes from the Digital Geometry theory [START_REF] Klette | Digital geometry : geometric methods for digital picture analysis[END_REF], where the discrete structure is the topological boundary of a subset of points in Z d+1 called a digital surface (an example of this object is pictured in Fig. 1). Such surfaces can be constructed from mathematical modeling or from boundaries of partitions in volumetric images. Indeed, digital objects naturally arise in many material sciences or medical imaging applications as tomographic volumetric acquisition devices usually generate regularly spaced data (e.g [START_REF] Hildebrand | Direct three-dimensional morphometric analysis of human cancellous bone: microstructural data from spine, femur, iliac crest, and calcaneus[END_REF][START_REF] Flin | An adaptive filtering method to evaluate normal vectors and surface areas of 3d objects. application to snow images from x-ray tomography[END_REF]).

Our goal here is to present a discretization of the Laplace-Beltrami operator on digital surfaces which satisfies strong consistency (i.e. pointwise convergence) with respect to the Laplace-Beltrami operator on the underlying manifold when the digital surface is the boundary of the digitization of a continuous object. As we demonstrate in our experiments, previous works fail to efficiently estimate the Laplace-Beltrami operator on these specific surfaces. The main obstacle is the fact that normal vectors to these surfaces do not converge to the normal vectors of the underlying manifold, whatever the sampling rate.

We adapt the operator of Belkin et al. [START_REF] Belkin | Discrete laplace operator on meshed surfaces[END_REF] to our specific data. The method uses an accurate estimation of areas associated with digital surface elements. This estimation is achieved through a multigrid convergent digital normal estimator of Coeurjolly et al. [START_REF] Coeurjolly | Multigrid convergent principal curvature estimators in digital geometry[END_REF]. This paper is a direct followup on [START_REF] Caissard | Heat kernel Laplace-Beltrami operator on digital surfaces[END_REF] where we experimentally investigate applications such as heat diffusion or shape approximation through the eigenvectors decomposition. We show strong consistency of the discrete operator, and compare it experimentally with various other discretizations adapted on digital surfaces.

Discretization schemes overview. The Laplacian being a second order differential operator, a discrete calculus framework is required to define this operator on embedded combinatorial structures such as meshes or digital surfaces. The first elements of discrete calculus may be traced back to Regge calculus [START_REF] Regge | General relativity without coordinates[END_REF] for quantum physics, where discrete domains are modeled with adjacent tetrahedra and metrics are only determined by edge lengths. The discrete Laplacian has also been present in spectral analysis of graphs since the 1950s. Then, with the development of geometric acquisition devices and modeling techniques, interest grew toward a calculus working on meshes and more generally simplicial complexes. Early works include the widely studied cotangent formula [START_REF] Taubin | A signal processing approach to fair surface design[END_REF][START_REF] Taubin | Geometric signal processing on polygonal meshes[END_REF][START_REF] Desbrun | Implicit fairing of irregular meshes using diffusion and curvature flow[END_REF][START_REF] Polthier | Computational aspects of discrete minimal surfaces[END_REF][START_REF] Fujiwara | Eigenvalues of laplacians on a closed riemannian manifold and its nets[END_REF][START_REF] Mayer | Numerical solutions for the surface diffusion flow in three space dimensions 20[END_REF][START_REF] Meyer | Discrete Differential-Geometry Operators for Triangulated 2-Manifolds[END_REF][START_REF] Pinkall | Computing discrete minimal surfaces and their conjugates[END_REF] for various applications, which may be derived directly from standard finite element method (e.g. see [START_REF] Lévy | Spectral Mesh Processing[END_REF]).

Discrete exterior calculus was then developed in the computational mathematics and geometry processing community, with a particular focus on triangulated meshes. The "German school" of discrete calculus developed an exact 2D calculus which generalizes the cotangent Laplacian, and is based on (conforming and non-conforming) finite elements [START_REF] Polthier | Identifying vector field singularities using a discrete Hodge decomposition[END_REF], thus obtaining expected theoretical results such as Stokes' theorem and Hodge decomposition. Its applications range widely: exact integration allows accurate remeshing via L 2 projection, shape morphing by prescribing first-order data on the surface, etc. This theory provides a sound base for actual computation, with one important limitation: the necessity to only use triangles (and, furthermore, triangles with good aspect ratios, for positive Laplacian).

A more versatile expression of discrete exterior calculus comes with Hirani's thesis [START_REF] Hirani | Discrete exterior calculus[END_REF] and the monograph of Desbrun, Hirani, Leok and Marsden [START_REF] Desbrun | Discrete exterior calculus[END_REF]. Their primal-dual construction does not impose the use of triangular meshes. The discretization is not an approximation of the smooth calculus, but rather a discrete analog:

We do not prove that these definitions converge to the smooth counterparts. The definitions are chosen so as to make some important theorems like the generalized Stokes' theorem true by definition, to preserve naturality with respect to pullbacks, and to ensure that operators are local. [START_REF] Hirani | Discrete exterior calculus[END_REF][START_REF] Desbrun | Discrete exterior calculus[END_REF] Metrics play a role in musical operators (flat and sharp which convert vector field to k-forms and conversely) and Hodge stars. Note that discrete exterior calculus coincides with the cotangent scheme on triangular meshes when the Voronoi dual is used.

In parallel, another discrete calculus emerges in the image, graph, electric circuits and network analysis communities, summed up in Grady and Polimeni's book [START_REF] Grady | Discrete calculus: Applied analysis on graphs for computational science[END_REF]. Metrics are also incorporated, although without the relation with the ambient space. This feature was desired since people frequently wish to analyze data without any knowledge of an embedding. Authors then show how classical filtering procedures and (discrete versions of) energy models (e.g. Mumford-Shah, Total Variation) fit well within this framework.

A much-alike discrete calculus on "chainlets" appears in geometric measure theory, for the mathematical analysis of general compact shapes like fractals [START_REF] Harrison | Stokes' theorem for nonsmooth chains[END_REF][START_REF] Harrison | Flux across nonsmooth boundaries and fractal gauss/green/stokes' theorems[END_REF]. The exterior derivative, a Hodge star and Laplace-Beltrami are defined there for very general spaces. However computational aspects are unclear. We can also mention a complex analysis approach to discrete calculus for 2D digital surfaces [START_REF] Mercat | Discrete riemann surfaces and the ising model[END_REF][START_REF] Mercat | Discrete complex structure on surfel surfaces[END_REF] with applications to digital surface parametrization and texture mapping [START_REF] Cartade | Mesh parameterization with generalized discrete conformal maps[END_REF].

Operators for point clouds can be found in [START_REF] Belkin | Constructing Laplace Operator from Point Clouds in R d[END_REF] and more recently in [START_REF] Qin | Laplace-beltrami operator on point clouds based on anisotropic voronoi diagram[END_REF]. A discretization on polygonal surfaces was proposed by Alexa and Wardetzky [START_REF] Alexa | Discrete laplacians on general polygonal meshes[END_REF]. As digital surfaces being specific quadrangulated polygonal surface, such approach perfectly fits with our data. However, we show in the experiments such polygonal Laplace-Beltrami operator gives inconsistent results. Indeed, digital surface quads are axis-aligned quads and do not capture the metric of the underlying continuous object properly. Note that Alexa and Wardetzky polygonal operator matches with the cotangent one on triangular meshes. Other operators on polyhedral surfaces can be found in [START_REF] Hildebrandt | On approximation of the laplace-beltrami operator and the willmore energy of surfaces[END_REF][START_REF] Wardetzky | Discrete Laplace operators: No free lunch[END_REF][START_REF] Hildebrandt | On the convergence of metric and geometric properties of polyhedral surfaces[END_REF] (see discussion below).

Convergence and consistency of the operator. We clarify notions of convergences for operators. Suppose that you want to solve the equation Au = f where A : X → Y is a bounded linear operator between two Banach spaces and f ∈ Y is given.

Suppose also that you have an approximate A ε of A and f ε of f and that u ε is the solution of A ε u ε = f ε . (e.g. we can consider ε as the grid step for example). We say that the approximation scheme is convergent if

lim ε→0 ||u ε -u|| X = 0.
We say that the approximation scheme is consistent if

lim ε→0 ||A ε v -Av|| Y = 0, for each v ∈ X [34].
We focus here on consistency results of approximations of the Laplace-Beltrami. The choice of space X and Y determines various properties regarding solution of the equation. Let M be a compact manifold embedded in R d+1 with its topological border ∂ M of class C 2 . In order to apply the Laplace-Beltrami operator, we require, in this paper, X to be C 2 (∂ M), the space of twice differentiable functions acting on ∂ M. We also require Y to be C 0 (∂ M), the set of smooth functions acting on ∂ M. An important consequence which arises from classical analysis is the extreme value theorem (see [START_REF] Rudin | Principles of Mathematical Analysis. International series in pure and applied mathematics[END_REF] for example). This theorem states that if K is a compact set, and u : K → R is a continuous function (in the topological sense), then u is bounded and reaches its maximum and minumum in K. Therefore, as ∂ M is compact and each function of C 2 (∂ M) is continuous, our input functions are bounded. Their gradients, which are continuous, are also bounded. We chose the infinity Lebesgue norm L ∞ for the consistence (note that as our functions are bounded, the space is complete with this norm, thus a Banach space). By setting A = ∆ and A ε = ∆ ε , we say that an operator is strongly consistent (or pointwise convergent) whenever

lim ε→0 ||∆ ε v -∆ v|| L ∞ =lim ε→0 sup x∈∂ M |(∆ ε v)(x) -(∆ v)(x)| = 0, for all v ∈ C 2 (∂ M
). An operator called the Mesh Laplacian satisfying this property was proposed by Belkin et al. [START_REF] Belkin | Discrete laplace operator on meshed surfaces[END_REF] for interpolating triangular meshes and later extended to point clouds [START_REF] Belkin | Constructing Laplace Operator from Point Clouds in R d[END_REF]. Carl gave a discretization on Semi-Discrete surfaces and proves (among many things) the strong consistency of his operator [START_REF] Carl | A Laplace Operator on Semi-Discrete Surfaces[END_REF]. Another approach was proposed by Hildebrandt and Polthier [START_REF] Hildebrandt | On approximation of the laplace-beltrami operator and the willmore energy of surfaces[END_REF] and is valid on polyhedral surfaces. This seems to be the closest setting to our digital surfaces. Although, in our case, the projection function between the discrete surface and the underlying manifold is generally not bijective (see Section 5), while this is mandatory in their work. The strong consistency is often required when it comes to approximate curvatures, or Willmore energies [START_REF] Willmore | Riemannian Geometry[END_REF]. For such problems, Hildebrandt et al. [START_REF] Hildebrandt | Generalized shape operators on polyhedral surfaces[END_REF] derive a strongly consistent curvature estimators.

Other problems require only convergence of solutions of boundary value problems (e.g. Poisson's problems). In this case, we can relax the requirements for input functions. This setting is called the weak consistency and is related to the weak form of the Laplace-Beltrami operator (see Wardetzky's thesis for a proper definition [START_REF] Wardetzky | Discrete Laplace operators: No free lunch[END_REF]). For example, strong consistency for the cotangent operator holds for very specific meshes [START_REF] Xu | Convergence of discrete laplace-beltrami operators over surfaces[END_REF][START_REF] Xu | Discrete laplace-beltrami operators and their convergence[END_REF], but fails in general (counterexamples can be found in [START_REF] Xu | Convergence of discrete laplace-beltrami operators over surfaces[END_REF][START_REF] Hildebrandt | On the convergence of metric and geometric properties of polyhedral surfaces[END_REF]). Both convergence in the operator norm and weak consistency was established first by Dziuk [START_REF] Dziuk | Finite Elements for the Beltrami operator on arbitrary surfaces[END_REF] and later extended to polyhedral surfaces by Hildebrandt, Polthier and Wardetzky [START_REF] Hildebrandt | On the convergence of metric and geometric properties of polyhedral surfaces[END_REF]. In particular, they show an equivalence between normal convergence, metric convergence, convergence of area and convergence of the Laplace-Beltrami operators in the operator norm (see Theorem 2. of [START_REF] Hildebrandt | On the convergence of metric and geometric properties of polyhedral surfaces[END_REF]). They use these results to prove weak consistency of the operator (Theorem 4. of [START_REF] Hildebrandt | On the convergence of metric and geometric properties of polyhedral surfaces[END_REF]).

Finally, the spectrum of the mesh Laplacian converges thanks to Dey et al. [START_REF] Dey | Convergence, Stability, and Discrete Approximation of Laplace Spectra[END_REF]. As for the cotangent operator, convergence of eigenvalues has been established by Wardetzky [START_REF] Wardetzky | Discrete differential operators on polyhedral surfaces -convergence and approximation[END_REF], but convergence of the eigenvectors is still an open problem.

In all the paper, we use the term strong consistency (we may omit "strong") which corresponds to convergence in the L ∞ space (i.e. pointwise convergence).

Measure estimations. When it comes to operator discretizations, we need a way to compute approximate measures on the discrete surface. In a classical interpolating triangular mesh, a good approximation of the underlying smooth structure area is simply given by triangle areas whereas for digital surfaces naive measures (such as the quadrangular face area for 2D surface in 3D) give poor approximations. Although multigrid-convergent estimators of object global volume and area/perimeter have existed for a long time [START_REF] Ellis | Measurement of the lengths of digitized curved lines[END_REF][START_REF] Lenoir | Fast computation of the normal vector field of the surface of a 3-D discrete object[END_REF], local measure estimators (length of a 1-cell, area of a 2-cell, etc) have seen advances only in the past ten years. Parameter-free tangent and normal estimation along 2D digital curves were established by Lachaud et al. [START_REF] Lachaud | Fast, accurate and convergent tangent estimation on digital contours[END_REF][START_REF] De Vieilleville | Convex digital polygons, maximal digital straight segments and convergence of discrete geometric estimators[END_REF], using properties of maximal digital straight segments. Coeurjolly, Lachaud and Levallois have defined a digital variant of integral invariants [START_REF] Pottmann | Principal curvatures from the integral invariant viewpoint[END_REF][START_REF] Pottmann | Integral invariants for robust geometry processing[END_REF] which induces convergent estimation of the normal vector field along digital surfaces in arbitrary dimension as well as the whole curvature tensor [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]. Cuel, Lachaud and Thibert use a digital Voronoi Covariance Measure to show the convergence and stability of a normal estimator [START_REF] Cuel | Voronoi-based geometry estimator for 3d digital surfaces[END_REF]. Integration of normals was used in [START_REF] Lenoir | Des outils pour les surfaces discrètes : estimation d'invariants géométriques, préservation de la topologie, tracé de géodésiques et visualisation[END_REF] to estimate the perimeter of digital curves or the surface area of a digital surface. The fact that convergent normals implies convergent measures for subsets of codimension one in digital spaces was established by Lachaud and Thibert [START_REF] Lachaud | Properties of gauss digitized shapes and digital surface integration[END_REF]. Consequently, we can estimate convergent length of 1-cells in 2D and area of 2-cells in 3D, even locally.

Outline. We provide formal definitions for various discrete operators on triangular meshes in Section 2. We detail the definition of the heat kernel Laplace-Beltrami operator in Section 3 and provide hints on the convergence proof of Belkin et al. In Section 4 we derive formal concepts about Gauss digitization, cubical grid and digital surface. Then in Section 5 we give some properties of the projection map ξ that links the digital surface to its smooth counterpart: these are key tools when it comes to consistency of digital operators. We adapt in Section 6 our discrete Laplace operator from the one of Belkin et al. presented in Section 3. We prove the strong consistency in Section 7, using theorems from [START_REF] Belkin | Towards a theoretical foundation for laplacian-based manifold methods[END_REF] and [START_REF] Lachaud | Properties of gauss digitized shapes and digital surface integration[END_REF]. The main contributions of this paper are Theorem 5 and its associated proof. Finally, we show empirical consistency results and a comparison between the literature and our operator in Section 8.

Discretizations of the Laplace-Beltrami operator and their properties

We summarize various discretizations of the Laplace-Beltrami operator on triangular meshes. As stated in the introduction, let M be a compact manifold of dimension d + 1 with a smooth boundary. The linear operator ∆ :

C 2 (∂ M) → C 0 (∂ M) defined by ∆ u = div(grad u)
is called the Laplace-Beltrami operator (the sign of the operator is arbitrary, and one can find in the literature the alternative definition "-div(grad u)" for ∆ u). It is a bounded linear operator acting between Banach spaces.

Let Γ be a combinatorial structure (a triangular mesh for instance), V (Γ ) its set of vertices and F(Γ ) its faces. Let u : ∂ M → R be a twice differentiable function on ∂ M. We suppose that V (Γ ) is a sampling of M (i.e. V (Γ ) ⊂ ∂ M). In other words, the images of the function u are perfectly defined for all w ∈ V (Γ ).

The first simple discretization, coming from elementary calculus, is the graph Laplacian or combinatorial Laplacian acting on Γ [69]:

(L COMBI u)(w) := -deg(w)u(w) + ∑ p∈link 0 (w) u(p) , (1)
for all w ∈ V (Γ ) where link 0 (w) is the set of points in V (Γ ) adjacent to w and deg(w) is the degree of w in Γ . L COMBI is widely used in graph theory and machine learning for its nice properties [START_REF] Hein | From Graphs to Manifolds -Weak and Strong Pointwise Consistency of Graph Laplacians[END_REF].

Then, a similar approach, yet more complicated comes from the Discrete Exterior Calculus framework [START_REF] Hirani | Discrete exterior calculus[END_REF][START_REF] Desbrun | Discrete exterior calculus[END_REF]. Given an arbitrary embedded dual structure of Γ , the operator is written as a weighted double finite difference:

(L DEC u)(w) := 1 | w| ∑ p∈link 0 (w) | e wp | |e wp | (u(p) -u(w)), (2)
where is the Hodge-duality star operator acting on discrete forms (see [START_REF] Hirani | Discrete exterior calculus[END_REF]) and | • | the measure of a k-cell. As illustrated in Fig. 2, | e wp | would be the length of the segment orthogonal to e wp (e wp being the 1-cell corresponding to the edge between vertices w and p). If we set all measures to one, L DEC coincides with L COMBI . By fixing the dual of Γ to be the Voronoi diagram of its vertices and by computing the measures as Euclidean lengths and areas of this dual complex, the DEC operator coincides exactly with the famous cotan Laplacian [START_REF] Taubin | A signal processing approach to fair surface design[END_REF][START_REF] Taubin | Geometric signal processing on polygonal meshes[END_REF][START_REF] Desbrun | Implicit fairing of irregular meshes using diffusion and curvature flow[END_REF][START_REF] Polthier | Computational aspects of discrete minimal surfaces[END_REF][START_REF] Fujiwara | Eigenvalues of laplacians on a closed riemannian manifold and its nets[END_REF][START_REF] Mayer | Numerical solutions for the surface diffusion flow in three space dimensions 20[END_REF][START_REF] Meyer | Discrete Differential-Geometry Operators for Triangulated 2-Manifolds[END_REF][START_REF] Pinkall | Computing discrete minimal surfaces and their conjugates[END_REF]:

(L COT u)(w) := 1 2A w ∑ p∈link 0 (w) cot(α wp ) + cot(β wp ) (u(p) -u(w)) , (3)
where A w is one third of the area of all incident triangles to the vertex w, α wp and β wp are the angles opposing the corresponding edge e wp (see Fig. 2). The matrix representation of L COT , namely L COT is given by

L COT := D -1 Q, with D i,i := 1 3 ∑ t∈St(v i ) |t|
, where St(i) is the set of incident triangles to a vertex v i and |t| is the area of a triangle and

Q i, j := 1 2 (cot(α i j ) + cot(β i j )), Q i,i := -∑ j Q i, j .
As mentioned in the introduction, both L COT and L DEC are not strongly consistent in general [START_REF] Xu | Discrete laplace-beltrami operators and their convergence[END_REF][START_REF] Xu | Convergence analysis of a discretization scheme for gaussian curvature over triangular surfaces[END_REF][START_REF] Hildebrandt | On the convergence of metric and geometric properties of polyhedral surfaces[END_REF]]. Yet, because of the convergence of their weak form, they suffice for many geometry processing applications such as geodesic computation [START_REF] Crane | Geodesics in heat: a new approach to computing distance based on heat flow[END_REF], spectral processing [START_REF] Vallet | Spectral geometry processing with manifold harmonics[END_REF], etc. Apart from convergence behavior, these operators are fast to compute, and the resulting linear operator (i.e. the matrix) is sparse. Hildebrandt and Polthier in [START_REF] Hildebrandt | On approximation of the laplace-beltrami operator and the willmore energy of surfaces[END_REF] proposed a strongly consistent discretization of the operator. The idea is to test the cotan operator against a family of "r-local" functions:

ϕ r,v i := φr,v i || φr,v i || , φr,v i (v j ) := max{1 - ||v i -v j || R 3 r , 0}, for some r ∈ R * + . Then the matrix representation of the oper- ator L R-LOC is L R-LOC := ΦL COT , (4)
with Φ i, j := ϕ v i (v j ). The operator can be viewed as a convolution between r-local functions and L COT . Finally, we briefly talk about the discretization of Alexa and Wardetzky [START_REF] Alexa | Discrete laplacians on general polygonal meshes[END_REF]. Their operator, named L QUAD here, is defined on polyhedral surfaces. The trick is to properly define the adjoint operator d * on such structure by computing inner products on 0-forms and 1-forms (see the paper for the details). The operator coincides with L COT on triangular meshes and can be seen as a generalization of it (even though it is not its initial purpose). For L COT the area of integration A w is one third the area of all triangles incident on vertex w in green. For L DEC , we represent the Voronoi dual structure in orange. The dual of the edge e wp is in blue.

Heat kernel Laplace-Beltrami operator on triangular meshes

We detail the definition of the mesh Laplacian from [START_REF] Belkin | Discrete laplace operator on meshed surfaces[END_REF]. Although the fact that the Laplacian solves the heat equation is not new and has been studied for quite a while in differential geometry [START_REF] Rosenberg | The Laplacian on a Riemannian Manifold[END_REF], probability theory and quantum mechanics (as it is an "Euclidean" version of the Schrödinger equation), the discretization comes from the work of Belkin et al. [START_REF] Belkin | Discrete laplace operator on meshed surfaces[END_REF] who defined it on triangular meshes. Studies have also been made in [START_REF] Crane | Geodesics in heat: a new approach to computing distance based on heat flow[END_REF] in the context of geodesic distance approximation.

Let g : ∂ M × (0, T ) → R be a time-dependent function which solves the partial differential equation called the heat equation:

∆ g(x,t) = ∂ ∂t g(x,t), (5)
with initial condition u = g(•, 0) : ∂ M → R which is the initial temperature distribution. An exact solution is:

y∈∂ M p(t, x, y)u(y)dy, (6)
where

p ∈ C ∞ (R + × ∂ M × ∂ M) is the heat kernel [56].
There is a wide range of studies on the behavior of p when t tends to zero (called small-time asymptotic of diffusion process). Early work includes the famous Varadhan formula [START_REF] Varadhan | On the behavior of the fundamental solution of the heat equation with variable coefficients[END_REF] on closed manifolds with or without borders later extended by Molchanov [START_REF] Molchanov | Diffusion processes and riemannian geometry[END_REF] on a wider class of shapes:

p(t, x, y) ∼ t→0 e -d(x,y) 2 4t (4πt) d 2
where d(•, •) corresponds to the intrinsic geodesic distance. This approximation is not robust in practice and very sensitive to both geodesic distance approximation and numerical errors [START_REF] Crane | Geodesics in heat: a new approach to computing distance based on heat flow[END_REF]. Fortunately we know from Belkin et al. [START_REF] Belkin | Towards a theoretical foundation for laplacian-based manifold methods[END_REF] that in small-time asymptotic, the geodesic distance can be approximated by the Euclidean distance:

p(t, x, y) ∼ t→0 p(t, x, y) := e -||x-y|| 2 4t (4πt) d 2
, which leads to the following approximated solution of the heat equation:

g(x,t) = y∈∂ M p(t, x, y)u(y)dy. (7)
By injecting Eq.( 7) in the heat equation Eq.( 5) we obtain:

∆ g(x,t) = ∂ ∂t y∈∂ M p(t, x, y)u(y)dy. (8)
Using a finite difference on t, and the basic property that the integral of the heat kernel must be 1:

∆ g(x,t) = lim t→0 1 t y∈∂ M p(t, x, y)u(y)dy -u(x) (9a) = lim t→0 1 t y∈∂ M p(t, x, y)(u(y) -u(x))dy. (9b)
The previous equation can be seen as a convolution between differences of u and a time dependent Gaussian. Note that the derivation holds for any approximations of the heat kernel p. Following these derivations, the mesh Laplace operator [START_REF] Belkin | Discrete laplace operator on meshed surfaces[END_REF] on Γ (an interpolation of a 2D surface in 3D) is:

(L MESH u)(w) := 1 4πt 2 ∑ f ∈F(Γ ) A f 3 ∑ p∈V ( f ) e -||p-w|| 2 4t (u(p) -u(w)), (10)
where A f is the area associated with the face f . In [START_REF] Belkin | Discrete laplace operator on meshed surfaces[END_REF] authors show that L MESH converges towards the real Laplace-Beltrami operator ∆ as the triangulation interpolates denser and denser the manifold ∂ M. The proof involves the definition of an intermediate object called the functional Laplace operator that we recall in Definition 1.

Definition 1 (Functional Laplace operator [START_REF] Belkin | Discrete laplace operator on meshed surfaces[END_REF]) Given a point x ∈ ∂ M and a function u ∈ C 2 (∂ M), the functional Laplace operator is defined as follows:

(L h u)(x) := 1 t(4πt) d 2 y∈∂ M e -||y-x|| 2 4t (u(y) -u(x))dy,
where t h is a function in R * + tending to zero as h tends to zero.

Later, we will need Theorem 1 in our own convergence proof. It shows that the functional approximation converges toward the real Laplace operator on manifolds in small-time asymptotic.

Theorem 1 (Functional convergence, Lemma 5. of [START_REF] Belkin | Towards a theoretical foundation for laplacian-based manifold methods[END_REF])

Given a point x ∈ ∂ M, a function u ∈ C 2 (∂ M) : lim h→0 (L h u)(x) = (∆ u)(x).
In [START_REF] Belkin | Towards a theoretical foundation for laplacian-based manifold methods[END_REF], Belkin et al. show that for a particular family of triangulations interpolating the continuous manifold, their operator L MESH tends toward L h hence toward the real Laplace operator.

Digital surfaces and digital curves

Definitions of digital structures can be found in [START_REF] Lachaud | Properties of gauss digitized shapes and digital surface integration[END_REF][START_REF] Klette | Digital geometry : geometric methods for digital picture analysis[END_REF]. Topological aspects are described in [START_REF] Herman | Geometry of Digital Spaces[END_REF]. We consider as before a d + 1-manifold with a smooth rectifiable boundary embedded in R d+1 . We recall the definition of the Gauss Digitization process, which makes the link between M and its digital approximation: Definition 2 (Gauss digitization) Let h > 0 be the sampling grid step. The Gauss Digitization of a compact shape M ⊂ R d+1 is defined as D h (M) := M ∩ (hZ) d+1 where d is the dimension.

The digitization process has therefore a very simple scheme: it considers the discrete points of the infinite regular grid with grid step h and keeps only the ones inside the shape (see Fig. 3). We call D h (M) (or Z when we want an arbitrary object not derived from a Gauss digitization) a digital set. It is a subset of Z d+1 scaled by h by definition. In the next sections, we use two extra objects to represent the boundary of a digital approximation: the h-boundary (Fig. 3, middle) and its decomposition into cells of dimension d (Fig. 3,right). Let us first define the boundary of a digital set Z: for every point z in (hZ) d+1 (called digital points), we denote the d + 1-dimensional axis-aligned closed cube centered on z as Q h z and refer to it as an h-cube (their side length is h). We define then the h-cube embedding of a digital set Z as

Q h [Z] := ∪ z∈Z Q h z .
Definition 3 (h-boundary) The h-boundary of M, denoted by ∂ h M, is the topological boundary of the h-cube embedding of the Gauss digitization of M:

∂ h M := ∂ z∈D h (M) Q h z .
This set is represented in blue in Fig. 3. Now we construct the cubical grid associated to a digital set Z. We construct such set by Cartesian product of segments of dimension one as in [START_REF] Lachaud | Properties of gauss digitized shapes and digital surface integration[END_REF]. More precisely, we assign coordinates in ( h 2 Z) d+1 to each cell of the space. For each t ∈ h 2 Z, we associate the set

I h (t) such that if t ∈ hZ then I h (t) := [t -h 2 ;t + h 2 ] otherwise I h (t) := {t}. Now, if z ∈ ( h 2 Z) d+1 then I h (z) := I h (z 1 ) × • • • × I h (z d+1 )
, where z i is the i-th coordinate of z. The primal cubical grid is defined as follows:

Definition 4 (Primal cubical grid) The set F h := {I h (z)} z∈( h 2 Z) d+1
tiles the Euclidean space R d+1 with hypercubes and its faces. It is called the primal cubical grid at step h. Elements of F h are called cells. The set of elements of dimension d is denoted by F d h . As mentioned in the definition, a cubical grid tiles the entire space: for example F 2 h is the set of squares centered on the digital points of the grid (in green in Fig. 3). Therefore, when we want to select all the elements of a boundary (for example a digital curve in red on the right of Fig. 3), we take the intersection between the cubical grid F d h and the h-boundary

∂ h M (or ∂ Q h [Z]
when we have an arbitrary set of digital points).

Relationship between ∂ M and ∂ h M

We summarize properties described in [START_REF] Lachaud | Properties of gauss digitized shapes and digital surface integration[END_REF]. Associated proofs can be found in [START_REF] Lachaud | Non-Euclidean spaces and image analysis : Riemannian and discrete deformable models, discrete topology and geometry[END_REF][START_REF] Lachaud | Properties of gauss digitized shapes and digital surface integration[END_REF]. Topological or geometric inference regarding ∂ h M can be studied using a functional approach of the distance function to a compact set A and the related projection map. If A ⊂ R d+1 , the distance function δ A is the function on R d+1 such that

δ A (x) := inf{||x -a|| : a ∈ A}.
The R-offset of A, denoted by A R is the set of points whose distance to A is less than R. The medial axis Med(A) ⊂ R d+1 of A is the set of points with at least two closest points on A. The reach [START_REF] Federer | Curvature measures[END_REF]

D h (M), its h-cube embedding Q h [D h (M)
] in yellow and its h-boundary ∂ h M. The cubical grid F h 2 is displayed in green on the right, and also in the related "cubical border" F h 1 ∩ ∂ h M in red. The topological border of the h-cube embedding is used to push integral from ∂ M in the continuous setting to the discrete setting. Then, the cubical grid is used to split the integral on elements of various dimensions d, thus approximating the continuous sum by a discrete one.

Theorem 2 (Theorem 1 of [START_REF] Lachaud | Properties of gauss digitized shapes and digital surface integration[END_REF]) Let M be a compact domain of R d+1 such that the reach of ∂ M is greater than R. Then, for any digitization step

0 < h < 2R/ √ d + 1, the Hausdorff distance between sets ∂ M and ∂ h M is less than √ d + 1h/2. More precisely: ∀x ∈ ∂ M, ∃y ∈ ∂ h M,        ||x -y|| ≤ √ d + 1 2 h
and y ∈ n(x,

√ d + 1 2 h), ∀y ∈ ∂ h M, ||y -ξ (y)|| ≤ √ d + 1 2 h.
where n(x, a) is the segment of length 2a centered on x and aligned with the normal vector to ∂ M at x.

When studying the topology of ∂ h M through ξ , it has been shown that this function is not always bijective: more precisely, it is surjective everywhere, but non-injective on some subset of ∂ h M:

mult(∂ M) := {x ∈ ∂ M, s.t. ∃y 1 , y 2 ∈ ∂ h M, y 1 = y 2 , ξ (y 1 ) = ξ (y 2 ) = x}.
Fortunately, we know that the size of mult h (∂ M) is bounded by a quantity in O(h) (see Theorem 3. of [START_REF] Lachaud | Properties of gauss digitized shapes and digital surface integration[END_REF]). We define the digital surface integration as follows:

Definition 6 (Digital surface integration) Let Z be a digital set and h the grid step. Let f : R d+1 → R be an integrable function and n be a digital normal estimator. We define the digital surface integral by

DI h ( f , Z, n) = ∑ r∈F d h ∩∂ Q h [Z] h d f (ṙ)µ(r),
where ṙ is the centroid of the d-cell r and µ(r

) = | n(ṙ) • n(ṙ)|
the estimated area of a surfel r with n the trivial normal to the d-cell r.

The continuous sum is approximated by a discrete one over elements of dimension d. Given a cell r, we value the function on its centroid, and use an area approximation µ given by the scalar product between an estimated normal of r and the elementary normal orthogonal to r. This estimated area is called the measure of a cell; it is the area of the projected cell r onto the tangent plane induced by the estimated normal n (see Fig. 4), which has been used for a long time [START_REF] Ellis | Measurement of the lengths of digitized curved lines[END_REF][START_REF] Lenoir | Fast computation of the normal vector field of the surface of a 3-D discrete object[END_REF]. Normal vectors are estimated using the estimator presented in [START_REF] Coeurjolly | Multigrid convergent principal curvature estimators in digital geometry[END_REF][START_REF] Levallois | Parameter-free and multigrid convergent digital curvature estimators[END_REF], which has the multigrid convergence property. Note that summing this measure for each cell of the surface leads to an estimation of the global area of the shape boundary, which itself has a multigrid convergence property [START_REF] Lachaud | Properties of gauss digitized shapes and digital surface integration[END_REF]. This Fig. 4 Black dots represent the pointels, black segments are the linels. The measure µ(s) (in orange) of a surfel s (in green) is the area of the projection of s onto the tangent plane induced by the estimated normal.

cell measure is a key ingredient of the digital formalization of the integral leading to both theoretical and experimental multigrid convergence.

When taking Z to be the Gauss digitization of a compact shape with positive reach in Definition 6, theoretical convergence is given by T heorem 3.

Theorem 3 (Theorem 4. of [START_REF] Lachaud | Properties of gauss digitized shapes and digital surface integration[END_REF]) Let M be a compact domain whose boundary has positive reach R. For h ≤ R √ d+1 , the digital integral is multigrid convergent toward the integral over ∂ M. More precisely, for any integrable function f : R d+1 → R, one gets

∂ M f (x)dx -DI h ( f , D h (M), n) ≤ 2 d+3 (d + 1) 3 2 Area(∂ M) Lip( f ) √ d + 1 h + || f || ∞ • || n( ċ) -n( ċ)|| est , where || f || ∞ := max x∈R d+1 | f (x)| and Lip( f ) := max x =y | f (x) -f (y)|/||x -y|| 2 .
It involves the convergence speed of the normal estimator, but also bounds on the input function related to its || • || ∞ norm and its Lipschitz constant. Note that in our case, the positive reach is a consequence of the compactness of M and the smoothness of ∂ M. Steps for proving convergence involve for example showing that the integral of a quantity over mult h (∂ M) is negligible and then computing bounds on the remaining integral using various properties of the function ξ to link ∂ M and ∂ h M. We now carry on with the formal definition of our digital Laplace operator.

Heat kernel based Laplace-Beltrami operator on digital surfaces

We adapt the formulation of Belkin et al. on digital surfaces. In the continuous heat kernel formulation, the parameter t must tend to zero. On digital surfaces, we set t as a function of the grid step h, denoted t h , that tends to zero as h tends to zero. Section 8 clarifies such function. As stated in Theorem 2, the h-boundary ∂ h M is an O(h)-Hausdorff approximation of ∂ M (whereas the triangulated surface Γ is a sampling of ∂ M which is a stronger assumption). As a consequence, we need to map the smooth function u defined on ∂ M to ∂ h M: Definition 7 (Tube extension of u) Given a smooth function u on ∂ M, we define the extension ũ of u to R d+1 as ũ :

x →      u(x) if x ∈ ∂ M (u • ξ )(x) if x ∈ R d+1 \ ∂ M and δ ∂ M (x) ≤ √ d+1 2 h 0 otherwise
√ d+1 2 h-offset of ∂ M (namely (∂ M) √ d+1
2 h , that is the interior and the boundary of the tube, represented in yellow and blue on the figure), and has value 0 everywhere else.

where ξ is the map defined as before and δ ∂ M is the distance function. An illustration of this definition can be found in Fig. 5.

Applying the discretization scheme defined in Definition 6 to Eq.(9b), we derive a definition for our digital Laplace-Beltrami operator in Definition 8. Motivation for choosing Eq.(9b) over Eq.(9a) is theoretical: the continuous Lipschitz property of u (which is inherited from the bounded gradient property) will be applied to |u(y)u(x)| in our proof.

Definition 8 (Digital Laplace-Beltrami operator) Let Z be a digital set and h the grid step. Let f be some function defined at least in ∂ h M. The digital Laplace-Beltrami operator is:

(L h f )(s) := 1 t h (4πt h ) d 2 ∑ r∈F d h ∩∂ Q h [Z] e - ||ṙ-ṡ|| 2 4t h [ f (ṙ) -f (ṡ)]µ(r),
where ṙ (resp. ṡ) is the centroid of the surfel r (resp. s), µ(r) is equal to the dot product between an estimated normal and the trivial normal orthogonal to the surfel s and t h is a function of h tending to zero as h tends to zero, which will be specified later.

Strong consistency of the Digital Laplace-Beltrami operator

In the sequel, t h is defined as h α , for some α > 0. The positiveness of α is given by Theorem 1 as t h must tend toward zero as h tends toward zero. We also assume that || n( ċ)n( ċ)|| est , the error on the normal vector estimation, to be in O(h β ). The speed of convergence (i.e. the value of β ) depends on the estimator [START_REF] Coeurjolly | Multigrid Convergence of Discrete Geometric Estimators[END_REF]. For example the convergence speed of the integral invariant normal estimator [START_REF] Coeurjolly | Multigrid convergent principal curvature estimators in digital geometry[END_REF][START_REF] Levallois | Parameter-free and multigrid convergent digital curvature estimators[END_REF] is O(h 3 ). We prove the strong consistency of our operator when considering the digital set Z to be D h (M), the Gauss discretization of M. Let s be a surfel in

F d h ∩ ∂ h M. We show that (∆ u)(ξ (ṡ)) -(L h ũ)(ṡ) (11)
tends toward zero as the grid step h tends toward zero. First, let us use the triangle inequality in Eq.( 11) to highlight the important steps of the proof:

(∆ u)(ξ (ṡ)) -(L h ũ)(ṡ) ≤ (∆ u)(ξ (ṡ)) -(L h u)(ξ (ṡ)) (Q1) + (L h u)(ξ (ṡ)) -(L h ũ)(ṡ) (Q2) + (L h ũ)(ṡ) -(L h ũ)(ṡ) (Q3)
The main contribution of this article is stated in both Lemma 1 and Lemma 2. Lemma 1 from Section 7.1 gives a bound on (Q2). Lemma 2 from Section 7.2 gives a bound on (Q3) and is a combination of Lemma 3 and Theorem 3. The consistency of (Q1) is given by Theorem 1. The overall consistency is given in Section 7.3 through Theorem 5. We do not provide a convergence speed for the entire result as the convergence speed of (Q1) is unknown (Theorem 1 just indicates convergence).

Bound on (Q2)

Lemma 1 Let s ∈ F d h ∩ ∂ h M, a function u ∈ C 2 (∂ M) and its extension ũ from Definition 7. For t h = h α , 0 < α ≤ 2 2+d
and h ≤ h max with h max the minimum between Diam(∂ M), K 3 (d, α, Diam(∂ M)) (see Eq.( 15) for an explicit value) and R/ √ d + 1, we have

|(L h u)(ξ (ṡ)) -(L h ũ)(ṡ)| ≤ Area(∂ M) ||∇u|| ∞ K 1 (d)h 1-α(1+ d 2 ) + K 2 (d)h 2-α 3+d 2 with K 1 (d) := √ d + 1 2 d-1 eπ d 2 and K 2 (d) := 3(d + 1) 2 d+ 5 2 √ eπ d 2
.

Proof Using Definition 1, we have:

(L h u)(ξ (ṡ)) -(L h ũ)(ṡ) = 1 t h (4πt h ) d 2 y∈∂ M e - ||y-ξ (ṡ)|| 2 4t h [u(y) -u(ξ (ṡ))] - y∈∂ M e - ||y-ṡ|| 2 4t h [ ũ(y) -ũ(ṡ)] . Since [u(y) -u(ξ (ṡ))] = [ ũ(y) -ũ(ṡ)],
we factorize by [u(y)u(ξ (ṡ))] to get:

|(L h u)(ξ (ṡ)) -(L h ũ)(ṡ)| = 1 t h (4πt h ) d 2 y∈∂ M e - ||y-ṡ|| 2 4t h -e - ||y-ξ (ṡ)|| 2 4t h [u(y) -u(ξ (ṡ))] dy . ( 13
)
As the infinity norm of the gradient of u is bounded (from the extreme value theorem), we know that u is Lipschitz continuous with constant equals to ||∇u|| ∞ . Therefore Eq.( 13) becomes: Let us first find the upper bound. We use the squared triangle inequality in (ṡ, ξ (ṡ), y) (see Fig. 6):

|(L h u)(ξ (ṡ)) -(L h ũ)(ṡ)| = 1 t h (4πt h ) d 2 y∈∂ M e - ||y-ṡ|| 2 4t h -e - ||y-ξ (ṡ)|| 2 4t h [u(y) -u(ξ (ṡ))] dy ≤ 1 t h (4πt h ) d 2 y∈∂ M e - ||y-ṡ|| 2 4t h -e - ||y-ξ (ṡ)|| 2 4t h • |u(y) -u(ξ (ṡ))| dy ≤ ||∇u|| ∞ t h (4πt h )
||y -ṡ|| 2 ≥ ||y -ξ (ṡ)|| 2 -||ṡ -ξ (ṡ)|| 2 -2 ||y -ṡ|| ||ṡ -ξ (ṡ)||.
Using the fact that ||y -ṡ|| ≤ ||yξ (ṡ)|| + ||ṡξ (ṡ)|| we have:

||y -ṡ|| 2 ≥ ||y -ξ (ṡ)|| 2 -||ṡ -ξ (ṡ)|| 2 -2 (||y -ξ (ṡ)|| + ||ṡ -ξ (ṡ)||) ||ṡ -ξ (ṡ)|| = ||y -ξ (ṡ)|| 2 -3||ṡ -ξ (ṡ)|| 2 -2 ||y -ξ (ṡ)|| ||ṡ -ξ (ṡ)||.
Then we apply the Hausdorff property of ∂ h M with respect to ∂ M to the term ||ṡξ (ṡ)|| (Theorem 2):

||y -ṡ|| 2 ≥ ||y -ξ (ṡ)|| 2 - 3(d + 1) 4 h 2 -||y -ξ (ṡ)|| √ d + 1h.
h 2 + ||y-ξ (ṡ)|| √ d+1 4t h h -1 .
For 0 ≤ x ≤ 2.51286 the inequality e x/2 -1 ≤ x is true (see Lemma 5 in the appendix). We apply it to the last equation:

b ≤ e - ||y-ξ (ṡ)|| 2 4t h 3(d + 1) 8t h h 2 + ||y -ξ (ṡ)|| √ d + 1 2t h h =: c.
The bound c is valid when 0 ≤ x ≤ 2.51286. Replacing x by 3(d+1) 8 h 2-α + ||y-ξ (ṡ)|| √ d+1 2 h 1-α and supposing that h ≤ Diam(∂ M) (which is reasonable in our context) we have:

3(d + 1) 8 h 2-α + ||y -ξ (ṡ)|| √ d + 1 2 h 1-α ≤ 2.51286 ⇐= 3(d + 1) 8 h 2-α + Diam(∂ M) √ d + 1 2 h 1-α ≤ 2.51286 ⇐= h 2-α 3(d + 1) 8 + √ d + 1 2 ≤ 2.51286 ⇐= h ≤ exp 3.00086 -log(3(d + 1) + 4 √ d + 1) 2 -α ) := K 3 (d, α). ( 15
)
The next step is to find a negative minor bound on b. We use a triangle inequality in (ṡ, ξ (ṡ), y) (see Fig. 6) and again the Hausdorff property of ∂ h M:

||y -ṡ|| 2 ≤ ||y -ξ (ṡ)|| 2 + ||ṡ -ξ (ṡ)|| 2 + 2 ||y -ξ (ṡ)|| ||ṡ -ξ (ṡ)|| ≤ ||y -ξ (ṡ)|| 2 + d + 1 4 h 2 + ||y -ξ (ṡ)|| √ d + 1 h.
Using the same derivation as for the upper bound we have

b ≥ e - ||y-ξ (ṡ)|| 2 4t h e -d+1 16t h h 2 - ||y-ξ (ṡ)|| √ d+1 4t h h -1 =: a.
Using Lemma 3 we have |b| ≤ max{|a|, c}. We now bound the absolute value of a:

|a| = e - ||y-ξ (ṡ)|| 2 4t h e -d+1 16t h h 2 - ||y-ξ (ṡ)|| √ d+1 4t h h -1 = e - ||y-ξ (ṡ)|| 2 4t h 1 -e -d+1 16t h h 2 - ||y-ξ (ṡ)|| √ d+1 4t h h . For 0 ≤ x, 1 -e -x ≤ x holds which leads to |a| ≤ e - ||y-ξ (ṡ)|| 2 4t h d + 1 16t h h 2 + ||y -ξ (ṡ)|| √ d + 1 4t h h .
We see that max{|a|, c} = c and using Lemma 3:

|b| = e - ||y-ṡ|| 2 4t h -e - ||y-ξ (ṡ)|| 2 4t h ≤ e - ||y-ξ (ṡ)|| 2 4t h 3(d + 1) 8t h h 2 + ||y -ξ (ṡ)|| √ d + 1 2t h h . ( 16
)
Injecting Eq.( 16) in Eq.( 14) we have

(L h u)(ξ (ṡ)) -(L h ũ)(ṡ) ≤ ||∇u|| ∞ t h (4πt h ) d 2 • 3(d + 1) 8t h h 2 y∈∂ M ||y -ξ (ṡ)||e - ||y-ξ (ṡ)|| 2 4t h dy + ||∇u|| ∞ t h (4πt h ) d 2 • √ d + 1 2t h h y∈∂ M ||y -ξ (ṡ)|| 2 e - ||y-ξ (ṡ)|| 2 4t h dy = 3(d + 1)||∇u|| ∞ 8(4πt h ) d 2 h 2 y∈∂ M ||y -ξ (ṡ)|| t 2 h e - ||y-ξ (ṡ)|| 2 4t h dy + ||∇u|| ∞ √ d + 1 2(4πt h ) d 2 h y∈∂ M ||y -ξ (ṡ)|| 2 t 2 h e - ||y-ξ (ṡ)|| 2 4t h dy
To bound the first integral, we use the inequality xe

-x 2 ≤ 1/ √ 2e for all x ∈ R. Putting x = ||y-ξ (ṡ)|| 2 √ t h we have ||y -ξ (ṡ)|| t 2 h e - ||y-ξ (ṡ)|| 2 4t h ≤ √ 2 √ e • t 3 2 h .
For the second one, we know that for all x ∈ R, x 2 e -x Continuing:

|(L h u)(ξ (ṡ)) -(L h ũ)(ṡ)| ≤   3 √ 2(d + 1)||∇u|| ∞ 8 √ e(4π) d 2 • h 2 t 3+d 2 h + 2 ||∇u|| ∞ √ d + 1 e(4π) d 2 h t 1+ d 2 h   y∈∂ M dy =   3 √ 2(d + 1)||∇u|| ∞ 8 √ e(4π) d 2 • h 2 t 3+d 2 h + 2 ||∇u|| ∞ √ d + 1 e(4π) d 2 h t 1+ d 2 h   Area(∂ M).
We replace t h by h α :

|(L h u)(ξ (ṡ)) -(L h ũ)(ṡ)| ≤ Area(∂ M) 3 √ 2(d + 1)||∇u|| ∞ 8 √ e(4π) d 2 • h 2-α 3+d 2 + 2 ||∇u|| ∞ √ d + 1 e(4π) d 2 h 1-α(1+ d 2 ) = Area(∂ M) ||∇u|| ∞ 3 • (d + 1) 2 d+ 5 2 √ eπ d 2 h 2-α 3+d 2 + √ d + 1 d d-1 eπ d 2 h 1-α(1+ d 2 ) .
The convergence holds when h's exponents are positive, that is 2α 3+d 2 > 0 and 1α(1 + d 2 ) > 0. Simple analysis leads to α ≤ 2/(2 + d).

Bound on (Q3)

This section proves the following lemma: Lemma 2 Let the normal estimator convergence speed be in O(h β ) and let t h = h α . For h

≤ h 0 = R/ √ d + 1, L h is strongly consistent with L h if 0 < α < min 2 d + 2 , 2β d + 1 with a convergence speed in C h 1-α( d 2 +1) + O(h β -α 1+d 2 ),
where

C := 48(d + 1) 2 π d 2 Area(∂ M)||∇u|| ∞ .
We only explicit the constant for the first term of the convergence as the second term is related to an arbitrary normal estimator. We introduce the function g h ṡ : R d+1 → R defined as

∀x ∈ R d+1 , g h ṡ (x) := 1 t h (4πt h ) d 2 e - ||x-ṡ|| 2 4t h ( ũ(x) -ũ(ṡ)). ( 17
)
Lemma 2 consists in showing the convergence of the digital approximation of the integral of g h ṡ over ∂ M. More precisely, we want:

lim h→0 x∈∂ M g h ṡ (x)dx -∑ r∈F d h ∩∂ h M g h ṡ (ṙ)µ(r) = 0.
Using Theorem 3 we have

x∈∂ M g h ṡ (x)dx -∑ r∈F d h ∩∂ h M g h ṡ (ṙ)µ(r) ≤ 2 d+3 (d + 1) 3 2 Area(∂ M) Lip(g h ṡ ) √ d + 1 h + ||g h ṡ || ∞ O(h β ) . (18)
The challenge now is to carefully bound from above Lip(g h ṡ ) and ||g h ṡ || ∞ , which depends on h. We first show an upper bound on Lip( ũ) using properties of ξ in Section 7.2.1. Next, we show bounds on Lip(g h ṡ ) and ||g h ṡ || ∞ within the tube extension of ∂ M in Section 7.2.2. Finally, by using the definition of g h ṡ we extend the two previous bounds to R d+1 which proves Lemma 2.

Bound on Lip( ũ)

As ũ is defined to be the composition between u ∈ C 2 (∂ M) and the projection function ξ , we state first Theorem 4 showing a bound on the Lipschitz constant associated to ξ .

Theorem 4 (Proposition 1 of [START_REF] Lachaud | Properties of gauss digitized shapes and digital surface integration[END_REF] and theorem 4.8 of [START_REF] Federer | Curvature measures[END_REF]) Let A be a compact set with positive reach. Then for every p ∈ A and every ι ∈ [0, 1[, the projection ξ A is 1 1-ι -Lipschitz in the ball centered on p with radius ι × reach(A).

To find an upper bound on Lip( ũ) we bound the gradient of ũ. To do so, we need ũ to be differentiable. Unfortunately, derivatives of ũ are defined everywhere except on the boundary of the tube extension of ∂ M. Therefore, we restrict for now the analysis within the tube extension (i.e. the offset)

T := (∂ M) √ d+1 2 h
(see Fig. 5). We write for an arbitrary function f

Lip T f := max x,y∈T ,x =y | f (x) -f (y)|/||x -y|| 2 .
We know that the Lipschitz constant is bounded by the maximal (vector) norm of the gradient:

Lip T ( ũ) ≤ max x∈T (∇ ũ)(x) ∞ = max x∈T (∇u • ξ )(x) ∞ .
Using the chain rule property of the gradient, we have:

Lip T ( ũ) ≤ max x∈T (Jξ (x)) T (∇u)(ξ (x)) ∞ , (19)
where Jξ (x) is the Jacobian of ξ at point x and (∇u)(ξ (x)) is the application of the gradient of u to the point ξ (x). For h ≤ R/ √ d + 1 (with R the reach of ∂ M) we know from Theorem 4 that ξ is 2-Lipschitz in a ball of radius R/2 (here we chose ι to be 1/2 which gives an upper bound on h, and the 2-Lipschitz property). Therefore, the transposed Jacobian is bounded by 2 (as each of the derivatives are bounded) which leads to:

Lip T ( ũ) ≤ 2 max x∈T (∇u)(ξ (x)) ∞ (20)
Using the fact that ξ is surjective everywhere, Eq.( 20) becomes

Lip T ( ũ) ≤ 2 max x∈T (∇u)(ξ (x)) ∞ = 2 ∇u ∞ . (21)
Indeed, as every points of ∂ M has a pre-image in (∂ M)

√ d+1 2 h
with respect to the function ξ the infinity norm reaches the same value.

Bounds for Lip

T (g h ṡ ) and ||g h ṡ || T
We use the following shorthand notation for this section

||∇ ũ|| T := max x∈T (∇ ũ)(x) ∞ .
We show the following Lemma:

Lemma 3 For h ≤ h 0 = R/ √ d + 1 we have ||g h ṡ || T ≤ 2 • ||∇u|| ∞ (4π) d 2 • h -α 1+d 2 and Lip T (g h ṡ ) ≤ 6 • ||∇u|| ∞ (4π) d 2 • h -α(1+ d 2 )
.

Proof If we consider Eq.( 17) and by using the Lipschitz property of ũ we have

||g h ṡ || T = max x∈T    e - ||x-ṡ|| 2 4t h t h (4πt h ) d 2 • ( ũ(x) -ũ(ṡ))    ≤ max x∈T    e - ||x-ṡ|| 2 4t h t h (4πt h ) d 2 • Lip T ( ũ) • ||x -ṡ|| 2    .
We know that for all y ∈ R, ye -y 2 ≤ 1/2. Putting y = ||x-ṡ||

2 √ t h we have ||x -ṡ|| t h e - ||x-ṡ|| 2 4t h ≤ 1 √ t h .
Using Eq.( 21) and then replacing t h by h α gives

||g h ṡ || T ≤ 2 • ||∇u|| ∞ t 1 2 h (4πt h ) d 2 = 2 • ||∇u|| ∞ (4π) d 2 • h -α 1+d 2 .
Next, since Lip T (g h ṡ ) ≤ ||∇g h ṡ || T we first compute the derivative of g h ṡ and then its upper bound. We have:

∇g h ṡ (x) = ∇ ũ(x) t h (4πt h ) d 2 e - ||x-ṡ|| 2 4t h + ũ(x) -ũ(ṡ) t h (4πt h ) d 2 ∇e - ||x-ṡ|| 2 4t h .
Using elementary calculus, we know that ∇e

- ||x-ṡ|| 2 4t h = - ||x -ṡ|| 2t h e - ||x-ṡ|| 2 4t h .
Therefore we have

∇g h ṡ (x) = ∇ ũ(x) t h (4πt h ) d 2 e - ||x-ṡ|| 2 4t h - ||x -ṡ|| [ ũ(x) -ũ(ṡ)] 2t 2 h (4πt h ) d 2 e - ||x-ṡ|| 2 4t h .
We use the triangle inequality to continue:

∇g h ṡ T ≤ ∇ ũ T t h (4πt h ) d 2 e - x-ṡ 2 4t h + max x∈T | ũ(x) -ũ(ṡ)| x -ṡ 2 2t 2 h (4πt h ) d 2 e - ||x-ṡ|| 2 4t h ≤ ||∇ ũ|| T t h (4πt h ) d 2 e - ||x-ṡ|| 2 4t h + Lip T ( ũ) • max x∈T ||x -ṡ|| 2 2t 2 h (4πt h ) d 2 e - ||x-ṡ|| 2 4t h .
As first mentioned in Lemma 1, elementary calculus shows that ∀y ∈ R, y 2 e -y 2 ≤ 1 (with use here a weaker bound for factorization purpose). Taking y to be ||x-ṡ||

2 √ t h we have ||x -ṡ|| 2 t 2 h e - ||x-ṡ|| 2 4t h ≤ 4 t h (22)
which gives the final bound on ||∇g h ṡ || T :

∇g h ṡ T ≤ ||∇ ũ|| T t h (4πt h ) d 2 + Lip T ( ũ) • max x∈T ||x -ṡ|| 2 2t 2 h (4πt h ) d 2 e - ||x-ṡ|| 2 4t h ≤ ||∇ ũ|| T t h (4πt h ) d 2 + 2 Lip T ( ũ) t h (4πt h ) d 2
(Using Eq.( 22))

≤ ||∇ ũ|| T t h (4πt h ) d 2 + 4||∇u|| ∞ t h (4πt h ) d 2
(Using Eq.( 21)).

Combining Eq.( 19) and Eq.( 21) we have

||∇ ũ|| T ≤ 2 ||∇u|| ∞ .
Injecting this quantity into the last equation gives:

∇g h ṡ T ≤ 2||∇u|| ∞ t h (4πt h ) d 2 + 4||∇u|| ∞ t h (4πt h ) d 2 = 6||∇u|| ∞ (4π) d 2 • h -α(1+ d 2 ) (putting t h = h α ),
which proves Lemma 3.

Conclusion on (Q3) convergence

Injecting Lemma 3 into Eq.( 18) we have:

x∈∂ M g h ṡ (x)dx -∑ r∈F d h ∩∂ h M g h ṡ (ṙ)µ(r) ≤ 2 d+3 (d + 1) 3 2 Area(∂ M) Lip(g h ṡ ) √ d + 1 h + ||g h ṡ || ∞ O(h β ) ≤ 2 d+3 (d + 1) 3 2 Area(∂ M) 6 • ||∇u|| ∞ (4π) d 2 • h -α(1+ d 2 ) √ d + 1 h + 2 • ||∇u|| ∞ (4π) d 2 • h -α 1+d 2 O(h β ) = 48(d + 1) 2 π d 2 Area(∂ M) ||∇u|| ∞ h 1-α( d 2 +1) + O(h β -α 1+d 2 ).
This result can be extended to the whole space (i.e. replacing T by ∞ in the norm) using the fact that Theorem 3 only needs bounded Lipschitz and the l ∞ error in T (as all proofs rely on bounds computed within this tube because of the Hausdorff property of ∂ h M). Therefore we write:

x∈∂ M g h ṡ (x)dx -∑ r∈F d h ∩∂ h M g h ṡ (ṙ)µ(r) ≤ 48(d + 1) 2 π d 2 Area(∂ M) ||∇u|| ∞ h 1-α( d 2 +1) + O(h β -α 1+d 2 ).
which proves Lemma 2.

Overall convergence result

Theorem 5 Let s be a surfel in F d h ∩ ∂ h M, a function u ∈ C 2 (∂ M) and its extension ũ from Definition 7. Let t h = h α and let the convergence speed of the normal estimator be in O(h β ). Let h 0 be the minimum between Diam(∂ M), R/ √ d + 1 and K 3 (d, α, Diam(∂ M)) (where K 3 is a constant defined in Eq.( 15)). For 0 < h ≤ h 0 we have

lim h→0 (∆ u)(ξ (ṡ)) -(L h ũ)(ṡ) = 0 if 0 < α < min 2 d+2 , 2β d+1 .
Proof We remind the reader that the condition α > 0 is given by Theorem 1 as t h must tends toward zero as h tends toward zero. Then Lemma 1 gives the following condition for the consistency:

α < 2 d + 2 ,
and Lemma 2 the same condition plus

α < 2β d + 1 .
Combining these conditions on α, Theorem 1 from Belkin et al. for the convergence of Eq.(Q1), Lemma 1 for the bound on Eq.(Q2) and Lemma 2 for the bound on Eq.(Q3), the theorem holds. Hence our digital operator is strongly consistent with the smooth Laplace-Beltrami operator on manifolds. The upper bound h 0 for h is given by Lemma 1 and Lemma 2.

Experiments

We only investigate in this section the empirical consistency property under the l ∞ norm. Associated geometry processing applications can be found in [START_REF] Caissard | Heat kernel Laplace-Beltrami operator on digital surfaces[END_REF]. We consider a unit ball S 3 and three different smooth functions u : ∂ S 3 → R, namely z, x 2 and e x (see. Fig. 8). Let θ be the azimuth angle, and φ the polar angle. The spherical Laplacian is then:

∆ ∂ S 3 u(θ , φ ) = 1 sin 2 φ ∂ 2 u ∂ θ 2 + 1 sin φ ∂ ∂ φ sin φ ∂ u ∂ φ . ( 23
)
We compute the Gauss digitization D h (S 3 ) of the ball for decreasing grid steps h. We take a ball of Euclidean radius one. For h = 0.1, the digital surface ∂ h S 3 has 1902 surfels, whereas for h = 0.01 it contains 188502 surfels. Since the elements of ∂ h S 3 do not interpolate the sphere, u is extended to ũ as defined in Definition 7. We use the normal vector estimator described in [START_REF] Coeurjolly | Multigrid convergent principal curvature estimators in digital geometry[END_REF] to compute µ the measure of the surfels. All tests are made using the DGtal library [START_REF]DGtal: Digital geometry tools and algorithms library[END_REF] written in C++. We first investigate convergence speed for various α values. Then we compare our operator with ones adapted to digital surface.

Consistency results for various α

We evaluate the consistency property of our Laplace-Beltrami operator. We display graphs in Fig. 8 for various parameters t h : h (in red ),

h 2 3 (in blue ), h 1 3
(in green ), h 1 6 (in purple ) and h 1 12 (in orange ). Convergence speeds are summed up in Table .1. Our operator L h can be seen as a convolution between a Gaussian of standard deviation σ = √ 2t h and differences of functions. As the discretization becomes finer, the standard deviation σ of the Gaussian in the convolution decreases and the number of points within it increases.

We observe that, as the exponent α of t h increases, the l ∞error decreases. Moreover, for both x 2 and exp(x) functions, the speed of convergence increases alongside the value of α. Although theoretical convergence speed is achieved for t h = h 1 3 , the empirical behavior is the opposite of the one given by Lemma 1 and Lemma 2 where as α tends toward 0 the convergence speed increases. We strongly think that this difference between the theory and the application is related with Theorem 1 where the theoretical convergence speed is not explicit. Furthermore, the sphere is a very specific shape and may also bias the result.

Comparison against other discretizations

We compare in Fig. 9 convergence speeds of various operators. We compute L COMBI (from Eq.( 1)), L QUAD and L h directly on ∂ h S 3 ; L COT (from Eq.( 3)), L MESH (from Eq.( 10)) and L R-LOC (from Eq.( 4)) on the associated marching cubes triangulation. Since the vertices of this mesh coincide with the centroids of the surfels of ∂ h S 3 , all these operators are evaluated at the same points. The parameter t for L MESH should depend on the sampling of the triangulated mesh [START_REF] Belkin | Discrete laplace operator on meshed surfaces[END_REF].

In our case, as the mesh comes from our digital surface, we use the parameter t h from L h in L MESH and chose t h = h 1 3 . As for the parameter r for L R-LOC we set it to be h For comparison, in order to mimic the setting of [START_REF] Belkin | Discrete laplace operator on meshed surfaces[END_REF], we have also considered the Laplacian L P MESH , which corresponds to L MESH when the vertices of the marching cubes are projected onto the sphere. In our framework, this operator is the gold standard as we perfectly know the vertex positions. Finally, we also compute L P COT , L P R-LOC and L P QUAD on the projection of the marching cubes. The various discrete spheres are depicted in Fig. 7. First, as theoretically expected, L COMBI ( ), L COT ( ) and L P COT ( ) are not strongly consistent on our setting on both the marching cubes and its projection. The polygonal Laplace of Alexa et al. on our surface (resp. projected surface) has the same behavior as L COT (resp. L P COT ). As for the cotan operator, the poor approximation of the tangent space through trivial normals leads to non-convergent areas, and thus a diverging operator on the sphere. As for the Table 1 Summary of convergence speed for various values of t h and different functions ũ. The model h γ has been fitted to the maximal error in the Laplacian evaluation (the same data set is used in Fig. 8). The table shows the values of the model parameter γ with respect to the input parameter α and the function ũ. The higher the values, the speeder the convergence. 1.0611 0.0861 0.0852 projection, even though the area converges, many triangles are ill-formed (as a result of the projection) thus breaking the consistency. Non-consistency is also observed for L MESH ( ) but with lower errors. On the opposite, L P MESH ( ) shows pointwise convergence, as expected in [START_REF] Belkin | Discrete laplace operator on meshed surfaces[END_REF]. Convergence speeds are close between L h and L P MESH , although there is a 10 2 gap for the u(x, y, z) = z function. L R-LOC shows non-consistency behaviors. Indeed, our setting does not fit in the theorem of [START_REF] Hildebrandt | On approximation of the laplace-beltrami operator and the willmore energy of surfaces[END_REF]: our projection function ξ is not bijective in general, and the surface normals of the marching cubes are not convergent breaking the weak consistency of L COT (see [START_REF] Wardetzky | Discrete Laplace operators: No free lunch[END_REF]) which is required in the strong consistency proof of L R-LOC . As for L P R-LOC , the projected marching cubes has the properties required by Theorem 7. of [START_REF] Hildebrandt | On approximation of the laplace-beltrami operator and the willmore energy of surfaces[END_REF], but the shape regularity ρ (which corresponds to the aspect ratio of the triangles) diverges on the projection, thus breaking the constant C from their theorem.

α ũ(x, y, z) = z ũ(x, y, z) = x 2 ũ(x, y, z) = e x 1 0
u(x, y, z) = e x h t h = h t h = h 2 3 t h = h 1 3 t h = h 1 6 t h = h 1 12
Fig. 8 Consistency results between ∆ and L h for the l ∞ norm are shown for various functions on the unit sphere ∂ S 3 . We use either t h = h in red,

t h = h 2 3 in blue, t h = h 1 3 in green, t h = h

Mean curvature estimation

We motivate here our strong consistency setting by computing the mean curvature through our discretization. It is known that the mean curvature vector of the underlying surface can be directly computed using the operator:

HN = ∆ I
where I is the embedding of the structure (i.e. its real coordinates) and N the unit normal. We use this relation to compute the mean curvature on a Goursat shape (see Fig. 10) of equation 0.03(x 4 + y 4 + z 4 ) -2(x 2 + y 2 + z 2 ) -8 = 0 and compute the l ∞ error between the estimated curvature and the real one. Two Laplace-Beltrami operators are used: L h and L R-LOC lead to two mean curvature estimators called H L h and H L R-LOC respectively. We also compare our method with integral invariant [START_REF] Coeurjolly | Multigrid convergent principal curvature estimators in digital geometry[END_REF] Fig. 11 Graph of l ∞ error for the mean curvature estimation. Estimation using Jet Fitting is in green ( ), using integral invariant in orange (

), using L R-LOC in magenta ( ) and using L h in purple (

).

Conclusion and future works

We have adapted the discretization of the Laplace-Beltrami operator proposed by Belkin et al. [START_REF] Belkin | Discrete laplace operator on meshed surfaces[END_REF] to our digital surface.

We have proved strong consistency (i.e. pointwise convergence or convergence in the l ∞ norm) of our operator in Theorem 5 and gave convergence speed for the functional Laplacian approximation in Lemma 1 and Lemma 2. We have given associated empirical tests for various values of the exponent α of t h = h α . We also compared our approach with existing discretizations to confirm that none of them achieves pointwise convergence. A first immediate future work would be to compute the convergence speed of the functional Laplace operator approximation of Belkin et al. in Theorem 1. Next, we would like to reduce the complexity of the algorithm (which is O(n 2 ) if n is the number of surfels as we compute a convolution between a Gaussian function and differences of functions). A natural way is to look at cuts of the Gaussian function. In fact it is known that almost all the mass under a Gaussian is contained within a few multiple of σ (the standard deviation) typically two or three times σ . We have empirical results showing consistency of the Laplace operator when we cut the Gaussian but the theoretical proof is still an open problem. In addition to reducing the computational cost of the convolution, cutting the Gaussian implies a sparse matrix representation of the Laplace operator when considered as a linear operator. Sparsity is an interesting property in many geometry processing applications using linear system solvers or eigen decompo-sition of operators containing the Laplacian. We would also like to investigate an adaptation of [START_REF] Alexa | Discrete laplacians on general polygonal meshes[END_REF] to our digital surfaces, where we could use the normal estimation to compute discrete inner products on k-forms and obtain better results than the straightforward implementation. where W -1 is the lower branch of Lambert W -function (also called omega function or product logarithms).

Proof We use the Lambert W -function to prove this lemma.

In-depth study of this object can be found in the book of Corless, Gonnet, Hare and Knuth [START_REF] Corless | On the lambertw function[END_REF]. The function is defined as the multivalued function W that satisfies z = W (z)e W (z)

for z ∈ C. It is equivalently the inverse function of f (w) = we w . The graph of the Lambert W -function in the real numbers is drawn in Fig. 12. The function has two real branches W 0 and W -1 in the interval -1/e < x < 0 which join at x = -1/e. This means that the equation x = we w has two solutions in this interval (one per branch). We will use both branches: W -1 which is decreasing in its interval, and W 0 which is increasing in this interval. We also use the identity W (xe x ) = x. We do a proof by equivalence of inequalities:

Fig. 1 A

 1 Fig. 1 A digital surface of dimension two embedded in R 3 .

Fig. 2

 2 Fig.2Illustration of L DEC (top), and L COT (down) on triangular meshes. For L COT the area of integration A w is one third the area of all triangles incident on vertex w in green. For L DEC , we represent the Voronoi dual structure in orange. The dual of the edge e wp is in blue.

  of A, denoted by reach(A), is inf{δ A (y) : y ∈ Med(A)}. Definition 5 (The projection map) The projection map onto a compact set A is the map ξ A : R d+1 \ Med(A) → A that maps any points x of R d+1 \Med(A) to its unique closest point on A. We denote by ξ := ξ ∂ M the projection onto ∂ M. First, Theorem 2 states the Hausdorff stability between ∂ h M and ∂ M: the distance between those two is bounded by the grid step h. In other words, given a point y ∈ ∂ h M, there always exists a point x ∈ ∂ M within a ball of radius √ d+1 2 h.

Fig. 3

 3 Fig.3Illustration of the notations used in this paper. A smooth shape M, its Gauss digitizationD h (M), its h-cube embedding Q h [D h (M)] in yellow and its h-boundary ∂ h M. The cubical grid F h2 is displayed in green on the right, and also in the related "cubical border" F h 1 ∩ ∂ h M in red. The topological border of the h-cube embedding is used to push integral from ∂ M in the continuous setting to the discrete setting. Then, the cubical grid is used to split the integral on elements of various dimensions d, thus approximating the continuous sum by a discrete one.

Fig. 5

 5 Fig. 5 Illustration of the tube definition of ũ on ∂ M. Function ũ is equal to u when x ∈ ∂ M (the black curve on the figure), and is equal to u • ξ when x lies in the

  (ṡ)|| 2 4t h and bound |b|. As b can be either negative or positive, we first derive a negative minor bound and a positive upper bound on it and conclude using Lemma 3 in the appendix.

Fig. 6

 6 Fig.6Illustration of the projection function ξ . Surfel centroids ṡ (in yellow) are mapped using ξ to the blue dots. The proof of Lemma 1 lies in comparing ||y -ṡ|| and ||yξ (ṡ)||. We use triangle inequalities in (ṡ, ξ (ṡ), y) for the proof.

Fig. 7

 7 Fig. 7 Discrete settings used for the comparison between various Laplace-Beltrami operators. The unit ball S 3 is discretized using the Gauss digitization process. The left image is the digital surface ∂ h S 3 where we compute L h and L COMBI . The image in the center is the marchingcubes associated with the digital surface where we compute L MESH , L R-LOC and L COT . Finally, the right image represents the projection of the marching-cubes on ∂ S 3 where L P MESH , L P R-LOC and L P COT are computed.

1 3

 1 as advised by the authors.

1 6 in purple and t h = h 1 12Fig. 9 1 3

 1191 Fig. 9 Multigrid convergence graphs for various functions on ∂ S 3 , the unit sphere. l ∞ is displayed for L COT , L P COT , L COMBI , L MESH , L P MESH , L R-LOC , L P R-LOC , L QUAD , L P QUAD and L h . Parameters t h and r are equal to h 1 3 for all five convolution operators.

Fig. 10

 10 Fig.10An illustration of the mean curvature values (from -0.11 in violet to 0.29 in yellow) on a Goursat shape using our discretization of the Laplace-Beltrami operator (H L h ). We used h = 0.17, t h = 1.02 and ball of radius 5 for the Integral Invariant normal estimator.

Appendices Lemma 4 Lemma 5 Fig. 12

 4512 Fig.12Plot of the two principal branches of the Lambert W -function. W 0 is in orange, and W -1 in green. The two branches join at the point (-1/e, -1).

  Next we divide by 4t h , apply the function e -x to the inequality

	e -	||y-ṡ|| 2 4t h	≤ e -	||y-ξ (ṡ)|| 2 4t h	e	3(d+1) 16t h	h 2 +	||y-ξ (ṡ)|| √ d+1 4t h	h
	and subtract e -	||y-ξ (ṡ)|| 2 4t h	on each side:
	b ≤ e -	||y-ξ (ṡ)|| 2 4t h		e	3(d+1) 16t h

  (H II ) and Monge Form via Jet Fitting [9] (H JET ). For all three H JET , H II and H L h , we set the convolution radius to h 1 3 . Results can be found in Fig. 11. H L R-LOC does not converge because L R-LOC does not either. The other three H JET , H II and H L h converge. H L h is slightly better than H II whereas both are better than H JET .

	H JET	H II
	H L R-LOC	H L h
	1
	error (px)
	0.1
	0.1	1
	h
	u(x, y, z) = z

This work has been partly funded by C O M E D I C ANR-15-CE40-0006 research grant. We would like to thank the anonymous reviewers for their detailed comments and suggestions for the manuscript.


Milena Tsvetkova
email: milenaic@uni-sofia.bg

FEAR FROM THE READERS

Keywords: reading studies, altered reader, right to read, copyright, free culture, intellectual freedom, libraries, UNESCO, Bulgaria

The reading bear the great cultural responsibility not to allow the mass communication to serve the individual sovereignty contra-adaptive and anti-manipulative under the pressure of the mass and the unifications. The reader cultivates in himself a tendency to counteracting, anarchism, self-confidence, optimism and information superiority. These are vested functions of reading as "dangerous" mental process. Referring Police Unit for Combating Organized Crime's operation "blocking the pirate book distribution in Internet through www.chitanka.info".

INTRODUCTION

Referring Police Unit for Combating Organized Crime's operation "blocking the pirate book distribution in Internet through www.chitanka.info" on June 22, 2010. On the 22 th of June 2010 Bulgarian Police Unit for Combating Organized Crime (GDBOP) hit the Internet-based network for distributing contents chitanka.info (in Bulgarian "chitanka" means an alphabet book). The police was alarmed by multiple sources including The Bulgarian Book Association, The Union of Translators in Bulgaria, and different publishing houses. The action was claimed a success in the fight against Internet piracy.

RIGHT AND DUTY TO READ

When we discuss the right of the author of reading's objects -my position is peremptory. I vote for free books and free reading. I vote for total elimination of readers' barriers -as good as of the market, the language, or the psychology. I vote against the juridical monopoly of "rights" holders and mercantile craftsmen. I vote against the market censorship over reading.

It seems too schizophrenic to block the access to texts to be read in time when all burst in tears for the "fading" reading. The most important thing is to read. Whether playing games with the guards of the texts, or passing over the system. The author's right is an ideal right and therefore impossible to be annulled.

Nevertheless, the author is not alone "on the saddle" -he also have been reading, listening, using ideas, copying, imitating, or stealing from others. "Author" is only The First One but who is he, is another topic. As Mark Twain reminds us, only Adam was happy because he knew everything said is said by him first. And a very important thing -the knowledge and the information are a priori free! We pay for them when the planet is invaded by economy, market, and money. This is the truth about our whole written civilization.

I am the author of several books, over 200 articles, and an immeasurable number of unpublished actual lectures I hand out to all my students. However I am a lobbyist for the Reader! He is the master in the circle of written knowledge and books. He is the client, he is the payer. He is the evaluator and the judge. He is able to put a publisher or a translator, an author or a book on pedestal. He also has the power the destroy them or to throw them in oblivion.

I want to say a thing to those who speculate with the mantra "Copyright Law" and to all the clappers to the police blockade of the site chitanka.info who jump with cow delight in the eyes and satisfied but thoughtless heads.

CONTENT VS BOOK

Law of Copyright and Connected Rights permits free usage, reproduction of already printed works and their direct or indirect multiplying in any mode and in any form. And this exactly is what the web-based virtual database of written works does. And no more.

Therefore every Internet on-and offline site for reading or downloading texts is a library. Not a bargain bookstall.

Naturally, every virtual library distributes contents, not books. In its core the content is digital. According to all the major laws of physics, without the hard wrapper it flows, spills, runs into our PCs. This is its job.

The wolfed craftsmen in the book field should try to calculate separately Content and their paper product known as a book. There we'll see whom to pay while reading Files. By films and songs the "downloaded" or perceived with eyes and ears work is identical to both copy and original but the book is a different story. We perceive the content from a paper copy in one way but from the display the way is different. Everyone settles his/her own account which way to chose. Because the file with the content of a paper book is even not a copy of this book.

Of course, every library -physical or digital -functions not as a commercial object but as an organization with ideal goals and exactly to provide the mankind with reproducted and brought to everyone of us readings for free usage.

Of course, the libraries are the only warrantor of Intellectual Freedom since the writing and reading exist. This is to be read in all the international manifestos of librarian associations. Who could forget that intellectual freedom is a freedom not only of writing, translating, publishing but most of access to information and reading?

Of course, every library -physical or digital -is an obstacle for publishers' and distributors' pockets and authors' and their heirs' purses. Since the world turns, the libraries legally multiply reading objects or texts infusing them into readers' heads. Or allowing people to copy them through hand-writing. Or allowing photocopies. Or giving the readers the freedom to distribute the only copy borrowed for home usage among neighbors and friends so their eyes could "steal" its content too.

There is no a living author who could precisely count his/her readers.

The reader is countless. The reader is Legion. The reader is the one who induces life into the book. Without eyes reading it the book is only a brick.

All the segments in the process of producing books exit only because of the reader.

Without him all the authors, editors, translators, publishers, distributors, and librarians are only mummies.

ALTERED READER

The main motive for police to hit virtual reading-rooms and libraries after being pressed by mercantile sources is transparent and discouraging: an inadmissible for the 21 st century ignorance of the altered reader.

The altered reader wants free books. He/she has evolved as a product of the "free culture" and still hopes there'll be someone here to believe in the American experience with the incredible commercial potential of the free book first published in Internet and afterwards on paper and against payment. According to the young readers the consumption of the written words is considered as a free operation but only the depository-type of persons and institutions calculates the reading in their investments (though they don't mean that exactly). Curious are the thoughts of Jean-François Barbier-Bouvet, an expert in reading, who grasps the "free culture" as a threat for the written words. He admits the dropping away of the "psychological" book price (which we are ready to pay to have it in our possession) though with the virus "free of charge", Internet uncontrollably infects the whole culture of deposits based on the printed media. The book appears expensive, grasped not absolutely (the book price objectively evolves slower than life's rise in price) but compared to other information sources more and more depending on the "free culture". The French researcher's anxiety is consequent from the space "loneliness" of the book where we have to pay for equipment not for its usage, for bearers not for the message, for access not for the content. The economic model of the book and printing (although they'll pass over to co-virtuality soon and fully) would hardly be able to make an exception.

The adepts of publishing industry barely continue believing that the book will keep them fed as an "island of the real value" in an empire of entroping prices and free of charge issues. The media industry responded in a hurry to this reader's caprice too -with free newspapers and magazines, gift books put into press packages, and a palette of crossmedia magic. The governments of Canada and France wholehearted fund the free books distribution -as PR actions or as a tool for stimulating the reading in general. Even the altered reader gets used to help him-/herself with this caprice -he/she established and multiplies the movement for free reading Bookcrossing. As a product of this actual claim appear also the initiatives for free reading in informal networks Bookboxes, Bookrings, Bookrays. The alternative countercultural Bookcrossing movement calls its participants on "sharing" (an operation in the free culture) their disposable readings with other people.

The altered reader is confident that this exactly is the way for the book to fulfill itself as a media, concept, thinking -among more people, independent on their solvency. Not just as an interior element at home, a prestigious decoration in the office, as an archive unite in the library or an object of commodity-monetary relations.

READERS ARE THE MOST DANGEROUS CREATURES

However the deeper motive to block the access to texts is the fear of the reading and well-read human being. The executors of such actions don't realize this. But it's not their fault. They carry out their duty in the next cult of justice and blindly sacrifice for the idol of skinflints and mercantilists.

But why would anyone fear of the readers? First, the reading is an asocial activity -it's committed in loneliness, realized invisibly (in the cranium) and therefore uncontrollably. The reader is almost out of sight for any sociological approach in the circumstance that the profound genuine reading is an auto-communication process. And the comfort environment for auto-communication is the "social isolation". Furthermore, the reading represents a loss of social time and intellectual energy of the society thus it is socially ineffective. UNESCO's report on The International Book Year in 1972 states that in the Third world the reading is still considered as a shameful, criminal, and reproachful act. The Afghanistan community bans the act of individual reading; the separating of an individual from the contiguous surrounding group is considered as a suspicious act and perceived as a threat of integrity. The reading alone is almost a sin similar to onanism. Nowadays the most liberal societies experience similar fears, too. So the reading person is lost for the predatory manipulators from the medium outside.

Here follows the third motive for fear and concealed hatred for the readers: the reading is an impeccable filter against rhetorical and audio-visual manipulations. In the act of reading the individual finds time to cultivate his/her critical judgment, to operate information destructively, to dissect opinions, to compare, juxtapose, expose.

Forth, the reading is a sabotage activity. A brain practice sabotaging the stereotypes and the authorities, all the known and unknown totalitarian concepts. The reading individual cultivates his/her inclination for counteraction, anarchism, self-confidence, optimism, and informational superiority.

The man reads not only in order to differ from the others, to be original, to be always on the top. He/she reads in order to overcome the boundaries of the crowds formed in front of the TV screens, to blow up the oneness of the mass, to tear the simultaneousness and transiency of the TV show once and for all left in the past. The reading is also the best method to avoid standardization and vulgarization (both inclined to the crowdism) typical for our high technologized century. The book as well as the reading are highly responsible to the culture not to allow under pressure of massovization and unifications the mass communications to serve the individual sovereignty counter adaptive and anti manipulative. These again are immanent functions of the reading as a "dangerous" brain process.

And fifth, in this planet only the reading man possesses the powers to changehimself, his environment, state, the world. The connoisseurs of neurophysiology of the reading know that the reader develops not really his/her "gray matter" but their cortex, meaning to free him-/herself from animal dependencies and behavior models genetically determined and tamed in a community (herd). As the brain evolves from inside out, so the reader's personality (after just has finished the last page) develops and changes his/her environment from inside out -absolutely decentralized, unsuspected and invisibly for "the guards from the tower" as well surprisingly -somewhere form the bosom, aside and up, similar to an earthquake's epicenter. The readers are disliked and dangerous because they change their opinion too (only the cow doesn't change it for the lack of information). Thus they make a dumping of the official theses and views.

CONCLUSION

If the world should remain like it is now, the readers are the most dangerous creatures. Arrest us all. But again you'll forget the most important player -the book. The book is a very sly media -it doesn't stand hindrances, chains and limitations. It is loaded with the breaking-through power of the light. Especially when it's wanted.