CCoupled quasi-harmonic bases
A. Kovnatsky , M. M. Bronstein , A. M. Bronstein ,K. Glashoff , and R. Kimmel Faculty of Informatics, Universit`a della Svizzera Italiana School of Electrical Engineering, Tel Aviv University Department of Mathematics, Hamburg University Department of Computer Science, Technion
October 2, 2012
Abstract
The use of Laplacian eigenbases has been shown to be fruitful in manycomputer graphics applications. Today, state-of-the-art approaches toshape analysis, synthesis, and correspondence rely on these natural har-monic bases that allow using classical tools from harmonic analysis onmanifolds. However, many applications involving multiple shapes are ob-stacled by the fact that Laplacian eigenbases computed independently ondifferent shapes are often incompatible with each other. In this paper, wepropose the construction of common approximate eigenbases for multipleshapes using approximate joint diagonalization algorithms. We illustratethe benefits of the proposed approach on tasks from shape editing, posetransfer, correspondence, and similarity.
It is well-established that the eigenfunctions of the Laplace-Beltrami operator( manifold harmonics ) of a 3D shape modeled as a 2-manifold play the roleof the Fourier basis in the Euclidean space [Tau95, L´ev06]. Methods basedon the Laplace-Beltrami operator have been used in a wide range of applica-tions, including remeshing [Kob97, NISA06], parametrization [FH05], compres-sion [KG00], recognition [RWP05, Rus07], and clustering. Many methods incomputer graphics and geometry processing draw inspiration from the worldof physics, finding analogies between physical processes such as heat diffusion[CL06] or wave propagation and the geometric properties of the shape [SOG09].Several papers have studied consistent discretizations of the Laplace-Beltramioperator [PP93, MDSB03, WMKG08]The influential paper of Taubin [Tau95] drew the analogy between the clas-sical signal processing theory and manifold harmonics, showing that standardtools in signal processing such as analysis and synthesis of signals can be carried1 a r X i v : . [ c s . C V ] S e p ψ ψ ψ ψ φ φ φ φ φ ˆ ψ ˆ ψ ˆ ψ ˆ ψ ˆ ψ ˆ φ ˆ φ ˆ φ ˆ φ ˆ φ Figure 1:
Top: Laplace-Beltrami eigenfunctions { φ i } and { ψ i } of cat and lion shapes:besides trivial sign flips in eigenfunctions 2 and 3, higher frequency eigenfunctions(4-6) manifest quite different behaviors. Bottom: coupled basis function { ˆ φ i } , { ˆ ψ i } obtained using approximate joint diagonalization procedure proposed in this paperbehave consistently. Hot and cold colors represent positive and negative values, re-spectively. out on manifolds. This idea was extended in [KR05] and later in [L´ev06, LZ09],who showed a practical framework for shape filtering and editing using the man-ifold harmonics transform. In [OBCS + compatible with each other.However, this assumption is often unrealistic. First, eigenfunctions are onlydefined up to sign flips for shapes having simple spectra (i.e., no multiplicityof eigenvalues). In this case every vector in the eigenspace is called an eigen-vector. In the more general case, eigenfunctions corresponding to an eigenvaluewith non-trivial multiplicity are not defined at all, and one can select an ar-bitrary orthonormal basis spanning each such an eigenspace. Second, due tonumerical instabilities, the ordering of the eigenfunctions, especially those rep-resenting higher frequencies, is not repeatable across shapes. Finally, harmonicbases computed independently on different shapes can be expected to be rea-sonably compatible only when the shapes are approximately isometric, sinceisometries preserve the eigenfunctions of the Laplace-Beltrami operator. Whenthis assumption is violated, it is generally impossible to expect that the n -thharmonic of one shape will correspond to the n -th harmonic of another shape.2hese drawbacks limit the use of harmonic bases in simultaneous shape analysisand processing to approximately isometric shapes, they do not allow to use highfrequencies, and usually require some intervention to order the eigenfunctionsor solve sign ambiguities. Contributions.
In this paper, we propose a general framework allowing toextend the notion of harmonic bases by finding a common (approximate) eigen-basis of multiple Laplacians. Numerically, this problem is posed as approximatejoint diagonalization of several matrices. Such methods have received limited at-tention in the numerical mathematics community [BGBM93] and have been em-ployed mostly in blind source separation applications [CS93,CS96,Yer02,Zie05].Philosophically similar ideas have been used in the machine learning communityfor multi-view spectral clustering [KRDI11]. To the best of our knowledge, thisis the first time they are applied to problems in computer graphics. We show afew example of applications of such coupled quasi-harmonic bases in Section 5.
Let us be given a shape modeled as a compact two-dimensional manifold X .Given a smooth scalar field f on X , the negative divergence of the gradient ofa scalar field, ∆ f = − div grad f , is called the Laplace-Beltrami operator of f and can be considered as a generalization of the standard notion of the Laplaceoperator to manifolds [Tau95, LZ09].Since the Laplace-Beltrami operator is a positive self-adjoint operator, it ad-mits an eigendecomposition with non-negative eigenvalues λ and correspondingorthonormal eigenfunctions φ , ∆ φ = λφ (1)where orthonormality is understood in the sense of the local inner product in-duced by the Riemannian metric on the manifold. Furthermore, due to theassumption that our manifold is compact, the spectrum is discrete, 0 = λ <λ < · · · . In physics, (1) is known as the Helmholtz equation representing thespatial component of the wave equation. Thinking of our shape as of a vibrat-ing membrane, the eigenfunctions φ i can be interpreted as natural vibrationmodes of the membrane, while the λ i ’s assume the meaning of the correspond-ing vibration frequencies. The eigenbasis of the Laplace-Beltrami operator isfrequently referred to as the harmonic basis of the manifold, and the functions φ i as manifold harmonics .In the discrete setting, we represent the manifold X as a triangular meshbuilt upon the vertex set { xxx , . . . , xxx n } . A function f on the manifold is rep-resented by the vector fff = ( f (xxx ) , . . . , f (xxx n )) T of its samples. A common ap-proach to discretizing manifold harmonics is by first constructing a discreteLaplace-Beltrami operator on the mesh, represented as an n × n matrix, fol-lowed by its eigendecomposition. A popular discretization is the cotangentscheme [MDSB03], resulting in the generalized eigenvalue problem WWWΦΦΦ = DDDΦΦΦΛΛΛ,3r equivalently, the eigenvalue problem LLLΦΦΦ = ΦΦΦΛΛΛ, where LLL = DDD − WWW,DDD = diag( s , . . . , s n ) / , (2) w ij = (cid:26) (cot( α ij ) + cot( β ij )) / i (cid:54) = j ; (cid:80) k (cid:54) = i w ik i = j, (3) s i > i , and α ij , β ij are the two angles opposite to the edge between vertices i and j inthe two triangles sharing the edge [WBH +
07, RWP06]. Here, ΦΦΦ = ( φφφ , . . . , φφφ n )is the matrix of eigenvectors arranged as columns, discretizing the eigenfunc-tions of the Laplace-Beltrami operator at the sampled points. Numerically, theaforementioned eigendecomposition problem can be posed as the minimizationmin ΦΦΦ off(ΦΦΦ T WWW X ΦΦΦ) s . t . ΦΦΦ T DDD X ΦΦΦ = III , (4)of the sum of squared off-diagonal elements, off(AAA) = (cid:80) i (cid:54) = j a ij [CS96]. Severalclassical algorithms for finding eigenvectors such as the Jacobi method in facttry to reduce the off-diagonal values in an iterative way.Laplace-Beltrami eigenbases are equivalent to Fourier bases on Euclideandomains, and allow to represent square-integrable functions on the manifold aslinear combinations of eigenfunctions, akin to Fourier analysis. In particular,solutions of PDEs on non-Euclidean domains can be expressed in the Laplacianeigenbasis, giving rise to numerous efficient methods for computing e.g. localdescriptors based on fundamental solutions of heat and wave equations [SOG09,ASC11], isometric embeddings of shapes [BN03,Rus07], diffusion metrics [CL06],shape correspondence and similarity [RWP05, BBK +
10, OMMG10, DK10]Unlike the Euclidean plane where the Fourier basis is fixed, in non-Euclideanspaces, the harmonic basis depends on the domain on which it is defined. Manyapplications working with several shapes (such as shape matching or pose trans-fer) rely on the fact that harmonic basis defined on two or more different shapesare consistent and behave in a similar way [L´ev06]. While experimentally it isknown that often low-frequency harmonics have similar behavior (finding pro-trusions in shapes, a fact often employed for shape segmentation [Reu10]), thereis no theoretical guarantee whatsoever of such behavior. Theoretically, consis-tent behavior of eigenfunctions can be guaranteed only in very restrictive set-tings of isometric shapes with simple spectrum [OBCS + Functional correspondence.
Ovsjanikov et al. [OBCS +
12] proposed anelegant way to avoid direct representation of correspondences as maps betweenshapes using a functional representation. The authors noted that when twoshapes X and Y are related by a bijective correspondence t : X → Y , thenfor any real function f : X → R , one can construct a corresponding function g : Y → R as g = f ◦ t − . In other words, the correspondence t uniquely definesa mapping between two function spaces T : F ( X, R ) → F ( Y, R ), where F ( X, R )denotes the space of real functions on X . Furthermore, such a mapping is linear .4aplacian eigenbases, T ( φ i ) = (cid:80) j> c ij ψ i Common bases, T ( ˆ φ i ) = (cid:80) j> c ij ˆ ψ i Figure 2:
Matrix CCC of coefficients expressing a given correspondence between twoposes of an elephant (left) and elephant and horse (right) in functional representationof [OBCS + Equipping X and Y with harmonic bases, { φ i } i ≥ and { ψ j } j ≥ , respectively,one can represent a function f : X → R using the set of (generalized) Fouriercoefficients { a i } i ≥ as f = (cid:80) i ≥ a i φ i . Then, translating the representation intothe other harmonic basis, one obtains a simple representation of the correspon-5ence between the shapes T ( f ) = (cid:88) i,j ≥ a i c ij ψ j , (5)where c ij are Fourier coefficients of the basis functions of X expressed in thebasis of Y , defined as T ( φ i ) = (cid:80) j ≥ c ij ψ j . The correspondence can be thus byapproximated using k basis functions and encoded by a k × k matrix CCC = ( c ij )of these coefficients, referred to as the functional matrix in [OBCS + t : X → Y is trans-lated into a simpler task of finding CCC from a set of correspondence constraints.This matrix has a diagonal structure if the harmonic bases are compatible, anassumption crucial for the efficient computation of the correspondence. How-ever, the authors report that in practice the elements of CCC spread off the diagonalwith the increase of the frequency due to the lack of perfect compatibility of theharmonic bases. XXX = (cid:80) i ≥ (cid:104) XXX , φ i (cid:105) φ i YYY = (cid:80) i ≥ (cid:104) YYY , ψ i (cid:105) ψ i ZZZ = (cid:80) i =1 (cid:104) XXX , φ i (cid:105) ψ i + (cid:80) i> (cid:104) YYY , ψ i (cid:105) ψ i ZZZ = (cid:80) i =1 (cid:104) XXX , ˆ φ i (cid:105) ˆ ψ i + (cid:80) i> (cid:104) YYY , ˆ ψ i (cid:105) ˆ ψ i Figure 3:
Pose transfer from horse (leftmost) to camel shape (second from left) bysubstituting the first 6 Fourier coefficients in the decomposition of extrinsic coordinatesof the shape in the Laplacian eigenbasis as done in [L´ev06] (third from left) and jointapproximate eigenbasis (rightmost). Coupling between the basis function is crucial forthis approach to work.
Pose transfer.
L´evy [L´ev06] proposed a pose transfer approach basedon the Fourier decomposition of the manifold embedding coordinates. Giventwo shapes X and Y embedded in R with the corresponding harmonic bases { φ i } i ≥ and { ψ i } i ≥ , respectively, one first obtains the Fourier decompositionsof the embeddings XXX = (cid:88) i ≥ aaa i φ i , YYY = (cid:88) i ≥ bbb i ψ i (6)(we denote by XXX and YYY the Euclidean embeddings of manifolds X and Y , andby aaa i , bbb i the three-dimensional vectors of the Fourier coefficients correspondingto each embedding coordinate). Next, a new shape Z is composed according toZZZ = n (cid:88) j =1 aaa i ψ i + (cid:88) i>n bbb i ψ i , (7)6ith the first n low frequency coefficients taken from X , and higher frequenciestaken from Y . This transfers the “layout” (pose) of the shape X to the shape Y while preserving the geometric details of Y . This method works when theaaa i ’s and the bbb i ’s are expressed in the same “language”, i.e., when the Laplacianeigenfunctions behave consistently in X and Y . Let us be given two shapes
X, Y with the corresponding Laplacians ∆ X , ∆ Y . If X and Y are related by an isometry t : X → Y and have simple Laplacianspectrum (no eigenvalues with multiplicity greater than 1), the eigenfunctionsare defined up to a sign flip, ψ i = ± φ i ◦ t − . If some eigenvalue λ i = . . . = λ i + p has multiplicity p + 1, the individual eigenvectors are not defined, but ratherthe subspaces they span: span { ψ i , . . . , ψ i + p } = span { φ i , . . . , φ i + p } ◦ t − . Moregenerally, if the shapes are not isometric, the behavior of their eigenfunctionscan differ dramatically (Figure 1, top).This poses severe limitations on applications we mentioned in the previ-ous section. Figure 2 (middle) shows the coefficient matrix CCC defined in (5),representing the functional correspondence between two shapes. For near-isometric shapes (two deformations of an elephant, Figure 2, middle left), since ψ i ≈ ± φ i ◦ t − , the coefficients c ij ≈ ± δ ij , and thus the matrix CCC is nearly diag-onal. However, when trying to express correspondence between non-isometricshapes (elephant and horse, Figure 2, middle right), the Laplacian eigenfunc-tions manifest a very different behavior breaking this diagonality.The same problem is observed when we try to use the technique of L´evy[L´ev06] for pose transfer by expressing the embedding coordinates of the shapein the respective Laplacian eigenbasis and substituting the low-frequency coef-ficients from another shape. L´evy disclaims that his method works “providedthat the eigenfunctions that correspond to the lower frequencies match” [L´ev06].However, such a consistent behavior is not guaranteed at all; Figure 3 (thirdfrom left) shows how the pose transfer breaks when the eigenfunction are incon-sistent.The main idea of this paper is to try to find bases ˆ φ i , ˆ ψ i that approximatelydiagonalize the respective Laplacians (∆ X ˆ φ ≈ λ X ˆ φ , ∆ Y ˆ ψ ≈ λ Y ˆ ψ ) and arecoupled ( ˆ ψ i ≈ ˆ φ i ◦ t − ). Such new coupled bases , while being nearly harmonic,make the basis functions consistent across shapes and cure the problems weoutlined above (Figure 1, bottom). Figure 2 (bottom) shows the functionalcorrespondence represented in the coupled bases { ˆ φ i } i ≥ , { ˆ ψ i } i ≥ . Due to thecoupling, the basis functions behave consistently resulting in almost perfectlydiagonal matrices CCC even when the shapes are highly non-isometric (bottomright). Likewise, Figure 3 (rightmost) shows that pose transfer using Fouriercoefficients in the coupled bases works correctly for shapes with inconsistentLaplacian eigenfunctions. We consider the case of a pair of shapes for the mere sake of simplicity. Extension to acollection of more shapes is straightforward. pproximate joint diagonalization. We assume that the shapes aresampled at n X , n Y points, and their Laplacians are discretized as matricesWWW X , WWW Y and DDD X , DDD Y of size n X × n X and n Y × n Y respectively, as definedin (3). We further assume that we know the correspondence between l points x j , . . . , x j l and y j (cid:48) , . . . , y j (cid:48) l (as we show in Section 5, very few and very roughcorrespondences are required). We represent these correspondences by an l × n X matrix PPP with elements p i,j i = 1 and zero elsewhere; the l × n Y matrix QQQ isdefined accordingly as q i,j (cid:48) i = 1 zero elsewhere.The problem of joint approximate diagonalization (JD) can be formulatedas the coupling of two problems (4),min ˆΦΦΦ , ˆΨΨΨ off( ˆΦΦΦ T WWW X ˆΦΦΦ) + off( ˆΨΨΨ T WWW Y ˆΨΨΨ) + µ (cid:107) PPP ˆΦΦΦ − QQQ ˆΨΨΨ (cid:107) s . t . ˆΦΦΦ T DDD X ˆΦΦΦ = III , ˆΨΨΨ T DDD Y ˆΨˆΨˆΨ = III (8)where off denotes some off-diagonality penalty, e.g., the sum of the squaredoff-diagonal elements as defined in Section 2. The parameter µ determines thecoupling strength. For µ = 0, the problem becomes uncoupled and boils downto individual diagonalization (4) of the two Laplacians. The joint eigenvectorsin this setting coincide with the eigenvectors of the Laplacians: ˆΦΦΦ = ΦΦΦ andˆΨΨΨ = ΨΨΨ.We can parametrize the joint basis functions as linear combinations of theLaplacian eigenvectors, ˆΦΦΦ = ΦΦΦAAA and ˆΨΨΨ = ΨΨΨBBB, where AAA and BBB are matrices ofcombination coefficients of size n X × n X and n Y × n Y , respectively. Noticingthat ˆΦΦΦ T WWW X ˆΦΦΦ = AAA T ΦΦΦ T WWW X ΦΦΦAAA = AAA T ΛΛΛ X AAA, and same way, ˆΨΨΨ T WWW Y ˆΨΨΨ = BBB T ΛΛΛ Y BBB,we transform problem (8) intomin
AAA , BBB off(AAA T ΛΛΛ X AAA) + off(BBB T ΛΛΛ Y BBB) + µ (cid:107) PPPΦΦΦAAA − QQQΨΨΨBBB (cid:107) s . t . AAA T AAA = III , BBB T BBB = III (9)Since AAA and BBB are orthonormal, they act as isometries in the respective eigenspacesof the Laplacians of X and Y . We can thus think geometrically of problem (9)as an attempt to rotate and reflect the eigenbases ¯ΦΦΦ and ¯ΨΨΨ such that they alignin the best way (in the least squares sense) at corresponding points, while stillapproximately diagonalizing the Laplacians.Since in many applications we are not interested in the entire eigenbasisbut in the first k eigenvectors, this formulation is especially convenient, as itallows us to express the first k joint eigenvectors as a linear combination of k (cid:48) eigenvectors (we provide a justification of this assumption in Appendix A), thushaving the matrices AAA and BBB of size k (cid:48) × k :min AAA , BBB off(AAA T ¯ΛΛΛ X AAA) + off(BBB T ¯ΛΛΛ Y BBB) + µ (cid:107) PPP ¯ΦΦΦAAA − QQQ ¯ΨΨΨBBB (cid:107) s . t . AAA T AAA = III , BBB T BBB = III , (10)where ¯ΦΦΦ = ( φφφ , . . . , φφφ k (cid:48) ) and ¯ΛΛΛ X = diag( λ Xk , . . . , λ Xk (cid:48) ); matrices ¯ΨΨΨ , ¯ΛΛΛ Y are de-fined accordingly. Typically, k, k (cid:48) (cid:28) n X , n Y , and thus the problem is much8maller than the full joint diagonalization (8). Note that the coupling term pro-vides lk constraints, so in order not to over-determine the problem, we shouldhave 2 k (cid:48) > l . Typical values used in our experiments were k (cid:48) ∼ l ∼
15, and k ∼ Off-diagonality penalty.
It is important to note that in problem (10)the use of the sum of squared off-diagonal elements as the off-penalty does notproduce an ordered set of approximate joint eigenvectors. This issue can besolved by using an alternative penalty, (cid:107) ˆΦΦΦ T LLL X ˆΦΦΦ − ¯ΛΛΛ X (cid:107) = (cid:107) AAA T ¯ΛΛΛ X AAA − ¯ΛΛΛ X (cid:107) , (11)which is similar to the sum of squared off-diagonal elements but also includesthe difference of the diagonal elements. The penalty for the shape Y is definedin the same way. Coupling term.
It is possible to change the coupling strength at differ-ent frequencies, by using a more generic form of the coupling term | (PPP ¯ΦΦΦAAA − QQQ ¯ΨΨΨBBB)VVV (cid:107) , where VVV = diag( v , . . . , v k ) is a diagonal matrix of weights v i de-creasing with frequency. Such frequency-dependent coupling can be useful inapplications like pose transfer or functional correspondence we are interestedin strong coupling at low frequencies and can afford weaker coupling of higherones. Procrustes problem.
Interestingly, in the limit case µ → ∞ where wecan ignore the off-diagonal penalties, problem (10) becomesmin AAA , BBB (cid:107)
PPP ¯ΦΦΦAAA − QQQ ¯ΨΨΨBBB (cid:107) s . t . AAA T AAA = III , BBB T BBB = III , (12)Using the invariance of the Frobenius norm under orthogonal transformation,we can rewrite problem (12) as an orthogonal Procrustes problem min ΩΩΩ (cid:107)
PPP ¯ΦΦΦ − QQQ ¯ΨΨΨΩΩΩ (cid:107) s . t . ΩΩΩ T ΩΩΩ = III , (13)where ΩΩΩ = BBBAAA T . The problem has an analytic solution ΩΩΩ = SSSRRR T , where¯ΦΦΦ T PPP T QQQ ¯ΨΨΨ = SSSΣΣΣRRR T is the singular value decomposition of the matrix ¯ΦΦΦ T PPP T QQQ ¯ΨΨΨwith left- and right singular vectors SSS , RRR [Sch66]. Then, AAA = SSS and BBB = RRR.
Problem (10) is a non-linear optimization problem with orthogonality con-straints. In our experiments, we used the first-order constrained minimizationalgorithm implemented in MATLAB Optimization Toolbox. We provide belowthe gradients of our cost function.
Gradient of the off-diagonality penalty is given by ∇ AAA (cid:88) i (AAA T ¯ΛΛΛ X AAA) ii = 4(¯ΛΛΛ X AAAAAA T ¯ΛΛΛ X AAA − OOO ◦ AAA¯ΛΛΛ X ) , (14)9here OOO is a matrix of equal columns ooo = diag(AAA T ¯ΛΛΛ X AAA) and ◦ denotes element-wise product of matrices (see derivation in Appendix B). Gradient of the alter-native penalty (11) is derived similarly as ∇ AAA (cid:107)
AAA T ¯ΛΛΛ X AAA − ¯ΛΛΛ X (cid:107) = 4(¯ΛΛΛ X AAAAAA T ¯ΛΛΛ X AAA − ¯ΛΛΛ X AAA¯ΛΛΛ X ) . (15) Gradient of the coupling term w.r.t. to AAA is given by ∇ AAA (cid:107)
PPP ¯ΦΦΦAAA − QQQ ¯ΨΨΨBBB (cid:107) = 2 ¯ΦΦΦ T PPP T (PPP ¯ΦΦΦAAA − QQQ ¯ΨΨΨBBB); (16)the gradient w.r.t. BBB is obtained in the same way.
Initialization.
Assuming the Laplacians are nearly jointly diagonalizableand have simple spectrum, their joint eigenvectors will be equal to the harmonicbasis functions up to sign flips. Thus, in this case AAA and BBB are diagonal matricesof ±
1. This was found to be a reasonable initialization to the joint diagonaliza-tion procedure. We set AAA = III, and then solve sign flip by setting the elementsof BBB to b ij = +1 i = j, (cid:107) PPP φφφ i − QQQ ψψψ i (cid:107) ≤ (cid:107) PPP φφφ i + QQQ ψψψ i (cid:107) ; − i = j, (cid:107) PPP φφφ i − QQQ ψψψ i (cid:107) > (cid:107) PPP φφφ i + QQQ ψψψ i (cid:107) ;0 else . (17)Initialized this way, the problem starts with decoupled but diagonalizing bases,and the optimization tries to improve the coupling. Band-wise computation.
In applications requiring the computation ofmany approximate joint eigenvectors ( k (cid:29) k × k (cid:48) matrices AAA , BBB, we can split the eigenvectors into non-overlappingbands of size k (cid:48)(cid:48) , and solve k/k (cid:48)(cid:48) problems (10) for k (cid:48)(cid:48) × k (cid:48)(cid:48) matrices AAA , BBB (herewe assume that k is a multiple of k (cid:48)(cid:48) and that k (cid:48)(cid:48) > l ). For the i th band, we use¯ΦΦΦ = ( φφφ ( i − k (cid:48)(cid:48) , . . . , φφφ ik (cid:48)(cid:48) ) and ¯ΛΛΛ X = diag( λ X ( i − k (cid:48)(cid:48) , . . . , λ Xik (cid:48)(cid:48) ) and ¯ΨΨΨ , ¯ΛΛΛ Y definedaccordingly in (10). In this section, we show additional examples of coupled bases construction, aswell as some potential applications of the proposed approach. We used shapesfrom publicly available datasets [BBK08, SP04, SMKF04]. Mesh sizes variedwidely between 600 - 25K vertices. Discretization of the Laplace-Beltrami op-erator was done using the cotangent formula [MDSB03]. In all our examples,we constructed coupled bases solving the JD problem (10) with off-diagonalitypenalty (11), as described in Section 4. Typical time to compute 15 joint eigen-vectors was about 1 minute.Figure 4 (top) shows examples of joint diagonalization of Laplacians of dif-ferent shapes: near isometric (two poses of an elephant) and non-isometric (ele-phant and horse). We computed the first k = 20 joint approximate eigenvectors.The Laplacians are almost perfectly diagonalized by the obtained coupled bases10 .
02 0 .
018 0 .
44 0 . .
49 0 .
19 0 .
28 0 . Figure 4:
Examples of joint diagonalization of Laplacians of near-isometric shapes(two poses of an elephant, top left) and non-isometric shapes (elephant and horse,top right; four humanoids, bottom). Second and fourth rows show the approximatediagonalization of the respective Laplacians. Coupling was done using 40 points forthe elephant and horse shapes and 25 for the humanoid shapes. Numbers show theratio of the norm of the diagonal and off-diagonal values. in the case of near-isometric shapes; for non-isometric shapes, off-diagonal ele-ments are more prominent. Nevertheless, a clear diagonally-dominant structureis present. Figure 4 (bottom) shows an additional example of joint diagonaliza-tion of the Laplacian matrices of four humanoid shapes by first 15 elements ofthe coupled bases.
Sensitivity to correspondence error.
Figure 5 shows the sensitivityof the presented procedure to errors in correspondence used for coupling. Inthis experiment, we used coupling at 10 points that deviated from groundtruthcorrespondence by up to 15% of the geodesic diameter of the shape. Figure 6shows examples of the first elements of the coupled bases obtained in the lat-11% error .
065 0 .
3% error .
073 0 .
6% error .
092 0 .
15% error .
152 0 . Figure 5:
Sensitivity of joint diagonalization to errors in correspondence (in % ofgeodesic diameter of the shape) used in the coupling term. Shapes are shown withsimilar colors representing corresponding points. Correspondences between 10 pointsused for coupling are shown with lines. Numbers show the ratio of the norm of thediagonal and off-diagonal values. ˆ φ ˆ φ ˆ φ ˆ φ ˆ ψ ˆ ψ ˆ ψ ˆ ψ Figure 6:
Coupled bases elements of the human shapes from Figure 5 obtained us-ing 10 points with inaccurate correspondence (error of 15% geodesic diameter) forcoupling. ter setting. This experiment illustrates that very few roughly correspondingpoints are required for the coupling term in our problem, and that the proposedprocedure is robust to correspondence noise.
Shape correspondence.
Ovsjanikov et al. [OBCS +
12] computed cor-respondence between near-isometric shapes from a set of constraints on the12 × k functional matrix CCC (encoding the correspondence represented using thefirst k harmonic functions as described in Section 2). Given a set of functions f , . . . , f p on X and corresponding functions g , . . . , g q on Y (for example, in-dicator functions of stable regions) represented as n X × p and n Y × p matricesFFF = (fff , . . . , fff p ) and GGG = (ggg , . . . , ggg p ), respectively, CCC is recovered from a systemof pk equations with k variables,FFF T ΦΦΦ = GGG T ΨΨΨCCC . (18)Additional constrains stemming from the properties of the matrix CCC are alsoadded [OBCS + ≈ diag( c , . . . , c kk ), we can rewrite (18) as a system of pk equations with only k variables corresponding to the diagonal of CCC, diag(ggg T1 ˆΨΨΨ)...diag(ggg T p ˆΨΨΨ) c ... c kk = ˆΦΦΦ T fff ...ˆΦΦΦ T fff p . (19)Equation (19) allows to use significantly less data to fully determine the cor-respondence, and is also more computationally efficient. The method can beapplied iteratively, by first establishing some (rough) correspondence to pro-duce coupled bases, then finding CCC by solving (19) and using the computedcorrespondence to refine the coupling and recompute the bases.Figure 7 shows an example of finding functional correspondence betweennon-isometric shapes of human and gorilla. As functions { f i } pi =1 and { g i } pi =1 we used binary indicator functions of p = 25 regions detected using the MSERalgorithm [LBB11]. We compared the method described in [OBCS +
12] for com-puting CCC by solving the system (18) in the Laplace-Beltrami eigenbases (Fig-ure 7, middle) and the diagonal-only approximation (19) in the coupled bases. Inboth cases, we used k = 20 first basis vectors. Then, point-wise correspondencewas obtained from CCC using the ICP-like approach proposed in [OBCS + Simultaneous mesh editing.
Rong et al. [RCG08] proposed an approachfor mesh editing based on elastic energy minimization. Given a shape withembedding coordinates XXX, the method attempts to find a deformation field dddproducing a new shape XXX (cid:48) = XXX + ddd, providing a set of user-defined n (cid:48) anchorpoints for which the displacement is known (w.l.o.g. assuming to be the first n (cid:48) points, ddd i = ddd (cid:48) i for i = 1 , . . . , n (cid:48) ), as a solution of the system of equations (cid:20) k b LLL X − k c LLL X MMMMMM T (cid:21) (cid:20) ddd γγγ (cid:21) = (cid:20) (cid:48) (cid:21) , (20)where MMM = (III , T is an n × n (cid:48) identity matrix, γγγ are unknown Lagrange multipli-ers corresponding to the constraints on anchor points, and k b , k c are parameters13 ..... Figure 7: Functional correspondence computation. Left: some of the MSERregions. Second and third columns: correspondence and matrix CCC computedusing the method of [OBCS +
12] in Laplace-Beltrami eigenbases and diagonal-only approach in coupled bases, respectively. Fourth column: groundtruth cor-respondence matrix CCC represented in coupled bases.trading off between resistance to bending and stretching, respectively [RCG08] .The system of equations can be expressed in the frequency domain using k (cid:28) n first harmonic basis functions, (cid:20) ¯ΦΦΦ T ( k b LLL X − k c LLL X ) ¯ΦΦΦ ¯ΦΦΦ T MMMMMM T ¯ΦΦΦ 000 (cid:21) (cid:20) αααγγγ (cid:21) = (cid:20) (cid:48) (cid:21) (21)where ααα = ¯ΦΦΦ T ddd are the k Fourier coefficients. The desired deformation field isobtained by solving the system of equations for ααα and transforming it to thespatial domain ddd = ¯ΦΦΦ ααα (for details, the reader is referred to [RCG08]).Using coupled bases, it is possible to easily extend this approach to simul-taneous editing of multiple shapes, solving the system (21) with the coupledbasis ˆΦΦΦ in place of ¯ΦΦΦ, and applying the deformation to the second mesh us-ing ddd = ˆΨΨΨ ααα . Figure 8 exemplifies this idea, showing how a deformation of thecat shape is automatically transferred to the lion shape, which accurately andnaturally repeats the cat pose.
Shape similarity.
The diagonalization quality of the Laplacians can beused as a criterion for shape similarity, with isometric shapes resulting in idealdiagonalization. With this approach, it is possible to compare two shapes froma small number of inaccurate correspondences provided for coupling in the jointdiagonalization problem. It is possible to compare more than two shapes at once,with the diagonalization quality acting like some “variance” of the collection ofshapes.Figure 9 shows the similarity matrix between 25 shapes belonging to 8 dif-ferent classes. Each shape is present with 3-4 near-isometric deformation. Weused 25 points for coupling; dissimilarity of a pair of shapes was computed by14igure 8:
Simultaneous shape editing in the frequency domain using the approachof [RCG08]. Top: editing the cat shape (green points denote the anchor vertices usedas constraints in problem (20)). Bottom: the same pose is automatically transferredto the lion shape using coupled basis. jointly diagonalizing the respective Laplacians and then computing the averageratio of the norms of the diagonal and off-diagonal elements of both matrices.
We showed several formulations of numerically efficient joint approximate diag-onalization algorithms for the construction of coupled bases of the Laplacians ofmultiple shapes. Such quasi-harmonic bases allow to extend many shape anal-ysis and synthesis tasks to cases where the standard harmonic bases computedon each shape separately cease being compatible. The proposed constructioncan be used as an alternative to the standard harmonic bases. A particularlypromising direction is non-rigid shape matching in the functional correspon-dence representation. The sparse structure of the functional matrix CCC has notbeen yet taken advantage of for the computation of correspondence in a properway, and naturally calls for sparse modeling methods successfully used in signalprocessing. However, the correctness of this model largely depends on the basisin which CCC is represented, and standard Laplace-Beltrami eigenbasis performsquite poorly when one deviates from the isometry assumptions. In follow-upstudies, we intend to explore the relation between sparse models and joint ap-proximate diagonalization in shape correspondence problems.
Appendix A - Perturbation analysis of JD
We analyze here the solution of the basic problem (8). For simplicity of anal-ysis, we assume that n X = n Y = l , the points are ordered as i k = j k = k , and µ → ∞ . This setting of the joint diagonalization problem has been consideredin numerical mathematics [BGBM93] and blind source separation [CS96] liter-ature. In this case, the matrices are of equal size, have row- and column-wisecorrespondence, and a single basis ( ˆΦΦΦ = ˆΨΨΨ) is searched.15igure 9: Shape similarity using joint diagonalization. Darker colors represent moresimilar shapes. One can clearly distinguish blocks of isometric shapes. Also, twoclasses of two- and four-legged shapes (marked with green and blue) are visible. Smallfigures show representative shapes from each class.
We further assume that LLL X = ΦΦΦΛΛΛΦΦΦ T has a simple τ -separated spectrum(i.e., | λ i − λ j | ≥ τ ). If Y is a near-isometric deformation of X , its Laplacian canbe described as a perturbation LLL Y = ΦΦΦΛΛΛΦΦΦ T + (cid:15) RRR of LLL X . Ignoring permutationof eigenfunctions and sign flips, the joint approximate eigenbasis can be writtenas the first-order perturbation [Car95]ˆ φφφ i ≈ φφφ i + (cid:15) (cid:88) j (cid:54) = i α ij φφφ j , (22)where α ij = φφφ T i RRR φφφ j / λ j − λ i ). We can bound these coefficients by the spectralnorm | α ij | ≤ (cid:107) RRR (cid:107) / τ = α max . We can conclude that the first k joint approxi-mate eigenvectors ˆ φφφ , . . . , ˆ φφφ k can be well represented as linear combinations of φφφ , . . . , φφφ k (cid:48) , with square error bounded by (cid:15) (cid:80) j>k (cid:48) | α ij | ≤ (cid:15) ( n X − k (cid:48) − α max .This result also justifies the use of band-wise computation discussed in Section4. To relate this bound to the geometry of the shape, let us assume that X and Y have the same connectivity and that the angles β, ¯ β of triangles in the16wo meshes satisfy θ ≤ β, ¯ β ≤ π − θ (all triangles are at least θ fat), the areaelements are at least s, ¯ s ≥ s , and each vertex is connected to at most ν vertices.We assume that the mesh Y is obtained as a deformation of mesh Y changingthe angles by | ¯ β − β | ≤ ∆ θ and the area elements by 1 − δ ≤ ¯ s/s ≤ δ . We areinterested in a bound on the spectral norm (cid:107) LLL X − LLL Y (cid:107) = (cid:15) (cid:107) RRR (cid:107) expressed interms of parameters ∆ θ, δ (strength of non-isometric deformation) and constants θ , s . From trigonometric identities, we have for i (cid:54) = j | ¯ w ij − w ij | ≤ (cid:12)(cid:12)(cid:12)(cid:12) sin( ¯ β − β )sin β sin ¯ β (cid:12)(cid:12)(cid:12)(cid:12) ≤ sin ∆ θ sin θ ≤ ∆ θ sin θ , (23)and | ¯ w ii − w ii | ≤ ν ∆ θ sin θ Applying the triangle inequality, for i (cid:54) = j we get | ¯ l ij − l ij | = 3 | ¯ w ij ¯ s − i − w ij s − i | ≤ | ¯ w ij − w ij | s − i + 3 | ¯ w ij || ¯ s − i − s − i | . Pluggingin the bounds on | ¯ w ij − w ij | and ¯ s i , we get | ¯ l ij − l ij | ≤ s (cid:18) ∆ θ sin θ + cot( θ ) δ (cid:19) ≤ s (cid:18) ∆ θ sin θ + δθ (cid:19) ≤ θ + δ ) s sin θ , from which it follows by norm inequality (cid:107) LLL X − LLL Y (cid:107) ≤ n / X (cid:107) LLL X − LLL Y (cid:107) ≤ νn / x s sin θ (∆ θ + δ ) , where (cid:15) = ∆ θ + δ is the degree of shape deformation (“lack of isometry”). Appendix B - Off-diagonality penalty gradient
Let us rewrite the penalty as off(AAA) = (cid:107)
AAA T ¯ΛΛΛ X AAA (cid:107) − (cid:80) i (AAA T ¯ΛΛΛ X AAA) ii . The firstterm is bi-quadratic and its gradient derivation is trivial, so here we derive thegradient of the second term only. We have (cid:80) i (AAA T ¯ΛΛΛ X AAA) ii = (cid:80) i (cid:0)(cid:80) k a ki λ k (cid:1) ;differentiating in a coordinate-wise manner, ∂∂a pq (cid:88) i (cid:32)(cid:88) k a ki λ k (cid:33) = 2 (cid:88) i (cid:32)(cid:88) k a ki λ k (cid:33) a pq λ p δ iq = 4 (cid:88) k a kq λ k a pq λ p . (24)Observing that AAAΛΛΛ = ( a pq λ p ) and denoting by o pq = ( (cid:80) k a kq λ k ) the elementsof the equal-columns matrix OOO, we get the expression 4OOO ◦ AAAΛΛΛ.
References [ASC11] M. Aubry, U. Schlickewei, and D. Cremers. The wave kernel signa-ture: a quantum mechanical approach to shape analysis. In
Proc.Workshop on Dynamic Shape Capture and Analysis , 2011.17BBK08] A. M. Bronstein, M. M. Bronstein, and R. Kimmel.
Numericalgeometry of non-rigid shapes . Springer-Verlag New York Inc, 2008.[BBK +
10] A. M. Bronstein, M. M. Bronstein, R. Kimmel, M. Mahmoudi, andG. Sapiro. A Gromov-Hausdorff framework with diffusion geometryfor topologically-robust non-rigid shape matching. 89(2-3):266–286,2010.[BGBM93] A. Bunse-Gerstner, R. Byers, and V. Mehrmann. Numerical meth-ods for simultaneous diagonalization.
SIAM J. Matrix Anal. Appl. ,14(4):927–949, 1993.[BN03] M. Belkin and P. Niyogi. Laplacian eigenmaps for dimensionalityreduction and data representation.
Neural Computation , 13:1373–1396, 2003.[Car95] J. F. Cardoso. Perturbation of joint diagonalizers. 1995.[CL06] R. R. Coifman and S. Lafon. Diffusion maps.
Applied and Compu-tational Harmonic Analysis , 21:5–30, July 2006.[CS93] J.-F. Cardoso and A. Souloumiac. Blind beamforming for non-gaussian signals.
Radar and Signal Processing , 140(6):362–370,1993.[CS96] J.-F. Cardoso and A. Souloumiac. Jacobi angles for simultaneousdiagonalization.
SIAM J. Matrix Anal. Appl , 17:161–164, 1996.[DK10] A. Dubrovina and R. Kimmel. Matching shapes by eigendecompo-sition of the laplace-beltrami operator. 2010.[FH05] M.S. Floater and K. Hormann. Surface parameterization: a tutorialand survey.
Advances in Multiresolution for Geometric Modelling ,1, 2005.[KG00] Z. Karni and C. Gotsman. Spectral compression of mesh geome-try. In
Proc. Conf. Computer Graphics and Interactive Techniques ,pages 279–286, 2000.[Kob97] L. Kobbelt. Discrete fairing. In
In Proceedings of the Seventh IMAConference on the Mathematics of Surfaces , pages 101–131, 1997.[KR05] B. Moon Kim and J. Rossignac. Geofilter: Geometric selectionof mesh filter parameters.
Comput. Graph. Forum , 24(3):295–302,2005.[KRDI11] A. Kumar, P. Rai, and H. Daum´e III. Co-regularized multi-viewspectral clustering. In
Proc. NIPS , 2011.18LBB11] R. Litman, A. M. Bronstein, and M. M. Bronstein. Diffusion-geometric maximally stable component detection in deformableshapes. CG , 35(3):549 – 560, 2011.[L´ev06] B. L´evy. Laplace-Beltrami eigenfunctions towards an algorithmthat “understands” geometry. In Proc. SMA , 2006.[LZ09] B. L´evy and R. H. Zhang. Spectral geometry processing. In
ACMSIGGRAPH ASIA Course Notes , 2009.[MDSB03] M. Meyer, M. Desbrun, P. Schroder, and A. H. Barr. Discretedifferential-geometry operators for triangulated 2-manifolds.
Visu-alization and Mathematics III , pages 35–57, 2003.[NISA06] A. Nealen, T. Igarashi, O. Sorkine, and M. Alexa. Laplacian meshoptimization. In
Proc. Conf. Computer Graphics and InteractiveTechniques , pages 381–389, 2006.[OBCS +
12] M. Ovsjanikov, M. Ben-Chen, J. Solomon, A. Butscher, andL. Guibas. Functional maps: A flexible representation of mapsbetween shapes.
TOG , 31(4), 2012.[OMMG10] M. Ovsjanikov, Q. M´erigot, F. M´emoli, and L. Guibas. One pointisometric matching with the heat kernel. In
Computer GraphicsForum , volume 29, pages 1555–1564, 2010.[PP93] U. Pinkall and K. Polthier. Computing discrete minimal surfacesand their conjugates.
Experimental mathematics , 2(1):15–36, 1993.[RCG08] G. Rong, Y. Cao, and X. Guo. Spectral mesh deformation.
TheVisual Computer , 24(7):787–796, 2008.[Reu10] M. Reuter. Hierarchical shape segmentation and registration viatopological features of Laplace-Beltrami eigenfunctions.
IJCV ,89(2):287–308, 2010.[Rus07] R. M. Rustamov. Laplace-beltrami eigenfunctions for deformationinavriant shape representation. In
Proc. of SGP , pages 225–233,2007.[RWP05] M. Reuter, F. E. Wolter, and N. Peinecke. Laplace-spectra as finger-prints for shape matching. In
Proc. ACM Symp. Solid and PhysicalModeling , pages 101–106, 2005.[RWP06] M. Reuter, F.-E. Wolter, and N. Peinecke. Laplace-Beltrami spec-tra as “shape-DNA” of surfaces and solids.
Computer Aided Design ,38:342–366, 2006.[Sch66] P.H. Sch¨onemann. A generalized solution of the orthogonal pro-crustes problem.
Psychometrika , 31(1):1–10, 1966.19SMKF04] P. Shilane, P. Min, M. Kazhdan, and T. Funkhouser. The princetonshape benchmark. In
Proc. SMI , 2004.[SOG09] J. Sun, M. Ovsjanikov, and L. J. Guibas. A concise and provablyinformative multi-scale signature based on heat diffusion. In
Proc.SGP , 2009.[SP04] R.W. Sumner and J. Popovi´c. Deformation transfer for trianglemeshes. In
TOG , volume 23, pages 399–405, 2004.[Tau95] G. Taubin.
A signal processing approach to fair surface design .ACM, 1995.[WBH +
07] M. Wardetzky, M. Bergou, D. Harmon, D. Zorin, and E. Grinspun.Discrete quadratic curvature energies.
CAGD , 24:499–518, 2007.[WMKG08] M. Wardetzky, S. Mathur, F. K¨alberer, and E. Grinspun. DiscreteLaplace operators: no free lunch. In
Conf. Computer Graphics andInteractive Techniques , 2008.[Yer02] A. Yeredor. Non-orthogonal joint diagonalization in the least-squares sense with application in blind source separation.
Trans.Signal Proc. , 50(7):1545 –1553, 2002.[Zie05] A. Ziehe.