1 Introduction
The field of shape analysis has been evolving rapidly during the last decades. The constant increase in computing power allowed image and shape understanding algorithms to efficiently handle difficult problems that could not have been practically addressed in the past. A large set of theoretical tools from metric and differential geometry, and spectral analysis has been imported and translated into action within the shape understanding arena. Among the numerous ways of analyzing shapes, a common one is to embed them into a different space where they can be processed more efficiently.
1.1 Related efforts
[Elad and Kimmel (2003)] introduced a method for analyzing surfaces based on embedding the intrinsic geometry of a given shape into a Euclidean space, extending previous efforts of [Schwartz et al. (1989), Zigelman et al. (2002), Grossmann et al. (2002)]. Their key idea was to consider a shape as a metric space, whose metric structure is defined by geodesic distances between pairs of points on the shape. Two nonrigid shapes are compared by first having their respective geometric structures mapped into a lowdimensional Euclidean space using multidimensional scaling (MDS) [Cox and Cox (2008)], and then comparing rigidly the resulting images, also called canonical forms.
[Mémoli and Sapiro (2005)] proposed a metric framework for nonrigid shape comparison based on the GromovHausdorff distance that was suggested by Gromov as a theoretical tool to quantify disimilarity between metric spaces. Using the GromovHausdorff formalism, the distance between two shapes is defined by matching pairwise distances on the shapes. However, the GromovHausdorff distance is difficult to compute when treated in a straightforward manner. To overcome this difficulty [Bronstein et al. (2006a), Bronstein et al. (2006b)] proposed an efficient numerical solver based on a continuous optimization problem, known as Generalized MDS (GMDS). Recently, other relaxation schemes have been proposed, see for example [Chen and Koltun (2015), Aflalo et al. (2016)].
In the past decade, the LaplaceBeltrami operator (LBO) – the extension of the Laplacian to nonEuclidean manifolds, has become growingly popular. Its properties have been well studied in differential geometry and it was used extensively in computer graphics. The LBO can be found in countless applications such as mesh filtering [Vallet and Levy (2008)], mesh compression [Karni and Gotsman (2000)], shape retrieval [Bronstein et al. (2011)], to name just a few. It has been widely used in shape matching where several approaches treat the correspondence problem by comparing isometric invariant pointwise descriptors between the two shapes. For example, the Global Point Signature (GPS) [Rustamov (2007)], the Heat Kernel Signature (HKS) [Sun et al. (2009)] and the Wave Kernel Signature (WKS) [Aubry et al. (2011)]
, all use the eigenfunctions and eigenvalues of the LBO to compute local shape descriptors. Matching only signatures at a small set of points, the correspondence between the points on the two shapes can be found. These points can serve as anchors and interpolated for the entire shape
[Ovsjanikov et al. (2010)] where refinement of the basis can be performed to produce precise dense correspondence [Ovsjanikov et al. (2012), Pokrass et al. (2013), Shtern and Kimmel (2014a)].Recently, learning based approaches [Litman and Bronstein (2014), Wei et al. (2016), Boscaini et al. (2016)] have also become highly popular in the shape matching arena.
The use of the basis defined by the LBO is in many senses a natural choice for surfaces analysis. It was chosen in the functional map framework [Ovsjanikov et al. (2012)] because of its compactness, stability, and invariance to isometries. Subsequently, it was proven to be optimal [Aflalo et al. (2015)] for representing smooth functions on the surface. In an attempt to overcome the topological sensitivity of the LBO and the nonlocal support of its eigenfunctions, compressed eigenfunctions have been adapted from mathematical physics to shape analysis [Neumann et al. (2014), Bronstein et al. (2016)]. Here, we try to find a richer family of basis functions that are based on intrinsic properties that can go beyond the geometry of the shape. Exploring a similar goal, [Kovnatsky et al. (2011)] combined geometric and photometric information within a unified metric for shape retrieval. [Iglesias and Kimmel (2012)] used artificial surface textures on shapes to define elliptic operators that give birth to a new family of diffusion distances. Along the same line of thought, [Hildebrandt et al. (2012)] designed a new family of eigenvibrations using extrinsic curvatures and deformation energies.
We suggest to further explore those ideas and construct the socalled potential operator that is added to the Laplace Beltrami operator. Here, a designed perturbation to the Laplacian permits a supervised control of the vibrational modes on the manifold.
1.2 Contributions
The main contribution of this paper is the exploration of the Hamiltonian operator on manifolds. We study spectral properties of the operator and the impact of an additional potential function to the Laplacian for shape analysis applications. The properties of the Hamiltonian allow it to be efficiently utilized by many spectralbased methods. The potential part can lead to a more descriptive operator when treated as a truncated basis generator. Modulated harmonics on the surface are obtained by treating different regions of interest as different values of the potential. We show that by simply plugging the resulting basis into existing spectral shape analysis pipelines could improve their performance.
The rest of the paper is organized as follows: in Section 2
, we propose to study the Hamiltonian on manifolds from the variational calculus point of view with motivation from quantum mechanics. We prove optimality properties of its eigenspace, characterize the associated diffusion process, the resulting nodal sets, introduce a discretization method and analyze the robustness of the operator.
In Section 3, we propose a global optimization framework for variational problems involving the basis defined by the Hamiltonian. We provide an approach for computing derivatives with respect to the potential based on eigenvectors perturbation theory. We demonstrate the effectiveness of the framework on the task of data representation.
In Section 4, we review recent improvement of the computation of the compressed modes [Ozolins et al. (2013)] that make use of the decomposition of the Hamiltonian. Finally, in Section 5 we present properties of the proposed basis that make it a better alternative for the task of shape matching where priors can be inserted through the potential in order to improve performance.
2 Hamiltonian operator
2.1 The Laplace Beltrami Operator
Consider a parametrized surface
with a metric tensor
. The space of squareintegrable functions on is denoted by with the standard inner product , where is the area element induced by the Riemannian metric . The Laplace Beltrami Operator acting on a scalar function is defined as(1) 
where is the determinant of the metric matrix and is the inverse metric. If
is a domain in the Euclidean plane, the metric matrix is generally the identity matrix and the LBO reduces to the wellknown Laplacian
(2) 
The LBO is selfadjoint and thus admits a spectral decomposition , where and , such that,
(3)  
with the Kronecker delta. In case has boundary, we add homogeneous Neumann boundary condition
(4) 
where
is the normal vector to the boundary
.The LBO eigendecomposition can be extracted from the Euler Lagrange solution to Dirichlet energy minimization
(5)  
s.t. 
Here, each ordered eigenfunction composing the basis on the manifold corresponds to the function with the smallest possible energy that is orthogonal to all the previous ones. Therefore, the LBO eigenfunctions can be seen as an extension of the Fourier harmonics in Euclidean spaces to manifolds and are often referred to as Manifold Harmonics.
2.2 Hamiltonian
A Hamiltonian operator on a manifold acting on a scalar function , is an elliptic operator of the form
(6) 
where is a realvalued scalar function. It plays a fundamental role in the field of quantum mechanics appearing in the famous Schrödinger equation that describes the wave motion of a particle with mass under potential ,
(7) 
where is the Planck’s constant and represents the wave function of the particle such that
is interpreted as the probability distribution of finding the particle at a given position
at time .The Schrödinger equation can be analyzed via perturbation theory by solving the spectral decomposition of the Hamiltonian
(8) 
also known as the timeindependent Schrödinger equation, where is the eigenenergy of a particle at the stationary eigenstate .
Since the potential is a diagonal operator, the Hamiltonian is selfadjoint as a sum of two selfadjoint operators and its eigenfunctions form a complete orthonormal basis on the manifold . As a generalization of the regular Laplacian, its spectral theory can be derived almost straightforwardly from that of the latter. Classical examples of the influence of potential functions in a onedimensional Euclidean domain are depicted in Figure 1.
2.3 Variational principle
Let us consider the following variational problem
(9)  
s.t. 
whose the EulerLagrange equation defines the eigendecomposition of the Hamiltonian defined in (8).
The basis defined by the Hamiltonian operator corresponds to the orthogonal harmonics modulated by the potential function. The potential defines the tradeoff between the orientation and the compactness of the basis and its global support. Larger values of the potential will enforce smooth solutions that concentrate on the low potential regions, while smaller ones will give solutions that better minimize the total energy at the expense of more extended wave functions.
2.4 Finite step potential
The timeindependent Schrödinger equation can yield a rather complicated problem to solve analytically, even in one dimension. Let us consider a system with an ideal step potential in one dimension [Griffiths (2005)]. We need to solve the differential equation , with denoting the energy of the particle, and the Heaviside function with step of magnitude , at point , given by
(10) 
The step divides the space in two constantpotential regions. At the zero potential region, the particle is free to move and the harmonic solutions are known. In the high potential region, on the other hand, for , the solution is a decaying exponentially, meaning that the particle cannot pass the potential barrier and is reflected according to classical physics. If , the solution is also harmonic, which means there is a probability for the particle to penetrate into the effective potential region with a different energy than that of particles in the zero potential region. We illustrate this effect in Figure 2 by numerically computing the eigenvectors of the LBO and the Hamiltonian with a potential defined on a human body surface.
Therefore, the potential energy can be tuned to enforce localization of the basis at the expense of loss of smoothness.
Let and be the spectral decompositions of the Laplacian, and the Hamiltonian, respectively. Then, everywhere on the manifold implies that the eigenvalues satisfy
.
Proof. According to the CourantFischer minmax theorem, we have
(11)  
(13)  
(14) 
Similarly,
(16)  
(17) 
Since the family of eigenvalues of the Helmholtz equation (3) consist of a diverging sequence ( [Weyl
et al. (1950)]), there exists an such that
and the tradeoff between localcompact and global support of the basis elements can be controlled by the potential energy. Then, we can estimate the magnitude of the potential required in order to allow for oscillations outside the regions where the potential vanishes.
Given a scalar we can define the Hamiltonian as
(18) 
where controls the resistance to diffusion induced by the potential. Let and be the th eigenvalue of the LBO and Hamiltonian, respectively. We seek for a constant such that so the particle can penetrate the high potential region. Considering the potential as small perturbation of the Laplacian, up to first order, the eigenenergies are defined as . In order to contain the basis support at most until the th eigenfunction, must satisfy
(19) 
According to its potential energy, the basis can then provide supervised multiresolution analysis on the manifold by containing the first eigenfunctions and allow global analysis for the following.
2.5 Optimality of the Hamiltonian eigenspace
Let us consider a function . We define the representation residual function as
(20)  
Defining , we know that
(21)  
Thus, from (20) and (21) we obtain
(22) 
Recall that for we return to the LBO case. Among the numerous reasons that motivated the selection of the Laplacian for shape analysis, a major one is its efficiency in representing functions with bounded gradient magnitude. This result was subsequently proved to be optimal for representing functions with bounded gradient magnitude over surfaces in [Aflalo et al. (2015)], which says that there exists no other basis with better representation error for all possible functions.
In case of the Hamiltonian, the Dirichlet energy is coupled with the potential energy. Thus the Hamiltonian operator advocates measuring smoothness differently for different regions of the domain where smoothness remains a less important factor than avoiding vibrations in high potential areas. This is a useful property to exploit in different shape analysis scenarios.
Next, we show that the Hamiltonian is optimal in approximating functions with both bounded gradient and low values in high potential areas. Let . There is no integer and no sequence of linearly independent functions in such that
(23) 
The proof of Theorem 2.5 is given in the Appendix.
2.6 Diffusion process
Let us be given a Riemannian manifold . The heat equation governing the diffusion process on is defined as
(24) 
with appropriate boundary conditions. A natural extension to the new operator with a potential , can be written as
(25) 
The solutions of (24) and (25) have the form [Iglesias and Kimmel (2012)]
(26) 
that represents the diffusion in time of heat on the manifold with potential , where . We refer to as the heat kernel. A standard proof is given in the Appendix.
According to the FeynmanKac formula [Simon (2005)], the solution of the diffusion process is expressed in terms of the Wiener process,
(27) 
In the Laplacian case, the initial value is carried over random paths in time, while the expected value of the stochastic process is equal to the solution . For , the diffusion spreads according to the potential on the manifold, when the transported value is modulated exponentially by the potential , diffusing anisotropically to low potential regions, as shown in Figure 3.
2.7 Nodal sets
An interesting property of the Laplacian is the relation between its eigenfunctions, the number of connected nodal (zero) sets, and the number of complementary regions they define. Given an eigenfunction , a nodal set is defined as the set of points at which the eigenfunction values are zero. That is,
(28) 
The Nodal Theorem [Courant and Hilbert (1966)] states that the th eigenfunction of the LBO can split to at most connected subdomains. In other words, the zero set of the th eigenfunction can separate the manifold into at most connected components.
Proposition 1
Given the selfadjoint Hamiltonian operator on , with arbitrary boundary conditions; if its eigenfunctions are ordered according to increasing eigenvalues, then, the nodal set of the th eigenfunction divides the domain into no more than connected subdomains.
The proof is essentially the same as that of the Laplacian case. See [Courant and Hilbert (1966)] Vol.1 Sec. VI.6 for a proof.
As shown in Figure 2, the Hamiltonian eigenfunctions are tuned by the potential. Thus, shape segmentation can be obtained by separating the surface according to the induced nodal sets as described in [Levy (2006)]. Given a potential defined on the surface, meaningful segmentation can be induced by the nodal domains of the resulting eigenfunctions, as presented in Figure 4.
2.8 Discretization
In the discrete setting, we consider a triangular mesh in with the associated space of functions that are continuous and linear in every triangle. According to the Finite Element Method (FEM) [Dziuk (1988)], the solution of the Hamiltonian eigenvalue problem (8) can be computed by imposing that the equation is satisfied in a weak sense, that is,
(29) 
where denote the Lagrange basis of piecewiselinear hatfunctions on , that take the value one at a vertex and vanishes at all the other vertices. represents the eigenenergy of the Hamiltonian. On the mesh, the bilinear form is evaluated by splitting the integrals into a sum over the triangles of by
(30) 
Since the Hamiltonian is a linear operator we have
(31) 
The matrix representation of and with respect to the Lagrange basis are well known [Pinkall and Polthier (1993)] and define the stiffness matrix and the mass matrix with the entries
(32) 
Thus,
(33)  
where the last equality is obtained by representing the potential function as a diagonal matrix according to the Lagrange basis functions. The discretization of the eigenvalue problem (8) is defined by finding all pairs such that
(34) 
Efficient solution methods can be found in [Vallet and Levy (2008)]. Among the possible explicit representations of the matrices and , we use here the cotangent formula [Pinkall and Polthier (1993), Meyer et al. (2003)] where the stiffness matrix is defined as
(35) 
with , where is the set of edges of the triangulated surface interpreted as a graph and denote the angles and of the triangles sharing the edge . The mass matrix is replaced by a diagonal lumped mass matrix of the area of local mixed Voronoi cells about each vertex . The manifold inner product is discretized as . Since only modifies the diagonal of , our operator remains a sparse matrix with the same effective entries, and thus, there is no increase in the computational cost of the generalized eigendecomposition compared to that of the LBO.
2.9 Robustness to noise
As a generalization of the Laplacian, the Hamiltonian exhibits similar robustness to noise. Consider the Hamiltonian matrix with the potential. Then, the perturbed Hamiltonian has the form . Let us define and . Based on perturbation theory, and up to secondorder corrections, the th eigenfunction of has the form
(36) 
with and being respectively the
th eigenfunction and eigenvalue of the unperturbed Hamiltonian. Assuming uniformly distributed random noise on the mesh, the eigenfunctions of the regular Laplacian may present smaller distortion to noise than the Hamiltonian since the perturbation is amplified by area and potential distortions. Still, in case of potential with small values the distortion is insignificant. In Figure
5, we present the original surface and its noisy version in which vertex positions have been corrupted by additive Gaussian noise with of the mean edge length. The potential is also modified by adding a Gaussian noise withof the initial variance of the potential.
The construction of the Laplacian depends crucially on the mesh connectivity making it sensitive to topological noise such as holes and part removal that can be found in many depth acquisition scenarios. The compact support of the basis elements of the Hamiltonian makes it robust to noise compared to the basis elements that are generated by the Laplacian. We illustrate the robustness property in Figure 5 where of the surface area was removed due to topological noise in the form of small holes.
3 Optimization of the potential
One natural problem emerging when working with the Hamiltonian is the ability to define an optimal potential function for a specific task. The choice of the potential is application dependent but can be represented through minimization problem generically defined as
(37)  
s.t. 
where denotes the data term depending on the data matrix and the vector defining the diagonal potential matrix. Regularization terms can be further be added. If the analytical solution remains complex, a common approach is to minimize the goal function with an optimization algorithm involving the gradient of the goal function with respect to the potential. In this section we propose an optimization framework based on perturbation theory of the eigenvectors where optimal potential is obtained. To that end, we need to derive the gradient for a given objective .
Here we will consider the problem of data representation using the discrete basis of the Hamiltonian referred to as representing the eigenvectors of the Hamiltonian such that . The discretized minimization problem is defined as
(38)  
s.t. 
with . The objective defines the representation error of the data in the subspace spanned by the columns of and in the sense of the Frobenius norm . For a general orthonormal matrix
, the problem is equivalent to Principal Component Analysis (PCA). We can straightforwardly obtain that
(39)  
Thus, the differential of the loss function
with respect to is obtained by(40)  
It remains to derive the differential of the Hamiltonian eigenvectors. Let us consider the full matrix of eigenvectors , the diagonal matrix of eigenenergies and the discrete Hamiltonian operator . The eigenvalue decomposition problem is given by . Thus, the differential of the spectral decomposition problem is given by
(41) 
Multiplying by on the left side and denoting with , we have
(42)  
since . We readily obtain that the off diagonal elements of the matrix can be defined by
(43) 
Here represents the th column of the matrix of eigenvectors. The diagonal elements of are defined by the following
(44)  
The diagonal elements are then defined by . Since second order elements are negligible, we have . We obtain that , with denoting the Hadamard product and the matrix defined as
(45) 
The selection of the first eigenvectors are obtained by multiplying by the truncated identity matrix . The differential is now known and can be plugged into (40) in order to extract , that is:
(46)  
The passage in the fourth line stems from the equivalence . Since , we obtain finally
(47) 
Two problems arise from the suggested scheme. First, the high computational cost of a full (sparse) matrix diagonalization. Second, the matrix remains undefined when eigenvectors have nontrivial multiplicities. The first problem can be relaxed by approximating the matrix with less eigenvectors. This is especially justified for distant indices, where the eigenenergies are well separated and the corresponding elements of matrix become negligible. Also, the data can be projected onto the LBO basis so the solution complexity remains constant with the size of the mesh. Even if the second problem has been treated in [Van Der Aa et al. (2007)], it seems that lack of smoothness at isolated points is not critical for computation and convergence by resorting to a subgradient approach. The alternative opted for here is to stabilize the matrix in order to avoid exploding gradients. We use the approximation
(48) 
where the sign function is not vanishing.
In its geometric setting, the solution remains similar and is obtained by defining a new basis , such that , coupled with the consistently discretized Frobenius norm . In the following experiments we allowed negative potential for performance consideration only, since the potential is defined over the whole codomain . Also, for physically interpretable solutions we enforced positive potential by using quadratic function . The extension of the derivation is straightforward but decreased the performance since it is more restrictive.
As a toy experiment, we propose to find the best potential for the representation of a function in the one dimensional Euclidean domain. Given a function , we seek for the best potential minimizing . We compare in Figure 6 the reconstruction performance on a one dimensional linear function with the Laplacian and the Hamiltonian built from the optimized potential.
In Figure 7, we propose to reconstruct the matrix of coordinates of a mesh so the data matrix is defined by . The experiments were implemented using the quasiNewton method with initial zero potential.
An important application related data representation is spectral mesh compression. [Karni and Gotsman (2000)] proposed to project the coordinates functions of the mesh onto the LBO eigenfunctions in order to encode the mesh geometry via the first coefficients only. Since most of the function energy is generally contained in the first coefficients, the reconstruction distortion is low, up to fine details related to higher frequencies. Since matrix decomposition is an expensive operation, they suggested to segment the shape into smaller parts that can be processed separately. By sending the mesh topology (triangles) separately, the combinatorial graph Laplacian is built on the decoder side and the signal can be reconstructed with the received coefficients. We suggest to apply this idea to our basis which potential is obtained by the proposed optimization framework. However, one major drawback is that we need to encode the potential as well as the coefficients. Also, some methods use the ordering of the vertices in order to encode information [Touma and Gotsman (1998)]. Here we suggest to reorder the vertices such that the vertex with the smallest potential is be assigned the index and the vertex with the largest potential is be assigned the index . By using a fixed potential defined as
(49) 
the decoder simply applies in order to obtain the Hamiltonian basis. Here and are the regression coefficients minimizing that are also encoded. We present compression results in Figure 8.
4 Compressed Manifold Modes
[Ozolins et al. (2013)] proposed a a novel method to create a set of localized eigenfunctions in Euclidean domains. To that end, they modified the construction of standard differential operators by adding an regularization term to the variational leading to the decomposition of the operator. The resulting eigenfunctions were called compressed modes and were shown to be compactly supported [Brezis (1974)]. [Neumann et al. (2014)] extended this construction to manifolds, suggesting the following discrete regularization problem
(50)  
s.t. 
with the parameter that controlling the localization of the basis. Proposed solutions require the use of expensive and unstable optimization techniques ([Neumann et al. (2014), Kovnatsky et al. (2016)]) based on ADMM and proximal operators.
The latter optimization problem (50) can be written as an Hamiltonian eigendecomposition problem [Bronstein et al. (2016)]
(51)  
s.t. 
where is the diagonal matrix operator defining the potential that corresponds to the th eigenvector that localizes the support of in lowpotential areas. The potential is defined iteratively using a reweighted least squares scheme
(52) 
ensuring that the minimizers of (50) and (51) coincide. Interestingly, the potential here is defined as a function of the eigenfunction, namely . The potential and the resulting eigenstate are then intrinsically linked, meaning that the potential is influenced by the state of the particle itself. Consequently, a perturbation of the potential enforces perturbation of the eigenfunction and vice versa until reaching steady state.
We formulate the compressed manifold modes problem as
(53)  
s.t. 
with and where is a sufficiently large constant such that the third term guarantees that the th mode is orthogonal to the previously computed modes , . Observe that albeit nonconvex, the problem has a closed form global solution, that is the smallest generalized eigenvector of
(54) 
with
For small number of compressed modes, is a low rank matrix and finding the smallest generalized eigenvector can be solved efficiently since the involved matrix is the sum of a sparse and a lowrank matrix.
Several numerical eigendecomposition implementations use the Arnoldi iteration algorithm. In our matrix decomposition problem, the core operation is the multiplication by the inverse of the matrix with a vector, operation that cannot be solved straightforwardly. Also, shifting the maximum eigenvalue using power method is too unstable since it depend on gap of the first eigenvalues, generally tight. In our configuration the Woodbury identity
can be used to compute efficiently the vector multiplication with the inverse of the matrix as a cascade of sparse and lowrank systems. Unlike solutions of the inconsistently discretized problem 50, the basis obtained with the proposed Hamiltonian framework is more robust under various discretizations and can be computed at a fraction of the computational cost (Figure 9).
5 Shape matching
The task of matching pairs of shapes lies at the core of many shape analysis tasks and plays a central role in operations such as 3D alignment and shape reconstruction. While rigid shape matching has been well studied in the literature, nonrigid correspondence remains a difficult task even for nearly isometric surfaces. When dealing with rigid objects, it is sufficient to find the rotation and translation that aligns one shape to the other [Chen and Medioni (1992)]
. Therefore, the rigid matching problem amounts to determining only six degrees of freedom. At the other end, nonrigid matching generally requires dealing with many more degrees of freedom. Since the LBO is invariant to isometric deformations, it has been used extensively to aid the solution of correspondence problem. Several properties of the Hamiltonian operator make it a better choice for this task compared to its zeropotential particular case that is the LBO.
Invariance. The LaplaceBeltrami Operator is defined in terms of the metric tensor which is invariant to isometries. For a potential function defined intrinsically, the resulting Hamiltonian is also isometryinvariant.
Compactness. Compactness means that scalar functions on a shape should be well approximated by using only a small number of basis elements. From Theorem 2.5 and as a generalization of the Laplacian, the global support and compactness hold for a bounded (low) potential.
Descriptiveness. The LBO eigenvalues are related to frequency. Similarly, eigenenergies of the Hamiltonian relate to the number of oscillations on the manifold. Theorem 2.4 demonstrates that the modes corresponding to small eigenvalues of the Hamiltonian defined with a positive potential, encapsulate higher frequencies, even when localized, compared to the modes of the regular LBO. At the other end, highly oscillating eigenfunctions can be used to represent fine details of the shape that can be crucial for shape matching. Also, the potential enforces different oscillations in different regions on the manifold, allowing for better discrimination of similar areas and disambiguation of intrinsic symmetries with asymmetric potential.
Stability. Deformations of nonrigid shapes and articulated objects can stretch the surface. In such cases, the LBO eigendecomposition of the two shapes will be different. We could compensate for such local metric distortions by carefully designing a potential. Assigning high potential to strongly distorted regions would lead to lower values of the eigenfunctions in those areas (9). Such a potential will reduce the discrepancy between corresponding eigenfunctions at least for the lower eigenergies, as shown via the functional maps representations [Ovsjanikov et al. (2012)] in Figure 10. Let define and , the area at vertex on mesh and on the second mesh respectively and a bijection between two (discretized) surfaces and . Then, we define the potential at vertex as
(55) 
Among the few stable intrinsic invariants that can be extracted from the geometry, we will use the stable first eigenfunctions of the LBO and geodesic distances. Additional non necessarily intrinsic information such as photometric properties or even extrinsic shape properties such as principal curvatures [Hildebrandt et al. (2012)] can also be integrated into the potential field.
5.1 Experimental Evaluation
We tested the proposed basis and compared its matching performances to that the LBO basis as applied to pairs of triangulated meshes of shapes from the TOSCA dataset [Bronstein et al. (2008)] and the SCAPE dataset [Anguelov et al. (2004)]. The TOSCA data set contains densely sampled synthetic human and animal surfaces, divided into several classes with given groundtruth pointtopoint correspondences between the shapes within each class. The SCAPE data set contains scans of real human bodies in different poses. The evaluation method used is described in [Kim et al. (2011)] where the distortion curves describe the percentage of surface points falling within a relative geodesic distance from what is assumed to be their true locations. Symmetries were not allowed in all evaluations. Note that we assume that the sign ambiguity of the first eigenfunctions generating the potential is resolved [Shtern and Kimmel (2014b)].
Figure 11 compares the two operators by matching diffusion kernel descriptors derived from the corresponding eigenfunctions. The diffusion on the shape using the Hamiltonian as the diffusion operator is more descriptive than regular diffusion that cannot resolve the symmetries. Also, it would be natural to compute the WKS signature when the Schrödinger equation is governed by a given effective potential. As intrinsic positive potential we use the normalized sum of the four first nontrivial eigenfunctions of the LBO on each shape, adding a constant of minimal value in order to obtain a nonnegative potential. This way only the intrinsic unstable geometry of the shape is involved in defining the Hamiltonian operator.
In case we know which regions are prone to elastic distortions, like joints and stretchable skin in articulated objects, we could suppress the effect of those regions in our matching procedures by using an appropriate potential as a selective mask. Figure 12, compares the operator with and without potential by matching the spectral signatures computed by the framework of [Shtern and Kimmel (2015)]. The potential we used is the local area distortion when comparing the meshes of two corresponding objects, as in (55). The descriptiveness of the potential and the localization of the harmonics lead to more accurate matching results.
To investigate the performances of the Hamiltonian with photometric textures used as potential, we present in Figure 13 the results of different signatures matching with a dalmatian texture defined for the ”Dogs” shapes from the TOSCA data set.
Iterative refinement of functional representations have been proven to be powerful in shape matching [Ovsjanikov et al. (2012)]. Given an initial partial or dense map, it tries to recover iteratively dense and accurate matching between two given shapes. Here we use a similar refinement framework dubbed as Iterative Closest Spectral Kernel Maps (ICSKM) [Shtern and Kimmel (2014a)] for performance comparison between the two bases. Figure 14 compares the regular ICSKM algorithm working with the Laplacian eigenspace and the Hamiltonian method when we provided one, two, or three landmark points, that were randomly selected from the groundtruth mapping. The potential used in these examples is the geodesic distance from the landmark points. Note that again we use only the geometry of the shapes in order to refine the match between them using the new basis.
[scale=0.3]initPointsGraphs
\put(13.0,120.0){\scriptsize{{ }}} 
6 Conclusion
A classical operator was adopted from the field of quantum mechanics and adapted to shape analysis problems. Functional and spectral properties of the Hamiltonian operator were presented and compared to the popular Laplacian operator often used in many shape analysis procedures. A general optimization method for solving variational problems involving the Hamiltonian operator have been proposed and employed to the task of mesh compression. Features and texture properties can be incorporated into the new operator to obtain a descriptive and stable basis that provides a powerful domain of operation for shape matching. Various directions for future research include exploration of the operator on other shape analysis tasks such as partial shape matching where occluded areas could be refined via the potential. This work has been supported by Grant agreement No. 267414 of the European Community’s FP7ERC program.
References
 Aflalo et al. (2015) Aflalo, Y., Brezis, H., and Kimmel, R. 2015. On the optimality of shape and data representation in the spectral domain. SIAM J. Imaging Sciences 8, 2, 1141–1160.

Aflalo
et al. (2016)
Aflalo, Y., Dubrovina, A., and Kimmel, R. 2016.
Spectral generalized multidimensional scaling.
International Journal of Computer Vision (IJCV)
118, 3, 380–392.  Anguelov et al. (2004) Anguelov, D., Srinivasan, P., Pang, H.C., Koller, D., Thrun, S., and Davis, J. 2004. The correlated correspondence algorithm for unsupervised registration of nonrigid surfaces. In Proceedings of the 17th International Conference on Neural Information Processing Systems. NIPS’04. MIT Press, Cambridge, MA, USA, 33–40.
 Aubry et al. (2011) Aubry, M., Schlickewei, U., and Cremers, D. 2011. The wave kernel signature: A quantum mechanical approach to shape analysis. In Computer Vision Workshops (ICCV Workshops), 2011 IEEE International Conference on. 1626–1633.
 Boscaini et al. (2016) Boscaini, D., Masci, J., Rodolà, E., and Bronstein, M. M. 2016. Learning shape correspondence with anisotropic convolutional neural networks. Tech. Rep. arXiv:1605.06437.
 Brezis (1974) Brezis, H. 1974. Solutions with compact support of variational inequalities. Russian Mathematical Surveys 29, 2, 103–108.

Brezis (2010)
Brezis, H. 2010.
Functional analysis, Sobolev spaces and partial differential equations
. Springer Science & Business Media.  Bronstein et al. (2008) Bronstein, A., Bronstein, M., and Kimmel, R. 2008. Numerical geometry of nonrigid shapes, 1 ed. Springer Publishing Company, Incorporated.
 Bronstein et al. (2016) Bronstein, A., Choukroun, Y., Kimmel, R., and Sela, M. 2016. ”Consistent discretization and minimization of the L1 norm on manifolds”. In 3D Vision (3DV), 2016 4nd International Conference on. IEEE.
 Bronstein et al. (2011) Bronstein, A. M., Bronstein, M. M., Guibas, L. J., and Ovsjanikov, M. 2011. Shape google: Geometric words and expressions for invariant shape retrieval. ACM Trans. Graph. 30, 1 (Feb.), 1:1–1:20.
 Bronstein et al. (2006a) Bronstein, A. M., Bronstein, M. M., and Kimmel, R. 2006a. Efficient computation of isometryinvariant distances between surfaces. SIAM J. Scientific Computing 28, 5, 1812–1836.
 Bronstein et al. (2006b) Bronstein, A. M., Bronstein, M. M., and Kimmel, R. 2006b. Generalized multidimensional scaling: A framework for isometryinvariant partial surface matching. Proceedings of the National Academy of Sciences 103, 5, 1168–1172.
 Chen and Koltun (2015) Chen, Q. and Koltun, V. 2015. Robust nonrigid registration by convex optimization. In 2015 IEEE International Conference on Computer Vision (ICCV). 2039–2047.
 Chen and Medioni (1992) Chen, Y. and Medioni, G. 1992. Object modelling by registration of multiple range images. Image Vision Comput. 10, 3 (Apr.), 145–155.
 Courant and Hilbert (1966) Courant, R. and Hilbert, D. 1966. Methods of mathematical physics. Vol. 1. CUP Archive.
 Cox and Cox (2008) Cox, A. M. A. and Cox, F. T. 2008. Multidimensional Scaling. Springer Berlin Heidelberg, Berlin, Heidelberg, 315–347.
 Dziuk (1988) Dziuk, G. 1988. Finite elements for the beltrami operator on arbitrary surfaces. In Partial differential equations and calculus of variations. Springer, 142–155.
 Elad and Kimmel (2003) Elad, A. and Kimmel, R. 2003. On bending invariant signatures for surfaces. IEEE Trans. Pattern Anal. Mach. Intell. 25, 10 (Oct.), 1285–1295.
 Griffiths (2005) Griffiths, D. 2005. Introduction to Quantum Mechanics. Pearson international edition. Pearson Prentice Hall.
 Grossmann et al. (2002) Grossmann, R., Kiryati, N., and Kimmel, R. 2002. Computational surface flattening: a voxelbased approach. IEEE Transactions on Pattern Analysis and Machine Intelligence 24, 4 (Apr), 433–441.
 Hildebrandt et al. (2012) Hildebrandt, K., Schulz, C., von Tycowicz, C., and Polthier, K. 2012. Modal shape analysis beyond laplacian. Computer Aided Geometric Design 29, 5, 204–218.
 Iglesias and Kimmel (2012) Iglesias, J. A. and Kimmel, R. 2012. Schrödinger Diffusion for Shape Analysis with Texture. Springer Berlin Heidelberg, Berlin, Heidelberg, 123–132.
 Karni and Gotsman (2000) Karni, Z. and Gotsman, C. 2000. Spectral compression of mesh geometry. In Proceedings of the 27th Annual Conference on Computer Graphics and Interactive Techniques. SIGGRAPH ’00. ACM Press/AddisonWesley Publishing Co., New York, NY, USA, 279–286.
 Kim et al. (2011) Kim, V. G., Lipman, Y., and Funkhouser, T. 2011. Blended intrinsic maps. ACM Trans. Graph. 30, 4 (July), 79:1–79:12.
 Kovnatsky et al. (2011) Kovnatsky, A., Bronstein, M. M., Bronstein, A. M., and Kimmel, R. 2011. Photometric heat kernel signatures. In Scale Space and Variational Methods in Computer Vision  Third International Conference, SSVM 2011, EinGedi, Israel, May 29  June 2, 2011, Revised Selected Papers. 616–627.
 Kovnatsky et al. (2016) Kovnatsky, A., Glashoff, K., and Bronstein, M. M. 2016. MADMM: a generic algorithm for nonsmooth optimization on manifolds. In European Conference on Computer Vision. Springer, 680–696.
 Levy (2006) Levy, B. 2006. Laplacebeltrami eigenfunctions towards an algorithm that ”understands” geometry. In Proceedings of the IEEE International Conference on Shape Modeling and Applications 2006. SMI ’06. IEEE Computer Society, Washington, DC, USA, 13–.
 Litman and Bronstein (2014) Litman, R. and Bronstein, A. M. 2014. Learning spectral descriptors for deformable shape correspondence. IEEE transactions on pattern analysis and machine intelligence 36, 1, 171–180.
 Mémoli and Sapiro (2005) Mémoli, F. and Sapiro, G. 2005. A theoretical and computational framework for isometry invariant recognition of point cloud data. Foundations of Computational Mathematics 5, 3, 313–347.
 Meyer et al. (2003) Meyer, M., Desbrun, M., Schröder, P., and Barr, A. H. 2003. Discrete differentialgeometry operators for triangulated 2manifolds. Springer Berlin Heidelberg, Berlin, Heidelberg, 35–57.
 Neumann et al. (2014) Neumann, T., Varanasi, K., Theobalt, C., Magnor, M., and Wacker, M. 2014. Compressed manifold modes for mesh processing. Comput. Graph. Forum 33, 5 (Aug.), 35–44.
 Ovsjanikov et al. (2012) Ovsjanikov, M., BenChen, M., Solomon, J., Butscher, A., and Guibas, L. 2012. Functional maps: A flexible representation of maps between shapes. ACM Trans. Graph. 31, 4 (July), 30:1–30:11.
 Ovsjanikov et al. (2010) Ovsjanikov, M., Mérigot, Q., Mémoli, F., and Guibas, L. 2010. One point isometric matching with the heat kernel. Computer Graphics Forum 29, 5, 1555–1564.
 Ozolins et al. (2013) Ozolins, V., Lai, R., Caflisch, R., and Osher, S. 2013. Compressed modes for variational problems in mathematics and physics. Proceedings of the National Academy of Sciences 110, 46, 18368–18373.
 Pinkall and Polthier (1993) Pinkall, U. and Polthier, K. 1993. Computing discrete minimal surfaces and their conjugates. Experiment. Math. 2, 1, 15–36.
 Pokrass et al. (2013) Pokrass, J., Bronstein, A. M., Bronstein, M. M., Sprechmann, P., and Sapiro, G. 2013. Sparse modeling of intrinsic correspondences. Computer Graphics Forum 32, 2pt4, 459–468.
 Rustamov (2007) Rustamov, R. M. 2007. Laplacebeltrami eigenfunctions for deformation invariant shape representation. In Proceedings of the Fifth Eurographics Symposium on Geometry Processing. SGP ’07. Eurographics Association, AirelaVille, Switzerland, Switzerland, 225–233.
 Schwartz et al. (1989) Schwartz, E. L., Shaw, A., and Wolfson, E. 1989. A numerical solution to the generalized mapmaker’s problem: Flattening nonconvex polyhedral surfaces. IEEE Trans. Pattern Anal. Mach. Intell. 11, 9 (Sept.), 1005–1008.
 Shtern and Kimmel (2014a) Shtern, A. and Kimmel, R. 2014a. Iterative closest spectral kernel maps. In 3D Vision (3DV), 2014 2nd International Conference on. Vol. 1. IEEE, 499–505.
 Shtern and Kimmel (2014b) Shtern, A. and Kimmel, R. 2014b. Matching the LBO eigenspace of nonrigid shapes via high order statistics. Axioms 3, 3, 300–319.
 Shtern and Kimmel (2015) Shtern, A. and Kimmel, R. 2015. Spectral gradient fields embedding for nonrigid shape matching. Computer Vision and Image Understanding 140, 21 – 29.
 Simon (2005) Simon, B. 2005. Functional integration and quantum physics. 2nd ed., 2nd ed. ed. Providence, RI: AMS Chelsea Publishing.
 Sun et al. (2009) Sun, J., Ovsjanikov, M., and Guibas, L. 2009. A concise and provably informative multiscale signature based on heat diffusion. In Proceedings of the Symposium on Geometry Processing. SGP ’09. Eurographics Association, AirelaVille, Switzerland, Switzerland, 1383–1392.
 Touma and Gotsman (1998) Touma, C. and Gotsman, C. 1998. Triangle mesh compression. PROC GRAPHICS INTERFACE. pp. 2634. 1998.
 Vallet and Levy (2008) Vallet, B. and Levy, B. 2008. Spectral Geometry Processing with Manifold Harmonics. Computer Graphics Forum.
 Van Der Aa et al. (2007) Van Der Aa, N., Ter Morsche, H., and Mattheij, R. 2007. Computation of eigenvalue and eigenvector derivatives for a general complexvalued eigensystem. Electronic Journal of Linear Algebra 16, 1, 300–314.

Wei
et al. (2016)
Wei, L., Huang, Q., Ceylan, D., Vouga, E., and
Li, H. 2016.
Dense human body correspondences using convolutional networks.
In
Computer Vision and Pattern Recognition (CVPR)
.  Weyl et al. (1950) Weyl, H. et al. 1950. Ramifications, old and new, of the eigenvalue problem. Bulletin of the American Mathematical Society 56, 2, 115–139.
 Zigelman et al. (2002) Zigelman, G., Kimmel, R., and Kiryati, N. 2002. Texture mapping using surface flattening via multidimensional scaling. IEEE Transactions on Visualization and Computer Graphics 8, 2 (Apr), 198–207.
Appendix A Proof of Theorem 2.
Let us be given the Hamiltonian operator .
Recall the Courant–Fischer minmax principle; see also [Aflalo et al. (2015)] and [Brezis (2010)] Problems 37 and 49. We have for every ,
(56) 
That is, the min is taken over a linear subspace with is the Sobolev space
of codimension and the max is taken over all such subspaces.
Set , so that is a subspace of codimension
.
By 56 we have that for all
(57) 
and thus
(58) 
On the other hand, by 56,
(59) 
Appendix B Diffusion kernel of the Hamiltonian
In order to solve the diffusion equation, we first need to find the fundamental solution kernel to the Dirichlet problem that yields the heat equation
(60) 
Recall that for we return to the regular LBO diffusion case.
Suppose that has a eigendecomposition
. In that case, we can write
(61) 
and from the linearity of we have
(62) 
Since