1 Introduction

Objectives. In geometry processing, Partial Differential Equations (PDEs) containing Laplace-Beltrami operator arise in many applications such as surface fairing, mesh smoothing, mesh parametrization, remeshing, mesh compression, feature extraction or shape matching (see [16] for an extensive survey). On digital surfaces, few digital Laplace-Beltrami operators have been proposed and none has been evaluated in terms of multigrid convergence (convergence of the operator toward the continuous one in digitization of smooth manifolds on grid with decreasing gridstep).

Contributions. In this article, we propose a discrete Laplace-Beltrami operator on digital surfaces (boundaries of subsets of \(\mathbb {Z}^2\) embedded in 3D). This new operator adapts Belkin et al. [2] on our specific data. The method uses an accurate estimation of areas associated with digital surface elements. This estimation is achieved through a convergent digital normal estimator described in [5]. We show experimental convergence of our operator but also that none of the existing approaches adapted to digital surfaces achieves such convergence. Finally, we illustrate the interest of the discretized Laplacian on digital surface geometry processing.

Related works. The Laplace-Beltrami operator being a second order differential operator (divergence of the function gradient, see Sect. 2) a discrete calculus framework is required to define such operator on embedded combinatorial structures such as meshes or digital surfaces. First works on discrete calculus may be found in the Regge Calculus [21] for quantum physics, where tetrahedra in combination with edge lengths are used. Works on geometric acquisition devices and models drove studies toward calculus working on meshes and more generally on simplicial complexes. Early works include a definition of the Laplace-Beltrami operator using the classical cotangent formula [19] for solving the problem of minimal surfaces, which is an analog of the standard finite element method [16]. Exact calculus generalizing the cotangent discretization in 2D based on finite elements [20] emerged from the German school but with a restriction to triangular complexes.

In a more generic discrete calculus perspective, the Discrete Exterior Calculus (DEC) framework was then developed in the computational mathematics and geometry processing community. Another more recent formulation of the DEC comes from Hirani’s thesis [11] and later by the monograph [8]. On triangular meshes, DEC based Laplace-Beltrami operator and the cotangent based one coincide.

In [10], authors show that under some strong assumptions, the cotangent Laplacian on a triangular mesh converges to the continuous one when the mesh interpolates a smooth manifold with increasing precision (with a continuous one-to-one map between the mesh and the manifold, which would not be the case on digital surfaces). The operator converges in the sense of distributions and the authors show that pointwise convergence in the \(l_2\) sense does not generally hold. On triangular meshes, Belkin et al. [2] have proposed a first Laplace-Beltrami operator that converges in the \(l_\infty \) uniform case. The digital Laplacian operator we propose is an extension to digital surfaces of such operator.

In digital geometry, many estimators of differential quantities have been proposed and there exists multigrid convergent estimators for many quantities such as length, tangent and curvature in 2D (see [6] for a complete survey), surface area [14], normal vectors and curvature tensor in dimension 3 [5]. A preliminary approach can be found in [4, 17]. It focuses on the conformal map computation of a digital surface, which is a related problem involving the definition of a Laplace-Beltrami operator. However, their definition is based on the cotangent formula, which lacks pointwise convergence. When designing a discrete version of the Laplacian, not all properties of the continuous one can be expected at the same time. In [24], entitled “Discrete Laplace operators: No free lunch.”, the authors have proposed a formal evaluation of such properties. We position our new digital Laplace-Beltrami operator with respect to this analysis.

Outline. After introducing mathematical definitions, we review the classical approaches to define a discrete Laplacian and compare their properties in Sect. 2. We then formalize our operator in Sect. 3. In Sect. 4, we experimentally evaluate our proposal in terms of multigrid convergence and geometry processing applications.

2 Discretizations of the Laplace-Beltrami Operator and Their Properties

We first describe various discretizations of the Laplace-Beltrami operator on triangular meshes. We then check several desired properties of Laplacian using [24] as a baseline for comparisons.

2.1 Preliminaries and Classical Discretizations on Triangular Meshes

Let M be a 2-smooth manifold, with or without boundary, embedded in \(\mathbb {R}^3\). The intrinsic smooth Laplace-Beltrami operator [22] is defined as:

$$\begin{aligned} \varDelta : C^2 (M)&\rightarrow C^2 (M) \nonumber \\ u&\mapsto \text {div} (\nabla \, u), \end{aligned}$$
(1)

where \(C^2\) is the set of functions which are twice differentiable and with the second derivative continuous (in the literature, alternative definitions may consider “\(- \text {div} (\nabla \, u)\)” for \(\varDelta \, u\)).

Let \(\varGamma \) be a combinatorial structure (a triangular mesh for instance), \(V(\varGamma )\) its set of vertices and \(F(\varGamma )\) its faces. Let \({u} : M \rightarrow \mathbb {R}\) be a twice differentiable function. We suppose that \(V(\varGamma )\) is a sampling of M (i.e. \(V(\varGamma ) \subset M\)). In other words, u(w) is perfectly defined for \(w\in V(\varGamma )\).

A first simple discretization only considers the combinatorial structure of \(\varGamma \). Such Laplacian is either called graph Laplacian or combinatorial Laplacian of \(\varGamma \) [25]:

$$\begin{aligned} (\mathcal {L}_{COMBI} \, u)(w):=- deg(w) u(w) + \sum _{p \in {\text {link}}_0(w)} u(p) \,, \end{aligned}$$
(2)

for all \(w\in V(\varGamma )\) where \({\text {link}}_0(w)\) is the set of points \(V(\varGamma )\) adjacent to w and deg(w) is the degree of w in \(\varGamma \).

A more complex approach can be defined using DEC operators [8, 11]. Using an arbitrarily embedded dual structure of \(\varGamma \), the Laplace operator can be shown to be a classical weighted double finite difference:

$$\begin{aligned} (\mathcal {L}_{DEC} \, u)(w) :=\frac{1}{ |\star w| }\sum _{p \in {\text {link}}_0(w)} \frac{|\star e_{wp}|}{|e_{wp}|}(u(p) - u(w)), \end{aligned}$$
(3)

where \(\star \) is the Hodge-duality star operator acting on discrete forms (see [11]), and \(|\cdot |\) the measure of a k-cell. As illustrated in Fig. 1, \(|\star e_{wp}|\) would be the length of the segment orthogonal to \(e_{wp}\). If we set all measures to one, \(\mathcal {L}_{DEC}\) coincides with \(\mathcal {L}_{COMBI}\).

By fixing the dual of \(\varGamma \) to be the Voronoi diagram of its vertices and by computing the measures as Euclidean lengths and areas of such dual complex, the DEC operator coincides exactly with the famous cotan Laplacian [19]:

$$\begin{aligned} (\mathcal {L}_{COT} \, u)(w) :=\frac{1}{2 A_w} \sum _{p {\text {link}}_0(w)} \left( \cot (\alpha _{wp}) + \cot (\beta _{wp})\right) \left( u(p) - u(w)\right) , \end{aligned}$$
(4)

where \(A_w\) is one third of the area of all incident triangles to vertex w, \(\alpha _{wp}\) and \(\beta _{wp}\) are the angles opposing the corresponding edge \(e_{wp}\) (see Fig. 1).

Fig. 1.
figure 1

Illustration of \(\mathcal {L}_{DEC}\) (left), and \(\mathcal {L}_{COT}\) (right) on triangular meshes. For \(\mathcal {L}_{COT}\) the area of integration \(A_w\) is one third the area of all triangles incident on vertex w in green. For \(\mathcal {L}_{DEC}\), the dual structure is in orange and the dual of the edge \(e_{wp}\) is in blue. (Color figure on line)

Full size image

Finally, we detail the definition of the mesh Laplacian from [2]. Let \(g : M \times (0, T) \rightarrow \mathbb {R}\) be a time-dependent function which solves the partial differential equation called the heat equation:

$$\begin{aligned} \varDelta g(x, t) = \frac{\partial }{\partial t}g(x, t), \end{aligned}$$
(5)

with initial condition \(g_0 = g(\cdot , 0) : M \rightarrow \mathbb {R}\) which is the initial temperature distribution. An exact solution [22] is:

$$\begin{aligned} g(x, t) = \int _{y\in M} p(t, x, y) g_0(y) dy, \end{aligned}$$
(6)

where \(p \in C^\infty (R^+ \times M \times M)\) is called the heat kernel. The construction of the heat kernel involves complex technics [22]. Fortunately, there exists many approximations of p as t tends toward 0 (called small time asymptotics). Early work includes the famous Varadhan formula [23] and later extensions on a wider class of shapes [18]. Recently, an approximation using the ambient metric (i.e. the \(L_2\) norm) have been proposed by Belkin et al. It is known to converge toward the real heat kernel for small t (see Lemma 5 of [1]):

$$\begin{aligned} g(x, t) \underset{t\rightarrow 0}{\sim } \frac{1}{4 \pi t} \int _{y\in M} e^{ - \frac{||x - y||^2}{4t} }g_0(y)dy. \end{aligned}$$
(7)

Injecting Eq. (7) into Eq. (5), applying a finite time difference and knowing that the integral of e over M is one:

$$\begin{aligned} \varDelta g(x, t) = \lim \limits _{t \rightarrow 0} \frac{1}{4 \pi t^2} \int _{y\in M} e^{-\frac{||x-y||^2}{4t}}(g_0(y) - g_0(x))dy. \end{aligned}$$
(8)

The previous equation can be seen as a convolution between differences of g and a time dependent Gaussian. Then, the mesh Laplace operator [2] on \(\varGamma \) is:

$$\begin{aligned} (\mathcal {L}_{MESH} \, u)(w) :=\frac{1}{4 \pi t^2} \sum _{f\in F(\varGamma )} \frac{A_f}{3} \sum _{p \in V(f)} e^{- \frac{||p-w||^2}{4t}}(u(p) - u(w)), \end{aligned}$$
(9)

where \(A_f\) is the area associated to the face f.

2.2 Desired Properties of a Discrete Laplacian

As discussed in [24], all properties of the continuous Laplace-Beltrami operator may not be preserved when discretizing it. We consider the discrete Laplacian as a linear operator acting on values \(\mathbf {u}:=\{u_p\}\) on \(V(\varGamma )\) (represented as a vector in \(\mathbb {R}^{|V(\varGamma )|}\)). Such operator can thus be denoted as a matrix \(\mathbf {L}\) with components \(l_{ij}\). Hence, \(\mathbf {v}:=\mathbf {L}\mathbf {u}\) would be the resulting Laplacian of \(\mathbf {u}\). Expected properties of the discrete Laplace-Beltrami operators are:  

Symmetry (SYM).:

\(\ell _{ij} = \ell _{ji}\) for \(0\le i <|V(\varGamma )|\). This is very useful when it comes to solve linear systems as solvers are usually more performant when the matrix is symmetric.

Locality (LOC).:

\(\ell _{ij} \ne 0\) if and only if vertices i and j share a common edge. The locality property gives very sparse matrices decreasing drastically memory consumption. It also opens a panel of very fast linear system solvers.

Linear Precision (LIN).:

\(\mathbf {L}\mathbf {u}= 0\) whenever \(\mathbf {u}\) is a linear function restricted to a plane.

Positive Weights (POS).:

\(\ell _{ij} \ge 0\) for \(i\ne j\). Furthermore, for each vertex i, there exists a vertex j such that \(\ell _{ij}>0\).

Positive Semi-Definiteness (PSD).:

The matrix is symmetric positive semi-definite regarding the standard inner product, and has a one-dimensional kernel. (SYM) and (POS) imply (PSD), but (PSD) does not imply (POS). This property ensures that the basis generated by the eigenvectors of L is orthogonal and that the eigenvalues are real.

Dirichlet Convergence (CON).:

When considering a sequence of meshes \(\{\varGamma _i\}\) converging to M in a given sense as \(i\rightarrow \infty \), we want the Laplacian sequence \(\mathbf {L}_i\) to converge to \(\varDelta \) with respect to the discrete Dirichlet problem. The convergence is mandatory when we seek approximate solutions of partial differential equations.

  The (CON) property requires a formal definition of the sequence \(\{M_i\}\). In addition to [24], we add the following property  

Pointwise Convergence (PCON).:

For a given sequence meshes \(\{\varGamma _i\}\), we want the associated Laplacians \(\{\mathbf {L}_i\}\) to converge to \(\varDelta \) in a pointwise \(l_2\) or \(l_\infty \) sense. This notion of convergence is stronger than (CON) which is implied by (PCON).

  For the Laplacian on digital surfaces (e.g. \(L_h^\star \) we define in Sect. 3), we consider a sequence of combinatorial structures defined as the boundaries of the Gauss digitizations of M (see next section). We thus consider the (DPCON) property as the pointwise convergence of the operator for multigrid digital surfaces.

Table 1 compares classical Laplacian discretizations with respect to these properties. The (DPCON) property has been evaluated experimentally in Sect. 4. \(\mathcal {L}_{COMBI} \) being purely combinatorial, the same operator can be used for both meshes and digital surfaces. For \(\mathcal {L}_{COT}\) and \(\mathcal {L}_{MESH}\), we have considered the marching-cubes representation of the digital surfaces (more precisely, a continuous analog of the dual surface [13]). Properties of \(L_h^\star \) are discussed in the next section.

Table 1. Properties of various Laplacians. See text for the description of each property. means that the property is valid, whereas means not. N.A. stands for Not Applicable and ? means that we do not know whether it is true or not. Finally exp. holds for experimental convergence with no known theoretical proof.
Full size table

3 New Laplace-Beltrami Operator on Digital Surfaces

In the discretization schemes discussed above, the points of the combinatorial structure interpolate the underlying manifold. In order to define a discrete Laplace-Beltrami operator on digital surfaces, a challenge is to work with a combinatorial structure that only approximates the underlying manifold in a Hausdorff sense. First, we extend Eq. (9) to digital surfaces. Proper definitions of a digital surface can be found in [12, 14]. Let us recall the Gauss digitization process:

Definition 1

(Gauss digitization). Let \(h > 0\) be the sampling grid step. The Gauss Digitization of an Euclidean shape \(X \subset \mathbb {R}^d\) is defined as \(\mathtt {D}_h(X) := X \cap (h\mathbb {Z})^d\) where d is the dimension.

For smooth object X in dimension 3 with boundary \(M:=\partial X\), the digital surface is defined as the topological boundary of \(\mathtt {D}_h(X)\), denoted \(\partial _h X\) (see [14] for details). More precisely, the digital surface has a cellular representation in a Cartesian cubical grid and is composed of points of dimension 0 (pointels, \(\mathbb {E}^0\)), straight segments of dimension 1 (linels, \(\mathbb {E}^1\)) and squares of dimension 2 (surfels, \(\mathbb {E}^2\)).

The topological boundary \(\partial _h X\) is an O(h)-Hausdorff approximation of M [14]. As a consequence, we need to map the smooth function u defined on M to \(\partial _h X\):

Definition 2

(Extension of u to \(\partial _h X\) ). Given a smooth function u on M, we define the extension \(\tilde{u}\) of u to \(\partial _h X\) as

$$\begin{aligned} \tilde{u}(\mathbf {s}) :=u(\xi ( \dot{\mathbf {s}} ))\,, \end{aligned}$$

where \(\dot{\mathbf {s}}\) is the centroid of the surfel \(\mathbf {s}\in \mathbb {E}^2\), and \(\xi \) is the map that projects a point of \(\partial _h X\) onto the closest point of M.

We show below how to adapt the definition of [2], recalled in Eq. (9), to digital surfaces. We chose this approach because \(\mathcal {L}_{MESH}\) has an interesting pointwise convergence for triangular meshes. Our Laplace-Beltrami operator is thus defined as follows

Definition 3

(Digital Laplace-Beltrami operator). The digital Laplace-Beltrami operator is defined on \(\partial _h{X}\) as:

$$\begin{aligned} (L_h^\star \, \tilde{u})({\mathbf {s}}) :=\frac{1}{4 \pi t_h^2} \sum _{\mathbf {r}\in \mathbb {E}^2} e^{- \frac{||\dot{\mathbf {r}} - \dot{\mathbf {s}}||^2}{4t_h}} \mu (\mathbf {r}) (\tilde{u}({{\mathbf {r}}}) - \tilde{u}({{\mathbf {s}}})), \end{aligned}$$
(10)

where the sum is taken over all surfels of \(\partial _h{X}\), \(\dot{\mathbf {r}}\) is the centroid of the surfel \(\mathbf {r}\), \(\mu (\mathbf {s})\) is equal to the dot product between an estimated normal and the trivial normal orthogonal to the surfel \(\mathbf {s}\) and \(t_h\) is a function of h tending to zero as h tends to zero.

The quantity \(\mu (\mathbf {s})\) is called the measure of the surfel \(\mathbf {s}\): it is the area of the projected surfel \(\mathbf {s}\) onto the tangent plane induced by the estimated normal. Normal vectors are estimated using the estimator presented in [5, 15] which has the multigrid convergence property. Note that summing \(\mu \) for each surfel of the surface leads to an estimation of the global area of the shape boundary, which itself has a multigrid convergence property [14]. This surfel measure is a key ingredient of the digital formalization of the operator leading to an experimental multigrid convergence and isotropic properties when used for diffusion (see Sect. 4).

If we index surfels in \(\mathbb {E}^2\), \(L_h^\star \) has a \(|\mathbb {E}^2|\times |\mathbb {E}^2|\) matrix representation \(\varvec{L}_h^\star \) defined as follows:

$$\begin{aligned} (\varvec{L}_h^\star )_{ij} = \left\{ \begin{array}{lcc} \displaystyle { \frac{1}{4 \pi t_h^2} e^{- \frac{||\dot{\varvec{s}}_j - \dot{\mathbf {s}}_i ||^2}{4t_h}} \mu (\varvec{s_j})} &{} \text { if } &{} i \ne j \\ \displaystyle { - \sum _{k \ne i} (\varvec{L}_h^\star )_{ik}} &{} \text { if } &{} i = j \end{array} \right. \end{aligned}$$
(11)

In other words, if \(\tilde{u}\) is constant and equal to 1, for a surfel \(\mathbf {s} \in \mathbb {E}^2\) of index i, we have \((\varvec{L}_h^\star \varvec{\tilde{u}})_i = (L_h^\star \tilde{u})(\dot{\varvec{\mathbf {s}}}_i)\).

For the properties listed in Table 1, we can first observe that Eq. (11) implies that we have property (POS) but not (SYM). As \(L_h^\star \)   performs a convolution on the complete surface with a Gaussian kernel, we do not have the (LOC) (similarly to \(\mathcal {L}_{MESH}\)). Although we do not provide a theoretical proof of (PSD) for \(L_h^\star \), eigenvalues of \(\varvec{L}_h^\star \), have always been positive through all experiments. The (PCON) property is not applicable to our framework. The pointwise convergence (DPCON) is observed experimentally and discussed in Sect. 4.1.

4 Experiments

4.1 Experimental Convergence

We first evaluate the multigrid convergence of our Laplace-Beltrami operator. We consider a unit sphere \(\mathbb {S}^2\) and three different smooth functions \(u : \mathbb {S}^2 \rightarrow \mathbb {R}\), namely z, \(x^2\) and \(e^x\) (see Fig. 2). Let \(\theta \) be the azimutal angle, and \(\phi \) the polar angle. The spherical Laplacian is then:

$$\begin{aligned} \varDelta _{\mathbb {S}^2} u(\theta , \phi ) = \frac{1}{\sin ^2 \phi } \frac{\partial ^2 u}{\partial \theta ^2} + \frac{1}{\sin \phi } \frac{\partial }{\partial \phi } \left( \sin \phi \frac{\partial u}{\partial \phi } \right) . \end{aligned}$$
(12)

We compute the Gauss digitization \(\mathtt {D}_{h}(\mathbb {S}^2)\) of the sphere for decreasing grid steps h. Since the elements of \(\partial _h \mathbb {S}^2\) do not interpolate the sphere, u is extended to \(\tilde{u}\) as defined in Sect. 3. We compute \(\mathcal {L}_{COMBI}\) and \(L_h^\star \) directly on \(\partial _h \mathbb {S}^2\), but \(\mathcal {L}_{COT}\) and \(\mathcal {L}_{MESH}\) on the associated marching-cubes triangulation. Since the vertices of this mesh coincide with the centroids of the surfels of \(\partial _h \mathbb {S}^2\), all these operators \(\mathcal {L}_{COMBI}\), \(\mathcal {L}_{MESH}\) and \(L_h^\star \) are evaluated at the same points.

For \(L_h^\star \), we use the normal vector estimator described in [5] to estimate the measure of the surfels. In addition, for both \(L_h^\star \) and \(\mathcal {L}_{MESH}\), the parameter \(t_h\) is set to \(0.1 \times h^{\frac{1}{3}}\). As the discretization becomes finer, the standard deviation \(\sigma := \sqrt{2t_h}\) of the Gaussian function decreases and the number of points within the standard deviation \(\sigma \) increases. The constant factor 0.1 is a scale term derived from the unit sphere that sets the kernel to 1 / 10 of the sphere.

For comparison, in order to mimic the setting of [2], we have also considered the Laplacian \(\mathcal {L}_{MESH} ^P\), which corresponds to \(\mathcal {L}_{MESH} \) when the vertices of the marching-cubes are projected onto the sphere. In our framework, this operator is the gold standard, since the position errors are corrected by the projection, which is usually unknown.

For all the above operators, we plot in Fig. 2 the \(l_2\) and \(l_\infty \) error between the computed Laplacians and the true spherical one against the grid step h. First, we observe that errors for \(\mathcal {L}_{COT}\) (in blue ) and \(\mathcal {L}_{COMBI}\) (in green ) are constant: clearly both operators are non-convergent. Non-convergence is also observed for \(\mathcal {L}_{MESH}\) (in red ) but with lower errors. On the opposite, \(\mathcal {L}_{MESH} ^P\) (in orange ) shows convergence behavior for both \(l_2\) and \(l_\infty \) error, as expected in [2]. Concerning \(L_h^\star \) (in purple ), experimental convergence holds for the three functions. For the periodic function z, the convergence speed is slower than \(\mathcal {L}_{MESH} ^P\) whereas for the non-linear functions, we can see that convergence speed is the same. Moreover, the \(l_2\) error for \(L_h^\star \) tends toward \(\mathcal {L}_{MESH} ^P\), and its \(l_\infty \) error is close to that of \(\mathcal {L}_{MESH} ^P\).

Fig. 2.
figure 2

Multigrid convergence graphs for various functions on \(\mathbb {S}^2\), the unit sphere. Both \(l_2\) error in plain line and \(l_\infty \) in dashed line are displayed for \(\mathcal {L}_{COMBI}\), \(\mathcal {L}_{MESH}\), \(\mathcal {L}_{MESH} ^P\), and \(L_h^\star \).

Full size image

4.2 Shape Approximation Using Eigenvectors Decomposition

In this section, we consider the spectral analysis framework to process shapes geometry [16]. Given a shape and its Laplace-Beltrami operator, we compute the eigenvalues and eigenvectors of the operator and project the geometry onto the eigenvector basis of the first k eigenvalues. More formally, given an operator \(\varvec{L}\), we denote by \(\varvec{e}_1, \varvec{e}_2, \dots , \varvec{e}_n\) its normalized eigenvectors and the matrix \(\varvec{E}\) whose columns are those eigenvectors. By \(\lambda _1, \lambda _2, \dots , \lambda _n\) we denote the associated increasing eigenvalues where n is the number of rows of \(\varvec{L}\). Given three input vectors \((\varvec{X},\varvec{Y},\varvec{Z})\) encoding the vertex positions in \(\mathbb {R}^3\), we can approximate the input shape using a fixed number k of eigenvectors:

$$\begin{aligned} \varvec{X}^{(k)}&= \varvec{E}^{(k)} (\varvec{E}^{(k)})^T \varvec{X},\\ \varvec{Y}^{(k)}&= \varvec{E}^{(k)} (\varvec{E}^{(k)})^T \varvec{Y}, \\ \varvec{Z}^{(k)}&= \varvec{E}^{(k)} (\varvec{E}^{(k)})^T \varvec{Z}, \end{aligned}$$

where \(\varvec{E}^{(k)}\) is a matrix of size \(n\times k\) containing the first k eigenvectors columnwise. We compute the eigen decomposition on \(\mathtt {D}_h(M)\) for \(\mathcal {L}_{COMBI}\), and \(L_h^\star \) in Fig. 3 on a bunny object (\(64^3\), 13236 eigenvectors). We illustrate the reconstruction for increasing number of eigenvectors k. For low frequencies (\(k \le 100\)), we observe that \(L_h^\star \) captures more geometrical details than \(\mathcal {L}_{COMBI}\). For \(k=100\), we clearly have a better approximation of both ears of the bunny shape. As k increases, both reconstructions converge to the original bunny shape (\(\varvec{E}^{(n)} (\varvec{E}^{(n)})^T=\varvec{I}\)). In Fig. 4, we show the first 20 eigenvectors on a axis aligned cube.

Fig. 3.
figure 3

Images of the reconstruction using an increasing number k of eigenvectors. (First row) using \(\mathcal {L}_{COMBI}\), (second row) with \(L_h^\star \) (\(r=6\) for [5] and \(t_h=3\)).

Full size image
Fig. 4.
figure 4

Eigenfunctions are displayed on a simple cube with faces aligned with the grid axes (with a red to blue colormap and zero-crossing in white). (Color figure on line)

Full size image

4.3 Heat Diffusion

In this section, we highlight an interesting isotropic property of \(L_h^\star \) compared to the combinatorial one. We compute a heat diffusion when the source is a Dirac in the center of a rotated cube face (heat diffusion is a preliminary step of [7] to estimate geodesics on a manifold). We can derive an expression of g(x, t) from Eq. (5):

$$\begin{aligned} (\varvec{I} - t\varvec{L})g(x, t) = g_0(x), \end{aligned}$$
(13)

where \(\varvec{I}\) is the identity matrix, and \(\varvec{L}\) is the Laplacian matrix associated to any discrete laplace operator. This diffusion only makes sense for small t and around the Dirac \(g_0(x)\). As the computed heat diffusion decreases exponentially, we display the absolute value of its \(\log \). For \(\mathcal {L}_{COMBI}\) (first column in Fig. 5), both small and large values of t lead anisotropic estimations of the intrinsic metric due to the staircase effect of the face (concentric rhombi or concentric ellipses depending on t, see similar discussion in [7]). When using Eq. (13) with our matrix based operator \(\varvec{L}_h^\star \) (second column), even if numerical instabilities occur far from the Dirac depending on t, the intrinsic metric is perfectly estimated (concentric circles). Note that both \(\mathcal {L}_{MESH}\) and \(L_h^\star \) already compute a diffusion without the approximation of Eq. (13). The third column shows the associated diffusion. Even if \(\mathcal {L}_{MESH}\) does not have multigrid convergence properties, on this specific object, both methods provide good isotropic behaviors.

Fig. 5.
figure 5

Heat diffusion on a cube aligned with \(\mathbb {R}^3\) axis. (First column) using \(\mathcal {L}_{COMBI}\), (second column) using \(L_h^\star \) and (third column, with \(t_h=4\)) using the diffusion computed through the ambient heat kernel. The rightmost picture shows staircases on the rotated cube.

Full size image

5 Conclusion and Future-Works

In this paper, we have investigated different discretization schemes of the Laplace-Beltrami operator on digital surfaces. The contribution is twofold: first, we have shown that classical schemes either do not asymptotically converge to the expected operator, or contain anisotropic artifacts when used for geometry processing tasks. Second, we have proposed a new Laplace-Beltrami operator that incorporates multigrid convergent surfel measures allowing us to have both an experimental multigrid convergence and isotropic properties on digital surfaces.

A natural future-work consists of focusing on the multigrid convergence proof of the operator. In dimension 2, a preliminary proof has been derived using digital integration results from [14] but it is still an open problem in dimension 3.