Nonisometric Surface Registration via Conformal Laplace-Beltrami Basis Pursuit

09/19/2018 ∙ by Stefan C. Schonsheck, et al. ∙ Rensselaer Polytechnic Institute 4

Surface registration is one of the most fundamental problems in geometry processing. Many approaches have been developed to tackle this problem in cases where the surfaces are nearly isometric. However, it is much more challenging to compute correspondence between surfaces which are intrinsically less similar. In this paper, we propose a variational model to align the Laplace-Beltrami (LB) eigensytems of two non-isometric genus zero shapes via conformal deformations. This method enables us compute to geometric meaningful point-to-point maps between non-isometric shapes. Our model is based on a novel basis pursuit scheme whereby we simultaneously compute a conformal deformation of a 'target shape' and its deformed LB eigensytem. We solve the model using an proximal alternating minimization algorithm hybridized with the augmented Lagrangian method which produces accurate correspondences given only a few landmark points. We also propose a reinitialization scheme to overcome some of the difficulties caused by the non-convexity of the variational problem. Intensive numerical experiments illustrate the effectiveness and robustness of the proposed method to handle non-isometric surfaces with large deformation with respect to both noise on the underlying manifolds and errors within the given landmarks.



There are no comments yet.


page 14

page 15

page 16

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

The computation of meaningful point-to-point mappings between pairs of manifolds lies at the heart of many shape analysis tasks. It is crucial to have robust methods to compute dense correspondences between two or more shapes in different applications including shape matching, label transfer, animation and recognition [1, 2, 3, 4, 5, 6]. In cases where shapes are very similar (isometric or nearly isometric), there are many approaches for computing such correspondences [7, 8, 9, 10, 11, 12, 13, 6, 14, 15]. However, it is still challenging to compute accurate correspondences when the shapes are far away from near isometry.

One of the key challenges in largely deformed non-isometric shape matching is that the intrinsic features of the two shapes are not similar enough for standard techniques to recognize their similarity. For example, when computing the correspondence between human faces, it is not particularly difficult to geometrically characterize the structure of a ‘nose’. However, similar techniques can not work well to compute a map between a horse and an elephant face since these two surfaces have many largely deformed local structures including the drastic difference between the trunk of the elephant and the nose of the horse. Because of this, it is crucial to develop new methods to adaptively characterize large deformations on surfaces.

The LB eigensystem is a ubiquitous tool for 3D shape analysis (see [16, 3, 17, 18, 19, 20, 21, 22, 10, 23, 24, 25, 14]

and references therein). It is invariant under isometric transformations and intrinsically characterizes the local and global geometry of manifolds through its eigensystem up to an isometry. In principle, the LB eigensystem reduces infinite dimensional nonlinear isomorphism ambiguities between two isometric shapes to a linear transformation group between two LB eigensystems. This linear transform is necessary due to the possible sign or sub-eigenspace ambiguity of LB eigensystems 

[13]. Additionally, similar shapes often have similar eigensystems which allows for joint analysis of similar shapes their spectral properties [6]. However, when the deformation between two shapes is far from an isometry, the large dissimilarity between LB eigensystems of two shapes is the major bottleneck to adapt existing spectral geometry approach to conduct registration.

A natural idea to extend spectral geometry methods to register non-isometric surfaces is to deform the metric of the target surface to the source one such that two surfaces share similar LB eigensystem after deformation. However, the amount of imposed deformation is essentially relied on the correspondence between these two surfaces. In this work, we proposed to simultaneously compute the deformation and learn features for registration. Mathematically, one way to characterize these deformations is through the conformal factor. It is well known that there exists a conformal mapping between any two genus zero surfaces [26]. Rather than reconstruct conformally deformed surfaces and related conformal map, we exploit a fundamental link between the conformal factor and the LB eigesystem by considering the conformally deformed LB eigensystem. This allows us to compute a new basis on the target surface to align the naturally defined LB eigensystem on the source surface. It leads to a variational method for non-isometric shape matching which enables us to overcome the natural ambiguities of the LB eigensystem, to align the bases of non-isometric shapes and to avoid the direct computation of conformal maps.

Numerically, we solve our model using a proximal alternating minimization (PAM) method [27] hybridized with the augmented Lagrangian method [28]. The method is iteratively composed of a curvilinear search method on orthogonality constrained manifold [29]

in one direction to compute the conformally deformed LB eigenfunctions and the BFGS 

[30] method for the other direction to compute the conformal factor. Theoretically, we can guarantee the local convergence of the proposed algorithm since the objective function and constraints satisfy the necessary Kurdyka-Lojasiewicz (KL) condition [27]. Numerical results on largely deformed test problems, including horse-to-elephant and Faust benchmark database [31], validate the effectiveness and robustness of our method.

Related Works.

A large number of 3D nonrigid shape matching approaches are based on analysis of the LB eigen system (see [32, 3, 17, 18, 20, 6, 12, 33, 14, 13] and reference therein). The LB eigensytem is intrinsic and invariant to isomorphism, and also characterizes the local and global geometry of a manifold. This makes it ideal for many shape processing tasks and many early works in the field involve directly comparing the LB spectrum of the shapes to determine how alike shapes are [32, 3, 17]. More recently, the general concept of functional maps [6] has played a central role in many new methods which have allowed for the formulation of accurate correspondence maps. This technique essentially reduces the non-linear transform between two shapes to a linear transform between their eigensystems. In general, these techniques work for well for isometric and near isometric cases, but can not produce satisfactory results when the LB eigensystems of shapes are very dissimilar. This occurs when the deformation between shapes is far from an isometry. To overcome this, the concept of coupled bases (also known as joint-diagonalization) was introduced for shape processing tasks in [12]. In this work the authors propose a variational model to define a shared basis for a pair of shapes which is ‘nearly harmonic’ on one shape and ’similar’ to the natural LB basis on the other. This joint optimization allows for much more accurate correspondence maps, but does not characterize the underlying deformations which lie at the heart of the non-isometric shape matching problem.

Conformal maps have been widely applied to various shape processing tasks in order to characterize these deformations [34, 35, 8, 36]. In one of the first works to combine spectral and deformation based approaches, [37] presents a scheme to find optimal conformal deformation to align two shapes in the embedded LB Space. Additionally, the authors present a general framework for computing LB eigensytems of conformaly deformed surfaces as well as several other imported related quantities. Continuing on this line of work in [38]

, the authors use the LB eigenvalues as a tool to guide conformal deformations. Using derivatives of the the LB eigenvalues, they compute optimal conformal metrics which approximate conformal and topological eigenvalues. In our work, we use the spectral coefficients of known features to guide the deformation, so rather than align the eigenvalues we align the eigenfunctions. This allows us to avoid the subspace ambiguity of the the LB eigensystem and computational errors in calculating high frequency eigenvalues.

Major Contributions.

We introduce a novel variational basis pursuit model for computing non-isometric shape correspondences via conformal deformation of the LB eigensystem. This model enhances spectral approaches from handling nearly isometric surface registration to tackling surfaces with large deformed metrics. It naturally combines the conformal deformation to the LB eigensystem and simultaneously computes surface deformations and LB eigenbasis which also automatically overcomes the ambiguities of LB eigensystems in surface registration. We also propose a numerical scheme to solve the variational model with local convergence guarantee. Additionally, we introduce a reinitializaiton scheme to help tackle local minima and improve the quality of the computed bases. This algorithm successfully handles non-isomorphic shape correspondence problems given only a few landmarks and is shown to be robust to noise and perturbations of landmarks.

The rest of this paper is organized as follows: In section 2, we review the theoretical background of conformal deformation of LB eigensystem and functional maps. After that, we propose the variational basis pursuit model for conformal deformation of the LB eigensystem in section 3. In section 4, we discretize the model and develop an optimization scheme based on PAM to solve the variational problem. Section 5 is further devoted to discuss a few details of the model and a reinitialization scheme to improve our numerical solver. In section 6, numerical results on several data sets are presented to show that the model accurately produces point-to-point mappings on non-isometric manifolds with large deformation given only a few landmark points. We also show that our approach is robust to both noise in the underlying manifolds and inaccuracies in the initial landmarks. Furthermore, we test the model to a benchmark data based to show its effectiveness. Lastly, we conclude the paper in Section 7.

2 Mathematical Background

In this section, we discuss the mathematical background of the proposed method. We first review a few key properties of the LB eigensystem of a Riemannian surface and discuss its conformal deformations with respect to deformations of the Riemannian surface metric [39, 26]. After this, we review the functional maps framework in [6] which will be closely related to our work.

2.1 Conformal deformation of LB eigensystem on Riemannian Surfaces

Given a closed Riemannian surface , its LB operator in a given local coordinate system, , is defined as [39, 26]:


where is the inverse of the metric matrix

and . The LB operator is self-adjoint and elliptic, therefore it has a discrete spectrum. We denote the eigenvalues of as with the corresponding eigenfunctions satisfying:


where is the area element on with respect to . It is well-known that forms an orthonormal basis for the real-valued, smooth function space on the manifold . This basis can be viewed as a generalization of the Fourier basis from flat space to a differentiable manifold. The LB eigensystem is invariant under both rigid and nonrigid isometric transformations, and it uniquely determines a manifold up to isometry [16].

In differential geometry, a conformal map is one which preserves angles locally. Formally, a conformal map preserves the first fundamental form up to a positive scaling factor. Given two manifolds and , a map is conformal if and only if the pullback with a positive function (written this way to emphasize positivity). A conformal deformation of a surface is a transformation which changes the local metric by a positive scaling factor. A well known result in conformal geometry is that there exists a conformal map between any two genus zero surfaces [26].

Given a closed surface with conformal deformation , the LB eigensystem of the deformed manifold can be viewed as a weighted LB eigensystem on the original surface . This simple fact intrinsically links the LB eigensystem of the deformed manifold to a weighed LB eigensystem on the original manifold. It allows us to compute the LB eigensytem of the conformally deformed manifold without explicitly reconstructing its embedding or coordinates. This also relates information about the local deformation and global eigensytem and later becomes the cornerstone of our approach. Formally, we have:

Proposition 1.

Let be a LB eigensystem of a conformally deformed surface , then is equivalent to the following weighted LB eigensystem on :


This is because

Hence the eigen problem: is equivalent to . In addition,

The problem of finding the LB eigensystem of a Riemannian manifold is equivalent to finding an orthonormal set of functions which have minimal harmonic energy on the surface. From the above proposition, the LB eigensystem of a conformally deformed manifold can be formulated as the following variational problem:


2.2 Functional Maps

Functional maps were introduced in [6] for isometric and nearly isometric shape correspondence. This method has been shown a very effective tool for various shape processing tasks [6, 12, 40]. Here we provide a basic overview of their framework. Consider Riemannian surfaces and , a smooth bijection induces a linear transformation between functional spaces of these two manifolds as:

Instead of computing surface map , the crucial idea of functional map is to compute the linear map between these two functional spaces. After that, the desired surface map can be encoded by considering images of indicator functions under .

Finding a functional map, , associated with a map is equivalent to finding the matrix representation of under a fixed orthonormal basis of and a fixed orthonormal basis of , respectively. Namely, if we write , then any two given corresponding functions and under can be represented using as:

Each entry of the matrix can be found by finding the coefficient of expressed in the coordinate system, i.e. . In practice, one can use two finite sets of orthonormal functions to approximate and , thus the functional map can be approximated by a finite dimensional matrix. For instance, the first eigenfunctions of the LB eigensystem is one common choice of such a basis. Then, the problem of finding the transformation can be approximated by the problem of seeking a finite dimension matrix . As long as is computed, the desired map can be computed through operating on indicator functions.

3 Conformal LB Basis Pursuit for Nonisometric Surface Registration

In this section, we propose a LB basis pursuit model for non-isometric surface registration. On the target surface , the model simultaneously finds a conformal deformation and a conformally deformed LB eigensystem so that the coefficients of the corresponding feature functions expressed on the deformed LB eigensystem of are the same as the coefficients on the fixed source surface .

3.1 Variational PDE model

Given two non-isometric genus-0 closed Riemannian surfaces , we aim at finding a geometrically meaningful correspondence between these two surfaces. In the case that and are nearly isometric, there are many successful methods to constructing maps between and by comparing their isometric invariant features. Using spectral descriptors from solutions of the LB eigensystem on manifolds is a common way of constructing such descriptors [3, 17, 41, 42, 20]. As extensions, some other descriptors such as Heat kernel signature [19], wave kernel signature [10] and optimal spectral descriptors [43] have also been proposed in the literature. However, most of the existing methods consider the construction of descriptors for nearly isometric manifolds. Registration methods based on the existing LB spectral descriptors can not provide satisfactory results for constructing correspondence between two non-isometric surfaces as their eigensystems are possibly quite far apart.

We propose to overcome the limitation of the LB spectral descriptors for largely deformed non-isometric shape registration by considering a continuous deformation of the LB spectral descriptors. Intuitively, given two non-isometric shapes and , our idea is to deform the metric of such that the deformed surface is isometrically the same as . Then the LB spectral descriptors can be applied as in isometric shape matching. However, it is challenging to find an appropriate deformation as the accurate amount of deformation on each local region of depends exactly on an accurate correspondence which is precisely the problem we would like to solve.

To handle this challenge, we propose to simultaneously find an optimal correspondence and an optimal deformation. More specifically, by fixing the LB eigensystem of , we seek a map and a conformal factor such that the LB eigensystem of can be aligned to the LB eigensystem of via . This problem can be written as the following variational PDE problem:


where and .The first term measures the alignment of two bases as the correct correspondence should map one LB eigensystem to another one, and the second term solves the first LB eigenfunctions for the deformed manifold due to the variational problem (4). Existence of a solution to this variational problem (5) is guaranteed as any two genus-0 surfaces are conformally equivalent and the LB operator is invariant under isometric transformations.

Computationally, the numerical search for in the mapping space is usually very time consuming. Inspired by the idea of functional maps [6] and the coupled quasi-harmonic bases [12], we choose to represent in the functional space. Instead of finding directly, we look for a basis which is nearly harmonic on and represents the corresponding features with the same coefficients as does. More precisely, given a set of corresponding features on and on , such that if and are corresponding points on and , we can replace the direct measurement of the basis alignment term with a coefficient matching term. That is, instead of measuring the alignment of and via , we measure how closely the coefficients for in the computed basis match the coefficients for in the fixed LB basis . Formally, we measure the coefficient alignment by constructing a matrix of the coefficients in for in and for in so that the term represents the coefficient for the corresponding function in the basis and computing their difference under the Frobenius norm. With this in mind, we propose the following model:



In practice we use indicator functions for and , but heat signatures [19], wave kernel signatures [10], or any other corresponding functions will also work. Once is obtained, we can easily compute the functional map as


The main advantage of this model over previous existing methods for shape correspondence is that we are able to employ much more of the information encoded in the differential structures of and in our algorithm by combining the spectral descriptors and local deformations. This additional flexibility enables us to compute correspondences between largely deformed shapes. Information about the conformal deformation of the metric allows us to find a harmonic basis on the deformed shape, meanwhile information about the alignment of the functional spaces guides our calculation of the conformal deformation. Furthermore the additional constraint of the feature alignment overcomes ambiguity casued by the fact that there is no unique conformal deformation between any two genus zero surfaces. To the best of our knowledge, the link between the conformal factor and deformed LB basis has not been exploited in such a way. Previous works have used only the conformal factor [8, 11] or only the functional space [6, 12] as stand alone tools rather than in concert as we present here.

3.2 Regularization and Area Constraint

We add harmonic energy term to smooth the conformal deformation and regularize the problem. This can both increase the speed of the algorithm and improve the quality of the map, both in terms of the geodesic errors of the final correspondence, and the accuracy of the resulting conformal factor. This is particularly helpful to handle deformations between the shapes which are far from isometry and to reduce the required number of features. Rather than smooth the conformal factor directly, we instead add the harmonic energy of to our objective function. Using instead of allows for easier analytic computation of the derivatives and a more efficient algorithm. In cases where the deformations are likely to be highly localized, this term may be omitted.

Lastly, we add an area preservation constraint to our model. That is, we would like the final deformed shape to be of the same size as the one we are matching it to. To enforce this, we mandate that the deformed manifold have the same surface area as the original manifold. This eliminates any scaling ambiguity. Then the final version of our model can be stated as:


where and

4 Discretization and Numerical Algorithms

In this section, we describe a discretization of the proposed variational model (8) using on triangular representation of surfaces. After that, we design a numerical algorithm to solve the proposed model based on proximal alternating minimization method.

4.1 Discretization of the Model

The main method we use to discretize surfaces and differential operators is based on a finite element scheme similar to that developed in [3, 19, 44]. Let be a set of vertices sampled on the manifold . A surface can be discretized as a triple made of vertices (), connected by edges () which form triangular faces . We define the first ring of , the set of all triangles which contain as . For each edge connecting points and , we define the angles opposite as angles and .

We define a diagonal mass matrix, , a positive definite matrix with entries given by:

We use this simplified version, rather than the standard finite element discretization, for convince in order to avoid expensive factorizations later in our algorithm. We remark that the standard version can also be used in our algorithm at the cost of speed. The surface area can be approximated as . Similarly, a function with discretization , then we have the approximation . The stiffness matrix, S, is a symmetric positive semidefinite matrix given by:

where is a linear pyramid function which is at and zero elsewhere. The mass and stiffness matrices can be used to approximate the LB eigenvalue problem as: [45].

We remark that one can also work with point clouds representation instead of triangulated meshes. These definitions for the stiffness and mass matrices can be approximated by the point clouds method discussed in [46]. The only change we would need to make is to use only the diagonal entries of the version of the mass matrix proposed in their paper to populate the strictly diagonal version we employ here.

Suppose two surfaces are represented by triangular meshes with the same number of points111In fact, we do not need to require that the surfaces have the same number of points, but doing so for now will allow for more convenient notation.. We denote as the mass and stiffness matrices of and let be the first LB eigenfunctions of , and be feature functions. Similarly, we write as the mass and stiffness matrices of , as the first LB eigenfunctions of (under ), and as corresponding feature functions, ordered the same as in . We also write as the discretized conformal factor on and as a diagonal matrix.

Therefore, the discretized optimization model (6) can be written as:


Here is the identity matrix and . Since is symmetric positive definite and diagonal, we can easily calculate the matrix decomposition . If we also substitute , then (9) can be written as:


where Note that this parameterization of the problem moves the conformal factor out of the orthogonality constraint (and into ). We will soon see that, for any fixed , this will make the problem for easier to solve.

4.2 Numerical Optimization

The two variables and in (10) make the optimization problem different from orthogonality constrained problems solved by nonconvex alternating direction method of multipliers (ADMM) methods considered in [47, 48, 49, 50]. Rather than solve this problem directly for and simultaneously by directly minimizing (10), we employ a method based on the framework of proximal alternating minimization (PAM) method [27]. The idea of PAM is to update variables using a Gauss-Seidel method with proximal regularizations in the following way.

Let and . We also define indicator functions

Then it is clear that and are semi-algebraic functions as and are zero sets of polynomial functions [51]. Therefore, we write an equivalent form of (10) as


Using the PAM method, we have the following iterative scheme


Here is a step size parameter. Essentially, these proximal terms penalizes large step sizes in and prevents the algorithm from “jumping” between multiple local minimums. The addition of these proximity terms allows us to analyze our proposed method in the framework of the PAM algorithm [27]. It has been shown in [27, 51, 52] that such proximal terms can guarantee the solutions generated at each step converge to a critical point o bbbbbbbbf the objective function. Formally, we have the following convergence theorem in accordace with Theorem 9 in [27].

Theorem 1.

Let be the sequence produced by (12), then the following statement hold:

  1. .

  2. .

  3. converges to a critical point of .


To prove this, we show that our model obeys the conditions required for local convergence of PAM in [27]. To do so, we need:

(1) Each of term in the objective which contains only one primal variable is bounded below and lower semicontinous.

(2) Each term in the objective which contains both variables is and has a locally Lipschitz continuous gradient.

(3) The objective satisfies the Kurdyka-Lojasiewicz (KL) property.
It is immediately clear that that the first two properties are satisfied by our objective. Furthermore, it is known that all semi-algebraic functions have KL property [27, 51, 48]. Our objective is semi-algebraic so we can guarantee local convergence of the proposed optimization method. ∎

We use the augmented Lagrangian method to solve the constrained sub-optimization problem for in (12). For convenience, let’s write


Overall, we solve (10) in the following way by hybridizing PAM with the augmented Lagrangian method.


The subproblem for minimizing requires a different approach to be solved. The main challenge in minimizing the first sub-optimization problem is the nonconvex orthogonality constraints. Recently, several approaches have been developed to solve orthogonally constrained problems in feasible or infeasible ways [29, 47, 48, 49, 50]. For our implementation, we have chosen the feasible approach developed in [29] which uses a curvilinear method based on the Cayley transform together with Barzilai-Bowein step size line search. This method updates variables along a geodesic curve on the Stiefel manifold, a geometric description of the orthogonality. It preserves the orthogonality constraints and guarantees convergence to critical points in our scenario. More precisely, given a feasible starting point and the coordinate gradient at this point, the update scheme is as follows:


Here is a step size parameter chosen by the Barzilai-Bowein criteria developed in [53]. Although convergence to a global minimum is not guaranteed, this method has proven effective for our purposes and only requires the computation of the objective function and its coordinate gradient with respect to at each step provided by:


The subproblem for (as written in (14)), on the other hand is smooth and unconstrained. For our implementation, we use the well known quasi-Newton BFGS algorithm [30]. The gradient of objective function with respect to can be written as:


where denotes the diagonal of the matrix, signifies element-wise Hadamard product and is the inverse of diagonal matrix multiplied with itself.

4.3 Computation of Point-to-Point Map

One naive way to compute a point-to-point map is to find the functional map by using the final deformed manifold and its LB eigensystem with respect to the deforamtion. However, this may not work well because of the ambiguity of LB eigensystem. Additional effort is needed to handle possible ambiguity of LB eigensystem such as the method discussed in [13]. As an advantage of the proposed method, the resulting basis generated by the proposed algorithm (recovered as ) to will naturally correct ambiguities of LB eigensystem. This is similar to the method discussed in [12]. Thus, we can compute the functional map as . However, this method is still quite inefficient and may be sensitive to small errors in the resulting basis.

Instead, after we recover the final basis from our method, we can compute the point-to-point map between the two surfaces by comparing the values of each of the basis functions. This is essentially the same scheme presented in [6]

, but applied to our new basis. We use a KNN search (with

) to match rows of and . This requires a search of points in dimension, but is much more efficient and accurate than using the delta function approach described in the previous paragraph. Other methods used to refine functional maps such as [40] can be applied in this setting without changes.

We summarize our numerical method for nonisometric surface registration as Algorithm LABEL:alg:LBBasisPursuit1.


5 Discussion

In this section, we discus our choice of feature functions, as well as ways to overcome problems which may arise from the non-convexity of the proposed optimization problem. In addition, we present a novel way to jointly measure the quality of the correspondence and alignment of the bases without any prior knowledge about the ground truth of the point-to-point map.

5.1 Choice of Feature Functions

The simplest, and in many applications, most natural features to choose for and are indicator functions. Let be a set of points on and be a corresponding set on . We can view each and as a -function on and respectively to indicate these landmarks.

Another option is to use heat diffusion functions. Given a corresponding pair of points we can use delta functions to define an initial condition and solve the heat diffusion problem using the Crank-Nicholson scheme where is a step size parameter. By taking “snap shots” (solutions of the equation for various values) of at different time values we can generate multiple functions from a single corresponding pair. This choice allows for a multi-scale selection of features.

The wave kernel signature (WKS) has also been used for characterizing points on non-rigid three dimensional shapes [10]. These functions are defined as the solutions to the Schrodinger equation: at different points on the surface. Given two corresponding points we can solve the equation at each point and use these as our corresponding functions. However, the solutions to these equations are highly dependent on both local and global geometries of the manifold. Because of this, they are only suitable for shape correspondence when the shapes are very similar and, in general, do not work well for non-nearly-isometric problems. The same problem exists for heat diffusion features, however, in general heat diffusion tends to be much more stable with respect to local deformations.

5.2 Reinitialization Schemes

Although we have shown that the proposed PAM based optimization algorithm converges to a critical point of the objective function, it is still challenge to achieve global optimizers as the problem is non-convex. In practice, we have found that the numerical results can often be improved in terms of both accuracy and speed of computation by adding a simple reinitialization scheme to our algorithm. The motivation for the scheme comes from an observation that if we know the exact conformal deformation and the source surface has a simple eigensystem (no repeated eigenvalues), then the LB eigensystem of is the same as the LB eigensystem of up to a change in sign. With this in mind, we propose to reinitialize the problem by resetting to be the solution to weighed eigenproblem . Computationally, to avoid introducing ambiguities of LB eigensystem by calling standard eigen-solvers, we solve a discrete counterpart to (4) as based on the curvilinear search method discussed in Section 4.2 with the current eigensystem, , as an initial guess for this problem. By using as warm start for the eigenproblem we can avoid re-introducing sign or multiplicity ambiguities into the problem which our algorithm has already resolved.

If using heat diffusion, wavelet kernel signatures, or any other functions which are defined based on local geometry as the input feature functions, then we also need to recalculate these functions with respect to the conformally deformed metric. For example, if we are using heat diffusions, we can recompute the heat diffusion functions on the deformed manifold by multiplying the mass matrix by in the Crank-Nicholson scheme: , A similar re-computation technique can be applied if using WKS features.

With this re-initialization procedure, we propose a modified version of our numerical solver as Algorithm LABEL:alg:LBBasisPursuit2.


6 Numerical Experiments

In this section, we apply our algorithm to several problems. We begin by working on a typical non-isomorphic matching problem for a pair of shapes with a large deformation: a horse and an elephant. We preform tests showing the effectiveness of our approach given different amounts of landmark points, and demonstrate robustness with respect to noise both on the manifold and in the initial correspondences. Finally, we conduct experiments on the Faust benchmark data set [31]. All numerical experiments are implemented in MATLAB on a PC with a 32GB RAM and two 2.6GHz CPUs.

In all of our experiments, we use randomly chosen correspondence points to create indicator functions as the input features. The first 100 non-trivial LB eigenfunctions are chosen to calculate the coefficient matching term, as well as for computing the final correspondence. We set for all experiments, even though the data sets and experimental conditions are very different. This choice of and allows the coefficient matching terms and eigenfunction term to balance each other out, with the choice of still being large enough to preserve the area constraint. is chosen to be small so that the harmonic energy, which tends to be quite large, does not dominate the others. We observe that our algorithm is robust to different choices of the parameters.

6.1 A large deformation pair: Horse to Elephant

Figure 1: Left: visualization of point-to-point map. Right: normalized geodesic errors for various numbers of landmarks with and without reinization.

The first experiment is designed to test the effectiveness of the proposed method on a pair of shapes with large deformation. Each surface, a horse and an elephant, is represented by a mesh with 1200 points. One of the challenges in this pair is the large deformations in the sharp corner and elongated regions including ears, teeth, noses and tails on the horse and elephant surfaces. Those regions make the registration problem very difficult.To demonstrate the efficacy of our approach, we perform this experiment under several different conditions. Our algorithm produces excellent results given a sufficient number of landmarks, and it still finds reliable correspondences given limited landmarks. We also show that using our reinitialization scheme (Algorithm LABEL:alg:LBBasisPursuit2) produces a more accurate map than without this extra step (Algorithm LABEL:alg:LBBasisPursuit1).

Figure 2: Left: convergence curves of our method. The coefficient matching term measures: . The eigen problem is: and the harmonic energy measures: and the total energy is the entire model derived in (10). Right: resulting and exact conformal factors.

Figure 1 shows the results of using 100, 75, 50 and 25 known landmark points with and without our reinitialization scheme. To qualitatively measure the mapping quality, we calculate the normalized geodesic distance from the point on the target surface produced by the map to ground truth following the Princeton Benchmark method [11]. These distances are collected into a cumulative error on the right of Figure 1 where the -axis measures the percent of points whose distances are less than or equal to the -axis value. For example, in the case of known landmarks, our algorithm matches over of the points to exact correct point and more than within a error margin.

Figure 3: First two rows: The first 9 non-trivial natural LB eigenfunctions of manifolds. Third row: results from the proposed basis pursuit algorithm. Fourth row: ground truth.

The left of Figure 2 shows the convergence of the objective function and illustrates the effectiveness of the reinitialization step. We plot the three terms in the objective function separately as well as the overall objective. We typically observe that the convergence curves in the coefficient matching and total energy flatten quickly as the algorithm tends to a local minimizer. However, each reinitialization significantly reduces the objective function, thus effectively overcoming the non-convexity of the model. We further demonstrate the validity of our algorithm by examining the resulting conformal factor. In the right image of Figure 2, we show the conformal factor calculated by our algorithm as well as the ground truth. The ground truth conformal factor is calculated by using the ground truth point-to-point map to compare the area of the first ring structure around each point on the source and target surface. Here we plot where for better visualization. From this figure we can confirm that the conformal mapping our algorithm produces is very close to the true factor.

Figure 4: Top row: 9th, 11th and 44th natural LB eigenfunctions on source. Bottom row: results and ground truth.

Since the elephant and horse are dramatically different shapes, the large dissimilarity of their natural LB eigenfunctions (first two rows of Figure 3) cannot be expected to produce meaningful correspondence. However, our model overcomes this by capturing the conformal deformation between the surfaces. As a result, the basis computed for the horse (target surface) by our model is consistent with the LB eigenfunctions of the elephant (source surface). We further compare these results with the ground truth which is calculated through the push forward of the LB eigenfunction of the source to the target surface using the a priori map. Figure 4 highlights the consistency of the produced bases on several highly distorted regions. Specifically, we focus on each of the ears, the nose/trunk and the tails. From Figures 3 and 4, we can see that our approach produces a new basis on the target that aligns very closely to the natural LB basis on the sources manifold. This is the reason that accurate registration results can be obtained using the new basis.

6.2 Noisy Data

In this experiment, we demonstrate that our algorithm can handle noisy data. Since noise on the surfaces can be viewed as local deformations, our algorithm is automatically robust to geometric noise. Medical scans often have noise resulting from the imaging instruments and segmentations. Our model can solve registration problems for this type of data. To demonstrate this, we generate noisy data by adding noise along the normal of each point. Figure 5 shows the results of two experiments: a noisy elephant to an elephant and a noisy horse to an elephant. We observe that our algorithm still produces very accurate results despite this noise.

Figure 5: Left: point-to-point maps for noisy data. Right: normalized geodesic errors for noisy data.

6.3 Perturbed Landmarks

We also demonstrate the robustness of our algorithm to landmark perturbations. Working again on the horse and elephant, we test cases where the landmarks are perturbed to another vertex within the first ring. The magnitude of these perturbations depends on the uniformality and meshing of the surface. The left graph in Figure 6 shows the size of the perturbations of the landmarks points as well as the error in their final mapping. The right graph in Figure 6 compares the geodesic error of the for all points when 25%, 50% and 100% of the landmarks points are perturbed. From these tests we conclude that our method can successfully reduce the error introduced in the perturbed landmarks and still produce accurate maps in the presence of perturbations.

Figure 6: Left: Initial perturbations to landmarks and final error of landmarks. Right: Final registration geodesic errors for all points using perturbed landmarks.

6.4 Benchmark test using the Faust Dataset

In our next experiment, we test our algorithm on a larger dataset to demonstrate its effectiveness and robustness on a variety of shapes. The Faust dataset is a collection of 100 3D shapes composed of 10 real individuals in 10 distinct poses (Figure 7) [31]. Instead of testing all 9900 possible correspondences be each of the pairs, we select two smaller subsets of shapes to formulate to smaller test sets. For the first test, we randomly choose 100 pairs of shapes and compute the correspondences. In the second test, we choose l0 scans and ensure that each individual and each pose is represented exactly once in the test set and compute all the 90 correspondences. This selection criteria ensures that no pairs are from the same the pose or individual. The bottom left graph in Figure 7 shows the average error of the mappings for each of these tests. We see that our algorithm again computes very accurate correspondences for both tests. Furthermore, we see that the results for the harder test set are very close to the results for the first test set. This indicates that our approach can effectively handle non-isometric matching problems with large deformations.

Figure 7: Top: selected subjects from the Faust data set [31]. Bottom: geodesic errors for randomly selected and least isomorphic pairs.

7 Conclusions

In this work, we have developed a variation method for computing correspondence between pairs of largely deformed non-isometric manifolds. Our approach considers conformal deformation of the manifolds and combines with traditional LB spectral theory. This method naturally connects metric deformations to the spectrum of the manifold and therefore allows us to register manifolds with large deformations. Our approach simultaneously aligns the bases of the manifolds and computes a conformal deformation without having to explicitly reconstruct the deformed manifolds. We have also proposed an efficient, locally convergent method to solve this model based on the PAM framework. Finally, we have conducted intensive numerical experiments to demonstrate the effectiveness and robustness of our methods.


The research of S. Schonsheck and R. Lai is supported in part by NSF DMS–1522645 and an NSF Career Award DMS–1752934.


  • [1] Kaleem Siddiqi, Yves Bérubé Lauziere, Allen Tannenbaum, and Steven W Zucker. Area and length minimizing flows for shape segmentation. IEEE Transactions on Image Processing, 7(3):433–443, 1998.
  • [2] Vladislav Kraevoy and Alla Sheffer. Cross-parameterization and compatible remeshing of 3d models. ACM Transactions on Graphics (TOG), 23(3):861–869, 2004.
  • [3] Martin Reuter, Franz-Erich Wolter, and Niklas Peinecke. Laplace–beltrami spectra as ‘shape-dna’of surfaces and solids. Computer-Aided Design, 38(4):342–366, 2006.
  • [4] Tobias Heimann and Hans-Peter Meinzer. Statistical shape models for 3d medical image segmentation: a review. Medical image analysis, 13(4):543–563, 2009.
  • [5] Oliver Van Kaick, Hao Zhang, Ghassan Hamarneh, and Daniel Cohen-Or. A survey on shape correspondence. In Computer Graphics Forum, volume 30, pages 1681–1707. Wiley Online Library, 2011.
  • [6] Maks Ovsjanikov, Mirela Ben-Chen, Juston Solomon, Adrian Butscher, and Leonidas Guibas. Functional maps: a flexible representation of maps between shapes. ACM Transactions on Graphics (TOG), 31(4):30, 2012.
  • [7] Asi Elad and Ron Kimmel. On bending invariant signatures for surfaces. IEEE Transactions on pattern analysis and machine intelligence, 25(10):1285–1295, 2003.
  • [8] Xianfeng Gu, Yalin Wang, Tony F Chan, Paul M Thompson, and Shing-Tung Yau. Genus zero surface conformal mapping and its application to brain surface mapping. IEEE Transactions on Medical Imaging, 23(8):949–958, 2004.
  • [9] Alexander M Bronstein, Michael M Bronstein, and Ron Kimmel. Efficient computation of isometry-invariant distances between surfaces. SIAM Journal on Scientific Computing, 28(5):1812–1836, 2006.
  • [10] Mathieu Aubry, Ulrich Schlickewei, and Daniel Cremers. The wave kernel signature: A quantum mechanical approach to shape analysis. In Computer Vision Workshops (ICCV Workshops), 2011 IEEE International Conference on, pages 1626–1633. IEEE, 2011.
  • [11] Vladimir G Kim, Yaron Lipman, and Thomas Funkhouser. Blended intrinsic maps. In ACM Transactions on Graphics (TOG), volume 30, page 79. ACM, 2011.
  • [12] Artiom Kovnatsky, Michael M Bronstein, Alexander M Bronstein, Klaus Glashoff, and Ron Kimmel. Coupled quasi-harmonic bases. In Computer Graphics Forum, volume 32, pages 439–448. Wiley Online Library, 2013.
  • [13] Rongjie Lai and Hongkai Zhao. Multiscale nonrigid point cloud registration using robust sliced-wasserstein distance via laplace-beltrami eigenmap. SIAM Journal on Imaging Sciences, 10(2):449–483, 2017.
  • [14] Alon Shtern and Ron Kimmel. Spectral gradient fields embedding for nonrigid shape matching. Computer Vision and Image Understanding, 140:21–29, 2015.
  • [15] Alon Shtern, Matan Sela, and Ron Kimmel. Fast blended transformations for partial shape registration. Journal of Mathematical Imaging and Vision, pages 1–16, 2016.
  • [16] Pierre Bérard, Gérard Besson, and Sylvain Gallot. Embedding riemannian manifolds by their heat kernel. Geometric and Functional Analysis, 4(4):373–398, 1994.
  • [17] Bruno Levy. Laplace-beltrami eigenfunctions: Towards an algorithm that understands geometry. In IEEE International Conference on Shape Modeling and Applications, invited talk, 2006.
  • [18] Raif M. Rustamov. Laplace-beltrami eigenfunctions for deformation invariant shape representation. Eurographics Symposium on Geometry Processing, 2007.
  • [19] Jian Sun, Maks Ovsjanikov, and Leonidas Guibas. A concise and provably informative multi-scale signature based on heat diffusion. In Computer graphics forum, volume 28, pages 1383–1392. Wiley Online Library, 2009.
  • [20] Michael M. Bronstein and Iasonas Kokkinos. Scale-invariant heat kernel signatures for non-rigid shape recognition.

    IEEE Conference on Computer Vision and Pattern Recognition (CVPR)

    , pages 1704–1711, 2010.
  • [21] Rongjie Lai, Yonggang Shi, Kevin Scheibel, Scott Fears, Roger Woods, Arthur W. Toga, and Tony F. Chan. Metric-induced optimal embedding for intrinsic 3D shape analysis. Computer Vision and Pattern Recognition (CVPR), pages 2871–2878, 2010.
  • [22] Dan Raviv, Alexander M. Bronstein, Michael M. Bronstein, Ron Kimmel, and Nir Sochen. Affine-invariant diffusion geometry for the analysis of deformable 3d shapes. CVPR, pages 2361–2367, 2011.
  • [23] Rongjie Lai, Yonggang Shi, Nancy Sicotte, and Arthur W Toga. Automated corpus callosum extraction via laplace-beltrami nodal parcellation and intrinsic geodesic curvature flows on surfaces. In Computer Vision (ICCV), 2011 IEEE International Conference on, pages 2034–2040. IEEE, 2011.
  • [24] Yonggang Shi, Rongjie Lai, and Arthur W Toga.

    Cortical surface reconstruction via unified reeb analysis of geometric and topological outliers in magnetic resonance images.

    IEEE transactions on medical imaging, 32(3):511–530, 2013.
  • [25] Yonggang Shi, Rongjie Lai, Danny JJ Wang, Daniel Pelletier, David Mohr, Nancy Sicotte, and Arthur W Toga. Metric optimization for surface analysis in the laplace-beltrami embedding space. IEEE transactions on medical imaging, 33(7):1447–1463, 2014.
  • [26] Jürgen Jost. Riemannian geometry and geometric analysis. Springer Science & Business Media, 2008.
  • [27] Hédy Attouch, Jérôme Bolte, Patrick Redont, and Antoine Soubeyran. Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the kurdyka-łojasiewicz inequality. Mathematics of Operations Research, 35(2):438–457, 2010.
  • [28] Roland Glowinski and Patrick Le Tallec. Augmented Lagrangian and operator-splitting methods in nonlinear mechanics. SIAM, 1989.
  • [29] Zaiwen Wen and Wotao Yin. A feasible method for optimization with orthogonality constraints. Mathematical Programming, 142(1-2):397–434, 2013.
  • [30] Mokhtar S. Bazaraa, Hanif D. Sherali, and Chitharanjan M. Shetty. Nonlinear programming: theory and algorithms. John Wiley & Sons, 2013.
  • [31] Federica Bogo, Javier Romero, Matthew Loper, and Michael J Black. Faust: Dataset and evaluation for 3d mesh registration. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 3794–3801, 2014.
  • [32] Martin Reuter, Franz-Erich Wolter, and Niklas Peinecke. Laplace-spectra as fingerprints for shape matching. In Proceedings of the 2005 ACM symposium on Solid and physical modeling, pages 101–106. ACM, 2005.
  • [33] Dan Raviv and Ron Kimmel. Affine invariant non-rigid shape analysis. International Journal of Computer Vision, to appear, 2014.
  • [34] M. K. Hurdal, K. Stephenson, P.L. Bowers, D.W.L. Sumners, and D.A. Rottenberg. Coordinate systems for conformal cerebellar flat maps. NeuroImage, S467, 2000.
  • [35] S. Haker, S. Angenent, A. Tannenbaum, R. Kikinis, G. Sapiro, and M. Halle. Conformal surface parameterization for texture mapping. IEEE Trans. Visualization and Computer Graphics, 6(2):181—189, 2000.
  • [36] B. Springborn, P. Schröder, and U. Pinkall. Conformal equivalence of triangle meshes. ACM Transactions on Graphics (TOG) - Proceedings of ACM SIGGRAPH 2008, 27(3), 2008.
  • [37] Yonggang Shi, Rongjie Lai, Raja Gill, Daniel Pelletier, David Mohr, Nancy Sicotte, and Arthur W Toga. Conformal metric optimization on surface (cmos) for deformation and mapping in laplace-beltrami embedding space. In International Conference on Medical Image Computing and Computer-Assisted Intervention, pages 327–334. Springer, 2011.
  • [38] Chiu-Yen Kao, Rongjie Lai, and Braxton Osting. Maximization of laplace- beltrami eigenvalues on closed riemannian surfaces. ESAIM: Control, Optimisation and Calculus of Variations, 23(2):685–720, 2017.
  • [39] Isaac Chavel. Eigenvalues in Riemannian geometry, volume 115. Academic press, 1984.
  • [40] Emanuele Rodola, Michael Moeller, and Daniel Cremers. Point-wise map recovery and refinement from functional correspondence. arXiv preprint arXiv:1506.05603, 2015.
  • [41] Bruno Vallet and Bruno Levy. Spectral geometry processing with manifold harmonics. Computer Graphics Forum (Proceedings Eurographics), 2008.
  • [42] Yonggang Shi, Rongjie Lai, Shelia Krishna, Nancy Sicotte, Ivo Dinov, and Aurthor W. Toga. Anisotropic Laplace-Beltrami eigenmaps: bridging Reeb graphs and skeletons. In Computer Vision and Pattern Recognition Workshops, pages 1–7, 2008.
  • [43] Roee Litman and Alexander M Bronstein. Learning spectral descriptors for deformable shape correspondence. IEEE transactions on pattern analysis and machine intelligence, 36(1):171–180, 2014.
  • [44] Gerhard Dziuk and Charles M Elliott. Finite element methods for surface pdes. Acta Numerica, 22:289–396, 2013.
  • [45] Mark Meyer, Mathieu Desbrun, Peter Schröder, Alan H Barr, et al. Discrete differential-geometry operators for triangulated 2-manifolds. Visualization and mathematics, 3(2):52–58, 2002.
  • [46] Rongjie Lai, Jiang Liang, and Hongkai Zhao. A local mesh method for solving pdes on point clouds. Inverse Prob. and Imaging, 7(3):737–755, 2013.
  • [47] Rongjie Lai and Stanley Osher. A splitting method for orthogonality constrained problems. Journal of Scientific Computing, 58(2):431–449, 2014.
  • [48] Yanfei You Weiqiang Chen, Jui JI. An augmented lagrangian method for l1-regularized prolems with orthogonality constrains. SIAM Journal on Scientific Computing, 38(7), 2016.
  • [49] Yu Wang, Wotao Yin, and Jinshan Zeng. Global convergence of admm in nonconvex nonsmooth optimization. arXiv preprint arXiv:1511.06324, 2015.
  • [50] Artiom Kovnatsky, Klaus Glashoff, and Michael M Bronstein. Madmm: a generic algorithm for non-smooth optimization on manifolds. In European Conference on Computer Vision, pages 680–696. Springer, 2016.
  • [51] Hedy Attouch, Jérôme Bolte, and Benar Fux Svaiter. Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward–backward splitting, and regularized gauss–seidel methods. Mathematical Programming, 137(1-2):91–129, 2013.
  • [52] Jérôme Bolte, Shoham Sabach, and Marc Teboulle. Proximal alternating linearized minimization for nonconvex and nonsmooth problems. Mathematical Programming, 146(1-2):459–494, 2014.
  • [53] Jonathan Barzilai and Jonathan M Borwein. Two-point step size gradient methods. IMA journal of numerical analysis, 8(1):141–148, 1988.