Laplacian Spectral Basis Functions

TL;DR

This paper explores various Laplacian-based spectral basis functions for 3D shape analysis, comparing their properties and computational aspects to improve signal approximation and shape processing techniques.

Contribution

It introduces diffusion and Laplacian spectral basis functions, analyzing their properties and providing a comprehensive comparison with existing basis functions for shape analysis.

Findings

Laplacian spectral basis functions exhibit desirable invariance and locality properties.

Diffusion basis functions offer improved stability over traditional eigenfunctions.

Geometric metrics effectively characterize basis function properties.

Abstract

Representing a signal as a linear combination of a set of basis functions is central in a wide range of applications, such as approximation, de-noising, compression, shape correspondence and comparison. In this context, our paper addresses the main aspects of signal approximation, such as the definition, computation, and comparison of basis functions on arbitrary 3D shapes. Focusing on the class of basis functions induced by the Laplace-Beltrami operator and its spectrum, we introduce the diffusion and Laplacian spectral basis functions, which are then compared with the harmonic and Laplacian eigenfunctions. As main properties of these basis functions, which are commonly used for numerical geometry processing and shape analysis, we discuss the partition of the unity and non-negativity; the intrinsic definition and invariance with respect to shape transformations (e.g., translation, rotation, uniform scaling); the locality, smoothness, and orthogonality; the numerical stability with respect to the domain discretisation; the computational cost and storage overhead. Finally, we consider geometric metrics, such as the area, conformal, and kernel-based norms, for the comparison and characterisation of the main properties of the Laplacian basis functions.

Tables

Table 1. Table 1: Overview of the main properties (Sect.  1 ) of the basis functions discussed in this paper. The symbols  ∙ ∙ \bullet ,  ∘ \circ ,  ⊗ tensor-product \otimes indicate that the corresponding property is satisfied, is not satisfied, or might be satisfied according to additional conditions, which are detailed in Sect.  7 .

Property	Harm.	Hamil.	Lapl.	Lapl. compr.	Hamil.	Spectral	Diff.
	basis	basis	eigenf.	modes	eigenf.	basis	basis
Partition of unity	$\otimes$	$\otimes$	$\otimes$	$\otimes$	$\otimes$	$\otimes$	$\otimes$
Non-negativity	$\bullet$	$\circ$	$\circ$	$\circ$	$\circ$	$\circ$	$\bullet$
Intrinsic def.	$\bullet$	$\bullet$	$\bullet$	$\bullet$	$\bullet$	$\bullet$	$\bullet$
Locality	$\circ$	$\circ$	$\circ$	$\bullet$	$\circ$	$\otimes$	$\otimes$
Orthogonality	$\otimes$	$\otimes$	$\bullet$	$\bullet$	$\bullet$	$\otimes$	$\otimes$
Isometry-inv.	$\bullet$	$\circ$	$\bullet$	$\bullet$	$\circ$	$\bullet$	$\bullet$
Numer. stability	$\bullet$	$\bullet$	$\otimes$	$\otimes$	$\bullet$	$\bullet$	$\bullet$
Comput. cost $𝒪(\cdot )$	$n$	$n$	$kn\log n$	$kn\log n$	$kn\log n$	$rn\log n$	$rn\log n$
Storage overhead $𝒪(\cdot )$	$n$	$n$	$k{n}^2$	$kn$	$k{n}^2$	$n$	$n$

Equations

$$u(\mathbf{p})=(2\pi R{)}^{-1}{\int }_{\Gamma }u(s)\text{d}s=(\pi R^2{)}^{-1}{\int }_{\mathcal{B}}u(\mathbf{q})\text{d}\mathbf{q},$$

$$\left\{\begin{array}[]{ll}f=\sum_{n=0}^{+\infty}\alpha_{n}\phi_{n},&\alpha_{n}:=\langle f,\phi_{n}\rangle_{2};\\ \|f\|_{2}^{2}=\sum_{n=0}^{+\infty}\alpha_{n}^{2},&\|\nabla f\|_{2}^{2}=\sum_{n=0}^{+\infty}\alpha_{n}^{2}\lambda_{n}.\end{array}\right.$$

$$\Vert r_n{\Vert }_2^2\le \frac{\Vert \nabla f{\Vert }_2^2}{{\lambda }_{n+1}};$$

$${\Vert f-\sum _{i=0}^n⟨f,{\psi }_i{⟩}_2{\psi }_i\Vert }_2\le \alpha \frac{\Vert \nabla f{\Vert }_2}{{\lambda }_{n+1}},\forall f.$$

$$\Phi (f)=\sum _{n=0}^{+\infty }{\alpha }_n{\Delta }^{n}f=\sum _{n=0}^{+\infty }\phi ({\lambda }_n)⟨f,{\varphi }_n{⟩}_2{\varphi }_n,$$

$$K_{\phi }(\mathbf{p},\mathbf{q})=\sum _{n=0}^{+\infty }\phi ({\lambda }_n){\varphi }_n(\mathbf{p}){\varphi }_n(\mathbf{q}),$$

$$K_t(\mathbf{p},\cdot )=\sum _{n=0}^{+\infty }\exp (-{\lambda }_{n}t){\varphi }_n(\mathbf{p}){\varphi }_n.$$

$${\mathbf{K}}_{\phi }\mathbf{f}=\sum _{i=1}^n\phi ({\lambda }_i)⟨\mathbf{f},{\mathbf{x}}_i{⟩}_{\mathbf{B}}{\mathbf{x}}_i.$$

$$\begin{gathered} {\mathbf{K}}_{\phi }{\mathbf{e}}_i \\ \approx p_r(\tilde{\mathbf{L}}){\mathbf{e}}_i \\ =\sum _{j=1}^r{\alpha }_j{\left(\mathbf{I}+{\beta }_j\tilde{\mathbf{L}}\right)}^{-1}{\mathbf{e}}_i \\ =\sum _{j=1}^r{\alpha }_j{\mathbf{g}}_j, \end{gathered}$$

$$\left(\mathbf{B}+{\beta }_j\mathbf{L}\right){\mathbf{g}}_j=\mathbf{B}{\mathbf{e}}_i,j=1,\dots ,r.$$

$${\mathbf{K}}_{\phi }\mathbf{f}\approx {\alpha }_0\mathbf{f}+\sum _{i=1,\dots ,n}^{j=1,\dots ,r}{\alpha }_j\frac{1}{1+{\beta }_j{\lambda }_i}⟨\mathbf{f},{\mathbf{x}}_i{⟩}_{\mathbf{B}}{\mathbf{x}}_i.$$

$$\kappa_{2}(\mathbf{B}+\beta_{i}\mathbf{L})=\frac{\max_{j}\{1+\beta_{i}\lambda_{j}\}}{\min_{j}\{1+\beta_{i}\lambda_{j}\}}=\left\{\begin{array}[]{ll}1+\beta_{i}\lambda_{n}&\beta_{i}>0,\\ (1+\beta_{i}\lambda_{n})^{-1}&\beta_{i}<0,\end{array}\right.$$

$$\left\{\begin{array}[]{ll}G_{\Delta}(\mathbf{p},\mathbf{q})=\sum_{n=0}^{+\infty}\frac{1}{\lambda_{n}}\phi_{n}(\mathbf{p})\phi_{n}(\mathbf{q})&\textrm{harmonic oper. }(\mathcal{L}=\Delta);\\ K_{t}(\mathbf{p},\mathbf{q})=\sum_{n=0}^{+\infty}\exp(-t\lambda_{n})\phi_{n}(\mathbf{p})\phi_{n}(\mathbf{q})&\textrm{diffusion oper. }(\mathcal{L}=\exp(t\Delta));\\ G_{\varphi}(\mathbf{p},\mathbf{q})=\sum_{n=0}^{+\infty}\varphi(\lambda_{n})\phi_{n}(\mathbf{p})\phi_{n}(\mathbf{q})&\textrm{spectral oper. }(\mathcal{L}=\varphi^{{\dagger}}(\Delta)).\end{array}\right.$$

$$h_a({\varphi }_i,{\varphi }_j)={\delta }_{ij},h_c({\varphi }_i,{\varphi }_j)=⟨{\varphi }_i,\Delta {\varphi }_j{⟩}_2={\lambda }_j{\delta }_{ij}.$$

$$h_c(f,g)=⟨f,\Delta g{⟩}_2\approx {\mathbf{f}}^{\top }\mathbf{B}\tilde{\mathbf{L}}\mathbf{g}=⟨\mathbf{f},\mathbf{g}{⟩}_{\mathbf{L}}.$$

$$⟨f,g⟩:={\int }_{\mathcal{N}\times \mathcal{N}}K(\mathbf{p},\mathbf{q})f(\mathbf{p})g(\mathbf{q})\text{d}\mathbf{p}\text{d}\mathbf{q}.$$

$$⟨f,g⟩=\sum _{n=0}^{+\infty }\phi ({\lambda }_n)⟨f,{\varphi }_n{⟩}_2⟨g,{\varphi }_n{⟩}_2.$$

Full text

††journal: CAGD - New Progress in Geometry for Applications

1 Introduction

In Computer Graphics, scalar functions are ubiquitous to represent the values of a physical phenomenon (e.g., the heat that diffuses from one or more source points), a shape descriptor (e.g., curvature, spherical harmonics), a distance (e.g., geodesic, bi-harmonic, diffusion distance) from a set of source points, or the pixels of a surface texture.

In all these cases, representing the input function in terms of a basis allows us to address a large number of applications, such as multi-resolution representations by selecting a set of multi-scale basis functions (e.g., Laplacian eigenfunctions, diffusion basis functions), sparse representations by choosing a low number of basis functions in order to achieve a target approximation accuracy, or compression by quantising the representation coefficients. We can also address deformation by modifying the coefficients that express the geometry of the input surface in terms of geometry-driven or shape-intrinsic basis functions, smoothing by neglecting the coefficients associated with large Laplacian frequencies, and the definition of Laplacian spectral kernels and distances as a filtered combination of the Laplacian spectral eigenfunctions.

In this context, our paper focuses on the main aspects of signal approximation on arbitrary 3D shapes, such as the definition, computation, and comparison of basis functions for applications in numerical geometry processing and shape analysis (Fig. 1). To define the functional space associated with an input 3D shape, we select a set of its basis functions and then represent any signal as a linear combination of these functions. In this work, we focus on the class of bases induced by the Laplace-Beltrami operator and its spectrum; i.e., the harmonic basis functions (Sect. 3), as solution to the Laplace equation; the Laplacian eigenfunctions (Sect. 4) associated with the spectrum of the Laplace-Beltrami operator; and the Laplacian spectral basis functions (Sect. 5), which are defined by filtering the Laplacian spectrum and include the diffusion basis functions.

The diffusion and Laplacian spectral functions are novel classes of functions, which will be defined by further developing our recent results on the Laplacian spectral distances [44, 43, 42]. Then, these functions will be compared with the harmonic functions and the Laplacian eigenfunctions. As a main contribution with respect to the previous work, the diffusion and Laplacian spectral basis functions are intrinsic to the input shape and local; i.e., they have a compact support that can be tuned easily by selecting the diffusion scale and the decay of the filter function to zero, respectively.

Defining basis functions with a compact support is particularly important to localise their behaviour around feature points and to reduce their storage overhead in the discrete setting. In fact, given a discrete domain $\mathcal{M}$ with $n$ points, the space of functions on $\mathcal{M}$ has dimension $n$ and its basis functions are encoded as a $n\times n$ matrix. This matrix can be stored only if we work with compactly-supported or analytically defined (e.g., radial basis or polynomial) functions, or if we select a subset of $k$ basis functions, with $k$ much smaller than $n$.

For the comparison of functions on the same surface (Sect. 6), we review different metrics (e.g., area, conformal, kernel-based metrics), which measure differential properties of the input scalar function and geometric properties of the underlying domain. As main properties, we discuss the partition of the unity and non-negativity; the intrinsic definition and invariance with respect to shape transformations (e.g., translation, rotation, uniform scaling); the locality, smoothness, and orthogonality; the numerical stability with respect to the domain discretisation; the computational cost and storage overhead.

Our experiments (Sect. 7, 8) show that the diffusion and the Laplacian spectral functions provide a valid alternative to the harmonic functions and Laplacian eigenfunctions. In fact, the diffusion basis functions can be centred at any point of the input domain, have a compact support, and encode local/global shape details according to the values of the temporal parameter or to the decay of the selected filter to zero. Similar properties also apply to the Laplacian spectral basis functions induced by a filter with a strong (e.g., exponential) decay to zero.

2 Laplace-Beltrami operator

Let $\mathcal{N}$ be a smooth surface, possibly with boundary, equipped with a Riemannian metrics and let us consider the inner product $\langle f,g\rangle_{2}:=\int_{\mathcal{N}}f(\mathbf{p})g(\mathbf{p})\textrm{d}\mathbf{p}$ defined on the space $\mathcal{L}^{2}(\mathcal{N})$ of square integrable functions on $\mathcal{N}$ and the corresponding norm $\|\cdot\|_{2}$. The Laplace-Beltrami operator $\Delta:=-\textrm{div}(\textrm{grad})$ is self-adjoint $\langle\Delta f,g\rangle_{2}=\langle f,\Delta g\rangle_{2}$, and positive semi-definite; i.e., $\langle\Delta f,f\rangle_{2}\geq 0$, $\forall f,g$, [48]. Finally, the value $\Delta f(\mathbf{p})$ does not depend on $f(\mathbf{q})$, for any couple of distinct points $\mathbf{p}$, $\mathbf{q}$ (locality).

According to [43], we represent the Laplace-Beltrami operator on surface and volume meshes in a unified way as $\tilde{\mathbf{L}}:=\mathbf{B}^{-1}\mathbf{L}$, where the mass matrix $\mathbf{B}$ is sparse, symmetric, positive definite, and the stiffness matrix $\mathbf{L}$ is sparse, symmetric, positive semi-definite, and $\mathbf{L}\mathbf{1}=\mathbf{0}$. Analogously to the continuous case, the Laplacian matrix satisfies the following properties:

$\mathbf{B}$-adjointness: $\tilde{\mathbf{L}}$ is adjoint with respect to the $\mathbf{B}$-inner product $\langle\mathbf{f},\mathbf{g}\rangle_{\mathbf{B}}:=\mathbf{f}^{\top}\mathbf{B}\mathbf{g}$; i.e., $\langle\tilde{\mathbf{L}}\mathbf{f},\mathbf{g}\rangle_{\mathbf{B}}=\langle\mathbf{f},\tilde{\mathbf{L}}\mathbf{g}\rangle_{\mathbf{B}}=\mathbf{f}^{\top}\mathbf{L}\mathbf{g}$. If $\mathbf{B}:=\mathbf{I}$, then this property reduces to the symmetry of $\mathbf{L}$;

• 2.

positive semi-definiteness: $\langle\tilde{\mathbf{L}}\mathbf{f},\mathbf{f}\rangle_{\mathbf{B}}=\mathbf{f}^{\top}\mathbf{L}\mathbf{f}\geq 0$. In particular, the Laplacian eigenvalues are positive, with one null eigenvalue associated with a constant eigenvector;

• 3.

locality: assuming that $\mathbf{B}$ is diagonal, the value $(\tilde{\mathbf{L}}\mathbf{f})_{i}$ depends only on the $f$-values at $\mathbf{p}_{i}$ and its $1$-star neighbourhood.

On triangle meshes, for the stiffness matrix we can consider the linear FEM Laplacian weights [47, 54] that reduce to the Voronoi-cotangent weights [18], by lumping the FEM mass matrix, and extend the cotangent weights [45]. The mean-value weights [23] have been derived from the mean-value theorem for harmonic functions and are always positive. In [17], the weak formulation of the Laplacian eigenproblem is achieved by selecting a set of volumetric test functions, which are defined as $k\times k\times k$ B-splines (e.g., $k:=4$) and restricted to the input shape. For the anisotropic Laplacian [5], the entries of $\mathbf{L}$ are a variant of the cotangent weights (i.e., with respect to different angles) and the entries of the diagonal mass matrix are the areas of the Voronoi regions.

On polygonal meshes, the Laplacian discretisation in [3, 28] generalises the Laplacian matrix with cotangent weights to surface meshes with non-planar, non-convex faces. An approximation of the Laplace-Beltrami operator with point-wise convergence has been proposed in [11].

In [8, 9, 10, 12], the Laplace-Beltrami operator on a point set has been defined as the Gram matrix associated with the exponential kernel. Starting from this approach, a new discretisation [38] has been achieved through a finer approximation of the local geometry of the surface at each point through its Voronoi cell. The discretisation of the Laplace-Beltrami operator on volumes represented as a tetrahedral mesh has been addressed in [4, 36, 53].

Finally, in the paper examples we apply the linear FEM weights and the level-sets of a given function are associated with iso-values uniformly sampled in its range.

3 Harmonic basis

Given a set $\mathcal{S}:=\{\mathbf{p}_{i}\}_{i\in\mathcal{I}}$ of seed points (Sect. 7.1), we compute the harmonic basis functions $\mathcal{B}:=\{\psi_{i}\}_{i\in\mathcal{I}}$, as solution to the harmonic equation $\Delta\psi_{i}=0$, $\psi_{i}(\mathbf{p}_{j})=\delta_{ij}$, $i,j\in\mathcal{I}$. Each harmonic function, and in particular any harmonic basis function,

minimizes the Dirichlet energy $\mathcal{E}(u):=\int_{\mathcal{N}}\|\nabla u(\mathbf{p})\|_{2}^{2}\textrm{d}\mathbf{p}$;

• 2.

satisfies the locality property; i.e., if $\mathbf{p}$ and $\mathbf{q}$ are two distinct points, then $\Delta u(\mathbf{p})$ is not affected by the value of $u$ at $\mathbf{q}$;

• 3.

verifies the mean-value theorem

[TABLE]

where $\mathcal{B}\subseteq\mathcal{N}$ is a disc of center $\mathbf{p}$, radius $R$, and boundary $\Gamma$ (mean-value theorem);

• 4.

satisfies the strong maximum principle [48]. More precisely, if $\Omega$ is a connected open set, $u\in\mathcal{C}^{2}(\Omega)$, $u$ is harmonic and attains either a global minimum or maximum in $\Omega$, then $u$ is constant.

The maximum principle states that the values of a harmonic function in a compact domain are bounded by its maximum and minimum values on the boundary. In particular, the values of the harmonic basis function $\psi_{i}$ belong to the interval $[0,1]$ and its maximum is attained at the selected seed point $\mathbf{p}_{i}$.

The discrete harmonic basis associated with a seed point $\mathbf{p}_{i}$ is the solution to the linear system $\mathbf{L}^{\star}\mathbf{g}=\mathbf{e}_{i}$, where $\mathbf{L}^{\star}$ is achieved by replacing the $i$-th row of the stiffness matrix with the vector $\mathbf{e}_{i}$, where $\mathbf{e}_{i}^{(j)}:=\delta_{ij}$ (Fig. 2). Harmonic functions are efficiently computed in $\mathcal{O}(n)$ time with iterative solvers of sparse linear systems; their computation is stable for the mean-value weights while negative Voronoi cotangent weights might induce local undulations in the resulting harmonic function. Finally, the (globally-supported) geometry-aware basis [52] has been defined by minimising the Dirichlet energy, with interpolating and least-squares constraints.

Pre-conditioners and fast solvers of elliptic problems

In the context of the iso-geometric analysis of elliptic PDEs, previous work has proposed (i) additive multilevel preconditioners [15], whose spectral conditioning number is independent of the grid spacing, and (ii) preconditioning methods based on the solution of a Sylvester-like equation [50], which are robust with respect to the mesh size and the spline degree. In the context of Computer Graphics, specialised pre-conditioners for the Laplacian matrix [33] can be applied to further attenuate numerical instabilities.

As solvers, we can select (i) robust and efficient multi-grid methods for single-patch iso-geometric discretisations of elliptic problems [30], based on tensor product B-splines of maximum smoothness, and (ii) multi-grid methods [20] for the Galerkin iso-geometric analysis based on B-splines and with optimal convergence rate (i.e., bounded independently of the matrix size). As alternatives, we mention the fast solvers for large linear systems arising from the Galerkin approximation based on B-splines [19], which are linear with respect to the matrix size and have a convergence speed that is independent of the matrix size, the spline degree, and the dimensionality of the input problem.

For the analysis of the main spectral properties (e.g., non-singularity, conditioning, spectral distribution of the eigenvalue) of the stiffness matrix related to the iso-geometric analysis of elliptic problems, we refer the reader to [25].

4 Laplacian and Hamiltonian eigenbasis

Recalling that the Laplace-Beltrami operator is self-adjoint and positive semi-definite (Sect. 2), it has an orthonormal eigensystem $\mathcal{B}:=\{(\lambda_{n},\phi_{n})\}_{n=0}^{+\infty}$, $\Delta\phi_{n}=\lambda_{n}\phi_{n}$, with $0=\lambda_{0}\leq\lambda_{n}\leq\lambda_{n+1}$, in $\mathcal{L}^{2}(\mathcal{N})$ (Sect. 2). Any function $f$ in $\mathcal{L}^{2}(\mathcal{N})$ satisfies the following relations

[TABLE]

The Laplacian eigenfunctions are intrinsic to the input shape and those ones related to smaller eigenvalues correspond to smooth and slowly-varying functions. Increasing the eigenvalues, the corresponding eigenfunctions generally show rapid oscillations. If two shapes are isometric, then they have the same Laplacian spectrum (iso-spectral property); however, the viceversa does not hold [27, 56] and we cannot recover the metric of a given surface.

Hamiltonian basis and eigenbasis

A Hamiltonian operator $\mathcal{H}:\mathcal{L}^{2}(\mathbb{R})\rightarrow\mathbb{R}$ is defined as $\mathcal{H}f=-\Delta f+\mu Vf$, where $V:\mathcal{N}\rightarrow\mathbb{R}$ is a potential function and $\mu$ is a real parameter that controls the influence of the potential on the Hamiltonian operator and defines the trade-off between the local and global support of the corresponding eigenfunctions. Analogously to the harmonic basis and Laplacian eigenbasis, we can define the Hamiltonian basis, as the solution to the problem $\mathcal{H}\psi_{i}=0$, $\psi_{i}(\mathbf{p}_{j})=\delta_{ij}$, and the Hamiltonian eigenbasis $\mathcal{H}\psi_{n}=\sigma_{n}\psi_{n}$, $n\in\mathbb{N}$.

Optimality of the Laplacian eigenbasis

According to [2], the Laplacian eigenfunctions represent an optimal basis for the representation of signals with bounded gradient magnitude. In fact,

the spectral decomposition $f=\sum_{i=0}^{n}\langle f,\phi_{i}\rangle_{2}\phi_{i}$ is optimal in approximating functions with $\mathcal{L}^{2}$ bounded gradient magnitude; i.e., the residual error $r_{n}:=f-f_{n}$ is bounded as

[TABLE]

• 2.

the spectral decomposition is optimal in approximating functions with respect to the error estimates in (1). In fact, for any $0\leq\alpha<1$ there is no integer $n$ and no sequence $(\psi_{i})_{i=0}^{n}$ of linearly independent functions in $\mathcal{L}^{2}$ such that

[TABLE]

Discrete Laplacian eigenbasis

The generalised Laplacian eigensystem $(\lambda_{i},\mathbf{x}_{i})_{i=1}^{n}$, with $0=\lambda_{1}<\lambda_{2}\leq\cdots\leq\lambda_{n}$, satisfies the identity $\mathbf{L}\mathbf{x}_{i}=\lambda_{i}\mathbf{B}\mathbf{x}_{i}$ and the eigenvectors are orthonormal with respect to the $\mathbf{B}$-inner product; i.e., $\langle\mathbf{x}_{i},\mathbf{x}_{j}\rangle_{\mathbf{B}}=\mathbf{x}_{i}^{\top}\mathbf{B}\mathbf{x}_{j}=\delta_{ij}$. In particular, $\tilde{\mathbf{L}}=\mathbf{X}\Gamma\mathbf{X}^{\top}\mathbf{B}$ is the spectral decomposition of the Laplacian matrix, where $\mathbf{X}$ is the eigenvectors’ matrix and $\Gamma$ is the diagonal matrix of the eigenvalues. Fig. 3 shows the level-sets of three eigenfunctions. In [7, 41], Laplacian eigenfunctions with a compact support (compressed manifold modes) have been defined by solving an orthogonality-constrained optimisation problem with a $\mathcal{L}^{1}$ penalty term. In [32], the compact support of the Laplacian eigenfunctions is enforced with a $\mathcal{L}^{p}$, $0<p<1$, penalisation term; then, these eigenfunctions are computed by combining splitting strategies with an ADMM-based (Alternating Direction Method of Multipliers) iterative algorithm. Finally, the optimality of the Laplacian eigenfunctions for spectral geometry processing and for signal approximation has been addressed in [13] and in [1, 2], respectively.

Computation of the Laplacian eigenbasis

For the computation of the Laplacian eigenfunctions, numerical methods generally exploit the sparsity of the Laplacian matrix and reduce the high-dimensional eigenproblem to one of lower dimension, by applying a coarsening step. Main examples include the algebraic multi-grid method [22], Arnoldi iterations [35, 51], and the Nystrom method [24]. Even though the eigenvalues and eigenvectors are computed in super-linear time [54], this computational cost and the required $\mathcal{O}(n^{2})$ storage are expensive for densely sampled domains. Indeed, modifications of the Laplacian eigenproblem are applied to evaluate specific sub-parts of the Laplacian spectrum (e.g., shift, power, and inverse methods). Finally, the Laplacian eigenvectors are computed only for a small set of eigenvalues and do not provide a flexible alignment of the function behaviour to specific shape features.

Stability of the Laplacian eigenbasis

Theoretical results on the sensitivity of the Laplacian spectrum against geometry changes, irregular sampling density and connectivity have been presented in [29, 55]. Here, we briefly recall that the instability of the computation of the Laplacian eigenbasis is generally due to repeated or close eigenvalues, with respect to the numerical accuracy of the solver of the eigen-equation [43]. While repeated eigenvalues are quite rare and typically associated with symmetric shapes, numerically close or switched eigenvalues can be present in the spectrum, in spite of the regularity of the input discrete surface.

5 Laplacian spectral basis

We address the definition (Sect. 5.1) and computation (Sect. 5.2) of the Laplacian spectral basis by filtering the Laplacian spectrum. Then, we study the relations between the Laplacian spectral basis and the Green kernel (Sect. 5.3).

5.1 Spectral basis

Given a positive filter $\varphi:\mathbb{R}^{+}\rightarrow\mathbb{R}^{+}$, let us consider the power series $\varphi(s)=\sum_{n=0}^{+\infty}\alpha_{n}s^{n}$. Noting that $\Delta^{i}f=\sum_{n=0}^{+\infty}\lambda_{n}^{i}\langle f,\phi_{n}\rangle_{2}\phi_{n}$ and under a few hypotheses on the decay of the filter to zero [44], we define the spectral operator as

[TABLE]

which is linear, continuous, and $\Phi(f)=\langle K_{\varphi},f\rangle_{2}$, where

[TABLE]

is the spectral kernel. Then, the spectral basis function centred at a point $\mathbf{p}$ is defined as the action of the spectral operator on $\delta_{\mathbf{p}}$ or equivalently as the spectral kernel centred at $\mathbf{p}$; in fact, $\Phi(\delta_{\mathbf{p}})=K_{\varphi}(\mathbf{p},\cdot)$.

The filters $\varphi_{t}(s):=\exp{(-ist)}$, $\varphi(s):=s^{-k/2}$, $\varphi(s):=s^{-1/2}$ induce the wave [14, 6], poly-harmonic, and commute-time basis functions, respectively. Mexican hat wavelets [31] are generated by the filter $\varphi(s):=s^{1/2}\exp(-s^{2})$ and the filter function $\varphi(s):=\exp(is)$, $s\in[0,2\pi]$, defines the wave kernel. In $\varphi(s):=t^{k}s^{\alpha}\exp(-ts^{\alpha})$, the parameter $k$ scales the rate of diffusion and $\alpha$ controls the decay of the Laplacian eigenvalues to zero, similarly to random walks [46].

Diffusion basis

Selecting the filter $\varphi_{t}(s):=\exp(-ts)$, the spectral operator reduces to the diffusion operator $\Phi_{t}:=\exp(-t\Delta)$ and the diffusion basis function at $\mathbf{p}$ (Fig. 4, Fig. 5) is the heat kernel $K_{t}(\cdot,\mathbf{p})$ centred at $\mathbf{p}$, which solves the heat diffusion equation $(\partial_{t}+\Delta)H(\cdot,t)=0$, $H(\cdot,0)=\delta_{\mathbf{p}}$; i.e.,

[TABLE]

The spectral representation (4) shows the smoothing effect on the initial condition; as the scale increases, the component of $\delta_{\mathbf{p}}$ along the eigenfunctions associated with the larger Laplacian eigenvalue becomes null.

Assuming that $\Omega$ is an open and bounded set and that the solution $F(\cdot,t)$ to the heat equation is sufficiently smooth, we define the parabolic cylinder $\Omega_{T}:=\Omega\times(0,T]$ and the parabolic boundary $\Gamma_{T}:=\overline{\Omega}_{T}\backslash\Omega_{T}$. According to the strong maximum principle for the heat diffusion equation [34] (Ch. $2$), if $\Omega$ is connected and there exists a point $(\mathbf{p}_{0},t_{0})$ in $\Omega_{T}$ such that $H(\mathbf{p}_{0},t_{0})=\max_{\overline{\Omega_{T}}}H(\mathbf{p},t)$, then $H(\cdot,t)$ is constant in $\overline{\Omega_{T}}$. In particular, the diffusion basis functions are positive.

5.2 Computation of the spectral basis

Considering the generalised eigensystem $\mathbf{L}\mathbf{X}=\mathbf{B}\mathbf{X}\Gamma$, with orthonormal eigenvectors $\mathbf{X}^{\top}\mathbf{B}\mathbf{X}=\mathbf{I}$, the discretisation of the spectral operator (2) is $\mathbf{K}_{\varphi}=\varphi(\tilde{\mathbf{L}})=\mathbf{X}\varphi(\Gamma)\mathbf{X}^{\top}\mathbf{B}$, $\varphi(\Gamma):=\textrm{diag}(\varphi(\lambda_{i}))_{i=1}^{n}$; i.e.,

[TABLE]

In particular, $\tilde{\mathbf{L}}$ and $\mathbf{K}_{\varphi}$ have the same eigenvectors and $(\varphi(\lambda_{i}))_{i=1}^{n}$ are the (filtered) Laplacian eigenvalues of $\mathbf{K}_{\varphi}$. Assuming that $\mathbf{K}_{\varphi}:=\varphi(\tilde{\mathbf{L}})$ is invertible (i.e., the filter does not vanish at the Laplacian eigenvalues), the vectors in $\mathcal{B}:=\{\mathbf{K}_{\varphi}\mathbf{e}_{i}\}_{i=1}^{n}$ (i.e., $\mathbf{f}:=\mathbf{e}_{i}$ in Eq. (5)) are linearly independent; in fact, $\sum_{i=1}^{n}\alpha_{i}\mathbf{K}_{\varphi}\mathbf{e}_{i}=\mathbf{0}$ is verified if and only if $\mathbf{K}_{\varphi}\alpha=\mathbf{0}$; i.e., $\alpha=\mathbf{0}$. The coefficients $\alpha:=(\alpha_{i})_{i=1}^{n}$, which express $\mathbf{f}:=\sum_{i=1}^{n}\alpha_{i}\mathbf{K}_{\varphi}\mathbf{e}_{i}$ as a linear combination of the basis $\mathcal{B}$, are computed by solving the linear system $\mathbf{K}_{\varphi}\alpha=\mathbf{f}$. We notice that $\alpha=\mathbf{K}_{\varphi}^{{\dagger}}\mathbf{f}=\mathbf{X}\varphi^{{\dagger}}(\Gamma)\mathbf{X}^{\top}\mathbf{B}\mathbf{f}$, where $\varphi^{{\dagger}}(\Gamma)$ is the diagonal matrix whose non-null values are $(1/\varphi(\lambda_{i}))_{i=1}^{n}$.

Recalling that the computation of the Laplacian eigenpairs is numerically unstable in case of repeated eigenvalues [43], the filter function should be chosen in such a way that the filtered Laplacian matrix does not have additional (if any) repeated eigenvalues. This condition is generally satisfied by choosing an injective filter. The selection of periodic filters, the expensive cost of the computation of the Laplacian spectrum, and the sensitiveness of multiple Laplacian eigenvalues to surface discretisation motivate the definition of alternative approaches for the evaluation of the spectral basis functions. Among them, we discuss the truncated and spectrum-free approximations.

Truncated approximation

The computational limits for the evaluation of the whole Laplacian spectrum and the decay of the filtered coefficients are the main reasons behind the approximation of a spectral basis function as a truncated sum; i.e., applying the relation $\Phi_{k}\mathbf{e}_{i}=\sum_{j=1}^{k}\varphi(\lambda_{j})\langle\mathbf{e}_{i},\mathbf{x}_{j}\rangle_{\mathbf{B}}\mathbf{x}_{j}$, where $k$ is the selected number of eigenpairs. Even though the first $k$ Laplacian eigenpairs are computed in super-linear time [54], the evaluation of the whole Laplacian spectrum is unfeasible for storage and computational cost, which are quadratic in the number of surface samples. The selection of filters that are periodic or do not decrease to zero motivates the need of defining a spectrum-free computation of the corresponding kernels, which cannot be accurately approximated with the contribution of only a subpart of the Laplacian spectrum. Furthermore, the number of selected eigenpairs is heuristically adapted to the decay of the filter function and the approximation accuracy cannot be estimated without computing the whole spectrum.

Padé-Chebyshev (spectrum-free) approximation

Truncated spectral representations apply a low pass filter and the resulting reconstructed signal generally preserves only its global features. Local features, which are associated with a high frequency, are generally missing in the approximated signals and can be recovered by considering only a high number of Laplacian eigenpairs. However, the evaluation of a large part of the Laplacian spectrum is computationally unfeasible, requires a large memory overhead, and is numerically unstable in case of duplicated or numerically close eigenvalues [43].

To avoid these drawbacks, we introduce the spectrum-free evaluation of the basis functions, which is based on a rational approximation of the filter.

For an arbitrary filter, we consider the rational Padé-Chebyshev approximation $p_{r}(s)=\frac{a_{r}(s)}{b_{r}(s)}$ of $\varphi$ [26] (Ch. $11$) with respect to the $\mathcal{L}^{\infty}$ norm. Here, $a_{r}(\cdot)$ and $b_{r}(\cdot)$ are polynomials of degree equal to or lower than $r$. Let $p_{r}(s)=\sum_{i=1}^{r}\alpha_{i}(1+\beta_{i}s)^{-1}$ be the partial form of the Padé-Chebyshev approximation, where $(\alpha_{i})_{i=1}^{r}$ are the weights and $(\beta_{i})_{i=1}^{r}$ are the nodes of the $r$-point Gauss-Legendre quadrature rule [26] (Ch. $11$). The weights and nodes are precomputed for any degree of the rational polynomial [16]. Applying this approximation to the spectral kernel, we get that

[TABLE]

where $\mathbf{g}_{j}$ solves the symmetric and sparse linear system

[TABLE]

The Padé-Chebyshev approximation generally provides an accuracy higher than the polynomial approximation, as a matter of its uniform convergence to the filter. In the paper examples, we have selected $r:=5$ as degree of the rational polynomial approximation of the filter function.

We now derive the spectral representation of the Padè-Chebyshev approximation of the function $\mathbf{K}_{\varphi}\mathbf{f}$. Re-writing the solution to $(\mathbf{B}+\beta_{j}\mathbf{L})\mathbf{g}=\mathbf{B}\mathbf{f}$ (i.e., Eq. (6) with $\mathbf{f}=\mathbf{e}_{i}$) as $\mathbf{g}=\mathbf{X}\alpha$, we get that $\alpha=(\mathbf{I}+\beta_{j}\Gamma)^{{\dagger}}\mathbf{X}^{\top}\mathbf{B}\mathbf{f}$ and the Padè-Chebyshev approximation is

[TABLE]

Setting $\mathbf{f}:=\mathbf{e}_{i}$, we obtain that the spectral representation of the Padè-Chebyshev approximation of the spectral basis function $\mathbf{K}_{\varphi}\mathbf{e}_{i}$ centred at $\mathbf{p}_{i}$.

According to [40], the approximation of the matrix $\varphi(\tilde{\mathbf{L}})$ might be numerically unstable if $\|\tilde{\mathbf{L}}\|_{2}$ is large. From the bound $\|\mathbf{B}^{-1}\mathbf{L}\|_{2}\leq{\lambda_{\min}^{-1}(\mathbf{B})}\lambda_{\max}(\mathbf{L})$, a well-conditioned mass matrix $\mathbf{B}$ guarantees that $\|\mathbf{B}^{-1}\mathbf{L}\|_{2}$ is bounded. Noting that $(1+\beta_{i}\lambda_{j},\mathbf{x}_{j})_{j=1}^{n}$ are the eigenpairs of $\mathbf{B}+\beta_{i}\mathbf{L}$, the corresponding conditioning number is bounded as

[TABLE]

which is always greater than 1. The numerical solution to the linear system in Eq. (6) is generally stable; as reviewed in Sect. 3, pre-conditioners can be applied to further attenuate numerical instabilities.

The truncated spectral approximation (Fig. 6) is affected by small geometric undulations (especially at small scales), the use of heuristics for the selection of the number of Laplacian eigenpairs with respect to the target approximation accuracy, and the scale of features of the input shape. The spectrum-free computation generally provides better results in terms of smoothness, regularity, and accuracy of the computed spectral basis. Changing the filter function and its decay to zero allows us to define different basis functions with a different behaviour and support (Fig. 7), thus showing the flexibility and simplicity of our approach. For further comparison examples, we refer the reader to Sect. 7.

5.3 Green kernel and basis functions

To study the relation between the Green kernel and the basis previously defined, let us consider a linear differential operator $\mathcal{L}$ and the corresponding Green kernel $G:\mathcal{N}\times\mathcal{N}\rightarrow\mathbb{R}$ such that $\mathcal{L}G(\mathbf{p},\mathbf{q})=\delta(\mathbf{p}-\mathbf{q})$, where $\delta(\cdot)$ is the delta function. Then, the solution to the differential equation $\mathcal{L}u=f$ is expressed in terms of the Green kernel as $u(\mathbf{p})=\langle G(\mathbf{p},\cdot),f\rangle_{2}$. Noting that the eigensystem $(\lambda_{n},\phi_{n})_{n=0}^{+\infty}$ of the integral operator $\mathcal{K}_{G}f:=\langle G(\mathbf{p},\cdot),f\rangle_{2}$ induced by the Green kernel satisfies the relations $\mathcal{L}\phi_{n}=\lambda_{n}^{-1}\phi_{n}$, $\lambda_{n}\neq 0$, we have that the eigensystem of $\mathcal{L}$ is $(\lambda_{n}^{-1},\phi_{n})_{n=0}^{+\infty}$. Indeed, $\mathcal{L}$, $\mathcal{K}_{G}$ have the same eigenfunctions and reciprocal eigenvalues. In particular, the spectral representation of the Green kernel $G(\mathbf{p},\mathbf{q})=\sum_{n=0}^{+\infty}\lambda_{n}\phi_{n}(\mathbf{p})\phi_{n}(\mathbf{q})$ is uniquely defined by the spectrum of $\mathcal{L}$.

Combining the previous result with the fact that the spectral basis functions have been defined as solution to different differential equations that involve the Laplace-Beltrami operator, we get the following relations

[TABLE]

6 Comparison metrics for scalar functions

Since the level-sets and critical points (i.e., maxima, minima, saddles) of a scalar function $f:\mathcal{N}\rightarrow\mathbb{R}$ are independent of positive re-scalings of $f$, we can assume that the function values have been normalised in such a way that $\textrm{Image}(f)=[0,1]$ and compare two functions by evaluating the corresponding $\mathcal{L}_{\infty}$ or $\mathcal{L}_{2}$ norm. However, these norms do not measure the differential properties of the underlying functions and are better suited to evaluate the approximation accuracy. Indeed, in the following we disucss the area, conformal, and kernel-based metrics.

Area and conformal metric

To compare two functions $f,g:\mathcal{N}\rightarrow\mathbb{R}$ on the same surface $\mathcal{N}$, we consider the

area distortion: $h_{a}(f,g):=\int_{\mathcal{N}}f(\mathbf{p})g(\mathbf{p})\textrm{d}\mathbf{p}$;

• 2.

conformal distortion: $h_{c}(f,g):=\int_{\mathcal{N}}\langle\nabla f(\mathbf{p}),\nabla g(\mathbf{p})\rangle_{2}\textrm{d}\mathbf{p}$.

The area and conformal distortions are invariant with respect to area-preserving and conformal transformations, respectively. In fact, let $\mathcal{N}$, $\mathcal{Q}$ be 2 surfaces and $\mathcal{F}(\mathcal{Q})$ the space of scalar functions defined on $\mathcal{Q}$. Given a map $T:\mathcal{N}\rightarrow\mathcal{Q}$, we have that [49] for all $f,g\in\mathcal{F}(\mathcal{Q})$, $h_{a}^{\mathcal{Q}}(f,g)=h_{a}^{\mathcal{N}}(f\circ T,g\circ T)$ if and only if $T$ is locally area-preserving and $h_{c}^{\mathcal{Q}}(f,g)=h_{c}^{\mathcal{N}}(f\circ T,g\circ T)$ if and only if $T$ is conformal. For instance, for the Laplacian eigenbasis we have that

[TABLE]

The discrete area distortion is $h_{a}(f,g)=\langle\mathbf{f},\mathbf{g}\rangle_{\mathbf{B}}=\mathbf{f}^{\top}\mathbf{B}\mathbf{g}$, where $\mathbf{B}$ is the area-driven mass matrix (Sect. 2), and the discrete conformal distortion is

[TABLE]

Expressing the input functions as a linear combination of the Laplacian eigenbasis with coefficients $\alpha=\mathbf{X}^{\top}\mathbf{B}\mathbf{f}$, $\beta=\mathbf{X}^{\top}\mathbf{B}\mathbf{g}$, we have that $h_{a}(f,g)=\alpha^{\top}\beta$ and $h_{c}(f,g)=\alpha^{\top}\Gamma\beta$. For example, the area and conformal distortion of the delta functions at $\mathbf{p}_{i}$, $\mathbf{p}_{j}$ reduce to the entry $(i.j)$ of $\mathbf{B}$ and $\mathbf{L}$, respectively. For the Laplacian eigenvectors, we have that $h_{a}(\mathbf{x}_{i},\mathbf{x}_{j})=\delta_{ij}$ and $h_{c}(\mathbf{x}_{i},\mathbf{x}_{j})=\lambda_{j}\delta_{ij}$. Fig. 8 shows the area and conformal metrics of the diffusion basis functions at 100 randomly sampled seed points and at 4 scales.

Kernel-based metric

Any symmetric and positive kernel $K:\mathcal{N}\times\mathcal{N}\rightarrow\mathbb{R}$ (e.g., the Gaussian kernel) induces the following kernel-based inner product

[TABLE]

Combining the spectral representation (3) of $K_{\varphi}$ with the definition of the kernel-based inner product (7), we get that

[TABLE]

The discretisation of (7) is $\langle\mathbf{f},\mathbf{g}\rangle=\mathbf{f}^{\top}\mathbf{B}\mathbf{K}\mathbf{g}\\ =\langle\mathbf{f},\mathbf{K}\mathbf{g}\rangle_{\mathbf{B}}$, where $\mathbf{K}:=(K(\mathbf{p}_{i},\mathbf{p}_{j}))_{i,j=1}^{k}$ is the Gram matrix associated with the kernel $K(\cdot,\cdot)$. To guarantee the symmetry of this inner product, it is enough to assume that $\mathbf{K}$ is positive definite and  $\mathbf{B}$-adjoint; i.e., $\langle\mathbf{f},\mathbf{K}\mathbf{g}\rangle_{\mathbf{B}}=\langle\mathbf{K}\mathbf{f},\mathbf{g}\rangle_{\mathbf{B}}$. As main examples of kernels, we mention the filtered spectral kernel (e.g., the diffusion kernel), which is $\mathbf{B}$-adjoint, and the Laplacian matrix with cotangent weights (Sect. 2), which is symmetric.

7 Discussion and examples

We describe a simple criterion for the selection of the seed points of the spectral basis functions (Sect. 7.1); then, we discuss their main properties (Sect. 7.2) our experiments (Sect. 7.3) on 3D shapes.

7.1 Selection of the seed points

Since the diffusion basis functions at small scales and the spectral basis functions induced by filters with a strong decay to zero are compactly-supported, the coverage of the whole surface with their support depends on the selected seeds and on the support of each basis function. Increasing the number of seed points or the scale $t$, or reducing the filter decay to zero, all the vertices of the input mesh will belong to the support of one or more basis functions (Fig. 9). The seed points can be selected as high-curvature points, or by uniformly sampling the input 3D shape, or by applying a sampling or clustering method (e.g., the farthest point sampling, the principal component analysis).

In our approach, at the first step we select a seed point (e.g., a point of maximum curvature) and define $k$ (e.g., $k:=10$, in our experiments) seed points by applying the farthest point sampling method [21, 39]. Then, we compute the corresponding set $\mathcal{B}$ of basis functions and identify the points of the input surface that are not covered by the support of the basis functions in $\mathcal{B}$. These points are then clustered and the representative points of these clusters provide the seeds for the new basis functions, which will be included in $\mathcal{B}$. This process proceeds until all the vertices belong to the support of at least one basis function. At each iteration, the seed points are selected among the vertices of the input surface and are marked as visited.

7.2 Properties of the basis functions

For each of type of basis, we discuss the following properties

partition of the unity $\sum_{i}\psi_{i}(\mathbf{p})=1$ and non-negativity $\psi_{i}(\mathbf{p})\geq 0$;

• 2.

orthogonality with respect to an inner product on the space of scalar functions defined on the input shape;

• 3.

locality: the support $\textrm{supp}(\psi):=\overline{\{\mathbf{p}\in\mathcal{N}:\,\psi(\mathbf{p})\neq 0\}}$ of the basis function is compact;

• 4.

smoothness, intrinsic definition with respect to the domain metric and/or encoding of local geometric properties (e.g., curvature, geodesic distance), and invariance to isometric transformations (e.g., for shape analysis and comparison);

• 5.

sparsity of the representation: scalar functions are accurately approximated by a small number of basis functions;

• 6.

functional stability of the basis with respect to domain discretisation (e.g., connectivity, sampling, geometric/topological noise) and stability of the representation with respect to a basis under small shape perturbations;

• 7.

computational cost and storage overhead feasible for applications in geometry processing and shape analysis.

Table 1 summarises the properties discussed in Sect. 1 for each set of basis functions $(\psi_{i})_{i}$, defined with one the methods previously described. For those functions that might satisfy a given property under some additional hypothesis, indicated with the symbol $\otimes$, we notice that the partition of the unity can be always induced by normalising the basis functions as $\psi_{i}/\sum_{j}\psi_{j}$. However, this normalisation generally alters other properties, such as the general Lagrange property. The locality is satisfied by the diffusion basis function at small scales and by the spectral basis functions if the corresponding filter rapidly decays to zero.

The non-negativity is satisfied by the harmonic and diffusion bases, as a result of the maximum principle for the Laplace and heat diffusion equations. By definition, the compressed manifold modes are local and orthonormal; however, we cannot center these basis functions at a given seed point and easily control their support, as we can do with the diffusion basis functions. The orthogonality can be induced by applying the Gram-Schmidt method, at the cost of altering other properties of the input functions. Two diffusion basis functions without overlapping supports are always orthogonal; the orthogonality of a set of diffusion basis functions with overlapping supports can be achieved with the Gram-Schmidt method at the cost of a (generally) larger support of the orthonormal diffusion basis functions, as a matter of the linear combination of functions with a different support size during the orthonormalisation iterations. The numerical stability is generally satisfied by all those basis functions that are computed through the solution of a linear system (e.g., harmonic, Hamiltonian, spectral basis functions, etc.). Numerical instabilities are generally associated with the computation of the Laplacian and Hamiltonian eigenbasis associated with multiple or numerically close eigenvalues.

7.3 Examples

For an arbitrary 3D shape, we cannot compute the ground-truth diffusion basis functions at a given seed point through the analytic representation of its Laplacian eigenfunctions (e.g., as for the sphere and the cylinder). Recalling the upper bound to the accuracy of the Padé-Chebyshev approximation of the diffusion basis with respect to the selected degree of the rational polynomial, we can analyse the discrepancy of this approximation with respect to the truncated spectral approximation. In Fig. 10, we report the $\ell_{\infty}$ error ($y$-axis) between the Padé-Chebyshev approximation of the diffusion basis at different seed points and scales, and the truncated spectral approximation with a different number ($x$-axis) of Laplacian eigenpairs. At large scales (e.g., $t=10^{-1},1$), the truncated spectral approximation with a low number of eigenpairs (e.g., $k\approx 50$) provides an approximation of the solution to the heat equation that is close to the one computed by the Padé-Chebyshev approximation, with respect to the $\ell_{\infty}$-norm ($\epsilon_{\infty}=10^{-4}$).

The level-sets of these two approximations also show their different behaviour at small scales, while at larger scales it becomes analogous in terms of both the shape and distribution of the level-sets. At small scales (e.g., $t=10^{-2},10^{-3}$), the truncated spectral approximation is generally affected by small undulations far from the seed point and that do not disappear while increasing the number of Laplacian eigenpairs. Increasing the number of Laplacian eigenpairs makes these undulation more evident, as a matter of the small vibrations of the Laplacian eigenvectors associated with larger Laplacian eigenvalues, which are included in the truncated approximation. In this case, a larger number of eigenpairs is necessary to obtain a truncated spectral approximation close to the one computed with the Padé-Chebyshev rational polynomial.

The truncated spectral approximation of the diffusion basis is generally affected by the Gibbs phenomenon. This phenomenon is more evident at small cases, as the kernel values decrease fast to zero and are largely affected by small negative values (Fig. 11(e,f)). In fact, at small scales the values of the diffusion basis functions (or equivalently, the diffusion kernel) decrease fast to zero and the negative values are no more compensated by the Laplacian eigenvectors related to smaller eigenvalues, as they are not included in the approximation. For the Padé-Chebyshev approximation (Fig. 11(a-d)), the kernel values are positive at all the scales; in fact, we approximate the filter function without selecting a sub-part of the Laplacian spectrum.

We have also tested the robustness of the spectrum-free computation with respect to a different shape discretisation (Fig. 12), bordered surfaces (Fig. 13), geometric (Fig. 14) or topological (Fig. 15, Fig. 16) noise, and almost isometric deformations (Fig. 17). A higher resolution (Fig. 12) of the input surface improves the quality of the level-sets, which are always uniformly distributed, and an increase of the noise magnitude (Fig. 14) does not affect the shape and distribution of the level-sets. In case of topological noise (Fig. 16), the shape of the level-sets of the diffusion basis functions at large scales (e.g., $t=1$) changes only in a neighbourhood of the topological cut. At small/medium scales (e.g., $t=0.01$, $t=0.1$), the geometry of the surface around the cut has no influence on the diffusion basis functions, as a matter of their local supports. In our experiments, the analogous behaviour of the level-sets of the diffusion basis functions confirms the robustness of the Padé-Chebyshev of the approximation with respect to sampling and discretisation.

From the previous experiments, we conclude that the diffusion basis functions and the Laplacian spectral functions provide a valid alternative to the harmonic functions and Laplacian eigenfunctions. In fact, the diffusion basis functions can be centred at any point of the input domain, encode local/global shape details according to the values of the temporal parameter, and encode geometric information in terms of curvature values and geometric details.

8 Conclusions and future work

Representing a signal as a linear combination of a set of basis functions is used in a wide range of applications, such as approximation, de-noising, and compression. In this context, we have focused our attention on the main aspects of signal approximation, such as the definition, computation, and comparison of basis functions induced by the Laplace-Beltrami operator. In particular, we have introduced the diffusion and Laplacian spectral basis functions, which encode intrinsic shape properties, are multi-scale and sparse, efficiently computed, and provide an alternative to the harmonic functions and to the Laplacian eigenfunctions. As future work, we mention the design of filters that support user-driven constraints on the resulting spectral basis functions and the use of more general classes of differential operators, such as the Dirac operator [37] and elliptic operators.

Acknowledgements

We thank the Reviewers for their constructive comments, which helped us to improve the paper presentation and content, and the members of the Shapes$\&$Semantics Modelling Group at CNR-IMATI for helpful discussions. This work has been partially supported by the H2020 ERC Advanced Grant CHANGE (Project Agreement N. 694515). Shapes are courtesy of the AIM@SHAPE Repository, the TOSCA, SHREC 2010, and SHREC 2016 data sets.