Wavelet-based Heat Kernel Derivatives: Towards Informative Localized Shape Analysis
VVolume xx ( ), Number z, pp. 1–14
Wavelet-based Heat Kernel Derivatives:Towards Informative Localized Shape Analysis
M. Kirgo , , S. Melzi , , G. Patanè , E. Rodolà and M. Ovsjanikov LIX - École Polytechnique - IP Paris, EDF R&D, IMATI CNR, University of Rome La Sapienza
Abstract
In this paper, we propose a new construction for the Mexican hat wavelets on shapes with applications to partial shape matching.Our approach takes its main inspiration from the well-established methodology of diffusion wavelets. This novel constructionallows us to rapidly compute a multiscale family of Mexican hat wavelet functions, by approximating the derivative of the heatkernel. We demonstrate that this leads to a family of functions that inherit many attractive properties of the heat kernel (e.g.,local support, ability to recover isometries from a single point, efficient computation). Due to its natural ability to encode high-frequency details on a shape, the proposed method reconstructs and transfers δ -functions more accurately than the Laplace-Beltrami eigenfunction basis and other related bases. Finally, we apply our method to the challenging problems of partial andlarge-scale shape matching. An extensive comparison to the state-of-the-art shows that it is comparable in performance, whileboth simpler and much faster than competing approaches. CCS Concepts • Computing methodologies → Shape analysis; • Theory of computation → Computational geometry; • Mathematics ofcomputing → Functional analysis;
1. Introduction
In the last decade, advances in 3D shape analysis have seen theemergence of a class of methods falling under the umbrella of diffusion geometry . Based on the seminal work of Coifman andLafon [CL06], such approaches leverage the relation between thegeometry of the underlying space and the diffusion process de-fined on it, as encoded especially by the spectrum of the Laplace-Beltrami operator (LBO, for short). This general strategy hasbeen successfully exploited for the construction of point signa-tures [SOG09, GBAL09] and shape matching [OBCS ∗
12] amongother tasks. More recently, progress in this field has shifted towardsa more “local” notion of shape analysis [OLCO13, MRCB18],where descriptors are computed only on small and properly se-lected neighborhoods (Sect. 2). This choice is motivated by sev-eral relevant settings dealing with real-world 3D data, where theacquired shapes have missing subparts, due to self-occlusions, ora wildly different mesh connectivity. To date, however, combininginformative diffusion-based geometric techniques with robust lo-calized shape analysis has remained an elusive goal addressed byfew methods [OMMG10, HVG11, HQ12, OLCO13, MRCB18].In this paper, we propose an extension to the classical diffusion-based constructions by considering functions that are obtained astime derivatives of the heat kernel (Sect. 3). Such functions have lo-cal support, thus providing a natural tool for capturing multi-scaleshape properties. Furthermore, they inherit fundamental properties of the heat kernel [OMMG10], such as an efficient computation to-gether with the ability to recover isometries from a single point. Ourconstruction is also related to Mexican hat wavelets that we builddirectly on the surface while avoiding spectral approximations.From a functional standpoint, the resulting family of functionsforms an over-complete basis (a frame or, as we refer to below, a dictionary ) that provides a richer functional representation power,compared to standard LBO eigenfunctions or heat kernel functions.For example, delta functions supported at surface points are recon-structed more faithfully through our representation under a lowermemory budget (Sect. 5). This aspect has direct consequences inseveral applications, such as dense correspondence, function trans-fer across shapes, and partial shape matching (Fig. 1 and Sect. 6).Our contributions can be summarized (Sect. 7) as follows: • we introduce heat kernel derivatives as a novel tool for localizedshape analysis; • our proposed representation is compact and efficient to com-pute, while allowing an accurate representation of Mexican hatwavelets. Furthermore, it is demonstrably competitive with full-fledged algorithmic pipelines for partial shape correspondenceand similarity, at a fraction of the computational cost; • we compare our approach to popular Mexican hat wavelet for-mulations and prove that it achieves the best trade-off betweenefficiency and accuracy. submitted to COMPUTER GRAPHICS Forum (9/2020).
M. Kirgo, S. Melzi, G. Patanè, E. Rodolà and M. Ovsjanikov / Wavelet-based Heat Kernel Derivatives
Source wavelet family Target wavelet family small scale medium scale large scale small scale medium scale large scale point-to-point map conversion
Source
OursL-Beigenbasis
Figure 1:
Example functions from the proposed wavelet family on a pair of shapes, respectively a full model (left) and a partial near-isometry(center). Each function is represented at three scales (from left to right) and localized around two different samples (1 for each row). Therightmost column shows the point-to-point correspondence obtained using our wavelet family (top) and the standard Laplace-Beltrami (L-B)eigenbasis (bottom). Both maps are estimated on the same set of 13 landmarks, and visualized by color coding. Our construction is designedto be the preferred choice in the partial setting.
2. Related work
The definition of a compact and efficient representation of signalsis a fundamental task in geometry processing. By far, the mostcommon approach is to use the eigenfunctions of the Laplace-Beltrami operator, which are a natural extension of the Fourierbasis to surfaces [Tau95, Lév06, VL08]. In most settings, a trun-cated approximation consisting of the low frequency eigenfunc-tions is used to guarantee numerical robustness and computationalefficiency. The LBO eigenfunctions basis lies at the core of manyglobal and pointwise shape signatures, such as [Rus07, RWP06,SOG09,GBAL09,ASC11,MRCB16] and has been widely used forshape deformation [RCG08], segmentation [RBG ∗ ∗ et al. showed thatthe LBO eigenfunctions are the optimal for representing continu-ous functions with bounded variation, thus providing a theoreticaljustification for its versatility.The LBO eigenfunctions are also commonly used in the func-tional map framework [OBCS ∗ ∗ ∗ et al. define a compatible basis on shape collections by performinga simultaneous diagonalization of the LBO. The compressed man-ifold modes [NVT ∗
14, OLCO13, KGB16] provide a set of sparseand localised basis functions. Modifying the LBO, Choukroun etal. [CSBK17] define a Hamiltonian operator, whose eigenfunctionsare localized in those regions that correspond to the modification ofthe LBO. In [MRCB18], a similar solution is applied to define a ba-sis that is also orthogonal to a given set of functions.In [NMR ∗ et al. use polynomial combinations ofthe LBO eigenfunctions basis in conjunction with standard linearcombinations of functions to allow the representation of higher fre- quencies. In [MMM ∗ Local and multi-scale processing
More closely related to ourwork are multi-scale shape analysis methods [HPPLG11] with(i) local descriptors [Joh99, CJ97, BMP01, HSKK01, MHYS04]and (ii) diffusion geometry [SOG09, VBCG10, BK10, Pat13, PS13,Pat16]. These latter methods typically exploit the multi-scale natureof the heat kernel, which captures progressively larger neighbor-hoods of a given point while being able to characterize local geom-etry efficiently. However, the signatures based on heat diffusion canfail to capture important (e.g., medium frequency) shape details,which has led to other descriptors, such as the Wave Kernel Signa-ture [ASC11] and optimal spectral descriptors [Bro11, WVR ∗ discriminative power ofthe computed descriptors, wavelet-based techniques aim explic-itly to construct locality-aware functional families. With respectto the spectral graph wavelet signature [LH13] and the spectralgraph wavelet transform [HVG11], our approach does not rely onan eigen-decomposition and solves a small set of sparse linearsystems, which allows to capture local details and to operate oncomplex geometries. We provide an extensive comparison with themost closely related wavelet methods in [HVG11] (Sect. 3.2) andshow that our approach leads to a rich functional family that can becomputed more efficiently compared to [HVG11], while capturinglocal high frequency details, crucial for partial shape matching. Wavelets on surfaces
Finally, our work is inspired by theconstruction of wavelet-based functional families on trianglemeshes [Zho12, Ch. 4] based on subdivision [LDW97], diffu-sion [CM06], and eigendecomposition [HVG11]. While our workdoes not fit directly in this field, we base our construction on diffu-sion wavelets and specifically propose to consider the negative time submitted to COMPUTER GRAPHICS
Forum (9/2020). . Kirgo, S. Melzi, G. Patanè, E. Rodolà and M. Ovsjanikov / Wavelet-based Heat Kernel Derivatives derivative of the heat kernel to construct our multiscale functionalfamily. In the Euclidean domain, this time derivative (or equiva-lently second derivative in space) corresponds to the Mexican hatwavelet. Moreover, since it is constructed without relying on theLBO eigendecomposition [HQ12], it provides a very efficient andpowerful tool for local shape analysis. Our main goal is to construct a family of functions that is both lo-cal and provides a multi-scale description of the shape geometry ,analogously to wavelets in Euclidean domains. The most classicalapproach for generating a family of wavelet functions is via shift-ing and dilation (or scaling) of a generating function, referred toas the mother wavelet . Extending this approach to curved surfacesis challenging because shifting and dilation are not canonically de-fined on non-Euclidean domains. As a result, a large number of ap-proaches [CM06, HQ12] circumvent these challenges by replacingthese operations with those easier to mimic on surfaces.Our construction is based on the notion of diffusion wavelets ,which broadly exploit the link between diffusion and function dila-tion. As a way of motivation, consider a standard zero-mean Gaus-sian function on the real line: f ( x ) = ( σ √ π ) − exp ( − x / ( σ )) .If f is dilated and re-scaled by 1 / s , then we obtain another Gaus-sian s f ( x / s ) , whose standard deviation is multiplied by s . On theother hand, if we consider a diffusion process ∂ t f ( x , t ) = ∂ xx f ( x , t ) ,then its fundamental solution is given by the classical heat ker-nel f ( x , t ) = ( t π ) / exp ( − x / ( t )) . Recalling that the heat kernelsatisfies f ( x , σ / ) = f and noting that f ( x , s σ / ) is a Gaussianwith standard deviation s σ , we get that: f ( x , s σ / ) = s f ( x / s ) .This computation shows that in certain cases, dilation and scal-ing can be equivalently computed by solving the diffusion equa-tion starting with f . While the above computation is done with aGaussian function f , a similar result also holds for the Mexicanhat (Ricker) wavelet, which is defined as the negative second orderderivative of a Gaussian function.According to these observations, the key idea of diffusionwavelets [CM06, HQ12] is to replace dilation by diffusion , whichis particularly useful on curved surfaces. In fact, while defining di-lation is itself difficult, diffusion is well defined by replacing theLaplacian ∂ xx f ( x , t ) with the Laplace-Beltrami operator. Followingthis line of work, our main goal is to construct a multi-scale familyof functions that both have strong locality properties and exhibitgood approximation of other functions through linear combina-tions. Previous approaches have exploited these links either by us-ing a multi-scale family built directly from the heat kernel [CM06]or by operating in the spectral domain, through a truncated repre-sentation [HQ12]. Instead, we build a multi-scale functional fam-ily using the derivative (in time or, equivalently, in space) of theheat kernel and operate purely in the spatial domain by explicitlysimulating heat diffusion. This choice allows us to both avoid anexpensive eigen-decomposition necessary to approximate very lo-cal functions and to achieve better function reconstruction accu-racy, exploiting the multi-resolution properties of the Mexican hatwavelet. (a) (b)(c) (d) t max ∈ R ∗ s ∈ S n scales ∈ N ∗ Figure 2: (a-c) Illustration of our diffusion wavelets (orange) ona 1D manifold (in black), and corresponding “scaling functions"(red), which approximate the heat kernel evaluated at the samesample. (d) Parameters of our approach: largest diffusion time t max ,number of scales n scales , and samples s. The end of the support ofsuccessive wavelets (cid:110) ψ Ms , n (cid:111) n ∈ [ n scales ] is represented by blue linesand the light blue region is covered by the wavelet at t max . In the spirit of [HQ12], we define a wavelet-like family by means ofa diffusion process on a manifold M (Fig. 2(a-c), red curve). Theresulting family provides a dictionary of functions and the corre-sponding linear vector space spanned by them is used for the repre-sentation of a function in wavelet coefficients. Let u : M × R → R be the solution to the heat equation ∂ t u ( x , t ) = − ∆ u ( x , t ) , u ( x , ) = u ( x ) , t ∈ R + . (1)If the initial condition is defined at a single point (i.e. u ( x ) = δ y ( x ) with y ∈ M ), then the solution of Eq. (1) is the heat ker-nel K t ( x , y ) . K t provides a family of Gaussian-like functions onthe surface M , with increasing standard deviation (or increas-ing “scale“) as t grows. At a given scale, the negative first-orderderivative of such a function constitutes the diffusion Mexican hatwavelet. Equivalently, the Mexican hat wavelets ψ t ( x , y ) at scale t can be computed from the heat kernel K t ( x , y ) as follows: ψ t ( x , y ) = − ∂ t K t ( x , y ) = ∆ x K t ( x , y ) , (2)where ∆ x denotes the Laplace-Beltrami operator with respect to thepoint x .Given the Laplacian eigensystem { λ k , Φ k } + ∞ k = , an exact spectralformulation of the heat kernel in the continuous setting exists andis given by: K t ( x , y ) = ∞ ∑ k = exp ( − t λ k ) Φ k ( x ) Φ k ( y ) . (3)Therefore, the associated family of wavelets is defined as: ψ t ( x , y ) = ∞ ∑ k = λ k exp ( − t λ k ) Φ k ( x ) Φ k ( y ) . (4)In [HQ12], this property is used to define Mexican hat wavelets submitted to COMPUTER GRAPHICS Forum (9/2020).
M. Kirgo, S. Melzi, G. Patanè, E. Rodolà and M. Ovsjanikov / Wavelet-based Heat Kernel Derivatives
Algorithm 1
Computation of a dictionary of Mexican hat-likefunctions for a set of samples S . Ω M is the area of M . A M and W M designate the normalized area and cotangent weight ma-trices, computed using Neumann boundary conditions. ρ ∈ (
0; 1 ] isan adjustment ratio, and || . || the L norm w.r.t. A M . Input: set of samples S , number of scales n scales , maximal diffu-sion time t max , ratio ρ Output: Ψ S (multi-scale dictionary for all s ∈ S ) t ← ρ t max n scales √ Ω M Ψ S ← {} ; ψ MS ← A †M W M δ S ; ψ MS , ← ψ MS for n ← n scales do ψ MS , n ← ( A M + tW M ) † A M ψ MS , n − ; Ψ S ← { Ψ S , ψ MS , n } end for Normalize each column c of Ψ s with || c || Normalize each column c of Ψ s by max ( c ) − min ( c ) on M as a truncated version up to N =
300 LBO eigenpairs: ψ t ( x , y ) = N ∑ k = λ k exp ( − t λ k ) Φ k ( x ) Φ k ( y ) , (5)In this work, we construct a dictionary of localized functions,based on the same intuition but avoiding the eigen-decomposition,and instead solving the diffusion equation directly. Our dictionaryshares the following properties (Fig. 2d) with the spectral Mexicanhat diffusion wavelets: • it is based on the same defining relation (through derivativesin time or space) between the heat kernel and the Mexican hatwavelets as in the Euclidean setting; • our functions are located at a set of chosen sample positions S ,and the resulting dictionary provides a multi-scale representa-tion via a maximum diffusion time t max and a chosen number ofscales n scales .
3. Proposed approach
We introduce our construction of diffusion wavelets on discretesurfaces (Sect. 3.1). We compare their accuracy to other diffusionwavelet constructions (Sect. 3.2), their conversion to a point-to-point map (Sect. 3.3), and analyze their main properties (Sect. 3.4).In Sect. 4, we perform an in-depth empirical study of these proper-ties.We assume that shapes are represented as triangle meshes in thediscrete setting. We also assume that each shape M is endowedwith the Laplace-Beltrami Operator L M = A M† W M , where A M and W M are respectively the area and cotangent weight matrices[MDSB03] and A †M is the pseudo-inverse of A M . For convenience, in the following we consider the case with onesample location at vertex s . To build a dictionary of functions atvarious scales, we make three observations.1. Given the δ -function at s , denoted δ s , in the discrete setting, onecan compute a Mexican hat wavelet by applying the LBO to δ s . 2. Given a function f , one can compute a scaled version of f byapplying the diffusion operator D t to f . Additionally, the “scal-ing factor” is controlled by the diffusion time t .3. D t can be approximated precisely and efficiently via abackward-Euler scheme.Observation 1. follows from the relation between the heat ker-nel and the Mexican hat wavelet summarized in Eq. (2) and thefact that the Laplace-Beltrami and diffusion operators commute.In other words, computing the heat kernel K t ( s , x ) and then ap-plying the LBO to obtain ψ s ( x ) = ∆ K t ( s , x ) is equivalent to com-puting D t ∆δ s . This approach leads to an analogue of the “motherwavelet” and provides the means to save computational effort, sinceit avoids applying the Laplacian to each scale of the heat kernelindependently. In practice, the mother wavelet is obtained by com-puting ψ M s = A †M W M δ s , where A M and W M are computed usingNeumann boundary conditions and the vertex coordinates of M aredivided by Ω M which is the total area of M . By applying Obser-vation 2, for increasing diffusion times t to ψ M s , we obtain a setof multi-scale Mexican hat-like functions Ψ s = (cid:110) ψ M s , n (cid:111) n ∈ [ n scales ] .Finally, Observation 3 provides an efficient way to compute 2.In practice, we found that 10 −
50 Euler-steps allow to approxi-mate D t ψ M s better than a truncated spectral formulation (Sect. 3.2).Moreover, we directly use each intermediate function ψ M s , n ob-tained at the n -th backward Euler step as a function of our dictio-nary. In other words, the number of scales n scales represents thenumber of backward-Euler steps that we use to produce the Mex-ican hat wavelet associated to a diffusion time t max . Given a func-tion f on M , applying one backward-Euler step to approximatethe effect of the discretized diffusion operator D t amounts to com-puting the quantity f di f f = ( A M + tW M ) † A M f . Therefore, giventhe ( n − ) -st wavelet at a sample s , we compute the n -th waveletas : ψ M s , n = ( A M + tW M ) † A M ψ M s , n − .In our applications, we use a linear time sampling: t = ρ t max n scales √ Ω M . If two shapes N and M are involved in the com-putation, the ratio parameter ρ is set to ρ = √ Ω N √ Ω M . ρ adjusts thediffusion scales on the two shapes so that they relate well in prac-tice. If a single shape is involved, ρ =
1. This ratio is especiallyuseful in the case of partial shape matching, where N is the partialshape and M the full shape.Storing the set of sample locations S = (cid:110) s , ..., s |S| (cid:111) in thematrix δ S (the k th column of this matrix is δ s k ) allows us tocompute the wavelets at all sample locations in parallel: in-stead of computing a single mother wavelet ψ M s , we com-pute a set of mother wavelets, stored as the columns of a ma-trix ψ MS : ψ M s ← A †M W M δ S , which are propagated via backwardEuler steps to t max . This procedure (Algorithm 1) enables us tocompute the full dictionary Ψ MS = Ψ s , . . . , Ψ s |S| efficiently. Our approach has two main competitors: the Mexican hatwavelets [HQ12] and the spectral graph wavelets [HVG11]. Wealso compare our construction to an alternative approach using thewFEM diffusion operator [Pat13], which replaces the backward submitted to COMPUTER GRAPHICS
Forum (9/2020). . Kirgo, S. Melzi, G. Patanè, E. Rodolà and M. Ovsjanikov / Wavelet-based Heat Kernel Derivatives − Number of vertices C o m pu t a ti on ti m e ( s . ) Ours [HQ12]
Figure 3:
Comparison of the scalability of eigen-decomposition-based wavelet computation method [HQ12] and our approach. Thewavelets are computed at 10 sample locations, using 25 scales.
Euler approximation scheme in the generation of the Mexican hatwavelets. Note that we do not compare our method with the workof Coifman et al. [CM06], which defines orthogonal wavelets, butnot Mexican hat wavelets, using a diffusion operator. Moreover,their construction does not allow to select a set of samples fromwhich to compute the wavelet functions, whereas we rely on thisinformation. Finally, their method is performed via a full bottom-up approach, starting from all vertices of the considered shape, tolarge-scale orthogonal wavelet functions. This process uses a costlyiterative procedure that is not well suited to our applications that in-volve dense meshes.To illustrate the scalability of our approach compared to meth-ods leveraging an eigen-decomposition of the Laplace-Beltrami op-erator, we measure in Fig. 3 the time required to compute a dic-tionary at 10 sample locations and 25 scales for [HQ12] and ourmethod. We use 8 shapes from the SHREC’19 data set (connectiv-ity track) [MMR ∗ . L error to the ground truth Mexican hat wavelets(Fig. 4 (left), Table 1), (ii) L ∞ error to the ground truth Mexicanhat wavelets (Fig. 14 of the Appendix, Table 1), (iii) computationtime (Fig. 4 (right), Table 1). The first and second criteria provide away to assess how well the compared approaches approximate theground truth Mexican hat wavelet functions, while the third crite-rion measures the computational efficiency of the approaches. Wecompute the ground truth Mexican hat wavelet family in Eq. (5)with the complete Laplacian spectrum, and the intermediate dif-fusion times t GTn = log ( nt ) , where t is introduced in Algorithm 1.The evaluations are performed on all 100 shapes of the FAUST dataset (remeshed to shapes with 5K vertices), using the average errorat 10 sample locations, picked using farthest point sampling withrandom initialization.As shown in Figures 4 (as well as Fig. 3 and Table 1 and Fig. 14of the Appendix), our approach produces the most accurate approx-imation with respect to the ground truth and moreover is more com-putationally efficient than the competitors. Note that our approxi- -3 -6 -9 % t max L E rr o r
20 40 60 80 10005101520 T i m e ( s ) Ours [Pat13] [HVG11] [HQ12]
Figure 4:
Comparison between various Mexican hat diffusionwavelet definitions and our approach as L error to ground truthwavelets (left) and computation time (right), on the complete setof shapes of the FAUST data set (remeshed to shapes with 5Kvertices). See the averaged values in Table 1. Table 1:
Comparison of four Mexican hat wavelet formulations onthe FAUST data set (100 remeshed shapes with 5K vertices). The L and L ∞ norms are used as accuracy error measures. Ours [Pat13] [HQ12] [HVG11]Av. L . × − . × − . × − . × − Av. L ∞ . × − . × − . × − . × − Av. t.(s) . . × .
46 9 . mation of the ground truth diffusion wavelets is the best at smallscales in terms of both L and L ∞ norms. This is especially im-portant for partial matching where local details must be capturedfaithfully. Nevertheless, we note that [HQ12] provides the best rep-resentation of Mexican hat wavelets at large scales ( ≥
54 scales).
We use two approaches to perform point-wise signal recovery onshapes N and M , equipped with their dictionaries Ψ NS and Ψ MS .The choice of the approach is conditioned by the task we deal with. δ -function reconstruction on a single shape In the case of δ -function reconstruction on a shape M , the mapping T that asso-ciates to each vertex index of M its image provided by Ψ MS can becomputed as T = arg max rows Ψ MS ( Ψ MS ) † . (6)The left pseudo-inverse of Ψ MS , denoted ( Ψ MS ) † , is used here asthe representation in dictionary space of all δ -functions of M :the k -th column of Ψ MS † corresponds to the coordinates in thedictionary space of the δ -function located at vertex k . To con-vert back this dictionary representation in the basis of δ -functionson M , Ψ MS ( Ψ MS ) † ∈ R n M × n M is computed. The k -th col-umn of the resulting n M × n M matrix is the image of the δ -function at vertex k . Taking the argmax over the rows of this ma-trix provides the location of the δ -function according to the dic-tionary Ψ MS . Since Ψ MS is rank deficient, the computation ofits pseudo-inverse via Eq. (6) is unstable. To remedy this, we submitted to COMPUTER GRAPHICS Forum (9/2020).
M. Kirgo, S. Melzi, G. Patanè, E. Rodolà and M. Ovsjanikov / Wavelet-based Heat Kernel Derivatives use a Tikhonov (or ridge) regularization-like approach, that in-troduces the variable α in T = arg max rows Ψ MS α , which is the so-lution to argmin α || Ψ MS α − I M || + || Γα || , where Γ is an n M × ( |S| × n scales ) matrix, whose columns contain the values k ,with k ∈ [ n scales ] repeated |S| times each, and I M is the iden-tity matrix on M . Then, the solution satisfies the linear system (cid:104) Ψ MS ; Γ (cid:105) α = (cid:104) I M ; 0 n M × ( |S|× n scales ) (cid:105) , where the semi-colon nota-tion represents the column-wise concatenation of two matrices. δ -function transfer and shape matching for multiple shapes For shape matching (or δ -function transfer ) from a source shape M to a target shape N , we do not rely on the “spectrum” of the dic-tionaries. Instead, we use as a representation of a vertex k on ashape M the set of values taken by each function constituting Ψ NS .In other words, instead of using ( Ψ MS ) † as a representation in ourdictionary, we simply use ( Ψ MS ) (cid:62) . In this last representation, theembedding of the k -th vertex consists of the value taken by eachwavelet of the dictionary at the k -vertex. Note that we do not as-sume the dictionary to be an orthogonal family (which it is notin most cases). To recover the mapping T that associates to eachvertex index of M its image on N , we perform a nearest neigh-bor search T = NN-search rows (cid:16) Ψ NS , Ψ MS (cid:17) , i.e., compute for each rowof Ψ MS its nearest neighbor among the rows of Ψ NS . Our construction of the Mexican hat wavelets above inheritsmany attractive properties of the heat kernel, including isometry-invariance (due to invariance of the LBO), locality and its multi-scale nature. Moreover, as we demonstrate below, generically therelation to a single seed point through the Mexican hat wavelet ψ M t ( p , x ) is enough both to encode each point on the surface andto recover an isometry across a pair of shapes. Specifically, wecall a point p generic if it does not belong to any nodal set ofthe Laplace-Beltrami eigenfunctions, i.e. if φ i ( p ) (cid:54) = i . Asshown in [OMMG10], the set of generic points has full measure.Moreover, a surface is called generic if its Laplace-Beltrami eigen-values are non-repeating. It is well-known [BU83] that an infinites-imal perturbation to a metric of any surface makes it generic. Withthese definitions, the following theorem guarantees that the unique-ness properties of the heat kernel also apply to our wavelet familyconstruction. Theorem 1
Let M be a generic connected compact manifoldwithout boundary and p a generic point on M . For any twopoints x , y , x = y if and only if ψ M t ( p , x ) = ψ M t ( p , y ) for all t .If M and N are two generic connected compact manifolds and p ageneric point on M , then a map T : M → N where T ( p ) is genericis an isometry if and only if ψ M t ( p , x ) = ψ N t ( T ( p ) , T ( x )) for all t .The proof of Theorem 1 follows the same reasoning as the proof ofthe main theorem in [OMMG10]. For the sake of completeness, weprovide it in the Appendix (Section 8).Theorem 1 implies that generically every point x on a surfacecan be uniquely characterized by its relation to some fixed point pvia ψ t ( p , x ) . Furthermore, an isometry can be recovered given a G e o . e rr . ( × − ) Poisson diskFPS (Euclidean)FPS (geodesic)Random
Figure 5:
Mean geodesic matching error (geo. err.) on the completeset of shapes from the SHREC’16 Partial cuts data set using varioussampling strategies using scales. correspondence between a single pair of seed points, analogouslyto the heat kernel [OMMG10]. As we demonstrate below, how-ever, our wavelet-inspired construction provides a more informa-tive characterization in practice, while retaining the locality andmulti-scale nature of the heat kernel.
4. Experimental analysis
We now study different aspects of the proposed approach, such asthe sample placement, the number of scales, the computational ro-bustness, and the choice of t max . The details of all the datasets areprovided in Section 8.2 of the Appendix. Sample selection
First, we assess the effect of the sampling strat-egy on the outcome of our method. To this end, we compare far-thest point sampling (Euclidean and geodesic), Poisson disk sam-pling (computed using the gptoolbox Matlab package [J ∗ δ -function reconstruction error, from which we drawidentical conclusions. Sample robustness
Second, we verify the resilience of our ap-proach to noise in the sample placement. We consider a set of 10samples, among which we displace 1, 2, 3, 5 or all samples within ageodesic radius around the original sample location (noise radius).Six different scales are compared: 1, 2, 3, 5, 25 and 50 with noiseradii varying between 1 . × − and 1 . × − of the greatestgeodesic distance on the shape. Fig. 6 summarizes our results col-lected on the complete SHREC’16 Partial cuts data set. As a base-line, we display the matching error produced by using a dictionaryof wavelets from [HQ12] and heat kernel functions, both with 25scales. This experiment furthermore illustrates the representativepower of our approach in the case of partial shape matching. submitted to COMPUTER GRAPHICS Forum (9/2020). . Kirgo, S. Melzi, G. Patanè, E. Rodolà and M. Ovsjanikov / Wavelet-based Heat Kernel Derivatives /
10 noisy 2 /
10 noisy 3 /
10 noisy 5 /
10 noisy 10 /
10 noisy
Noise rad. ( × − ) G e o . e rr . ( × − ) Noise rad. ( × − ) G e o . e rr . ( × − ) Noise rad. ( × − ) G e o . e rr . ( × − ) Noise rad. ( × − ) G e o . e rr . ( × − ) Noise rad. ( × − ) G e o . e rr . ( × − ) Figure 6:
Geodesic matching error as a function of the noise level applied to the positioning of the samples on the complete SHREC’16Partial cuts data set. The level of noise is given as a noise radius (noise rad.), representing the geodesic disc centered around the originalsample, which the noisy sample is drawn from. The geodesic radius is expressed as a fraction of the maximum geodesic distance. In eachcolumn, an increasing number of samples are noisy (from left to right: , , , and , out of a total of samples). Table 2:
Mean geodesic error on 200 shape pairs of the FAUSTdata set and 212 pairs of the TOSCA Isometric data set (original,with edges flipped and remeshed to 5K vertices in both cases), us-ing scales and samples. Data set Mean geodesic errorFaust (original) 9 . × − Faust (remeshed) . × − Faust (edges flipped) 9 . × − TOSCA (original) . × − TOSCA (remeshed) 4 . × − TOSCA (edges flipped) 6 . × − Number of scales
Fig. 6 empirically indicates that choosing 25scales is a good trade-off between robustness to noise and computa-tional efficiency, especially when the sample position is inaccurate.
Robustness to topological changes
We verify the robustness ofour method to topological changes by comparing three versions ofthe FAUST and TOSCA Isometric data sets: (i) the original datasets, with shapes counting respectively 6890 and around 25K ver-tices, (ii) the data sets remeshed to shapes with close to 5K ver-tices and (iii) the original data sets with random edge flips appliedto 12 .
5% of the original edges. Table 2 demonstrates that our com-putation is robust to these changes leading to similar low error inall these scenarios.
Choice of t max The maximum diffusion time t max remains a freeparameter of our method. Throughout the experiments that wepresent, we choose to fix its value to 1, since it provides good re-sults on the data sets that we are using. However, selecting its valuedepending on the shape could allow to improve the quality of thematching, in particular in situations where the samples cannot beplaced regularly on the shape’s surface. To illustrate this aspect,we conducted an experiment on the humerus bones data set. Allshapes have been remeshed to count 1K vertices. Fig. 7 shows thatthe geodesic error varies substantially depending on t max . Its valueis the smallest for t max =
14. Furthermore, to highlight the repre- t max value G e o . e rr . ( × − ) Ours
Heat Kernel[HQ12]
SourceTarget
Figure 7:
Left: mean geodesic error on shape pairs of thehumerus bones data set, remeshed to 1k shapes.The minimummatching error of . × − is attained for t max = (red dot).Right: illustration of a matching estimated between two bones us-ing the best t max value. Corresponding points are depicted with thesame color. The target shape was rescaled by a factor of × . tomatch the size of the source shape. sentative power of our construction, we use the heat kernel and thediffusion wavelets of [HQ12] as a baseline. For all diffusion timeselected in this experiment, we outperform both approaches.
5. Experimental Comparisons
To illustrate the benefits of the proposed approach, we discuss self-(Sect. 5.1) and regular (Sect. 5.2) shape matching, and compare ourperformance to the heat kernel (Sect. 5.3).
One feature of our approach is that it provides a better representa-tion for δ -functions. With only a small number of sample points,we provide an approximation of δ -functions that is significantlymore accurate than traditional functional bases, such as the LBOeigenfunctions. To illustrate this aspect, we consider self-matching ,in which we evaluate the expressive power of our family in re-constructing δ functions, thereby matching the vertices of a shapeto itself. Fig. 8 presents the results obtained on all shapes of theSHREC’16 Partial cuts data set and all shapes of the TOSCA nonIsometric data set (remeshed version). The evaluation indicates submitted to COMPUTER GRAPHICS Forum (9/2020).
M. Kirgo, S. Melzi, G. Patanè, E. Rodolà and M. Ovsjanikov / Wavelet-based Heat Kernel Derivatives · − . .
15 0 . . Geodesic distance % C o rr e s pond e n ce s · − . .
15 0 . . Geodesic distance % C o rr e s pond e n ce s · − . .
15 0 . . Geodesic distance % C o rr e s pond e n ce s · − . .
15 0 . . Geodesic distance % C o rr e s pond e n ce s · − . .
15 0 . . Geodesic distance % C o rr e s pond e n ce s · − . .
15 0 . . Geodesic distance % C o rr e s pond e n ce s Ours
LBOBCMMLMH T O S CA n . I . S H R E C ’ Figure 8:
Geodesic error (self-matching) on all shapes of the SHREC’16 Partial cuts data set (top row) and on all shapes of the TOSCA non-isometric data set (category 8) (bottom row). Averaged values are reported respectively in Tables 3 and 4). Each column is with a differentnumber of sample points/non constant basis functions. Our functional dictionary recovers points on the surfaces much more accurately forthe same basis budget.
Table 3:
Average geodesic error for the SHREC’16 Partial cutsdata set (self-matching), corresponding to the top row of Fig. 8. . × − . × − . × − LMH 4 . × − . × − . × − CMM 7 . × − . × − . × − Ours 2 . × − . × − . × − the error in terms of geodesic radius, identically to the procedurein [KLF11] but using the same shape as source and target.To build our family of functions, we use 2, 4, or 6 samples,placed using Euclidean farthest point sampling, 25 scales, t max = δ -function reconstruction described in Sect. 3.3. Wecompare the performance of our dictionary to the LBO eigen-functions basis (LBOB), the Localized Manifold Harmonic ba-sis (LMH) [MRCB18] and the Compressed Manifold Modes(CMM) [NVT ∗ |S| + δ -function location is determined by taking the position atwhich its approximation in the basis considered is maximal. Ac-cording to the mean geodesic error for all approaches (Tables 3, 4),the proposed method outperforms significantly the other methods. In a more practical application, we study how well our family offunctions recovers δ -functions basis after being transferred froma source M to a target shape N via shape matching. This cor-responds to the scenario of extending a set of known seed pointcorrespondences to the entire shapes.For our family of functions, we employ the same setup as forthe δ -function reconstruction with a few adjustments: we use 3, 10 Table 4:
Average geodesic error (self-matching) for all shapesfrom the TOSCA non isometric data set (category ), correspondingto the bottom row of Fig. 8. . × − . × − . × − LMH 4 . × − . × − . × − CMM 7 . × − . × − . × − Ours 3 . × − . × − . × − or 20 samples and the transfer point-to-point conversion introducedin Sect. 3.3. The remaining parameters stay identical (25 scalesand t max = |S| sample points on source and target shapes.For global bases such as the LBO, we assume the ground truthcorrespondence between |S| first non-constant basis functions. Inthe latter setting, we follow the procedure used in [MRCB18] andleverage this known correspondence to build a ground truth func-tional map [OBCS ∗ C gt = φ N (cid:62) A N Π gt φ M , where Π gt is the ground truth point-to-point map. We then use this groundtruth functional map C gt to compute the dense point-to-point mapfollowing the the standard nearest-neighbor procedure [OCB ∗ δ -function reconstruction experiment and comparingagain against the LBOB, LMH and CMM bases. The evaluation,performed according to the standard protocol proposed in [KLF11],indicates the error in terms of geodesic radius. According to the av-erage values in Tables 5, 6, our dictionary outperforms the compet-ing bases by a substantial margin. submitted to COMPUTER GRAPHICS Forum (9/2020). . Kirgo, S. Melzi, G. Patanè, E. Rodolà and M. Ovsjanikov / Wavelet-based Heat Kernel Derivatives . . . . Geodesic distance % C o rr e s pond e n ce s . . . . Geodesic distance % C o rr e s pond e n ce s . . . . Geodesic distance % C o rr e s pond e n ce s . . . . Geodesic distance % C o rr e s pond e n ce s . . . . Geodesic distance % C o rr e s pond e n ce s . . . . Geodesic distance % C o rr e s pond e n ce s Ours
LBOBCMMLMH T O S CA n . I . S H R E C ’ Figure 9:
Matching geodesic error on all pairs of the SHREC’16 Partial cuts data set (“SHREC’16”, top row, see averaged values in Table 5)and shape pairs of the TOSCA non-isometric data set (“TOSCA n.I.”, bottom row, see averaged values in Table 6) using an increasingnumber of ground truth sample points/non constant basis functions correspondences. On all plots, the x-axis represents the normalizedgeodesic distance and the y-axis is the fraction of correspondences in percent.
Table 5:
Average geodesic error (partial shape matching) for allshape pairs of the SHREC’16 Partial cuts data set, correspondingto the bottom row of Fig. 9. . × − . × − . × − LMH 5 . × − . × − . × − CMM 7 . × − . × − . × − Ours 8 . × − . × − . × − The construction of our functions is closely related to those pro-vided by the heat kernel. While both function types share the abil-ity to characterize uniquely every point on a surface, our heat kernelderivatives are more informative in practice. To assess this practi-cal advantage, we conduct the following experiment on a set of 10pairs of the dog class from the TOSCA data set. Given an increasingnumber of samples, we compute for each pair of shapes the AUC(Area Under the Curve: the probability that a point is matched withan error less or equal to 0 .
25 in normalized geodesic distance) andthe mean geodesic error using the proposed family and the heatkernel, associated with the conversion of a point-to-point map.We setup our dictionary using the same parameters as for the δ -function transfer experiment, using ground-truth correspondencesbetween the sample points on the source and target shapes. Thequantitative and qualitative evaluation of this experiment is de-picted in Fig. 10. Relying on a diffusion process to define bothfamilies of functions, we emphasize that this experiment can beseen as an additional comparison to standard diffusion wavelets.This result highlights that heat kernel derivatives are more infor-mative than heat kernel functions. Table 6:
Average geodesic error (full shape matching) for shape pairs from the TOSCA non-isometric data set, correspondingto the bottom row of Fig. 9. . × − . × − . × − LMH 6 . × − . × − . × − CMM 7 . × − . × − . × − Ours 2 . × − . × − . × −
6. Application to Partial Shape Matching
As the main application of our method, we tackle the problem ofpartial shape matching, one of the challenging scenarios in non-rigid shape matching. We experiment on the SHREC’16 PartialCorrespondence benchmark [CRB ∗
16] (Sect. 6.2) and on a newset of partial shapes, namely FARM partial (Sect. 6.1). If not speci-fied, we adopt as sparse set of correspondences for our approachthe fully automatic result obtained with the pipeline proposedin [RBA ∗ We first evaluate our method on the FARM partial dataset. Thisdataset contains partiality and shapes undergoing non-isometric de-formations and extremely different connectivities. This makes this submitted to COMPUTER GRAPHICS
Forum (9/2020). M. Kirgo, S. Melzi, G. Patanè, E. Rodolà and M. Ovsjanikov / Wavelet-based Heat Kernel Derivatives
Ours basis
Heat Kernel basisSource
Ours
Heat Kernel1 sample 4 samples 8 samples number of samples AU C number of samples m ea ng e od e s i ce rr o r Ours
Heat Kernel
Figure 10: (1st, 2nd Rows) Comparison with heat kernels in point-to-point map conversion with the same number of scales and asmall number of point samples. Four scales of the heat diffusionfrom a sample on a pair of shapes; the colormap ranges from blue(negative) to red (positive) with values close to zero in white. (3rdRow) Qualitative comparison of the resulting maps for 1, 4, 8 sam-ples (left to right), using color correspondence to show the resultingpoint-to-point map between a source and a target shape. (4th Row)Performance of our approach compared to the heat kernel in termsof Area Under the Curve (AUC) and mean geodesic error. Resultsare averaged over 10 pairs of the dog class from the TOSCA dataset. The evaluation highlights the better performance of our repre-sentation over the heat kernel. dataset particularly challenging as many shape matching pipelinesare known to overfit to similar mesh connectivities. On the left ofFig. 11, we show a quantitative comparison on FARM partial tostate-of-the-art Partial Functional Maps [RCB ∗
17] (PFM) method.For a fair comparison, we additionally evaluate the performance ofPFM when it is initialized with the same sparse correspondencethat we exploit to generate our family of functions (PFM sparse).We also provide a qualitative illustration of the computed maps inFigure 12.Note that the state-of-the-art PFM does not perform well on thesechallenging pairs. In contrast, our method is robust, significantly simpler and more efficient, leading to a dramatic improvement inaccuracy.
We also evaluate our method In the evaluation on the SHREC’16Partial Cuts data set, where each partial shape is matched to the fullshape of the same shape category. Remark that this dataset con-tains shape pairs undergoing near-isometric deformations, whichare well-captured by the LBO basis.The quantitative evaluation is shown in Fig. 11 (middle andright). On the left, we compare our approach on the entire cuts setfrom SHREC’16 [CRB ∗
16] with all the methods that were consid-ered in the challenge. Our performance is comparable with partialfunctional maps [RCB ∗
17] (PFM) the state-of-the-art for partialmatching. The constrained optimization performed by PFM pro-duces more accurate correspondences because it is able to solve theinaccuracies contained in the initial sparse correspondence. Note,however, that due to the way the data set was produced, the shapepairs of this data set have similar connectivity, which is a knownfactor of overfitting for shape matching techniques.In Fig. 11 (right), we compare our approach to PFM when bothare initialized with 20 and 30 ground-truth correspondences only onthe cat class. As can be seen, if the sparse correspondences are cor-rect, our method is comparable to PFM and even better. We high-light that this is the only evaluation in which we use a ground-truthinitialization. According to the qualitative results in Fig. 13, ourperformance is comparable to PFM.
Computational Efficiency
When considering the computationalefficiency (in seconds), our method outperforms PFM by a sig-nificant margin. On the complete SHREC16 data set, PFM sparsetakes on average per shape pair, PFM takes , whileour method requires . Moreover, the sparse set of samplestakes on average 38.7s per shape to be computed. Therefore, mostof the computation overhead lies in this preprocessing step forour method. Once a set of sparse correspondences is available,we require an average computation time of per shape pair,which represents an improvement of 13x compared to PFM sparse( ).
7. Conclusions
In this work, we have proposed an extension to the basic diffu-sion (heat kernel) construction by considering its derivatives intime, or equivalently in space. The resulting family of diffusion-based Mexican hat wavelets is local and allows to find accuratepoint-to-point correspondences between shapes; this includes par-ticularly challenging settings such as partial shape correspondence,and matching between shapes with highly different triangulations.At the same time, the efficient use of diffusion-based methods al-lows to solve these difficult problems at a fraction of the compu-tational cost compared to other approaches. We further proved thatour functions inherit properties of the heat kernel map, such as theability of only requiring one sample point to recover an isometry.Our experiments on δ -function reconstruction and transfer indi-cate that our family can be thought of as an over-complete basis submitted to COMPUTER GRAPHICS Forum (9/2020). . Kirgo, S. Melzi, G. Patanè, E. Rodolà and M. Ovsjanikov / Wavelet-based Heat Kernel Derivatives Different connectivity Similar connectivity Similar connectivityFARM partial SHREC16 cuts SHREC16 cuts cat shapes
Geodesic error % C o rr e s pond e n ce s PFMPFM sparse
Ours
Geodesic error % C o rr e s pond e n ce s PFMPFM sparse
Ours
RFIMENGT 0 0.05 0.1 0.15 0.2 0.25020406080100
Geodesic error % C o rr e s pond e n ce s PFM 20PFM 30
Ours Ours Figure 11:
Quantitative comparison on the FARM partial data set (shapes with different connectivity) and SHREC’16 partial cut benchmark(composed of shapes with similar connectivity). In all plots, the x-axis is the mean geodesic distance to the ground truth. Abbreviationsused: PFM (partial functional maps), PFM sparse (PFM initialized with the same sparse correspondence used to compute our frame), RF,IM, EN, GT. For PFM and ours applied to SHREC’16 cuts on the cat shape, an additional number specifies the number of ground-truthcorrespondences that were used for initialization (20 or 30). that provides a richer functional representation power compared toLBO eigenfunctions, diffusion functions, and other bases. More-over, the application of wavelet-like functions on partial and large-scale shapes show promising results compared to state-of-the-artmethods, especially when taking into consideration its simplicity.The main limitation of our approach currently lies in the depen-dency on an initial sparse correspondence, which is assumed to beroughly accurate. Although further progress in deformable sparse matching would have a direct and positive impact on our method,we believe that this problem can be solved jointly within our framecalculation algorithm, and leave this challenge as an exciting direc-tion of future research.
Acknowledgments
We thank the Reviewers for their comments,which helped us to improve the technical and experimental parts ofthe article. Parts of this work were supported by the KAUST OSRAward No. CRG-2017-3426, the ERC Starting Grant No. 758800(EXPROTEA), and the ERC Starting Grant No. 802554 (SPEC-GEO).
References [ABK15] A
FLALO
Y., B
REZIS
H., K
IMMEL
R.: On the optimality ofshape and data representation in the spectral domain.
SIAM Journal onImaging Sciences 8 , 2 (2015), 1141–1160. 2[ASC11] A
UBRY
M., S
CHLICKEWEI
U., C
REMERS
D.: The wave kernelsignature: A quantum mechanical approach to shape analysis. In
Conf.on Computer Vision Workshops (2011), pp. 1626–1633. 2[BBK08] B
RONSTEIN
A. M., B
RONSTEIN
M. M., K
IMMEL
R.:
Numer-ical Geometry of Non-Rigid Shapes . Springer, New York, NY, 2008. 13,14[BK10] B
RONSTEIN
M. M., K
OKKINOS
I.: Scale-invariant heat kernelsignatures for non-rigid shape recognition. In
Conf. on Computer Visionand Pattern Recognition (2010), pp. 1704–1711. 2[BMP01] B
ELONGIE
S., M
ALIK
J., P
UZICHA
J.: Shape context: a newdescriptor for shape matching and object recognition. In
Advances inneural information processing systems (2001), pp. 831–837. 2 [BRLB14] B
OGO
F., R
OMERO
J., L
OPER
M., B
LACK
M. J.: FAUST:dataset and evaluation for 3D mesh registration. In
Conf. on ComputerVision and Pattern Recognition (2014). 13[Bro11] B
RONSTEIN
A. M.: Spectral descriptors for deformable shapes. arXiv preprint arXiv:1110.5015 (2011). 2[BU83] B
ANDO
S., U
RAKAWA
H.: Generic properties of the eigenvalueof the laplacian for compact riemannian manifolds.
Tohoku Mathemati-cal Journal, Second Series 35 , 2 (1983), 155–172. 6[CC78] C
ATMULL
E., C
LARK
J.: Recursively generated b-spline sur-faces on arbitrary topological meshes.
Computer-aided design 10 , 6(1978), 350–355. 5[CJ97] C
HUA
C. S., J
ARVIS
R.: Point signatures: A new representa-tion for 3D object recognition.
Intern. Journal of Computer Vision 25 , 1(1997), 63–85. 2[CL06] C
OIFMAN
R. R., L
AFON
S.: Diffusion maps.
Applied and Com-putational Harmonic Analysis 21 , 1 (2006), 5–30. 1[CM06] C
OIFMAN
R. R., M
AGGIONI
M.: Diffusion wavelets.
Appliedand Computational Harmonic Analysis 21 , 1 (2006), 53–94. 2, 3, 5[CRB ∗
16] C
OSMO
L., R
ODOLÀ
E., B
RONSTEIN
M., T
ORSELLO
A.,C
REMERS
D., S
AHILLIO ˘GLU
Y.: Partial matching of deformableshapes. In
Proc. of the Eurographics Workshop on 3D Object Retrieval (2016), pp. 61–67. 9, 10, 14[CSBK17] C
HOUKROUN
Y., S
HTERN
A., B
RONSTEIN
A., K
IMMEL
R.:Hamiltonian operator for spectral shape analysis. arXiv:1611.01990 (2017). 2[GBAL09] G
EBAL
K., B
ÆRENTZEN
J. A., A
NÆS
H., L
ARSEN
R.:Shape analysis using the auto diffusion function.
Computer GraphicsForum 28 , 5 (2009), 1405–1413. 1, 2[GM13] G
UNZ
P., M
ITTEROECKER
P.: Semilandmarks: a method forquantifying curves and surfaces.
Hystrix, the Italian journal of mammal-ogy 24 , 1 (2013), 103–109. 14[HPPLG11] H
EIDER
P., P
IERRE -P IERRE
A., L I R., G
RIMM
C.: Localshape descriptors, a survey and evaluation. In
Proc. Eurographics Conf.on 3D Object Retrieval (2011), pp. 49–56. 2[HQ12] H OU T., Q IN H.: Continuous and discrete mexican hat wavelettransforms on manifolds.
Graphical Models 74 , 4 (2012), 221–232. 1,3, 4, 5, 6, 7, 14[HSKK01] H
ILAGA
M., S
HINAGAWA
Y., K
OHMURA
T., K
UNII
T. L.: submitted to COMPUTER GRAPHICS
Forum (9/2020). M. Kirgo, S. Melzi, G. Patanè, E. Rodolà and M. Ovsjanikov / Wavelet-based Heat Kernel Derivatives M e s h e s O u r s PF M s p a r s e PF M Source (SMPL) Target from different data sets
Figure 12: (First row) Different mesh connectivity. (Second-fourthrows) Qualitative comparison on the FARM partial data set be-tween our approach and the PFM in is original version (PFM) andinitialized with the sparse correspondence that we adopt for thedefinition of our family of functions (PFM sparse).
Topology matching for fully automatic similarity estimation of 3Dshapes. In
Conf. on Computer Graphics and Interactive Techniques (2001), pp. 203–212. 2[HVG11] H
AMMOND
D. K., V
ANDERGHEYNST
P., G
RIBONVAL
R.:Wavelets on graphs via spectral graph theory.
Applied and Computa-tional Harmonic Analysis 30 , 2 (2011), 129–150. 1, 2, 4, 5, 14[J ∗
18] J
ACOBSON
A.,
ET AL .: gptoolbox: Geometry processing toolbox,2018. http://github.com/alecjacobson/gptoolbox . 6[Joh99] J
OHNSON
A. E.: Using spin images for efficient object recog-nition in cluttered 3D scenes.
ACM Trans. on Pattern Analysis and Ma-chine Intelligence 21 , 5 (1999), 433–449. 2[KBB ∗
13] K
OVNATSKY
A., B
RONSTEIN
M. M., B
RONSTEIN
A. M.,G
LASHOFF
K., K
IMMEL
R.: Coupled quasi-harmonic bases.
ComputerGraphics Forum 32 (2013), 439–448. 2[KGB16] K
OVNATSKY
A., G
LASHOFF
K., B
RONSTEIN
M. M.:MADMM: a generic algorithm for non-smooth optimization on mani-folds. In
Proc. of the European Conference on Computer Vision (2016),pp. 680–696. 2[KLF11] K IM V. G., L
IPMAN
Y., F
UNKHOUSER
T.: Blended IntrinsicMaps.
ACM Trans. on Graphics 30 , 4 (2011), 79. 8[LDW97] L
OUNSBERY
M., D E R OSE
T. D., W
ARREN
J.: Multiresolu-tion analysis for surfaces of arbitrary topological type.
ACM Trans. onGraphics 16 , 1 (1997), 34–73. 2[Lév06] L
ÉVY
B.: Laplace-Beltrami eigenfunctions towards an algo- Source
Ours
PFM sparse PFM
Figure 13:
Qualitative comparison on the SHREC’16 partial cutbenchmark for 5 classes (wolf, horse, centaur, dog, cat) betweenour approach, the PFM (original version), and the PFM initial-ized with the sparse correspondence that we adopt for the defini-tion of our family of functions (PFM sparse). The resulting point-to-point mapping is displayed through color correspondence. Ourapproach, despite its simplicity, is comparable to PFM. rithm that understands geometry. In
Proc. of Shape Modelling Inten-rational (2006), pp. 13–25. 2[LH13] L I C., H
AMZA
A. B.: A multiresolution descriptor for de-formable 3D shape retrieval.
The Visual Computer 29 , 6-8 (2013), 513–524. 2[LMR ∗
15] L
OPER
M., M
AHMOOD
N., R
OMERO
J., P
ONS -M OLL
G.,B
LACK
M. J.: SMPL: A skinned multi-person linear model.
ACM Trans.on Graphics 34 , 6 (2015), 248:1–248:16. 14[MDSB03] M
EYER
M., D
ESBRUN
M., S
CHRÖDER
P., B
ARR
A. H.:Discrete differential-geometry operators for triangulated 2-manifolds. In
Visualization and mathematics III . Springer, 2003, pp. 35–57. 4[Mel19] M
ELZI
S.: Sparse representation of step functions on manifolds.
Computers & Graphics 82 (2019), 117 – 128. 2[MHYS04] M
ANAY
S., H
ONG
B.-W., Y
EZZI
A. J., S
OATTO
S.: Inte-gral invariant signatures. In
European Conf. on Computer Vision (2004),Springer, pp. 87–99. 2[MMM ∗
20] M
ELZI
S., M
ARIN
R., M
USONI
P., B
ARDON
F., T
ARINI
M., C
ASTELLANI
U.: Intrinsic/extrinsic embedding for functionalremeshing of 3D shapes.
Computers & Graphics 88 (2020), 1 – 12.2[MMR ∗
19] M
ELZI
S., M
ARIN
R., R
ODOLÀ
E., C
ASTELLANI
U., R EN J., P
OULENARD
A., W
ONKA
P., O
VSJANIKOV
M.: SHREC 2019:Matching Humans with Different Connectivity. In
Eurographics Work-shop on 3D Object Retrieval (2019). 5, 14[MMRC18] M
ARIN
R., M
ELZI
S., R
ODOLÀ
E., C
ASTELLANI
U.:Farm: Functional automatic registration method for 3D human bodies,2018. 14 submitted to COMPUTER GRAPHICS
Forum (9/2020). . Kirgo, S. Melzi, G. Patanè, E. Rodolà and M. Ovsjanikov / Wavelet-based Heat Kernel Derivatives
ELZI
S., R
ODOLÀ
E., C
ASTELLANI
U., B
RONSTEIN
M.:Shape analysis with anisotropic windowed fourier transform. In
Proc.3DV (2016). 2[MRCB18] M
ELZI
S., R
ODOLÀ
E., C
ASTELLANI
U., B
RONSTEIN
M. M.: Localized manifold harmonics for spectral shape analysis.
Com-puter Graphics Forum 37 , 6 (2018), 20–34. 1, 2, 8[NMR ∗
18] N
OGNENG
D., M
ELZI
S., R
ODOLÀ
E., C
ASTELLANI
U.,B
RONSTEIN
M. M., O
VSJANIKOV
M.: Improved functional mappingsvia product preservation.
Computer Graphics Forum 37 , 2 (2018), 179–190. 2[NVT ∗
14] N
EUMANN
T., V
ARANASI
K., T
HEOBALT
C., M
AGNOR
M.,W
ACKER
M.: Compressed manifold modes for mesh processing.
Com-puter Graphics Forum 33 , 5 (2014), 35–44. 2, 8[OBCS ∗
12] O
VSJANIKOV
M., B EN -C HEN
M., S
OLOMON
J.,B
UTSCHER
A., G
UIBAS
L.: Functional maps: a flexible repre-sentation of maps between shapes.
ACM Trans. on Graphics 31 , 4(2012), 30:1–30:11. 1, 2, 8[OCB ∗
17] O
VSJANIKOV
M., C
ORMAN
E., B
RONSTEIN
M., R
ODOLÀ
E., B EN -C HEN
M., G
UIBAS
L., C
HAZAL
F., B
RONSTEIN
A.: Com-puting and processing correspondences with functional maps. In
SIG-GRAPH 2017 Courses . 2017. 8[OLCO13] O
ZOLIN , Š V., L AI R., C
AFLISCH
R., O
SHER
S.: Compressedmodes for variational problems in mathematics and physics.
Proc. of theNational Academy of Sciences 110 , 46 (2013), 18368–18373. 1, 2[OMMG10] O
VSJANIKOV
M., M
ÉRIGOT
Q., M
ÉMOLI
F., G
UIBAS
L.:One point isometric matching with the heat kernel.
Computer GraphicsForum 29 , 5 (2010), 1555–1564. 1, 6, 13[Pat13] P
ATANÉ
G.: wFEM heat kernel: Discretization and applicationsto shape analysis and retrieval.
Computer-Aided Geometric Design 30 , 3(2013), 276–295. 2, 4, 5, 14[Pat16] P
ATANÈ
G.: STAR - Laplacian spectral kernels and distances forgeometry processing and shape analysis.
Computer Graphics Forum 35 ,2 (2016), 599–624. 2[Pat18] P
ATANÈ
G.: Laplacian spectral basis functions.
Computer-AidedGeometric Design 65 (2018), 31 – 47. 2[PS13] P
ATANÉ
G., S
PAGNUOLO
M.: Heat diffusion kernel and distanceon surface meshes and point sets.
Computers & Graphics 37 , 6 (2013),676–686. 2[RBA ∗
12] R
ODOLÀ
E., B
RONSTEIN
A. M., A
LBARELLI
A., B
ERGA - MASCO
F., T
ORSELLO
A.: A game-theoretic approach to deformableshape matching. In
Conf. on Computer Vision and Pattern Recognition (2012), pp. 182–189. 9[RBG ∗
09] R
EUTER
M., B
IASOTTI
S., G
IORGI
D., P
ATANÈ
G., S
PAG - NUOLO
M.: Discrete Laplace-Beltrami operators for shape analysis andsegmentation.
Computers and Graphics 33 , 3 (2009), 381–390. 2[RCB ∗
17] R
ODOLÀ
E., C
OSMO
L., B
RONSTEIN
M., T
ORSELLO
A.,C
REMERS
D.: Partial functional correspondence.
Computer GraphicsForum 36 , 1 (2017), 222–236. 10[RCG08] R
ONG
G., C AO Y., G UO X.: Spectral mesh deformation.
TheVisual Computer 24 , 7 (2008), 787–796. 2[Rus07] R
USTAMOV
R. M.: Laplace-Beltrami eigenfunctions for defor-mation invariant shape representation. In
Proc. of Symposium on Geom-etry Processing (2007), pp. 225–233. 2[RWP06] R
EUTER
M., W
OLTER
F.-E., P
EINECKE
N.: Laplace–Beltrami spectra as ‘Shape-DNA’of surfaces and solids.
Computer-Aided Design 38 , 4 (2006), 342–366. 2[SK14] S
HTERN
A., K
IMMEL
R.: Matching the LBO eigenspace of non-rigid shapes via high order statistics.
Axioms 3 , 3 (2014), 300–319. 2[SOG09] S UN J., O
VSJANIKOV
M., G
UIBAS
L.: A concise and prov-ably informative multi-scale signature based on heat diffusion.
ComputerGraphics Forum 28 , 5 (2009), 1383–1392. 1, 2 [Tau95] T
AUBIN
G.: A signal processing approach to fair surface design.In
ACM SIGGRAPH (1995), pp. 351–358. 2[TSDS10] T
OMBARI
F., S
ALTI
S., D I S TEFANO
L.: Unique signaturesof histograms for local surface description. In
Proc. of the EuropeanConference on Computer Vision (2010), Springer, pp. 356–369. 9[VBCG10] V
AXMAN
A., B EN -C HEN
M., G
OTSMAN
C.: A multi-resolution approach to heat kernels on discrete surfaces.
ACM Trans.on Graphics 29 , 4 (2010), 121. 2[VL08] V
ALLET
B., L
ÉVY
B.: Spectral geometry processing with mani-fold harmonics.
Computer Graphics Forum 27 , 2 (2008), 251–260. 2[WVR ∗
14] W
INDHEUSER
T., V
ESTNER
M., R
ODOLÀ
E., T
RIEBEL
R.,C
REMERS
D.: Optimal intrinsic descriptors for non-rigid shape analysis.In
Proc. of the British Machine Vision Conf. (2014). 2[XZC18] X U Z., Z
HANG
Q., C
HENG
S.: Multilevel active registrationfor kinect human body scans: from low quality to high quality.
Multime-dia Systems 24 , 3 (2018), 257–270. 14[YYZ ∗
14] Y
ANG
Y., Y U Y., Z
HOU
Y., D U S., D
AVIS
J., Y
ANG
R.:Semantic parametric reshaping of human body models. In
Conf. on 3DVision (2014). 14[Zho12] Z
HONG
M.:
Harmonic Shape Analysis: From Fourier toWavelets . Master’s thesis, Stony Brook University, 2012. 2
8. Appendix8.1. Proof of Theorem 1
Both statements of the theorem follow directly from the spectralexpansion of ψ Mt ( p , x ) = ∑ i λ i exp ( − t λ i ) Φ i ( p ) Φ i ( x ) and the fol-lowing lemma, proved in [OMMG10] (Lemma 3.2, Remark 3.3). Lemma 1
Given two strictly increasing sequences λ i and µ i ofnon-negative numbers that tend to infinity, if a ( t ) = ∑ i a i exp ( − t λ i ) and b ( t ) = ∑ i b i exp ( − tµ i ) where a i , b i (cid:54) = a ( t ) = b ( t ) for all t if and only if λ i = µ i and a i = b i for all i .Applying this lemma to the spectral expansion of ψ while recallingthat only the first eigenvalue on a connected manifold is zero, weimmediately get the first statement of Theorem 1 (using the sameproof as of Theorem 3.1 in [OMMG10]). The second statementfollows from the same argument as Theorem 3.5 in [OMMG10]).Namely, by first applying this lemma to x = p , and y = q , whichimplies preservation of eigenvalues and second to other points onthe surface implying preservation of all but first eigenfunction. To-gether this implies that T preserves ψ if and only if it preserves theLaplace-Beltrami operator, which is equivalent to an isometry. The
FAUST training data set [BRLB14] consists of 100 humanshapes, with 10 different humans in 10 different poses. All shapeshave a consistent manifold mesh structure, with 6890 vertices.The
TOSCA data set [BBK08] is composed of 80 shapes of var-ious categories: 11 cats, 9 dogs, 3 wolves, 8 horses, 6 centaurs, 4gorillas, 12 female figures, and 2 different male figures, in 7 and 20poses. The mean vertex count is about 50K. If not explicitly men-tioned, the shapes of this data set are remeshed to count around 5Kvertices each. In the
TOSCA Isometric data set, we consider shapepairs within the same category (e.g., cats matched to cats), whereasin the
TOSCA non-Isometric data set, we consider matches betweenthe gorilla shapes and the two human categories (male and female). submitted to COMPUTER GRAPHICS
Forum (9/2020). M. Kirgo, S. Melzi, G. Patanè, E. Rodolà and M. Ovsjanikov / Wavelet-based Heat Kernel Derivatives -3 -6 -9 % t max L ∞ E rr o r Ours [Pat13] [HVG11] [HQ12]
Figure 14:
Comparison of the L ∞ error to the ground-truth diffu-sion wavelets for various Mexican hat wavelets. See Table 1. The
Humerus Bones data set is composed of a collection of 15humerus bone models of wild boars acquired using a 3D sensor.Each bone was scanned independently, with 24 consistent land-marks provided by experts in the field [GM13] on each shape. Theoriginal resolution of the shapes is around 25K vertices.The
SHREC’16 partial cuts correspondence bench-mark [CRB ∗ ∼
10K vertices and share a similar density.The
SHREC’19 connectivity track benchmark [MMR ∗
19] iscomposed of 44 shapes of humans. The complete benchmark con-sists in 430 shape pairs composed by meshes that represent de-formable human body shapes. Shapes belonging to these categoriesundergo changes in pose and identity. The meshes exhibit varia-tions of two different types: density (from 5K to 50K vertices) anddistribution (uniform, non-uniform). For each shape, the full SMPLmodel [LMR ∗
15] (6890 vertices) serves as our ground-truth.The
FARM partial data set is a collection of partial shapes thatwe extract from a subset of 5 meshes of the SHREC’19 con-nectivity track [MMR ∗ ∗ ∗
15] (6890 vertices). The ground-truth correspon-dence is extended to these partial shapes from the FARM registra-tion [MMRC18], which provides a ground truth dense correspon-dence between SMPL and each of the full shapes involved. & samplingScalability As a complement to Fig. 3, Table 7 displays the com-putation time required by [HQ12] and our approach for variousorder of magnitude. In all cases, the proposed method outper-forms [HQ12] by at least a factor 2. L ∞ error In addition to the L error displayed in Fig.4, we alsomeasured the L ∞ of the error, using the same experimental setup Table 7:
Computation time (in sec.) of [HQ12] compared to ourapproach on shapes (Fig. 3) from the SHREC’19 data set. Ours (s.) Improv.10 . . × . . . × . . . × . . . × .
120 20 40 60 80 1000246810 R ec on . e rr . ( × − ) Poisson diskFPS EuclideanFPS GeodesicRandom
Figure 15:
Mean reconstruction error (recon. err.) as a function ofthe number of samples ( samples) on all δ -functions for 26 shapesof the SHREC’19 connectivity track data set, using scales. as in Sect.3.2. The result, shown in Fig.14, is similar to what weobserve for the L error. Sampling strategy