Matching problems, seeking some useful correspondence between two shapes or, more generally, discrete metric spaces, are central in computer graphics and vision. Matching problems are often modeled as optimization of a quadratic energy over permutations. Global optimization and approximation of such problems is known to be NP-hard [Loiola et al. 2007].
A common strategy for dealing with the computational hardness of matching problems is replacing the original optimization problem with an easier, similar problem that can be solved globally and efficiently. Perhaps the two most common scalable relaxations for such problems are the spectral relaxation for energies with non-negative entries [Leordeanu and Hebert 2005, Feng et al. 2013], and the doubly stochastic (DS) relaxation for convex energies [Aflalo et al. 2015, Fiori and Sapiro 2015].
Our work is motivated by the recent work of [Kezurer et al. 2015] who proposed a semi-definite programming (SDP) relaxation which is provably stronger than both spectral and DS relaxations. The obtained relaxation was shown empirically to be extremely tight, achieving the global ground truth in most experiments presented. However, a major limitation was the computational cost of solving a semi-definite program with variables. Accordingly in this paper we pursue the following question:
Question: Is it possible to construct a relaxation which is stronger than the spectral and DS relaxations, without compromising efficiency?
We give an affirmative answer to this question and show that by correctly combining the spectral and DS relaxations in the spirit of [Fogel et al. 2013] we obtain a relaxation which is provably tighter than both, and is in fact in a suitable sense exactly the intersection of both relaxations. We name this relaxation DS+. Moreover, we observe that a refined spectral analysis leads to a significant improvement to this relaxation and a provably tighter quadratic program we name DS++. This relaxation enjoys the same scalability as DS and DS+ as all three are quadratic programs with variables and the same number of constraints. Additional time efficiency is provided by specialized solvers for the DS relaxation such as [Solomon et al. 2016]. We note that DS++ is still less tight than the final expensive and accurate relaxation of [Kezurer et al. 2015] yet strikes a balance between tightness and computational complexity. The hierarchy between the relaxations is illustrated in the inset and proven in section 4.
Since DS++ is a relaxation, it is not guaranteed to output an integer solution (i.e., a permutation). To obtain a feasible permutation we propose a homotopy-type method, in the spirit of [Ogier and Beyer 1990, Zaslavskiy et al. 2009]. This method continuously deforms the energy functional from convex to concave, is guaranteed to produce an integer-solution and in practice outperforms standard Euclidean projection techniques. Essentially it provides a strategy for finding a local minima for the original non-convex problem using a good initial guess obtained from the convex relaxation.
Our algorithm is very flexible and can be applied to both convex and non-convex energies (in contrast with DS), and to energies combining quadratic and linear terms (in contrast with the spectral relaxation, which also requires energies with non-negative entries). It can also be easily modified to allow for additional linear constraints, injective and partial matching, and solving quadratic optimization problems over the doubly stochastic matrices. We present experiments demonstrating the effectiveness of our method in comparison to random initializations of the non-convex problem, spectral, DS, and DS+ relaxations, as well as lifted linear-programming relaxations.
We have tested our algorithm on three applications: (i) non-rigid matching; (ii) image arrangements; and (iii) coarse-to-fine matching. Comparison to state-of-the-art algorithms for these applications shows that our algorithm produces favorable results in comparable speed.
Our contributions in this paper are threefold:
We identify the optimal initial convex and concave relaxation.
We show, both theoretically and experimentally that the proposed algorithm is more accurate than other popular contemporary methods. We believe that establishing a hierarchy between the various relaxation methods for quadratic matching is crucial both for applications, and for pushing forward the algorithmic state of the art, developing stronger optimization algorithms in the future.
Lastly, we build a simple end-to-end algorithm utilizing recent advances in optimization over the doubly-stochastic matrices to provide a scalable yet accurate algorithm for quadratic matching.
2 Previous work
Many works in computer vision and graphics model correspondence problems as quadratic optimization problems over permutation matrices. In many cases these problems emerge as discretizations of isometry-invariant distances between shapes[Mémoli and Sapiro 2005, Mémoli 2011] . We focus here on the different methods to approximately solve these computationally hard problems.
The spectral relaxation for correspondence problems in computer vision has been introduced in [Leordeanu and Hebert 2005] and has since become very popular in both computer vision and computer graphics, e.g., [De Aguiar et al. 2008, Liu et al. 2009, Feng et al. 2013, Shao et al. 2013]
. This method replaces the requirement for permutation matrices with a single constraint on the Frobenious norm of the matrices to obtain a maximal eigenvalue problem. It requires energies with positive entries to ensure the obtained solution is positive. This relaxation is scalable but is not a very tight approximation of the original problem. A related relaxation appears in[Rodola et al. 2012], where the variable is constrained to be non-negative with . This optimization problem is generally non-convex, but the authors suggest a method for locally minimizing this energy to obtain a sparse correspondence.
An alternative approach relaxes the set of permutations to its convex hull of doubly stochastic matrices [Schellewald et al. 2001]. When the quadratic objective is convex, this results in a convex optimization problem (quadratic program) which can be minimized globally, although the minimum may differ from the global minima of the original problem. [Solomon et al. 2012] argue for the usefulness of the fuzzy maps obtained from the relaxation. For example, for symmetric shapes fuzzy maps can encode all symmetries of the shape.
[Aflalo et al. 2015] shows that for the convex graph matching energy the DS relaxation is equivalent to the original problem for generic asymmetric and isomorphic graphs. These results are strengthened in [Fiori and Sapiro 2015]. However when noise is present the relaxations of the convex graph matching energy will generally not be equivalent to the original problem [Lyzinski et al. 2016] . Additionally, for concave energies the DS relaxation is always equivalent to the original problem [Gee and Prager 1994], since minima of concave energies are obtained at extreme points. The challenge for non-convex energies is that global optimization over DS matrices is not tractable.
To achieve good initialization for local minimization of such problems, [Ogier and Beyer 1990, Gee and Prager 1994, Zaslavskiy et al. 2009] suggest to minimize a sequence of energies which gradually vary from a convex energy to an equivalent concave energy . In this paper we adopt this strategy to obtain an integer solution, and improve upon it by identifying the optimal convex and concave energies from within the energies .
The authors of [Fogel et al. 2013, Fogel et al. 2015] show that the DS relaxation can be made more accurate by adding a concave penalty of the form to the objective. To ensure the objective remains convex they suggest to choose to be the minimial eigenvalue of the quadratic objective. We improve upon this choice by choosing to be the minimial eigenvalue over the doubly stochastic subspace, leading to a provably tighter relaxation. The practical advantage of our choice (DS++) versus Fogel’s choice (DS+) is significant in terms of the relaxation accuracy as demonstrated later on. The observation that this choice suffices to ensure convexity has been made in the convergence proof of the softassign algorithm [Rangarajan et al. 1997].
Optimization of DS relaxation
Specialized methods for minimization of linear energies over DS matrices [Kosowsky and Yuille 1994, Cuturi 2013, Benamou et al. 2015, Solomon et al. 2015] using entropic regularization and the Sinkhorn algorithm are considerably more efficient than standard linear program solvers for this class of problems. Motivated by this, [Rangarajan et al. 1996] propose an algorithm for globally minimizing quadratic energies over doubly stochastic matrices by iteratively minimizing regularized linear energies using Sinkhorn type algorithms. For the optimization in this paper we applied [Solomon et al. 2016] who offer a different algorithm for locally minimizing the Gromov-Wasserstein distance by iteratively solving regularized linear programs. The advantage of the latter algorithm over the former algorithm is its certified convergence to a critical point when applied to non-convex quadratic energies.
Other convex relaxations
Stronger relaxations than the DS relaxation can be obtained by lifting methods which add auxiliary variables representing quadratic monomials in the original variables. This enables adding additional convex constraints on the lifted variables. A disadvantage of these methods is the large number of variables which leads to poor scalability. [Kezurer et al. 2015] propose in an SDP relaxation in the spirit of [Zhao et al. 1998], which is shown to be stronger than both DS (for convex objective) and spectral relaxations, and in practice often achieves the global minimum of the original problem. However, it is only tractable for up to fifteen points. [Chen and Koltun 2015] use a lifted linear program relaxation in the spirit of [Werner 2007, Adams and Johnson 1994]. To deal with scalability issues they use Markov random field techniques [Kolmogorov 2006] to approximate the solution of their linear programming relaxation.
Several works aim at globally solving the quadratic assignment problem using combinatorial methods such as branch and bound. According to a recent survey [Loiola et al. 2007] these methods are not tractable for graphs with more than points. Branch and bound methods are also in need of convex relaxation to achieve lower bounds for the optimization problem. [Anstreicher and Brixius 2001] provide a quadratic programming relaxation for the quadratic assignment problem which provably achieves better lower bounds than a competing spectral relaxation using a method which combines spectral, linear, and DS relaxations. Improved lower bounds can be obtained using second order cone programming [Xia 2008] and semi-definite programming [Ding and Wolkowicz 2009] in variables. All the relaxations above use the specific structure of the quadratic assignment problem while our relaxation is applicable to general quadratic objectives which do not carry this structure and are very common in computer graphics. For example, most of the correspondence energies formulated below and considered in this paper cannot be formulated using the quadratic assignment energy.
Other approaches for shape matching
A similar approach to the quadratic optimization approach is the functional map method (e.g., [Ovsjanikov et al. 2012]) which solves a quadratic optimization problem over permutations and rotation matrices, typically using high-dimensional ICP provided with some reasonable initialization. Recently [Maron et al. 2016] proposed an SDP relaxation for this problem with considerably improved scalability with respect to standard SDP relaxations.
Supervised learning techniques have been successfully applied for matching specific classes of shapes in [Rodolà et al. 2014, Masci et al. 2015, Zuffi and Black 2015, Wei et al. 2016]. A different approach for matching near isometric shapes is searching for a mapping in the low dimensional space of conformal maps which contains the space of isometric maps [Lipman and Funkhouser 2009, Zeng et al. 2010, Kim et al. 2011]. More information on shape matching can be found in shape matching surveys such as [Van Kaick et al. 2011].
Quadratic optimization problems over the set of permutation matrices arise in many contexts. Our main motivating example is the problem of finding correspondences between two metric spaces (e.g., shapes) and which are related by a perfect or an approximate isometry. This problem can be modeled by uniformly sampling the spaces to obtain and , and then finding the permutation which minimizes an energy of the form
Here is some penalty on deviation from isometry: If the points correspond to the points (resp.), then the distances between the pair on the source shape and the pair on the target shape should be similar. Therefore we choose
where is some function penalizing for deviation from the set . Several different choices of exist in the literature.
The linear term is sometimes used to aid the correspondence task by encouraging correspondences between points with similar isometric-invariant descriptors.
Our goal is to solve quadratic optimization problems over the set of permutations as formulated in (1). Denoting the column stack of permutations
by the vector
leads to a more convenient phrasing of (1):
This optimization problem is non-convex for two reasons. The first is the non-convexity of (as a discrete set of matrices), and the second is that is often non-convex (if is not positive-definite). As global minimization of (3) is NP-hard [Loiola et al. 2007] we will be satisfied with obtaining a good approximation to the global solution of (3) using a scalable optimization algorithm. We do this by means of a convex relaxation coupled with a suitable projection algorithm for achieving integer solutions.
3.1 Convex relaxation
We formulate our convex relaxation by first considering a one-parameter family of equivalent formulations to (3): observe that for any permutation matrix we have that . It follows that all energies of the form
coincide when restricted to the set of permutations. Therefore, replacing the energy in (3) with provides a one-parameter family of equivalent formulations. For some choices of the energy in these formulations is convex, for example, for any , where is the minimal eigenvalue of .
For each such equivalent formulation we consider its doubly stochastic relaxation. That is, replacing the permutation constraint (3b) with its convex-hull, the set of doubly-stochastic matrices:
Our goal is to pick a relaxed formulation (i.e., choose an ) that provides the best lower bound to the global minimum of the original problem (3). For that end we need to consider values of that make convex and consequently turn (5) into a convex program that provide a lower bound to the global minimum of (3).
Among all the convex programs described above we would like to choose the one which provides the tightest lower bound. We will use the following simple lemma proved in Appendix A:
For all doubly stochastic we have when , and if and only if is a permutation.
An immediate conclusion from this lemma is that and so the best lower bound will be provided by the largest value of for which is convex.
See for example the inset illustrating the energy graphs for different values for a toy-example: in red - the graph of the original (non-convex) energy with ; in blue the energy with ; and in green . Note that the green graph lies above the blue graph and all graphs coincide on the corners (i.e., at the permutations). Since the higher the energy graph the better lower bound we achieve it is desirable to take the maximal that still provides a convex program in (5). In the inset the green and blue points indicate the solution of the respective relaxed problems; in this case the green point is much closer to the sought after solution, i.e., the lower-left corner.
To summarize the above discussion: the desired is the maximal value for which is a convex function. As noted above choosing in the spirit of [Fogel et al. 2013], leads to a convex problem which we denote by DS+. However this is in fact not the maximal value in general. To find the maximal we can utilize the fact that is constrained to the affine space defined by the constraints (5b): We parameterize the affine space as , where is some permutation, is any parameterization satisfying , and . Plugging this into provides a quadratic function in of the form
where is some affine function of . It follows that (5) will be convex iff is positive semi-definite. The largest possible fulfilling this condition is the minimal eigenvalue of which we denote by . Thus our convex relaxation which we name DS++ boils down to minimizing (5) with .
We now describe our method for projecting the solution of our relaxation onto the set of permutations. This method is inspired by the ”convex to concave” method from [Ogier and Beyer 1990, Gee and Prager 1994, Zaslavskiy et al. 2009], but also improves upon these works by identifying the correct interval on which the convex to concave procedure should be applied as we now describe.
Lemma 1 tells us that the global optimum of over the doubly stochastic matrices provides an increasingly better approximation of the global optimum of the original problem (3) as we keep increasing even beyond the convex regime, that is . In fact, it turns out that if is chosen large enough so that is strictly concave, then the global optima of (5) and the global optima of the original problem over permutations are identical. This is because the (local or global) minima of strictly concave functions on a compact convex set are always obtained at the extreme points of the set. In our case, the permutations are these extreme points.
This leads to a natural approach to approximate the global optimum of (3): Solve the above convex problem with and then start increasing until an integer solution is found. We choose a finite sequence , where and is strictly concave. We begin by solving (5) with which is exactly the convex relaxation described above and obtain a minimizer . We then iteratively locally minimize (5) with using as an initialization the previous solution . The reasoning behind this strategy is that when and are close a good solution for the latter should provide a good initialization for the former, so that at the end of the process we obtain a good initial guess for the minimization of , which is equivalent to the original integer program. We stress that although the obtained solution may only be a local minimum, it will necessarily be a permutation.
To ensure that is strictly concave we can choose any larger than , which analogously to is defined as the largest eigenvalue of . In practice we select which in the experiments we conducted is sufficient for obtaining integer solutions. We then took by uniformly sampling where unless stated otherwise we used ten samplings (). Throughout the paper we will use the term DS++ algorithm to refer to our complete method (relaxation+projection) and DS++ or DS++ relaxation to refer only to the relaxation component.
Figure 2 shows the correspondences (encoded in a specific row of ) obtained at different stages of the projection procedure when running our algorithm on the FAUST dataset [Bogo et al. 2014] as described below. The figure shows the correspondences obtained from optimizing for .
4 Comparison with other relaxations
The purpose of this section is to theoretically compare our relaxation with common competing relaxations. We prove
The DS++ relaxation is more accurate than DS+, which in turn is more accurate than the spectral and doubly stochastic relaxation.
The SDP relaxation of [Kezurer et al. 2015] is more accurate than all the relaxations mentioned above.
Our strategy for proving this claim is formulating all relaxations in a unified framework, using the SDP lifting technique in [Kezurer et al. 2015], that in turn readily enables comparison of the different relaxations.
The first step in constructing SDP relaxations is transforming the original problem (3) into an equivalent optimization problem in a higher dimension.The higher dimension problem is formulated over the set:
Using the identity
we obtain an equivalent formulation to (3):
SDP relaxations are constructed by relaxing the constraint using linear constraints on and the semi-definite constraint .
[Kezurer et al. 2015] showed that the spectral and doubly stochastic relaxations are equivalent to the following SDP relaxations:
We note that the spectral relaxation is applicable only when , and the DS relaxation is tractable only when the objective is convex, i.e., . The equivalence holds under these assumptions.
Given this new formulation of spectral and DS, an immediate method for improving both relaxations is considering the Intersection-SDP, obtained by enforcing the constraints from both and . The relaxation can be further improved by adding additional linear constraints on . This is the strategy followed by [Kezurer et al. 2015] to achieve their final tight relaxation which is presented in Eq. (9) in the appendix. The main limitation of this approach is its prohibitive computational price resulting from solving SDPs with variables, in strong contrast to the original formulation of spectral and DS that uses only variables (i.e., the permutation ). This naturally leads to the research question we posed in the introduction, which we can now state in more detail:
Question: Is it possible to construct an SDP relaxation which is stronger than and , and yet is equivalent to a tractable and scalable optimization problem with variables?
We answer this question affirmatively by showing that the Intersection-SDP is in fact equivalent to DS+. Additionally DS++ is equivalent to a stronger SDP relaxation which includes all constraints from the Intersection-SDP, as well as the following additional constraints: Let us write the linear equality constraints appearing in the definition of the DS matrices (i.e., (5b) ) in the form . Then any in particular satisfies and therefore also:
Adding these constraints to the Intersection-SDP we obtain
Theorem 1 now follows from:
We prove the lemma in the appendix.
5 Implementation details
Optimization of (5) can be done using general purpose non-convex solvers such as Matlab’s fmincon, or solvers for convex and non-convex quadratic programs. We opted for the recent method of Solomon et al., Justin that introduced a specialized scalable solver for local minimization of regularized quadratic functionals over the set of doubly stochastic matrices.
The algorithm of [Solomon et al. 2016] is based on an efficient algorithm for optimizing the KL divergence
is the column stack of a doubly stochastic matrixand is some fixed positive vector. The solution for the KKT equations of this problem can be obtained analytically for , up to scaling of the rows and columns, which is performed by the efficient Sinkhorn algorithm. See [Cuturi 2013] for more details.
The algorithm of [Solomon et al. 2016] minimizes quadratic functionals (where in our case ) over doubly stochastic matrices by iteratively optimizing KL-divergence problems. First the original quadratic functional is regularized by adding a barrier function keeping the entries of away from zero to obtain a new functional
The parameter is chosen to be some small positive number so that its effect on the functional is small. We then define so that
We then optimize iteratively: In iteration , is held fixed at its previous value , and an additional term is added penalizing large deviations of from . More precisely, is defined to be the minimizer of
where denotes entry-wise multiplication of vectors. For small enough values of , [Solomon et al. 2016] prove that the algorithm converges to a local minimun of .
In our implementation we use . We choose the smallest possible so that all entries of the argument of the exponent in the definition of are in . This choice is motivated by the requirement of choosing small coupled with the breakdown of matlab’s exponent function at around . Note that this choice requires to update at each iteration. We find that with this choice of the regularization term has little effect on the energy and we obtain final solutions which are close to being permutations. To achieve a perfect permutation we project the final solution using the projection. The projection is computed by minimizing a linear program as described, e.g., in [Zaslavskiy et al. 2009].
We compute and by solving two maximal magnitude eigenvalue problems: We first solve for the maximal magnitude eigenvalue of . If this eigenvalue is positive then it is equal to . We can then find by translating our matrix by to obtain a positive-definite matrix whose maximal eigenvalue is related to the minimal eigenvalue of the original matrix via .
If the solution of the first maximal magnitude problem is negative then this eigenvalue is , and we can use a process similar to the one described above to obtain .
Solving maximal magnitude eigenvalue problems requires repeated multiplication of vectors by the matrix , where and . If is sparse, computing can become a computational bottleneck. To avoid this problem, we note that has the same maximal eigenvalue as the matrix and so compute the maximal eigenvalue of the latter matrix. The advantage of this is that multiplication by the matrix can be computed efficiently:
Since is the orthogonal projection onto , we can use the identity where is the projection onto the orthogonal complement of . The orthogonal complement is of dimension and therefore where .
We solve the maximal magnitude eigenvalue problems using Matlab’s function eigs.
Our method can be applied with minor changes to injective matching. The input of injective matching is points sampled from the source shape and points sampled from the target shape , and the goal is to match the points from injectively to a subset of of size .
Matrices representing injective matching have entries in , and have a unique unit entry in each row, and at most one unit entry in each column. This set can be relaxed using the constraints:
We now add a row with positive entries to the variable matrix to obtain a matrix . The original matrix satisfies the injective constraints described above if satisfies
These constraints are identical to the constraint defining DS, up to the value of the marginals which have no affect on our algorithm. As a result we can solve injective matching problems without any modification of our framework.
The input of partial matching is points sampled from the source and target shape, and the goal is to match points from injectively to a subset of of size . We do not pursue this problem in this paper as we did not find partial matching necessary for our applications. However we believe our framework can be applied to such problems by adding a row and column to the matching matrix .
Adding linear constraints
Modeling of different matching problems can suggest adding additional linear constraints on that can be added directly to our optimization technique. Additional linear equality constraints further decrease the dimension of the affine space is constrained to and as a result make the interval smaller, leading to more accurate optimization. We note however that incorporating linear constraints into the optimization method of [Solomon et al. 2016] is not straightforward.
Upsampling refers to the task of interpolating correspondences betweenpoints sampled from source and target metric spaces to a match between a finer sampling of source points and target points. We suggest two strategies for this problem: Limited support interpolation and greedy interpolation.
Limited support interpolation uses the initially matched points to rule out correspondences between the finely sampled points. The method of ruling out correspondences is discussed in the Appendix. We enforce the obtained sparsity pattern by writing , where the first matrix is zero in all forbidden entries and the second is zero in all permissible entries. We then minimize the original energy only on the permissible entries, and add a quadratic penalty for the forbidden entries. That is, we minimize
choosing some large . The sparsity of enables minimizing this energy for large for which minimizing the original energy is intractable.
When are large we use greedy interpolation. We match each source point separately. We do this by optimizing over correspondences between source points and target points, where the points are the known points and the point . Since there are only such correspondences optimization can be performed globally by checking all possible correspondences.
Optimization over doubly stochastic matrices
Our main focus was on optimization problems over permutations. However in certain cases the requested output from the optimization algorithm may be a doubly stochastic matrix and not a permutation. When the energy is non convex this still remains a non-convex problem. For such optimization problems our method can be applied by taking samples from the interval , since minimization of (5) with is the problem to be solved while minimization of (5) with forces a permutation solution.
In this section we evaluate our algorithm and compare its performance with relevant state of the art algorithms. We ran all experiments on the 100 mesh pairs from the FAUST dataset [Bogo et al. 2014] which were used in the evaluation protocol of [Chen and Koltun 2015].
Comparison with [Solomon et al. 2016]
In figure 5(b) we compare our method for minimizing non-convex functionals with the local minimization algorithm of [Solomon et al. 2016]. Since this method is aimed at solving non-convex functionals over doubly-stochastic matrices, we run our algorithm using samples in as explained in Section 6. We sample points from each mesh using farthest point sampling [Eldar et al. 1997], and optimize the Gromov-Wasserstein (GW) functional advocated in [Solomon et al. 2016], which amounts to choosing from (2) to be . As local minimization depends on initialization we locally minimize times per mesh pair, using different random initializations. The initializations are obtained by randomly generating a positive matrix in
with uniform distribution, and projecting the result onto the doubly stochastic matrices using the Sinkhorn algorithm. As can be seen in the figure our algorithm, using only ten iterations, was more accurate than all the local minima found using random initializations. As a baseline for comparison we note that the difference in energy between randomly drawn permutations and our solution was around, while the difference in energy shown in the graph is around . Figure 6 visualizes the advantages of the fuzzy maps obtained by our algorithm in this experiment over the best of the random maps generated by [Solomon et al. 2016].
In figure 5(a) we examine how the result obtained from our projection method is influenced by the number of points sampled from . We compared the behavior of our relaxation with several different choices of as well as with the standard projection onto the set of permutations. As expected, our projection is always better that the projection, and the projection improves as the number of samples is increased.
Comparison with other relaxations
We compare our method with other relaxation based techniques. In figure 5 (c)-(d) we compare our relaxation with the spectral relaxation, the DS+ relaxation, and the SDP relaxation of [Kezurer et al. 2015]. In this experiment the energy we use is non-convex so DS is not applicable.
We choose the parameter . For the minimization we used all four relaxations, obtaining a lower bound for the optimal value, Figure 5 (c). We then projected the solutions obtained onto the set of permutations, thus obtaining an upper bound, Figure 5 (d). For methods other than the DS++ algorithm we used the projection. In all experiments the upper and lower bounds provided by the SDP relaxation of [Kezurer et al. 2015] were identical, thus proving that the SDP relaxation found the globally optimal solution. Additionally, in all experiments the upper bound and lower bound provided by our relaxation were superior to those provided by the spectral method, and our projection attained the global minimum in approximately of the experiments in contrast to obtained by the projection of the spectral method. The differences between the spectral relaxation and the stronger DS+ relaxation were found to be negligible.
In figure 5(e)-(f) we perform the same experiment, but now we minimize the convex graph matching functional from [Aflalo et al. 2015] for which the classical DS relaxation is applicable. Here again the ground truth is achieved by the SDP relaxation. Our relaxation can be seen to modestly improve the lower bound obtained by the classical DS relaxation, while our projection method substantially improves upon the standard projection.
We have tested our method for three applications: non-rigid shape matching, image arrangement, and coarse-to-fine matching.
We evaluated the performance of our algorithm for non-rigid matching on the FAUST dataset [Bogo et al. 2014]. We compared to [Chen and Koltun 2015] which demonstrated superb state of the art results on this dataset (for non learning-based methods). For a fair comparison we used an identical pipeline to [Chen and Koltun 2015], including their isometric energy modeling and extrinsic regularization term. We first use the DS++ algorithm to match points, then upsampled to using limited support interpolation and to using greedy interpolation, as described in Section 6; the final point resolution is as in [Chen and Koltun 2015].
Figure 7 depicts the results of the DS++ algorithm and [Chen and Koltun 2015]. As can be read from the graphs, our algorithm compares favorably in both the inter and intra class matching scenarios in terms of cumulative error distribution and average error. These results are consistent for both the convex relaxation part (top row) and the upsampled final map (middle row); The graphs show our results both with our upsampling as described above (denoted by DS++(1)) and the results of combining our relaxation with the upsampling of [Chen and Koltun 2015] (DS++(2)). We find DS++(1) to be better on the inter class, and DS++(2) is marginally better on the intra class. The error is calculated on a set of 52 ground truth points in each mesh as in [Maron et al. 2016]. Figure 1 (left), and 3 show typical examples of maps computed using the DS++ algorithm in this experiment.
|dataset||feature||improvement||rand average||Fried mean||our mean||functional||swaps?||grid size|
|SUN dataset||deep feature object||55.97%||0.433||0.345||0.295||Fried||no||14|
|SUN dataset||deep feature object||6.31%||0.433||0.300||0.292||Fried||yes||14|
|LFW||deep feature face||50.70%||0.422||0.355||0.320||Fried||no||14|
|LFW||deep feature face||2.81%||0.422||0.321||0.318||Fried||yes||14|
The task of arranging image collections in a grid has received increasing attention in recent years [Quadrianto et al. 2009, Strong and Gong 2014, Fried et al. 2015, Carrizosa et al. 2016]. Image arrangement is an instance of metric matching: the first metric space is the collection of images and a dissimilarity measure defined between pairs of images; and the second, target metric space is a 2D grid (generally, a graph) with its natural Euclidean metric.
[Fried et al. 2015] suggested an energy functional for generating image arrangements, which are represented by a permutation matrix . Their choice of energy functional was supported by a user study. This energy functional is:
where are the distance measures between images and grid points respectively, and is the unknown scale factor between the two metric spaces. [Fried et al. 2015] suggested a two step algorithm to approximate the minimizer of the above energy over the set of permutations: The first step is a dimensionality reduction, and the second is linear assignment to a grid according to Euclidean distances. Fried et al., demonstrated significant quantitative improvement over previous state of the art methods.
We perform image arrangement by using an alternative method for optimizing the energy (8). We fix so that the mean of and are the same, which leads to a quadratic matching energy which we optimize over permutations using the DS++ algorithm.
Table 1 summarizes quantitative comparison of the DS++ algorithm and [Fried et al. 2015] on a collection of different image sets and dissimilarity measures. Each row shows the mean energies over 100 experiments of Fried et al., DS++, and random assignments which provide a baseline for comparison; in each experiment we randomized a subset of images from the relevant set of images and generated image arrangements using the two methods. [Fried et al. 2015] suggested an optional post processing step in which random swaps are generated and applied in case they reduce the energy; our experiment measures the performance of both algorithms with and without random swaps. The first set of experiments tries to arrange random colors in a grid. The second set of experiments uses the mean image color for images form the SUN database [Xiao et al. 2010]
. The third set uses the last layer of a deep neural network trained for object recognition[Chatfield et al. 2014] as image features, again for images in the SUN dataset. The fourth set of experiments organizes randomly sampled images from the Labeled Faces in the Wild (LFW) dataset [Huang et al. 2007] according to similar deep features taken from the net trained by [Parkhi et al. 2015]. For the last experiment, we rendered three 3D models from the SHREC07 dataset [Giorgi et al. 2007] from various illumination directions and ordered them according to the raw distance between pairs of images.
Our algorithm outperformed [Fried et al. 2015] in all experiments (in some cases our algorithm achieved an improvement of more than 50%). Figures 1, 4 and 11 show some image arrangements from these experiments. Note, for example, how similar faces are clustered together in Figure 11 (a), and similar objects are clustered in Figure 1 (right). Further note how the image arrangement in Figure 11 (b) nicely recovered the two dimensional parameter of the lighting direction, where Figure 8 shows the input random lighting directions renderings of a 3D models.
We consider the problem of matching two shapes and using sparse correspondences specified by the user. User input can be especially helpful for highly non-isometric matching problems where semantic knowledge is often necessary for achieving high quality correspondences. The inset shows such example where three points (indicated by colored circles) are used to infer correspondences between a horse and a giraffe.
We assume the user supplied a sparse set of point correspondences, , , and the goal is to complete this set to a full correspondence set between the shapes and . Our general strategy is to use a linear term to enforce the user supplied constraints, and a quadratic term to encourage maps with low distortion.
For a quadratic term we propose a ”log-GW” functional. This functional amounts to choosing for the definition of in (2), where is a metric on defined by
This metric punishes for high relative distortion between , and thus is more suitable for our cause than the standard Euclidean metric used for the GW functional.
As a linear term we propose
The first summand from the left penalizes matchings which violate the known correspondences, while the second summand penalizes matchings which cause high distortion of distances to the user supplied points. The parameter controls the strength of the linear term. In our experiments we chose , where is the spectral norm of the quadratic form..
We applied the algorithm for coarse-to-fine matching on the SHREC dataset [Giorgi et al. 2007], using of the labeled ground truth points and evaluating the error of the obtained correspondence on the remaining points. The number of labeled points is class dependent and varies between and .
Representative results are shown in Figure 10. The graph on the left hand side of Figure 9 shows that our algorithm is able to use this minimal user supplied information to obtain significantly better results than those obtained by the unaided BIM algorithm [Kim et al. 2011]. The graph compares the algorithms in terms of cumulative error distribution over all 218 SHREC pairs for which BIM results are available.
The graph on the right hand side of Figure 9 shows our results outperform the algorithm presented in [Solomon et al. 2016] for matching aided by user supplied correspondences. Both algorithms were supplied with 6 ground truth points. We ran the algorithm of [Solomon et al. 2016] matching points (we take since [Solomon et al. 2016] does not support injective matching) and using the maximal-coordinate projection they chose to achieve a permutation solution. These results are denoted by (2) in the graph. However we find that better results are achieved when matching only points, and when using the projection. These results are denoted by (1) in the graph.
Typical running times of our optimization algorithm for the energy of [Chen and Koltun 2015] matching points takes 6 seconds; takes 26 seconds; and , points takes around 2 minutes (130 seconds). The precomputation of and with these parameters requires around 15 seconds, the projection requires 5 seconds, and the remaining time is required for our optimization algorithm.
Parameter values of (as well as ) were used in the image arrangement task from Section 8, and parameter values were used for our results on the FAUST dataset. For the latter application we also upsampled to using limited support interpolation and then upsampled to using greedy interpolation as described in Section 6.Limited support interpolation required seconds and greedy interpolation required seconds. The total time for this application is around minutes.
The efficiency of our algorithm significantly improves if the product of the quadratic form’s matrix with a vector can be computed efficiently. This is illustrated by the fact that optimization of the sparse functional we construct for the task of resolution improvement with takes similar time as optimization of the non-sparse functional of [Chen and Koltun 2015] with .
Another case where the product can be computed efficiently is the GW or log-GW energy. In both cases the product can be computed by multiplication of matrices of size (see [Solomon et al. 2016] for the derivation), thus using operations instead of the operations necessary for general . Using this energy, matching points to points takes only 12 seconds, matching points to points takes 22 seconds, matching to takes 82 seconds, and for we require around six minutes (368 seconds).
The efficiency of our algorithm depends linearly on . Minimizing the GW energy with using sample takes 12 seconds, approximately half of the time needed when using samples. These parameters were used for our results on matching with user input described in Section 8. For this task we also used greedy interpolation to obtain full maps between the shapes, which required an additional 18 seconds. Overall this application required around half a minute.
Our algorithm was implemented on Matlab. All running times were computed using an Intel i7-3970X CPU 3.50 GHz.
We have introduced the DS++ algorithm for approximating the minimizer of a general quadratic energy over the set of permutations. Our algorithms contains two components: (i) A quadratic program convex relaxation that is guaranteed to be better than the prevalent doubly stochastic and spectral relaxations; and (ii) A projection procedure that continuously changes the energy to recover a locally optimal permutation, using the convex relaxation as an initialization. We have used recent progress in optimal transport to build an efficient implementation to the algorithm.
The main limitation of our algorithm is that it does not achieve the global minima of the energies optimized. Partially this is unavoidable and due to the computational hardness of our problem. However the experimental results in Figure 5 show that accuracy can be improved by the SDP method of [Kezurer et al. 2015] which while computationally demanding, can still be solved in polynomial time. Our future goal is to search for relaxations whose accuracy is close to those of [Kezurer et al. 2015] but which are also fairly scalable. One concrete direction of research could be finding the ‘best’ quadratic programming relaxation in variables.
- [Adams and Johnson 1994] Adams, W. P., and Johnson, T. A. 1994. Improved linear programming-based lower bounds for the quadratic assignment problem. DIMACS series in discrete mathematics and theoretical computer science 16, 43–75.
- [Aflalo et al. 2015] Aflalo, Y., Bronstein, A., and Kimmel, R. 2015. On convex relaxation of graph isomorphism. Proceedings of the National Academy of Sciences 112, 10 (Mar.), 2942–2947.
- [Anstreicher and Brixius 2001] Anstreicher, K. M., and Brixius, N. W. 2001. A new bound for the quadratic assignment problem based on convex quadratic programming. Mathematical Programming 89, 3, 341–357.
- [Benamou et al. 2015] Benamou, J.-D., Carlier, G., Cuturi, M., Nenna, L., and Peyré, G. 2015. Iterative bregman projections for regularized transportation problems. SIAM Journal on Scientific Computing 37, 2, A1111–A1138.
[Bogo et al. 2014]
Bogo, F., Romero, J., Loper, M., and Black, M. J.
FAUST: Dataset and evaluation for 3D mesh registration.
Computer Vision and Pattern Recognition (CVPR), 2014 IEEE Conference on, IEEE, 3794–3801.
- [Carrizosa et al. 2016] Carrizosa, E., Guerrero, V., and Morales, D. R. 2016. Visualizing proportions and dissimilarities by space-filling maps: a large neighborhood search approach. Computers & Operations Research.
- [Chatfield et al. 2014] Chatfield, K., Simonyan, K., Vedaldi, A., and Zisserman, A. 2014. Return of the devil in the details: Delving deep into convolutional nets. arXiv preprint arXiv:1405.3531.
- [Chen and Koltun 2015] Chen, Q., and Koltun, V. 2015. Robust nonrigid registration by convex optimization. In Proceedings of the IEEE International Conference on Computer Vision, 2039–2047.
- [Cuturi 2013] Cuturi, M. 2013. Sinkhorn distances: Lightspeed computation of optimal transport. In Advances in Neural Information Processing Systems, 2292–2300.
- [De Aguiar et al. 2008] De Aguiar, E., Stoll, C., Theobalt, C., Ahmed, N., Seidel, H.-P., and Thrun, S. 2008. Performance capture from sparse multi-view video. In ACM Transactions on Graphics (TOG), vol. 27, ACM, 98.
- [Ding and Wolkowicz 2009] Ding, Y., and Wolkowicz, H. 2009. A low-dimensional semidefinite relaxation for the quadratic assignment problem. Mathematics of Operations Research 34, 4, 1008–1022.
- [Dym and Lipman 2016] Dym, N., and Lipman, Y. 2016. Exact recovery with symmetries for procrustes matching. arXiv preprint arXiv:1606.01548.
- [Eldar et al. 1997] Eldar, Y., Lindenbaum, M., Porat, M., and Zeevi, Y. Y. 1997. The farthest point strategy for progressive image sampling. Image Processing, IEEE Transactions on 6, 9, 1305–1315.
- [Feng et al. 2013] Feng, W., Huang, J., Ju, T., and Bao, H. 2013. Feature correspondences using morse smale complex. The Visual Computer 29, 1, 53–67.
- [Fiori and Sapiro 2015] Fiori, M., and Sapiro, G. 2015. On spectral properties for graph matching and graph isomorphism problems. Information and Inference 4, 1, 63–76.
- [Fogel et al. 2013] Fogel, F., Jenatton, R., Bach, F., and d’Aspremont, A. 2013. Convex relaxations for permutation problems. In Advances in Neural Information Processing Systems, 1016–1024.
- [Fogel et al. 2015] Fogel, F., Jenatton, R., Bach, F., and d’Aspremont, A. 2015. Convex relaxations for permutation problems. SIAM Journal on Matrix Analysis and Applications 36, 4, 1465–1488.
- [Fried et al. 2015] Fried, O., DiVerdi, S., Halber, M., Sizikova, E., and Finkelstein, A. 2015. Isomatch: Creating informative grid layouts. In Computer Graphics Forum, vol. 34, Wiley Online Library, 155–166.
- [Gee and Prager 1994] Gee, A. H., and Prager, R. W. 1994. Polyhedral combinatorics and neural networks. Neural computation 6, 1, 161–180.
- [Giorgi et al. 2007] Giorgi, D., Biasotti, S., and Paraboschi, L. 2007. Shape retrieval contest 2007: Watertight models track. SHREC competition 8.
- [Huang et al. 2007] Huang, G. B., Ramesh, M., Berg, T., and Learned-Miller, E. 2007. Labeled faces in the wild: A database for studying face recognition in unconstrained environments. Tech. rep., Technical Report 07-49, University of Massachusetts, Amherst.
- [Kezurer et al. 2015] Kezurer, I., Kovalsky, S. Z., Basri, R., and Lipman, Y. 2015. Tight relaxation of quadratic matching. Comput. Graph. Forum 34, 5 (Aug.), 115–128.
- [Kim et al. 2011] Kim, V. G., Lipman, Y., and Funkhouser, T. 2011. Blended intrinsic maps. In ACM Transactions on Graphics (TOG), vol. 30, ACM, 79.
- [Kolmogorov 2006] Kolmogorov, V. 2006. Convergent tree-reweighted message passing for energy minimization. IEEE transactions on pattern analysis and machine intelligence 28, 10, 1568–1583.
- [Kosowsky and Yuille 1994] Kosowsky, J., and Yuille, A. L. 1994. The invisible hand algorithm: Solving the assignment problem with statistical physics. Neural networks 7, 3, 477–490.
- [Leordeanu and Hebert 2005] Leordeanu, M., and Hebert, M. 2005. A spectral technique for correspondence problems using pairwise constraints. In Tenth IEEE International Conference on Computer Vision (ICCV’05) Volume 1, vol. 2, IEEE, 1482–1489.
- [Lipman and Funkhouser 2009] Lipman, Y., and Funkhouser, T. 2009. Möbius voting for surface correspondence. In ACM Transactions on Graphics (TOG), vol. 28, ACM, 72.
- [Liu et al. 2009] Liu, J., Luo, J., and Shah, M. 2009. Recognizing realistic actions from videos “in the wild”. In Computer Vision and Pattern Recognition, 2009. CVPR 2009. IEEE Conference on, IEEE, 1996–2003.
- [Loiola et al. 2007] Loiola, E. M., de Abreu, N. M. M., Boaventura-Netto, P. O., Hahn, P., and Querido, T. 2007. A survey for the quadratic assignment problem. European journal of operational research 176, 2, 657–690.
- [Lyzinski et al. 2016] Lyzinski, V., Fishkind, D. E., Fiori, M., Vogelstein, J. T., Priebe, C. E., and Sapiro, G. 2016. Graph matching: Relax at your own risk. IEEE Trans. Pattern Anal. Mach. Intell. 38, 1, 60–73.
- [Maron et al. 2016] Maron, H., Dym, N., Kezurer, I., Kovalsky, S., and Lipman, Y. 2016. Point registration via efficient convex relaxation. ACM Trans. Graph. 35, 4 (July), 73:1–73:12.
[Masci et al. 2015]
Masci, J., Boscaini, D., Bronstein, M., and Vandergheynst, P.
Geodesic convolutional neural networks on riemannian manifolds.In Proceedings of the IEEE international conference on computer vision workshops, 37–45.
- [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.
- [Mémoli 2011] Mémoli, F. 2011. Gromov–wasserstein distances and the metric approach to object matching. Foundations of computational mathematics 11, 4, 417–487.
- [Ogier and Beyer 1990] Ogier, R. G., and Beyer, D. 1990. Neural network solution to the link scheduling problem using convex relaxation. In Global Telecommunications Conference, 1990, and Exhibition.’Communications: Connecting the Future’, GLOBECOM’90., IEEE, IEEE, 1371–1376.
- [Ovsjanikov et al. 2012] Ovsjanikov, M., Ben-Chen, M., Solomon, J., Butscher, A., and Guibas, L. 2012. Functional maps: a flexible representation of maps between shapes. ACM Transactions on Graphics (TOG) 31, 4, 30.
- [Parkhi et al. 2015] Parkhi, O. M., Vedaldi, A., and Zisserman, A. 2015. Deep face recognition. In British Machine Vision Conference, vol. 1, 6.
- [Quadrianto et al. 2009] Quadrianto, N., Song, L., and Smola, A. J. 2009. Kernelized sorting. In Advances in neural information processing systems, 1289–1296.
- [Rangarajan et al. 1996] Rangarajan, A., Gold, S., and Mjolsness, E. 1996. A novel optimizing network architecture with applications. Neural Computation 8, 5, 1041–1060.
- [Rangarajan et al. 1997] Rangarajan, A., Yuille, A. L., Gold, S., and Mjolsness, E. 1997. A convergence proof for the softassign quadratic assignment algorithm. Advances in neural information processing systems, 620–626.
- [Rodola et al. 2012] Rodola, E., Bronstein, A. M., Albarelli, A., Bergamasco, F., and Torsello, A. 2012. A game-theoretic approach to deformable shape matching. In Computer Vision and Pattern Recognition (CVPR), 2012 IEEE Conference on, IEEE, 182–189.
[Rodolà et al. 2014]
Rodolà, E., Rota Bulo, S., Windheuser, T., Vestner, M., and Cremers,
Dense non-rigid shape correspondence using random forests.In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 4177–4184.
- [Schellewald et al. 2001] Schellewald, C., Roth, S., and Schnörr, C. 2001. Evaluation of convex optimization techniques for the weighted graph-matching problem in computer vision. In Joint Pattern Recognition Symposium, Springer, 361–368.
- [Shao et al. 2013] Shao, T., Li, W., Zhou, K., Xu, W., Guo, B., and Mitra, N. J. 2013. Interpreting concept sketches. ACM Transactions on Graphics (TOG) 32, 4, 56.
- [Solomon et al. 2012] Solomon, J., Nguyen, A., Butscher, A., Ben-Chen, M., and Guibas, L. 2012. Soft maps between surfaces. In Computer Graphics Forum, vol. 31, Wiley Online Library, 1617–1626.
- [Solomon et al. 2015] Solomon, J., De Goes, F., Peyré, G., Cuturi, M., Butscher, A., Nguyen, A., Du, T., and Guibas, L. 2015. Convolutional wasserstein distances: Efficient optimal transportation on geometric domains. ACM Transactions on Graphics (TOG) 34, 4, 66.
- [Solomon et al. 2016] Solomon, J., Peyré, G., Kim, V. G., and Sra, S. 2016. Entropic metric alignment for correspondence problems. ACM Trans. Graph. 35, 4 (July), 72:1–72:13.
- [Strong and Gong 2014] Strong, G., and Gong, M. 2014. Self-sorting map: An efficient algorithm for presenting multimedia data in structured layouts. IEEE Transactions on Multimedia 16, 4, 1045–1058.
- [Van Kaick et al. 2011] Van Kaick, O., Zhang, H., Hamarneh, G., and Cohen-Or, D. 2011. A survey on shape correspondence. In Computer Graphics Forum, vol. 30, Wiley Online Library, 1681–1707.
- [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).
- [Werner 2007] Werner, T. 2007. A linear programming approach to max-sum problem: A review. IEEE transactions on pattern analysis and machine intelligence 29, 7.
- [Xia 2008] Xia, Y. 2008. Second order cone programming relaxation for quadratic assignment problems. Optimization Methods & Software 23, 3, 441–449.
[Xiao et al. 2010]
Xiao, J., Hays, J., Ehinger, K. A., Oliva, A., and Torralba, A.
Sun database: Large-scale scene recognition from abbey to zoo.In Computer vision and pattern recognition (CVPR), 2010 IEEE conference on, IEEE, 3485–3492.
- [Zaslavskiy et al. 2009] Zaslavskiy, M., Bach, F., and Vert, J.-P. 2009. A path following algorithm for the graph matching problem. IEEE Transactions on Pattern Analysis and Machine Intelligence 31, 12, 2227–2242.
- [Zeng et al. 2010] Zeng, Y., Wang, C., Wang, Y., Gu, X., Samaras, D., and Paragios, N. 2010. Dense non-rigid surface registration using high-order graph matching. In Computer Vision and Pattern Recognition (CVPR), 2010 IEEE Conference on, IEEE, 382–389.
[Zhao et al. 1998]
Zhao, Q., Karisch, S. E., Rendl, F., and Wolkowicz, H.
Semidefinite programming relaxations for the quadratic assignment
Journal of Combinatorial Optimization 2, 1, 71–109.
- [Zuffi and Black 2015] Zuffi, S., and Black, M. J. 2015. The stitched puppet: A graphical model of 3D human shape and pose. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 3537–3546.
Appendix A Proofs
Proof of Lemma 1
The function is strictly convex and satisfies for all extreme points of . Therefore
with equality iff is a permutation.
Proof of Lemma 2
We omit the proof of the first part of the lemma since it is similar to, and somewhat easier than, the proof of the second part.
Let be the minimizer of over the doubly stochastic matrices and let
be the eigenvector ofof unit Euclidean norm corresponding to its minimal eigenvalue . Denote . We define , where we choose so that (6b) holds. This is possible since . Further note that also satisfies (6f) since , and (6e) since
where we used the fact that is a solution to the homogeneous linear equation . Finally the energy satisfies (ignoring the constant )
Now let be a minimizer of (6), we show that is a feasible solution of our relaxation with the same energy. In fact due to the previous claim it is sufficient to show that . The feasibility of is clear since it is already DS. Next, denote
then (ignoring the constant )
and so follows from
where the last inequality follows from the fact that the two matrices in square brackets are positive semi-definite due to the definition of and (6f).
We now prove the third part of the lemma:
Comparison with SDP relaxation
The SDP relaxation of [Kezurer et al. 2015] is
where is the entry replacing the quadratic monomial . We note this relaxation contains all constraints from the SDP relaxation (6) with the exception of (6e). It also contains the additional constraints (9f)-(9h) which do not appear in (6). Thus to show that [Kezurer et al. 2015] is tighter than our relaxation it is sufficient to show that (6e) is implied by the other constraints of [Kezurer et al. 2015]. We recall that (6e) represent all constraints obtained by multiplying linear equality constraints by a linear monomial.
For a quadratic polynomial
let us denote by the linearized polynomial
We will use the following property of SDP relaxations (see [Dym and Lipman 2016]): If a quadratic polynomial is of the form then
For any feasible we have .
If is satisfied for all feasible , then for any quadratic of the form we have .
Accordingly, it is sufficient to show that the squares of all the linear equality polynomials
satisfy . We obtain from
the proof that is identical.
Appendix B Sparsity pattern for improving matching resolution
We construct a sparsity pattern for the task of matching to using known correspondences .
For each we use the following procedure to determine which correspondence will be forbidden: We find the five matched points which are closest to and compute the geodesic distance of these points from . This gives us a feature vector . We then compute the geodesic distances of each of the points from the matched points corresponding to the five closets points to . This gives us feature vectors . For to be a viable match we require that be small. We therefore allow the top of the correspondences according to this criteria.
To symmetrize this process, we use the same procedure to find permissible matches for each , and then select as permissible all matches which were found permissible either when starting from or when starting from .