LLaplacian Spectral Basis Functions
G. PatanèCNR-IMATI, Italy Abstract
Representing a signal as a linear combination of a set of basis functions is central ina wide range of applications, such as approximation, de-noising, compression, shapecorrespondence and comparison. In this context, our paper addresses the main as-pects of signal approximation, such as the definition, computation, and comparison ofbasis functions on arbitrary 3D shapes. Focusing on the class of basis functions in-duced by the Laplace-Beltrami operator and its spectrum, we introduce the diffusionand Laplacian spectral basis functions, which are then compared with the harmonicand Laplacian eigenfunctions. As main properties of these basis functions, which arecommonly used for numerical geometry processing and shape analysis, we discuss thepartition of the unity and non-negativity; the intrinsic definition and invariance with re-spect to shape transformations (e.g., translation, rotation, uniform scaling); the locality,smoothness, and orthogonality; the numerical stability with respect to the domain dis-cretisation; the computational cost and storage overhead. Finally, we consider geomet-ric metrics, such as the area, conformal, and kernel-based norms, for the comparisonand characterisation of the main properties of the Laplacian basis functions.
1. Introduction
In Computer Graphics, scalar functions are ubiquitous to represent the values of aphysical phenomenon (e.g., the heat that diffuses from one or more source points), a ∗ G. Patané, Consiglio Nazionale delle Ricerche, Istituto di Matematica Applicata e Tecnologie Infor-matiche, Via De Marini 6, 16149 Genova, Italy
Email address: [email protected] (G. PatanèCNR-IMATI, Italy)
Preprint submitted to CAGD - New Progress in Geometry for Applications June 11, 2019 a r X i v : . [ c s . G R ] J un hape 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 surfacetexture.In all these cases, representing the input function in terms of a basis allows us toaddress a large number of applications, such as multi-resolution representations by se-lecting a set of multi-scale basis functions (e.g., Laplacian eigenfunctions, diffusionbasis functions), sparse representations by choosing a low number of basis functionsin order to achieve a target approximation accuracy, or compression by quantising therepresentation coefficients. We can also address deformation by modifying the co-efficients that express the geometry of the input surface in terms of geometry-drivenor shape-intrinsic basis functions, smoothing by neglecting the coefficients associatedwith large Laplacian frequencies, and the definition of Laplacian spectral kernels anddistances as a filtered combination of the Laplacian spectral eigenfunctions.In this context, our paper focuses on the main aspects of signal approximation onarbitrary 3D shapes, such as the definition, computation, and comparison of basis func-tions for applications in numerical geometry processing and shape analysis (Fig. 1). Todefine the functional space associated with an input 3D shape, we select a set of itsbasis functions and then represent any signal as a linear combination of these func-tions. In this work, we focus on the class of bases induced by the Laplace-Beltramioperator and its spectrum; i.e., the harmonic basis functions (Sect. 3), as solution to theLaplace equation; the
Laplacian eigenfunctions (Sect. 4) associated with the spectrumof the Laplace-Beltrami operator; and the
Laplacian spectral basis functions (Sect. 5),which are defined by filtering the Laplacian spectrum and include the diffusion basisfunctions .The diffusion and Laplacian spectral functions are novel classes of functions, whichwill be defined by further developing our recent results on the Laplacian spectral dis-tances [44, 43, 42]. Then, these functions will be compared with the harmonic func-tions and the Laplacian eigenfunctions. As a main contribution with respect to theprevious work, the diffusion and Laplacian spectral basis functions are intrinsic to theinput shape and local; i.e., they have a compact support that can be tuned easily byselecting the diffusion scale and the decay of the filter function to zero, respectively.2efining basis functions with a compact support is particularly important to lo-calise their behaviour around feature points and to reduce their storage overhead in thediscrete setting. In fact, given a discrete domain M with n points, the space of func-tions on M has dimension n and its basis functions are encoded as a n × n matrix. Thismatrix 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 differentmetrics (e.g., area, conformal, kernel-based metrics), which measure differential prop-erties 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 intrin-sic definition and invariance with respect to shape transformations (e.g., translation,rotation, uniform scaling); the locality, smoothness, and orthogonality; the numericalstability with respect to the domain discretisation; the computational cost and storageoverhead.Our experiments (Sect. 7, 8) show that the diffusion and the Laplacian spectralfunctions provide a valid alternative to the harmonic functions and Laplacian eigen-functions. In fact, the diffusion basis functions can be centred at any point of the inputdomain, have a compact support, and encode local/global shape details according to thevalues of the temporal parameter or to the decay of the selected filter to zero. Similarproperties also apply to the Laplacian spectral basis functions induced by a filter witha strong (e.g., exponential) decay to zero.
2. Laplace-Beltrami operator
Let N be a smooth surface, possibly with boundary, equipped with a Riemannianmetrics and let us consider the inner product (cid:104) f , g (cid:105) : = (cid:82) N f ( p ) g ( p ) d p defined on thespace L ( N ) of square integrable functions on N and the corresponding norm (cid:107) · (cid:107) .The Laplace-Beltrami operator ∆ : = − div ( grad ) is self-adjoint (cid:104) ∆ f , g (cid:105) = (cid:104) f , ∆ g (cid:105) ,and positive semi-definite ; i.e., (cid:104) ∆ f , f (cid:105) ≥ ∀ f , g , [48]. Finally, the value ∆ f ( p ) doesnot depend on f ( q ) , for any couple of distinct points p , q ( locality ).3 igure 1: Overview of the proposed approach for the definition of Laplacian spectral basis functions basedon the solution to the harmonic equation, the Laplacian eigenproblem, and the diffusion equation, and on thefiltering of the Laplacian spectrum. According to [43], we represent the Laplace-Beltrami operator on surface and vol-ume meshes in a unified way as ˜ L : = B − L , where the mass matrix B is sparse, sym-metric, positive definite, and the stiffness matrix L is sparse, symmetric, positive semi-definite, and L1 = . Analogously to the continuous case, the Laplacian matrix satisfiesthe following properties: • B - adjointness : ˜ L is adjoint with respect to the B -inner product (cid:104) f , g (cid:105) B : = f (cid:62) Bg ;i.e., (cid:104) ˜ Lf , g (cid:105) B = (cid:104) f , ˜ Lg (cid:105) B = f (cid:62) Lg . If B : = I , then this property reduces to the sym-metry of L ; • positive semi-definiteness : (cid:104) ˜ Lf , f (cid:105) B = f (cid:62) Lf ≥
0. In particular, the Laplacian eigen-values are positive, with one null eigenvalue associated with a constant eigenvec-tor; • locality : assuming that B is diagonal, the value ( ˜ Lf ) i depends only on the f -values at p i and its 1-star neighbourhood.On triangle meshes , for the stiffness matrix we can consider the linear FEM Laplacian4eights [47, 54] that reduce to the Voronoi-cotangent weights [18], by lumping theFEM 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 alwayspositive. In [17], the weak formulation of the Laplacian eigenproblem is achieved byselecting a set of volumetric test functions, which are defined as k × k × k B-splines(e.g., k : =
4) and restricted to the input shape. For the anisotropic Laplacian [5], theentries of 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 Lapla-cian matrix with cotangent weights to surface meshes with non-planar, non-convexfaces. An approximation of the Laplace-Beltrami operator with point-wise conver-gence has been proposed in [11].In [8, 9, 10, 12], the Laplace-Beltrami operator on a point set has been defined asthe Gram matrix associated with the exponential kernel. Starting from this approach,a new discretisation [38] has been achieved through a finer approximation of the localgeometry of the surface at each point through its Voronoi cell. The discretisation ofthe Laplace-Beltrami operator on volumes represented as a tetrahedral mesh has beenaddressed in [4, 36, 53].Finally, in the paper examples we apply the linear FEM weights and the level-setsof a given function are associated with iso-values uniformly sampled in its range.
3. Harmonic basis
Given a set S : = { p i } i ∈I of seed points (Sect. 7.1), we compute the harmonic basisfunctions B : = { ψ i } i ∈I , as solution to the harmonic equation ∆ψ i = ψ i ( p j ) = δ i j , i , j ∈ I . Each harmonic function, and in particular any harmonic basis function, • minimizes the Dirichlet energy E ( u ) : = (cid:82) N (cid:107)∇ u ( p ) (cid:107) d p ; • satisfies the locality property ; i.e., if p and q are two distinct points, then ∆ u ( p ) is not affected by the value of u at q ;5 igure 2: Distribution and shape of the level-sets of harmonic basis functions, centred at seed points. • verifies the mean-value theoremu ( p ) = ( π R ) − (cid:90) Γ u ( s ) d s = ( π R ) − (cid:90) B u ( q ) d q , where B ⊆ N is a disc of center p , radius R , and boundary Γ ( mean-value theo-rem ); • satisfies the strong maximum principle [48]. More precisely, if Ω is a connectedopen set, u ∈ C ( Ω ) , u is harmonic and attains either a global minimum or max-imum in Ω , then u is constant.The maximum principle states that the values of a harmonic function in a compactdomain are bounded by its maximum and minimum values on the boundary. In partic-ular, the values of the harmonic basis function ψ i belong to the interval [ , ] and itsmaximum is attained at the selected seed point p i .The discrete harmonic basis associated with a seed point p i is the solution to thelinear system L (cid:63) g = e i , where L (cid:63) is achieved by replacing the i -th row of the stiff-ness matrix with the vector e i , where e ( j ) i : = δ i j (Fig. 2). Harmonic functions are ef-ficiently computed in O ( n ) time with iterative solvers of sparse linear systems; theircomputation is stable for the mean-value weights while negative Voronoi cotangentweights might induce local undulations in the resulting harmonic function. Finally, the(globally-supported) geometry-aware basis [52] has been defined by minimising theDirichlet energy, with interpolating and least-squares constraints.6 re-conditioners and fast solvers of elliptic problems. In the context of the iso-geometricanalysis of elliptic PDEs, previous work has proposed (i) additive multilevel precondi-tioners [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 ofComputer Graphics, specialised pre-conditioners for the Laplacian matrix [33] can beapplied to further attenuate numerical instabilities.As solvers, we can select (i) robust and efficient multi-grid methods for single-patchiso-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., boundedindependently of the matrix size). As alternatives, we mention the fast solvers for largelinear systems arising from the Galerkin approximation based on B-splines [19], whichare linear with respect to the matrix size and have a convergence speed that is indepen-dent 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, condition-ing, 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 B : = { ( λ n , φ n ) } + ∞ n = , ∆φ n = λ n φ n , with0 = λ ≤ λ n ≤ λ n + , in L ( N ) (Sect. 2). Any function f in L ( N ) satisfies the fol-lowing relations f = ∑ + ∞ n = α n φ n , α n : = (cid:104) f , φ n (cid:105) ; (cid:107) f (cid:107) = ∑ + ∞ n = α n , (cid:107)∇ f (cid:107) = ∑ + ∞ n = α n λ n . The Laplacian eigenfunctions are intrinsic to the input shape and those ones relatedto smaller eigenvalues correspond to smooth and slowly-varying functions. Increasingthe eigenvalues, the corresponding eigenfunctions generally show rapid oscillations.If two shapes are isometric, then they have the same Laplacian spectrum ( iso-spectral φ φ Figure 3: Level-sets of three Laplacian eigenfunctions. property ); however, the viceversa does not hold [27, 56] and we cannot recover themetric of a given surface.
Hamiltonian basis and eigenbasis. A Hamiltonian operator H : L ( R ) → R is definedas H f = − ∆ f + µV f , where V : N → R is a potential function and µ is a real parame-ter that controls the influence of the potential on the Hamiltonian operator and definesthe trade-off between the local and global support of the corresponding eigenfunctions.Analogously to the harmonic basis and Laplacian eigenbasis, we can define the Hamil-tonian basis , as the solution to the problem H ψ i = ψ i ( p j ) = δ i j , and the Hamiltonianeigenbasis H ψ n = σ n ψ n , n ∈ N . Optimality of the Laplacian eigenbasis.
According to [2], the Laplacian eigenfunc-tions represent an optimal basis for the representation of signals with bounded gradientmagnitude. In fact, • the spectral decomposition f = ∑ ni = (cid:104) f , φ i (cid:105) φ i is optimal in approximating func-tions with L bounded gradient magnitude; i.e., the residual error r n : = f − f n isbounded as (cid:107) r n (cid:107) ≤ (cid:107)∇ f (cid:107) λ n + ; (1) • the spectral decomposition is optimal in approximating functions with respect tothe error estimates in (1). In fact, for any 0 ≤ α < n and no8a)(b)(c) Figure 4: (b,c) Level-sets of the diffusion basis functions at different scales and centered at different seedpoints in (a). The example confirms the locality, smoothness, and shape-awareness of the diffusion basisfunctions, computed with the spectrum-free method (Sect. 5.2). sequence ( ψ i ) ni = of linearly independent functions in L such that (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) f − n ∑ i = (cid:104) f , ψ i (cid:105) ψ i (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ≤ α (cid:107)∇ f (cid:107) λ n + , ∀ f . Discrete Laplacian eigenbasis.
The generalised Laplacian eigensystem ( λ i , x i ) ni = , with0 = λ < λ ≤ · · · ≤ λ n , satisfies the identity Lx i = λ i Bx i and the eigenvectors are or-thonormal with respect to the B -inner product; i.e., (cid:104) x i , x j (cid:105) B = x (cid:62) i Bx j = δ i j . In particu-lar, ˜ L = X Γ X (cid:62) B is the spectral decomposition of the Laplacian matrix, where X is theeigenvectors’ matrix and Γ is the diagonal matrix of the eigenvalues. Fig. 3 shows the9 = . t = . t = . Figure 5: Level-sets of diffusion basis functions at three scales and centred at the same seed point. level-sets of three eigenfunctions. In [7, 41], Laplacian eigenfunctions with a compactsupport ( compressed manifold modes ) have been defined by solving an orthogonality-constrained optimisation problem with a L penalty term. In [32], the compact sup-port of the Laplacian eigenfunctions is enforced with a L p , 0 < p <
1, penalisationterm; then, these eigenfunctions are computed by combining splitting strategies withan ADMM-based (
Alternating Direction Method of Multipliers ) iterative algorithm. Fi-nally, the optimality of the Laplacian eigenfunctions for spectral geometry processingand for signal approximation has been addressed in [13] and in [1, 2], respectively.
Computation of the Laplacian eigenbasis.
For the computation of the Laplacian eigen-functions, numerical methods generally exploit the sparsity of the Laplacian matrix andreduce the high-dimensional eigenproblem to one of lower dimension, by applying acoarsening step. Main examples include the algebraic multi-grid method [22], Arnoldiiterations [35, 51], and the Nystrom method [24]. Even though the eigenvalues andeigenvectors are computed in super-linear time [54], this computational cost and therequired O ( n ) storage are expensive for densely sampled domains. Indeed, modifi-cations of the Laplacian eigenproblem are applied to evaluate specific sub-parts of theLaplacian spectrum (e.g., shift, power, and inverse methods). Finally, the Laplacian10igenvectors are computed only for a small set of eigenvalues and do not provide aflexible alignment of the function behaviour to specific shape features. Stability of the Laplacian eigenbasis.
Theoretical results on the sensitivity of the Lapla-cian spectrum against geometry changes, irregular sampling density and connectivityhave been presented in [29, 55]. Here, we briefly recall that the instability of the com-putation 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]. Whilerepeated eigenvalues are quite rare and typically associated with symmetric shapes, nu-merically close or switched eigenvalues can be present in the spectrum, in spite of theregularity of the input discrete surface.
5. Laplacian spectral basis
We address the definition (Sect. 5.1) and computation (Sect. 5.2) of the Laplacianspectral basis by filtering the Laplacian spectrum. Then, we study the relations betweenthe Laplacian spectral basis and the Green kernel (Sect. 5.3).
Given a positive filter ϕ : R + → R + , let us consider the power series ϕ ( s ) = ∑ + ∞ n = α n s n .Noting that ∆ i f = ∑ + ∞ n = λ in (cid:104) f , φ n (cid:105) φ n and under a few hypotheses on the decay of thefilter to zero [44], we define the spectral operator as Φ ( f ) = + ∞ ∑ n = α n ∆ n f = + ∞ ∑ n = ϕ ( λ n ) (cid:104) f , φ n (cid:105) φ n , (2)which is linear, continuous, and Φ ( f ) = (cid:104) K ϕ , f (cid:105) , where K ϕ ( p , q ) = + ∞ ∑ n = ϕ ( λ n ) φ n ( p ) φ n ( q ) , (3)is the spectral kernel . Then, the spectral basis function centred at a point p is definedas the action of the spectral operator on δ p or equivalently as the spectral kernel centredat p ; in fact, Φ ( δ p ) = K ϕ ( p , · ) .The filters ϕ t ( s ) : = exp ( − ist ) , ϕ ( s ) : = s − k / , ϕ ( s ) : = s − / induce the wave [14,6], poly-harmonic , and commute-time basis functions , respectively. Mexican hat wavelets [31]11 runcated spectral approximation P.C. approximation (a) k =
100 (b) k = k = K (d) ε ∞ = − s = − s = − s = − Figure 6: Smoothness and locality of the diffusion basis functions at a seed point (red dot) and at differentscales, which have been computed with (a-c) the truncated spectral approximation (i.e., k Laplacian eigen-pairs, linear FEM Laplacian weights) and (d) the Padè-Chebyshev approximation. Since the input shapehas 2 K vertices, the spectral approximation (c) provides the ground-truth. are generated by the filter ϕ ( s ) : = s / exp ( − s ) and the filter function ϕ ( s ) : = exp ( is ) , s ∈ [ , π ] , defines the wave kernel. In ϕ ( s ) : = t k s α exp ( − ts α ) , the parameter k scalesthe rate of diffusion and α controls the decay of the Laplacian eigenvalues to zero,similarly to random walks [46]. Diffusion basis.
Selecting the filter ϕ t ( s ) : = exp ( − ts ) , the spectral operator reducesto the diffusion operator Φ t : = exp ( − t ∆ ) and the diffusion basis function at p (Fig. 4,Fig. 5) is the heat kernel K t ( · , p ) centred at p , which solves the heat diffusion equation ( s ) = ( + s ) − ϕ ( s ) = ( + s ) − ϕ ( s ) = + s ( + s ) ϕ t ( s ) = exp ( − ts ) t = . Figure 7: Level-sets of Laplacian spectral basis functions induced by different filters. ( ∂ t + ∆ ) H ( · , t ) = H ( · , ) = δ p ; i.e., K t ( p , · ) = + ∞ ∑ n = exp ( − λ n t ) φ n ( p ) φ n . (4)The spectral representation (4) shows the smoothing effect on the initial condition; asthe scale increases, the component of δ p along the eigenfunctions associated with thelarger Laplacian eigenvalue becomes null.Assuming that Ω is an open and bounded set and that the solution F ( · , t ) to the heatequation is sufficiently smooth, we define the parabolic cylinder Ω T : = Ω × ( , T ] andthe parabolic boundary Γ T : = Ω T \ Ω T . According to the strong maximum principle for the heat diffusion equation [34] (Ch. 2), if Ω is connected and there exists a point ( p , t ) in Ω T such that H ( p , t ) = max Ω T H ( p , t ) , then H ( · , t ) is constant in Ω T . Inparticular, the diffusion basis functions are positive. Considering the generalised eigensystem LX = BX Γ , with orthonormal eigenvec-tors X (cid:62) BX = I , the discretisation of the spectral operator (2) is K ϕ = ϕ ( ˜ L ) = X ϕ ( Γ ) X (cid:62) B , ϕ ( Γ ) : = diag ( ϕ ( λ i )) ni = ; i.e., K ϕ f = n ∑ i = ϕ ( λ i ) (cid:104) f , x i (cid:105) B x i . (5)In particular, ˜ L and K ϕ have the same eigenvectors and ( ϕ ( λ i )) ni = are the (filtered)Laplacian eigenvalues of K ϕ . Assuming that K ϕ : = ϕ ( ˜ L ) is invertible (i.e., the fil-ter does not vanish at the Laplacian eigenvalues), the vectors in B : = { K ϕ e i } ni = (i.e.,13 = − t = − t = − t = Area metricsConformal metrics
Figure 8: The pixel ( i , j ) of each 100 ×
100 image indicates the area and conformal metrics of the couples ( K t e i , K t e j ) , i , j = ,..., f : = e i in Eq. (5)) are linearly independent; in fact, ∑ ni = α i K ϕ e i = is verified if andonly if K ϕ α = ; i.e., α = . The coefficients α : = ( α i ) ni = , which express f : = ∑ ni = α i K ϕ e i as a linear combination of the basis B , are computed by solving the linear system K ϕ α = f . We notice that α = K † ϕ f = X ϕ † ( Γ ) X (cid:62) Bf , where ϕ † ( Γ ) is the diagonal matrixwhose non-null values are ( / ϕ ( λ i )) ni = .Recalling that the computation of the Laplacian eigenpairs is numerically unstablein case of repeated eigenvalues [43], the filter function should be chosen in such a waythat the filtered Laplacian matrix does not have additional (if any) repeated eigenvalues.14 igure 9: Percentage ( y -axis) of vertices belonging to the support of k ( x -axis) diffusion basis functionsat three different scales t = . , . , t = t = . k = k =
26, respectively) diffusion basis functions cover the whole surface. At small scales ( t = . This condition is generally satisfied by choosing an injective filter. The selection ofperiodic filters, the expensive cost of the computation of the Laplacian spectrum, andthe sensitiveness of multiple Laplacian eigenvalues to surface discretisation motivatethe 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 wholeLaplacian spectrum and the decay of the filtered coefficients are the main reasons be-hind the approximation of a spectral basis function as a truncated sum; i.e., applyingthe relation Φ k e i = ∑ kj = ϕ ( λ j ) (cid:104) e i , x j (cid:105) B x j , where k is the selected number of eigenpairs.Even though the first k Laplacian eigenpairs are computed in super-linear time [54], theevaluation of the whole Laplacian spectrum is unfeasible for storage and computationalcost, which are quadratic in the number of surface samples. The selection of filters thatare periodic or do not decrease to zero motivates the need of defining a spectrum-free15 = − t = − t = Figure 10: (cid:96) ∞ error ( y -axis) between the Padé-Chebyshev approximation of the diffusion basis functionsat different scales and the truncated spectral approximation with a different number ( x -axis) of Laplacianeigenpairs. computation of the corresponding kernels, which cannot be accurately approximatedwith the contribution of only a subpart of the Laplacian spectrum. Furthermore, thenumber of selected eigenpairs is heuristically adapted to the decay of the filter func-tion and the approximation accuracy cannot be estimated without computing the wholespectrum. Padé-Chebyshev (spectrum-free) approximation.
Truncated spectral representations ap-ply a low pass filter and the resulting reconstructed signal generally preserves only itsglobal features. Local features, which are associated with a high frequency, are gen-erally missing in the approximated signals and can be recovered by considering onlya high number of Laplacian eigenpairs. However, the evaluation of a large part of theLaplacian 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 basisfunctions, which is based on a rational approximation of the filter. For an arbitrary16a) t = − (b) t = − (c) t = − (d) t = − (e) t = − (f) t = − Figure 11: (a-d) Robustness of the Padé-Chebyshev approximation of the diffusion basis functions and (e,f)sensitiveness of truncated spectral approximation to the Gibbs phenomenon. At all scales (a-d), the functionvalues (red curve) computed with the Padé-Chebyshev approximation are positive; at large scales (e,f), thetruncated spectral approximation is affected by the Gibbs phenomenon, as represented by the part of the plotbelow the zero line (black curve). filter, we consider the rational Padé-Chebyshev approximation p r ( s ) = a r ( s ) b r ( s ) of ϕ [26](Ch. 11) with respect to the L ∞ norm. Here, a r ( · ) and b r ( · ) are polynomials of de-gree equal to or lower than r . Let p r ( s ) = ∑ ri = α i ( + β i s ) − be the partial form ofthe Padé-Chebyshev approximation, where ( α i ) ri = are the weights and ( β i ) ri = are thenodes of the r -point Gauss-Legendre quadrature rule [26] (Ch. 11). The weights andnodes are precomputed for any degree of the rational polynomial [16]. Applying thisapproximation to the spectral kernel, we get that K ϕ e i ≈ p r ( ˜ L ) e i = r ∑ j = α j (cid:0) I + β j ˜ L (cid:1) − e i = r ∑ j = α j g j , where g j solves the symmetric and sparse linear system ( B + β j L ) g j = Be i , j = , . . . , r . (6)The Padé-Chebyshev approximation generally provides an accuracy higher than thepolynomial approximation, as a matter of its uniform convergence to the filter. In the17aper examples, we have selected r : = K ϕ f . Re-writing the solution to ( B + β j L ) g = Bf (i.e., Eq. (6) with f = e i )as g = X α , we get that α = ( I + β j Γ ) † X (cid:62) Bf and the Padè-Chebyshev approximationis K ϕ f ≈ α f + j = ,..., r ∑ i = ,..., n α j + β j λ i (cid:104) f , x i (cid:105) B x i . Setting f : = e i , we obtain that the spectral representation of the Padè-Chebyshev ap-proximation of the spectral basis function K ϕ e i centred at p i .According to [40], the approximation of the matrix ϕ ( ˜ L ) might be numerically un-stable if (cid:107) ˜ L (cid:107) is large. From the bound (cid:107) B − L (cid:107) ≤ λ − ( B ) λ max ( L ) , a well-conditionedmass matrix B guarantees that (cid:107) B − L (cid:107) is bounded. Noting that ( + β i λ j , x j ) nj = arethe eigenpairs of B + β i L , the corresponding conditioning number is bounded as κ ( B + β i L ) = max j { + β i λ j } min j { + β i λ j } = + β i λ n β i > , ( + β i λ n ) − β i < , 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 furtherattenuate numerical instabilities.The truncated spectral approximation (Fig. 6) is affected by small geometric undu-lations (especially at small scales), the use of heuristics for the selection of the numberof Laplacian eigenpairs with respect to the target approximation accuracy, and the scaleof features of the input shape. The spectrum-free computation generally provides betterresults 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 basisfunctions with a different behaviour and support (Fig. 7), thus showing the flexibilityand simplicity of our approach. For further comparison examples, we refer the readerto Sect. 7. To study the relation between the Green kernel and the basis previously defined,let us consider a linear differential operator L and the corresponding Green kernel : N × N → R such that L G ( p , q ) = δ ( p − q ) , where δ ( · ) is the delta function. Then,the solution to the differential equation L u = f is expressed in terms of the Greenkernel as u ( p ) = (cid:104) G ( p , · ) , f (cid:105) . Noting that the eigensystem ( λ n , φ n ) + ∞ n = of the inte-gral operator K G f : = (cid:104) G ( p , · ) , f (cid:105) induced by the Green kernel satisfies the relations L φ n = λ − n φ n , λ n (cid:54) =
0, we have that the eigensystem of L is ( λ − n , φ n ) + ∞ n = . Indeed, L , K G have the same eigenfunctions and reciprocal eigenvalues. In particular, the spectralrepresentation of the Green kernel G ( p , q ) = ∑ + ∞ n = λ n φ n ( p ) φ n ( q ) is uniquely definedby the spectrum of L .Combining the previous result with the fact that the spectral basis functions havebeen defined as solution to different differential equations that involve the Laplace-Beltrami operator, we get the following relations G ∆ ( p , q ) = ∑ + ∞ n = λ n φ n ( p ) φ n ( q ) harmonic oper. ( L = ∆ ) ; K t ( p , q ) = ∑ + ∞ n = exp ( − t λ n ) φ n ( p ) φ n ( q ) diffusion oper. ( L = exp ( t ∆ )) ; G ϕ ( p , q ) = ∑ + ∞ n = ϕ ( λ n ) φ n ( p ) φ n ( q ) spectral oper. ( L = ϕ † ( ∆ )) .
6. Comparison metrics for scalar functions
Since the level-sets and critical points (i.e., maxima, minima, saddles) of a scalarfunction f : N → R are independent of positive re-scalings of f , we can assume thatthe function values have been normalised in such a way that Image ( f ) = [ , ] andcompare two functions by evaluating the corresponding L ∞ or L norm. However,these norms do not measure the differential properties of the underlying functions andare better suited to evaluate the approximation accuracy. Indeed, in the following wedisucss the area, conformal, and kernel-based metrics. Area and conformal metric.
To compare two functions f , g : N → R on the same sur-face N , we consider the • area distortion : h a ( f , g ) : = (cid:82) N f ( p ) g ( p ) d p ; • conformal distortion : h c ( f , g ) : = (cid:82) N (cid:104)∇ f ( p ) , ∇ g ( p ) (cid:105) d p .The area and conformal distortions are invariant with respect to area-preserving andconformal transformations, respectively. In fact, let N , Q be 2 surfaces and F ( Q ) the19 = n = n = K Figure 12: Robustness of the Padé-Chebyshev approximation of the diffusion basis functions with respect todifferent resolutions ( n vertices). space of scalar functions defined on Q . Given a map T : N → Q , we have that [49] forall f , g ∈ F ( Q ) , h Q a ( f , g ) = h N a ( f ◦ T , g ◦ T ) if and only if T is locally area-preservingand h Q c ( f , g ) = h N c ( f ◦ T , g ◦ T ) if and only if T is conformal. For instance, for theLaplacian eigenbasis we have that h a ( φ i , φ j ) = δ i j , h c ( φ i , φ j ) = (cid:104) φ i , ∆φ j (cid:105) = λ j δ i j . The discrete area distortion is h a ( f , g ) = (cid:104) f , g (cid:105) B = f (cid:62) Bg , where B is the area-drivenmass matrix (Sect. 2), and the discrete conformal distortion is h c ( f , g ) = (cid:104) f , ∆ g (cid:105) ≈ f (cid:62) B ˜ Lg = (cid:104) f , g (cid:105) L . Expressing the input functions as a linear combination of the Laplacian eigenbasis withcoefficients α = X (cid:62) Bf , β = X (cid:62) Bg , we have that h a ( f , g ) = α (cid:62) β and h c ( f , g ) = α (cid:62) Γβ .For example, the area and conformal distortion of the delta functions at p i , p j reduceto the entry ( i . j ) of B and L , respectively. For the Laplacian eigenvectors, we have that h a ( x i , x j ) = δ i j and h c ( x i , x j ) = λ j δ i j . Fig. 8 shows the area and conformal metrics ofthe diffusion basis functions at 100 randomly sampled seed points and at 4 scales. Kernel-based metric.
Any symmetric and positive kernel K : N × N → R (e.g., theGaussian kernel) induces the following kernel-based inner product (cid:104) f , g (cid:105) : = (cid:90) N ×N K ( p , q ) f ( p ) g ( q ) d p d q . (7)20 = . t = . t = . t = . Figure 13: Stability of the Padé-Chebyshev computation of the diffusion basis functions with respect toincomplete surfaces.
Combining the spectral representation (3) of K ϕ with the definition of the kernel-basedinner product (7), we get that (cid:104) f , g (cid:105) = + ∞ ∑ n = ϕ ( λ n ) (cid:104) f , φ n (cid:105) (cid:104) g , φ n (cid:105) . The discretisation of (7) is (cid:104) f , g (cid:105) = f (cid:62) BKg = (cid:104) f , Kg (cid:105) B , where K : = ( K ( p i , p j )) ki , j = isthe Gram matrix associated with the kernel K ( · , · ) . To guarantee the symmetry of thisinner product, it is enough to assume that K is positive definite and B -adjoint; i.e., (cid:104) f , Kg (cid:105) B = (cid:104) Kf , g (cid:105) B . As main examples of kernels, we mention the filtered spectralkernel (e.g., the diffusion kernel), which is B -adjoint, and the Laplacian matrix withcotangent weights (Sect. 2), which is symmetric.
7. Discussion and examples
We describe a simple criterion for the selection of the seed points of the spectralbasis functions (Sect. 7.1); then, we discuss their main properties (Sect. 7.2) our exper-iments (Sect. 7.3) on 3D shapes.
Since the diffusion basis functions at small scales and the spectral basis functionsinduced by filters with a strong decay to zero are compactly-supported, the coverage ofthe whole surface with their support depends on the selected seeds and on the supportof each basis function. Increasing the number of seed points or the scale t , or reducing21a) (b)(c) (d) Figure 14: Robustness of the computation of the diffusion basis functions centred at a seed point placed onthe head. Level-sets of the diffusion basis functions on (b-d) the noisy shapes (bottom part, red surfaces)have been plotted on (a) the initial shape. the filter decay to zero, all the vertices of the input mesh will belong to the support ofone or more basis functions (Fig. 9). The seed points can be selected as high-curvaturepoints, or by uniformly sampling the input 3D shape, or by applying a sampling orclustering 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 maximumcurvature) and define k (e.g., k : =
10, in our experiments) seed points by applyingthe farthest point sampling method [21, 39]. Then, we compute the correspondingset B of basis functions and identify the points of the input surface that are not coveredby the support of the basis functions in B . These points are then clustered and therepresentative points of these clusters provide the seeds for the new basis functions,which will be included in B . This process proceeds until all the vertices belong to thesupport of at least one basis function. At each iteration, the seed points are selectedamong the vertices of the input surface and are marked as visited.22 able 1: Overview of the main properties (Sect. 1) of the basis functions discussed in this paper. Thesymbols • , ◦ , ⊗ indicate that the corresponding property is satisfied, is not satisfied, or might be satisfiedaccording 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 ⊗ ⊗ ⊗ ⊗ ⊗ ⊗ ⊗
Non-negativity • ◦ ◦ ◦ ◦ ◦ •
Intrinsic def. • • • • • • •
Locality ◦ ◦ ◦ • ◦ ⊗ ⊗
Orthogonality ⊗ ⊗ • • • ⊗ ⊗
Isometry-inv. • ◦ • • ◦ • •
Numer. stability • • ⊗ ⊗ • • •
Comput. cost O ( · ) n n kn log n kn log n kn log n rn log n rn log n Storage overhead O ( · ) n n kn kn kn n n For each of type of basis, we discuss the following properties • partition of the unity ∑ i ψ i ( p ) = non-negativity ψ i ( p ) ≥ • orthogonality with respect to an inner product on the space of scalar functionsdefined on the input shape; • locality : the support supp ( ψ ) : = { p ∈ N : ψ ( p ) (cid:54) = } of the basis function iscompact; • smoothness , intrinsic definition with respect to the domain metric and/or encod-ing of local geometric properties (e.g., curvature, geodesic distance), and invari-ance to isometric transformations (e.g., for shape analysis and comparison); • sparsity of the representation : scalar functions are accurately approximated by asmall number of basis functions; • functional stability of the basis with respect to domain discretisation (e.g., con-nectivity, sampling, geometric/topological noise) and stability of the representa-tion with respect to a basis under small shape perturbations; • computational cost and storage overhead feasible for applications in geometryprocessing and shape analysis. 23able 1 summarises the properties discussed in Sect. 1 for each set of basis functions ( ψ i ) i , defined with one the methods previously described. For those functions thatmight satisfy a given property under some additional hypothesis, indicated with thesymbol ⊗ , we notice that the partition of the unity can be always induced by normal-ising the basis functions as ψ i / ∑ j ψ j . However, this normalisation generally altersother properties, such as the general Lagrange property. The locality is satisfied bythe diffusion basis function at small scales and by the spectral basis functions if thecorresponding filter rapidly decays to zero.The non-negativity is satisfied by the harmonic and diffusion bases, as a result ofthe maximum principle for the Laplace and heat diffusion equations. By definition,the compressed manifold modes are local and orthonormal; however, we cannot centerthese basis functions at a given seed point and easily control their support, as we can dowith the diffusion basis functions. The orthogonality can be induced by applying theGram-Schmidt method, at the cost of altering other properties of the input functions.Two diffusion basis functions without overlapping supports are always orthogonal; theorthogonality of a set of diffusion basis functions with overlapping supports can beachieved with the Gram-Schmidt method at the cost of a (generally) larger supportof the orthonormal diffusion basis functions, as a matter of the linear combination offunctions with a different support size during the orthonormalisation iterations. The numerical stability is generally satisfied by all those basis functions that are computedthrough the solution of a linear system (e.g., harmonic, Hamiltonian, spectral basisfunctions, etc.). Numerical instabilities are generally associated with the computationof the Laplacian and Hamiltonian eigenbasis associated with multiple or numericallyclose eigenvalues. For an arbitrary 3D shape, we cannot compute the ground-truth diffusion basisfunctions at a given seed point through the analytic representation of its Laplacianeigenfunctions (e.g., as for the sphere and the cylinder). Recalling the upper bound tothe accuracy of the Padé-Chebyshev approximation of the diffusion basis with respectto the selected degree of the rational polynomial, we can analyse the discrepancy of this24 = − t = − t = t = − t = − Figure 15: Robustness of the Padé-Chebyshev computation of the diffusion basis functions at different scaleswith respect to (fist, second row) different poses. (third row) Zoom-in. approximation with respect to the truncated spectral approximation. In Fig. 10, we re-port the (cid:96) ∞ error ( y -axis) between the Padé-Chebyshev approximation of the diffusionbasis at different seed points and scales, and the truncated spectral approximation witha different number ( x -axis) of Laplacian eigenpairs. At large scales (e.g., t = − , k ≈ (cid:96) ∞ -norm( ε ∞ = − ).The level-sets of these two approximations also show their different behaviour atsmall scales, while at larger scales it becomes analogous in terms of both the shapeand distribution of the level-sets. At small scales (e.g., t = − , − ), the truncated25a) (b) (c) t = .
01 (d) t = . t = t = .
01 (g) t = . t = Figure 16: Input 3D shape with (a) one topological handle, (b) sliced surface around the thumb, and zoom-in.At small/medium scales (e.g., t = . t = .
1) the behaviour of the diffusion basis functions on the (c, d)input and (f,g) sliced surface is consistent in terms of shape and distribution of the level-sets. At large scales(e.g., t = spectral approximation is generally affected by small undulations far from the seedpoint 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 largerLaplacian eigenvalues, which are included in the truncated approximation. In this case,a larger number of eigenpairs is necessary to obtain a truncated spectral approximationclose to the one computed with the Padé-Chebyshev rational polynomial.The truncated spectral approximation of the diffusion basis is generally affectedby the Gibbs phenomenon. This phenomenon is more evident at small cases, as thekernel values decrease fast to zero and are largely affected by small negative values26 igure 17: Robustness of the computation of the linear FEM diffusion basis functions with respect to shapedeformations. (Fig. 11(e,f)). In fact, at small scales the values of the diffusion basis functions (orequivalently, the diffusion kernel) decrease fast to zero and the negative values are nomore compensated by the Laplacian eigenvectors related to smaller eigenvalues, asthey 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 approximatethe filter function without selecting a sub-part of the Laplacian spectrum.We have also tested the robustness of the spectrum-free computation with respectto 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 ofthe level-sets, which are always uniformly distributed, and an increase of the noisemagnitude (Fig. 14) does not affect the shape and distribution of the level-sets. In caseof topological noise (Fig. 16), the shape of the level-sets of the diffusion basis functionsat large scales (e.g., t =
1) changes only in a neighbourhood of the topological cut. Atsmall/medium scales (e.g., t = . t = .
8. Conclusions and future work
Representing a signal as a linear combination of a set of basis functions is used ina wide range of applications, such as approximation, de-noising, and compression. Inthis 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 theLaplace-Beltrami operator. In particular, we have introduced the diffusion and Lapla-cian spectral basis functions, which encode intrinsic shape properties, are multi-scaleand sparse, efficiently computed, and provide an alternative to the harmonic functionsand to the Laplacian eigenfunctions. As future work, we mention the design of filtersthat support user-driven constraints on the resulting spectral basis functions and the useof more general classes of differential operators, such as the Dirac operator [37] andelliptic operators.
Acknowledgements
We thank the Reviewers for their constructive comments, which helped us to improve thepaper presentation and content, and the members of the Shapes&Semantics Modelling Groupat CNR-IMATI for helpful discussions. This work has been partially supported by the H2020ERC Advanced Grant CHANGE (Project Agreement N. 694515). Shapes are courtesy of theAIM@SHAPE Repository, the TOSCA, SHREC 2010, and SHREC 2016 data sets. eferencesReferences [1] Aflalo, Y., Brezis, H., Bruckstein, A., Kimmel, R., Sochen, N., 2016. Best basesfor signal spaces. Comptes Rendus Mathematique 354 (12), 1155 – 1167.[2] Aflalo, Y., Brezis, H., R.Kimmel, 2015. On the optimality of shape and data rep-resentation in the spectral domain. SIAM Journal Imaging Sciences 8 (2), 1141–1160.[3] Alexa, M., Wardetzky, M., 2011. Discrete Laplacians on general polygonalmeshes. ACM Trans. on Graphics 30 (4).[4] Alliez, P., Cohen-Steiner, D., Yvinec, M., Desbrun, M., 2005. Variational tetra-hedral meshing. ACM Trans. on Graphics 24 (3), 617–625.[5] Andreux, M., Rodola, E., Aubry, M., Cremers, D., 2014. Anisotropic Laplace-Beltrami operators for shape analysis. In: Sixth Workshop on Non-Rigid ShapeAnalysis and Deformable Image Alignment (NORDIA).[6] Aubry, M., Schlickewei, U., Cremers, D., 2011. The wave kernel signature: aquantum mechanical approach to shape analysis. In: IEEE Computer VisionWorkshops. pp. 1626–1633.[7] Barekat, F., Caflisch, R., Osher, S., 2017. On the support of compressed modes.SIAM Journal on Mathematical Analysis 49 (4), 2573–2590.[8] Belkin, M., Niyogi, P., 2003. Laplacian eigenmaps for dimensionality reductionand data representation. Neural Computations 15 (6), 1373–1396.[9] Belkin, M., Niyogi, P., 2006. Convergence of Laplacian eigenmaps. In: NeuralInformation Processing Systems. pp. 129–136.[10] Belkin, M., Niyogi, P., 2008. Towards a theoretical foundation for Laplacian-based manifold methods. Journal of Computer System Sciences 74 (8), 1289–1308. 2911] Belkin, M., Sun, J., Wang, Y., 2008. Discrete Laplace operator on meshed sur-faces. In: Proc. of the Twenty-fourth Annual Symp. on Computational Geometry.pp. 278–287.[12] Belkin, M., Sun, J., Wang, Y., 2009. Constructing Laplace operator from pointclouds in R d . ACM Press, Ch. 112, pp. 1031–1040.[13] Ben-Chen, M., Gotsman, C., 2005. On the optimality of spectral compression ofmesh data. ACM Trans. Graph. 24 (1), 60–80.[14] Bronstein, M., Bronstein, A., 2011. Shape recognition with spectral distances.IEEE Trans. on Pattern Analysis and Machine Intelligence 33 (5), 1065 –1071.[15] Buffa, A., Harbrecht, H., Kunoth, A., Sangalli, G., 2013. BPX-preconditioningfor isogeometric analysis. Computer Methods in Applied Mechanics and Engi-neering 265, 63 – 70.[16] Carpenter, A., Ruttan, A., Varga, R., 1984. Extended numerical computations onthe “1/9” conjecture in rational approximation theory. In: Rational Approxima-tion and Interpolation. Vol. 1105 of Lecture Notes in Mathematics. Springer, pp.383–411.[17] Chuang, M., Luo, L., Brown, B. J., Rusinkiewicz, S., Kazhdan, M., 2009. Esti-mating the Laplace-Beltrami operator by restricting 3D functions. In: Proc. of theSymp. on Geometry Processing. pp. 1475–1484.[18] Desbrun, M., Meyer, M., Schröder, P., Barr, A. H., 1999. Implicit fairing of irreg-ular meshes using diffusion and curvature flow. In: ACM Siggraph. pp. 317–324.[19] Donatelli, M., Garoni, C., Manni, C., Serra-Capizzano, S., Speleers, H., 2015.Robust and optimal multi-iterative techniques for IgA Galerkin linear systems.Computer Methods in Applied Mechanics and Engineering 284, 230 – 264, iso-geometric Analysis Special Issue.[20] Donatelli, M., Garoni, C., Manni, C., Serra-Capizzano, S., Speleers, H., 2017.Symbol-based multigrid methods for Galerkin B-spline isogeometric analysis.SIAM Journal on Numerical Analysis 55 (1), 31–62.3021] Eldar, Y., Lindenbaum, M., Porat, M., Zeevi, Y. Y., Sep. 1997. The farthest pointstrategy for progressive image sampling. Trans. Img. Proc. 6 (9), 1305–1315.[22] Falgout, R. D., Nov. 2006. An introduction to algebraic multigrid. Computing inScience and Engineering 8 (6), 24–33.[23] Floater, M. S., 2003. Analysis of curve reconstruction by meshless parameteriza-tion. Numerical Algorithms (32), 87–98.[24] Fowlkes, C., Belongie, S., Chung, F., Malik, J., 2004. Spectral grouping using theNystrom method. IEEE Trans. Pattern Analysis and Machine Intelligence 26 (2),214–225.[25] Garoni, C., Manni, C., Pelosi, F., Serra-Capizzano, S., Speleers, H., 2014. On thespectrum of stiffness matrices arising from isogeometric analysis. NumerischeMathematik 127 (4), 751–799.[26] Golub, G., VanLoan, G., 1989. Matrix Computations. John Hopkins UniversityPress, 2nd Edition.[27] Gordon, C. S., Szabo, Z. I., 06 2002. Isospectral deformations of negativelycurved Riemannian manifolds with boundary which are not locally isometric.Duke Math. J. 113 (2), 355–383.[28] Herholz, P., Kyprianidis, J. E., Alexa, M., 2015. Perfect Laplacians for polygonmeshes. Computer Graphics Forum 34 (5), 211–218.[29] Hildebrandt, K., Polthier, K., Wardetzky, M., 2006. On the convergence of metricand geometric properties of polyhedral surfaces. Geometria Dedicata, 89–112.[30] Hofreither, C., Takacs, S., 2017. Robust multigrid for isogeometric analysis basedon stable splittings of spline spaces. SIAM Journal on Numerical Analysis 55 (4),2004–2024.[31] Hou, T., Qin, H., 2012. Continuous and discrete Mexican hat wavelet transformson manifolds. Graphical Models 74 (4), 221–232.3132] Huska, M., Lazzaro, D., Morigi, S., Feb 2018. Shape partitioning via l pp