Matching LBO eigenspace of non-rigid shapes via high order statistics

A fundamental tool in shape analysis is the virtual embedding of the Riemannian manifold describing the geometry of a shape into Euclidean space. Several methods have been proposed to embed isometric shapes in flat domains while preserving distances measured on the manifold. Recently, attention has been given to embedding shapes into the eigenspace of the Lapalce-Beltrami operator. The Laplace-Beltrami eigenspace preserves the diffusion distance, and is invariant under isometric transformations. However, Laplace-Beltrami eigenfunctions computed independently for different shapes are often incompatible with each other. Applications involving multiple shapes, such as pointwise correspondence, would greatly benefit if their respective eigenfunctions were somehow matched. Here, we introduce a statistical approach for matching eigenfunctions. We consider the values of the eigenfunctions over the manifold as sampling of random variables, and try to match their multivariate distributions. Comparing distributions is done indirectly, using high order statistics. We show that the permutation and sign ambiguities of low order eigenfunctions, can be inferred by minimizing the difference of their third order moments. The sign ambiguities of antisymmetric eigenfunctions can be resolved by exploiting isometric invariant relations between the gradients of the eigenfunctions and the surface normal. We present experiments demonstrating the success of the proposed method applied to feature point correspondence.


Introduction
The embedding of nonrigid shapes into a Euclidean space is well established, and widely used by shape analysis applications. Usually, the mapping from the manifold to the Euclidean space preserves distances, that is the distance measured between two points on the manifold is approximated by the respective distance calculated in the Euclidean space. The embedding of multiple isometric shapes into the same common Euclidean space seems to be ideal for applications like pointwise correspondence and shape editing. A useful property of this common embedding would be if any corresponding points of different isometric shapes were mapped to nearby target points in the Euclidean space. If this property is fulfilled then the simultaneous processing of shapes in the target domain can be done in a straightforward manner.
Elad et al . [Ela03] used classical MDS embedding into the geodesic kernel eigenspace. The MDS dissimilarity measure was based on the geodesic distances computed by the fast marching procedure [Kim98]. Bérard et al . [Bér94] used the heat operator spectral decomposition to define a metric between two manifolds M and M . They embedded the two manifolds into their respective eigenspaces. They showed that the Gromov Haussdorf distance between the embedded manifolds d GH (M, M ) = 0 if and only if the Riemannian manifolds M and M are isometric. Lafon et al . [Coi05] defined the diffusion maps and showed that the embedding into the heat kernel eigenspace is isometry invariant, and preserves the diffusion metric. Rustamov [Rus07] introduced the Global Point Signature (GPS) embedding for deformation invariant shape representation.
Although the diffusion maps computed independently on isometric shapes have a nearly compatible eigenbasis, several inconsistencies arise: • Eigenfunctions are defined up to a sign.
• The order of the eigenfunctions, especially those representing higher frequencies, is not repeatable across shapes.
• The eigenvalues of the Laplace-Beltrami operator may have multiplicity greater than one, with several eigenfunctions corresponding to each such eigenvalue.
• It is generally impossible to expect that an eigenfunction with large eigenvalue of one shape will correspond to any eigenfunction of another shape.
• Intrinsic symmetries introduce self-ambiguity, adding complexity to the sign estimation challenge.
These drawbacks limit the use of diffusion maps in simultaneous shape analysis and processing, they do not allow using high frequencies, and usually require some intervention to order the eigenfunctions or solve sign ambiguities. In this paper we present a novel method for matching eigenfunctions that were independently calculated for two nearly isometric shapes. We rely on the fact that for low order eigenfunctions, inconsistencies are governed by a small number of discrete parameters characterized by the sign sequence and permutation vector. We estimate these parameters by matching statistical properties over the spectral domain. The matching of the corresponding eigenfunctions enables the use of diffusion maps for consistent embedding of multiple isometric shapes into a common Euclidean space.

Related Work
The problems of eigenfunctions permutation and sign ambiguity were previously addressed in the context of simultaneous shape processing. Several authors, among them Shapiro and Brady [Sha92], and Jain et al . [Jai07], proposed using either exhaustive search or greedy approach for the eigenvalue ordering and sign detection. Umeyama [Ume88] proposed using a combination of the absolute values of the eigenfunctions and an exhaustive search. Mateus et al . [Mat07] expressed the connection between the eigenfunctions of two shapes by an orthogonal matrix. They formulated the matching as a global optimization problem, optimizing over the space of orthogonal matrices, and solved it using the expectation minimization approach. Later, Mateus et al . [Mat08] and Knossow et al . [Kno09] suggested using histograms of eigenfunctions values to detect their ordering and signs. Dubrovina et al . [Dub11], suggested using a coarse matching based on absolute values of eigenfunctions together with geodesic distances measured on the two shapes.
Most of these methods do not reliably resolve eigenfunction permutation [Jai07,Kno09,Mat07,Sha92]. Some of the above algorithms are limited by high complexity and do not allow the matching of more than a few eigenfunctions [Ume88,Dub11]. None of these methods reliably estimate the sign sequence of antisymmetric eigenfunctions.
At the other end, Kovnatsky et al . [Kov13] proposed to avoid the matching problem by constructing a common approximate eigenbases for multiple shapes using approximate joint diagonalization algorithms. Yet, it relies on a prior knowledge of a set of corresponding feature points.
Finally, the algorithm proposed by Pokrass et al . [Pok13] mostly resembles our approach. They used sparse modelling to match the LBO eigenfunctions that spans the Wave Kernel Signature (WKS). Yet, that approach does not reliably infers the signs of the antisymmetric eigenfunctions.

Laplace-Beltrami Eigendecomposition
Let us be given a shape modeled as a compact two-dimensional manifold M . The divergence of the gradient of a function f over the manifold, is called the Laplace-Beltrami operator (LBO) of f and can be considered as a generalization of the standard notion of the Laplace operator to manifolds [Tau95,Lév10]. The Laplace-Beltrami operator is completely derived from the metric tensor G.
where g ij = (G −1 ) ij are the components of the inverse metric tensor. Since the operator −∆ G is a positive self-adjoint operator, it admits an eigendecomposition with non-negative eigenvalues λ i and corresponding orthonormal eigenfunctions φ i , where orthonormality is understood in the sense of the local inner product induced by the Riemannian metric on the manifold. Furthermore, due to the assumption that our manifold is compact, the spectrum is discrete. We can order the eigenvalues as follows 0 = λ 1 < λ 2 < · · · < λ i < · · · . The set of corresponding eigenfunctions given by {φ 1 , φ 2 , · · · , φ i , · · · } forms an orthonormal basis of functions defined on M .

Diffusion maps
The heat equation describes the distribution of heat in time. On a manifold M , the heat equation is governed by the Laplace-Beltrami operator ∆ G , The heat kernel K t (x, y) is the diffusion kernel of the heat operator e t∆ G (t > 0). It is a fundamental solution of the heat equation with point heat source at x (heat value at point y after time t). The heat kernel can be represented in the Laplace-Beltrami eigenbasis as whereλ i are the eigenvalues of the heat operator, λ i are the eigenvalues of the LBO, Using the heat kernel we can define the diffusion distance [Coi05] where da is the area element of M . The diffusion distance d M,t (x, y) can be computed by embedding the manifold into the infinite Euclidean space spanned by the LBO eigenbasis The diffusion map {Φ t } embeds the data into the finite N -dimension Euclidean space so that in this space, the Euclidean distance is equal to the diffusion distance up to a relative truncation error (9)

Multivariate distribution comparison
The distribution of N continuous random variables φ 1 , φ 2 , ..., φ N is directly represented by the probability density function f φ 1 ,φ 2 ,...,φ N (φ 1 , φ 2 , ..., φ N ). The direct estimation of the multivariate probability density function from data samples is hard to accomplish. Therefore, an indirect representation is often being utilized. The probability distribution can be indirectly specified (under mild conditions) in a number of different ways, the simplest of which is by its raw moments In order to compare the multivariate distributions of two sets of N random we can use this indirect representation, and compare the raw moments of the random variables. In practice, only a small set of the moments I can be used for measuring the difference between the distributions where ρ i 1 ,i 2 ,...,i N are the weights associated with each raw moment.

Problem formulation
Let us denote by X and Y the two shapes we would like to match. We represent the correspondence between X and Y by a bijective mapping ϕ : We wish to find embeddings of shape X and shape Y to the finite dimensional Euclidean space, such that the corresponding points x ∈ X and ϕ(x) ∈ Y will be mapped to nearby points in the embedded space. Because of the inconsistencies described in the introduction, the diffusion maps of shapes X and Y do not necessarily fulfill this property. Our task is to modify the diffusion map Φ Y (y) by a small number of parameters θ such that the new embeddingΦ Y θ (y) will match Φ X , i.e.
. For the N low eigenvalues the matching is characterized by the following parameters: • The respective signs of the eigenfunctions s : s i ∈ {+1, −1}.
We would like to find the parametersθ = {ŝ;π}, that create the matched embedding
• ψ q : M → R are the components of an external point signature.
• ∇ G is the gradient induced by the metric tensor G.
• E[z] = M zda M , where da M is the area element of the manifold M .
• n is the normal to the surface.
• × is the cross product in R 3 and · is the inner product in R 3 .
• The weighting parameter α determines the relative weight of the gradient cost functions.
In Appendix A we give full details of the discretization we have used to implement the matching algorithm.
The application specific parameters include: • N -The number of eigenfunctions to be matched.
• {w p } P p=1 -The P nonlinear weighting functions.
• {ψ q } Q q=1 -The external point signature of size Q.
• α -The relative weight of the gradient cost functions.
In Appendix B we give the details of the application specific parameters that were used in our experiments.
Next, we review the different terms of the cost function.

Resolving sign ambiguities and permutations
For now, let us limit our discussion to resolving the sign ambiguity s. If we had known the correspondence between the two shapes, the sign of the i th eigenfunction s i could be inferred by pointwise comparison and the expectation is taken over the manifold where da X is the area element of the shape X. Unfortunately, the correspondence is unknown. Hence, pointwise comparison cannot be used in a straightforward manner. We now make the analogy between the values of the eigenfunctions over the manifold and N random variables. We consider the vector of values of the diffusion map Φ X (x) at point x as a sample out of a multivariate distribution f Φ (φ 1 (x), φ 2 (x), ..., φ N (x)). We wish to match the multivariate distributions f Φ X and f Φ Y θ . As explained in Section 1.2.3, an indirect representation of the distribution is suitable for comparing multivariate distributions. Specifically, we shall use the raw moments.
By way of construction, the non-trivial eigenfunctions have zero mean and are orthonormal. Hence, the first and second moments carry no information. Accordingly, we must use higher order moments to match the distributions. We propose to use the third order moments over the manifold M

Resolving antisymmetric eigenfunctions
For shapes with intrinsic symmetries (see [Rav07]) some of the eigenfunctions have antisymmetric distributions. The distribution of the antisymmetric eigenfunctions is agnostic to sign change. Hence, the signs of the antisymmetric eigenfunctions cannot be resolved by the simple scheme described in section 2.2.1.
The gradient of the eigenfunctions ∇φ k could be exploited to resolve the sign ambiguity: • The gradient ∇f of an antisymmetric eigenfunction f is not antisymmetric.
• The gradient is a linear operator. Consequently ∇(−f ) = −∇f, ∀f . Therefore, we can farther expand the set of variables that are used in the calculation of the raw moments, by incorporating the gradient. The gradient vector is contained in the tangent plane. Thus, the cross product of the gradients of two eigenfunctions points either outward or inward of an orientable surface. Changing the sign of one eigenfunction will flip the direction of the cross product. We can use this property to define new functions ν i,j over the manifold where n is the outward pointing normal to the tangent plane. We shall use the joint moments of the eigenfunctions and their gradients We note that Equation (16) can be farther simplified by (∇φ j × n) can be computed only once for each φ j .

Raw moments over regions
Taking the expectation over the whole shape may be too crude, especially for detecting antisymmetric sign ambiguities. We can refine the minimization criterion by taking the expectation over different regions. Remember that the correspondence between the shapes is yet unknown, therefore, directly dividing the shape into corresponding regions is impossible. Indirectly dividing the shape to different regions is possible by using the eigenfunctions themselves. The eigenfunctions φ k , k ∈ {1...N } have respective low eigenvalues which means that they have a slow rate of change. Therefore, it is possible to define functions w p (|φ k |), p ∈ {1..P } in a way that will output high or low values at different regions. For example, we can define h(|φ k |) = 1 if |φ k | > TH and zero otherwise, where TH is a scalar threshold. The output of these functions automatically divides the two shapes in a similar manner, without the use of pointwise correspondence. Moreover, because the function w p (|φ k |) is symmetric, its output does not depend on the sign of the eigenfunction φ k . We conclude that we can use w p (|φ k |) to make a weighted average of the raw moments according to different regions

Pointwise signatures as side information
We can easily use other signatures (ψ 1 , ψ 2 , ..., ψ Q ) as side information to refine the minimization criterion. Specifically, we can use signatures that carry no inconsistencies among different shapes. In our experiments we used the Heat Kernel Signature (HKS) as an additional signature [Sun09]. We can use the joint moments of the diffusion maps and the additional signatures ψ q and compute the cross of the eigenfunctions gradient ∇φ i and the signature functions gradients ∇ψ q ν S i,q = (∇φ i × ∇ψ q ) · n,

Solving the minimization problem
The minimization of Equation (12) is a non-convex optimization problem. Yet, it only involves a small number of discrete parameters. Therefore, an exhaustive search is possible. In practice, we implemented the search in four steps: • Step 1 -An initialization of s 0 is determined by s i = sign(µ X i,i,i µ Y i,i,i ) and π 0 = [0, 1, ...N ].
We make an educated guess for the possible permutations, limiting the search for two permutation profiles: two consecutive eigenfunction switching (with possible sign change), i.e.
In this step all possible quadruple sign changes are checked, setting the permutation vector found in Step 2. If the cost function was decreased in Step 2 or Step 3, then return to Step 2. While finding the optimal sign sequence and permutation vector, we keep a list of all possible good sign sequences for the next step. • Step 4: The optimal sign sequenceŝ is found by comparing the entire cost function C(s, π) + C S (s, π) + α(C P ∇ (s, π) + C P,S ∇ (s, π)) for each sign sequence in the list created in Step 3.
We note that the computation of the moments can be done prior to the minimization algorithm. This calculation of the raw moments of shape X is independent of shape Y and can be performed for each shape before the matching procedure. The entire cost function given in Equation (12) can be computed from the raw moments and the parameters, without the use of the eigenfunctions themselves.

Results
We tested the proposed method on pairs of shapes represented by triangulated meshes from the TOSCA database [Bro08]. Figures 1, 3  For example, we can see in Figure 1 that the eigenfunction matching algorithm swapped φ Y 3 and φ Y 4 . It also correctly flipped the signs of φ Y 1 , φ Y 2 and φ Y 4 . We also notice that the matching algorithm was able to detect the correct signs of the antisymmetric eigenfunctions. For example, in Figure 3, the sign of the antisymmetric eigenfunction φ Y 2 was correctly flipped, while keeping the sign of φ Y 3 . We applied the matched eigenfunctions for detecting feature point correspondence between the two shapes. A selected number of feature points from the first shape were matched to the second one using a combination of two signatures: • The matched low order eigenfunctions that represent global strucure of the shapes.
• The Heat Kernel Signature (HKS) derivative (Equation (33)), that being a bandpass filter, expresses more local features. Figures 2, 4 and 6 show that the correspondences between feature points were found correctly. Notice that this approach was able to resolve the symmetries of the given shapes.

Conclusion
The Laplace Beltrami operator (LBO) provides us with a flat eigenspace in which surfaces could be represented as canonical forms in an isometric invariant manner. However, the order and directions (signs) of the axises in this Hilbert space do not have to correspond when two isometric surfaces are considered. In order to resolve such potential ambiguities we resorted to high order statistics of the eigenfunctions of the LBO and their interaction with the surface normal. It appears that these cross moments allow for ordered directional matching of the components of corresponding eigenspaces. We demonstrated that resolving the sign and order correspondence allows for shape matching in various scenarios. In the future we plan on extending the proposed framework to enable it to deal with more generic transformations, like the scale invariant metric introduced by Aflalo et al . [Afl11]. solved the generalized eigendecomposition problem, as suggested by Rustamov [Rus07]. where and A is a diagonal matrix. A ii equals the voroni area about vertex i.

A.2 Gradient
We assume that the function f is linear over the triangle with vertices V i , V j , V k with values f i , f j , f k at the vertices. We define the local coordinates (u, v) with coordinates (0, 0), (0, 1), (1, 0) at the vertices V i , V j , V k . Because f is assumed to be linear and which can be written as The Jacobian By the chain rule ∂f ∂(u,v) = ∂f ∂(x,y,z) ∂(x,y,z) ∂(u,v) and in matrix form (DF ) T = (∇f ) T J T , or equivalently By taking the pseudoinverse we get the discrete gradient operator over a triangle

B Application specific parameters
In our experiments we used the following specific parameters for the matching algorithm: • We matched N = 10 eigenfunctions.
• For generating the external pointwise signature ψ q , the Heat Kernel Signature (HKS) was used [Sun09]. In the approximation of the Heat Kernel Signature HKS t (x) = h i=1 e −λ i t φ 2 i (x), we used h = 30 eigenfunctions. We used a bandpass filter form of the HKS by taking the derivative of the Heat Kernel Signature. The HKS derivative was logarithmically sampled Q = 6 times at t = t q , q = 1, 2, ...Q, with t 1 = 1 50λ 1 and t Q = 1 λ 1 . ψ q were normalized according to the inner product over the manifold.
• The relative weight parameter α was set by balancing the influence of the terms of the cost function.