1 Introduction
Shape matching is a widely researched topic in computer vision and graphics, and a diverse range of techniques that tackle the problem of correspondence have been proposed during the years
[26]. Of particular interest is the case in which the input shapes are allowed to undergo nonrigid deformations, which are typically assumed to be approximately isometric. Recent advancements in this area include the seminal work of Ovsjanikov et al.[17], who proposed modeling functional correspondence between shapes; in this view, the focus shifts from studying pointwise mappings to the definition of a linear operator (the functional map) relating spaces of functions on the two shapes. The classical pointwise representation constitutes then a special case in which the functional map corresponds deltafunctions to deltafunctions. A major advantage of the functional representation lies in the linearity of the operator: the functional map admits a matrix representation which can be made compact under a proper choice of bases for the two functional spaces. In [17]the authors advocated the use of the Laplacian eigenfunctions as the natural basis for smooth functions on the respective shapes; with this choice, one is allowed to “truncate” the representation by using the first few basis functions and still obtain a good approximation to the underlying pointwise correspondence.
Followup work by several authors showed how to extend the framework to deal with nonisometric deformations [19, 11, 21, 12], partial similarity [20], shape exploration [23, 8] and image segmentation [27] among others. However, despite the success of these methods, there has been a general lack of interest on the inverse problem of accurately reconstructing a pointwise map from its functional representation – a common requirement in many practical applications.
The established approach, originally proposed in [17], operates by formulating the conversion problem as a nearestneighbor search in the embedded functional space; as we show in the following sections the approach works well if the functional map is accurate enough, with significant decrease in accuracy as the number of basis functions is reduced. The resulting pointtopoint map can be iteratively refined by following a simple procedure, but this can only be done under specific assumptions on the initial functional map. A similar refinement technique was recently applied for nearisometric partial matching [20], and in a correspondenceless setting for shape retrieval tasks [7]. Finally, the pointwise recovery problem was sidestepped in [12] by adopting a soft error criterion to evaluate the quality of functional maps without converting them to a pointwise counterpart.
Contribution.
In this paper, we consider the problem of accurate pointwise map recovery from a given functional map. The key contributions can be summarized as follows:

We provide the first rigorous analysis of the pointwise map recovery problem. In particular, we show how a simple modification to the baseline approach can lead to consistent improvements.

We introduce a simple probabilistic model for map recovery and refinement. Our model does not impose any assumption on the input functional map, as well as on the specific choice of a functional basis on the two shapes.

Our method is efficient, and significantly outperforms the existing method in both, the nearlyisometric and the interclass settings.
2 Background
We model shapes as compact connected twodimensional Riemannian manifolds (possibly with boundary) endowed with the standard measure induced by the volume form. Let denote the space of squareintegrable functions on , and let be the standard manifold inner product. The space ( features the symmetric LaplaceBeltrami operator (or Laplacian) , which provides us with all the tools of Fourier analysis on our manifold. In particular, this operator admits an eigendecomposition for
, with eigenvalues
and eigenfunctions forming an orthonormal basis on .Drawing an analogy with classical signal processing theory, the eigenfunctions are often referred to as manifold harmonics, and the associated eigenvalues assume the interpretation of frequencies [24]. Any function then admits a Fourier series expansion as
(1) 
Functional correspondence.
Let us be given two manifolds and , and let be a bijective mapping between them. Departing from the traditional pointcentric setting, Ovsjanikov et al.[17] introduced the notion of functional map between two shapes as the linear operator , mapping functions on to functions on via the composition
(2) 
The approach is a natural generalization to classical pointwise correspondence, which can be seen as the special case in which maps indicator functions to indicator functions.
Let and denote orthonormal bases on and respectively. The functional correspondence with respect to these bases can be expressed as follows, for some function :
(3)  
Thus, the action of
amounts to linearly transforming the expansion coefficients of
from basis onto basis . The transformation is encoded in the coefficients , providing a representation for as the (possibly infinite) matrix . Seeking a functional correspondence among the two shapes then amounts to solving for the unknown that better preserves certain mapping constraints [17] or manifesting regularity at different levels [21, 12].Basis truncation.
A natural choice for a basis on the two shapes is given by the eigenfunctions , of the respective Laplacians (harmonic basis). The basis functions are said to be compatible among the two manifolds if the equality holds (approximately) for all , which is the case with the manifold harmonics when the shapes are related by a nearisometry. If the deformation is far from isometric, compatible basis functions can still be computed adhoc for the two manifolds, based on a minimal set of coupling functions (e.g., a sparse set of pointtopoint matches) [11]. In particular, since such basis functions still exhibit a natural ordering in the Fourier sense, they are said to form a quasiharmonic basis.
Assuming to have stable and compact (namely, harmonic) bases at disposal, in [17] the authors proposed to truncate the series (3) at the first coefficients, resulting in a matrix approximating the full map in a compact way. The reduced representation greatly simplifies correspondencebased tasks (e.g., shape matching); at the same time, the truncation has the effect of ‘lowpass’ filtering, thus producing smooth correspondences. In many applications, however, it is desirable to reconstruct the pointtopoint mapping induced by the functional map. Thus, the interest shifts to the inverse problem of recovering the bijection from its functional representation .
3 Pointwise map recovery
Let us now consider the discretized problem and assume that and are represented by a triangular mesh with nodes each. Let matrices contain the first
basis functions for the two shapes respectively, represented as column vectors. Note that, due to truncation, we have that
, but , and similarly for . We assume that the bijection is known. For the sake of simplicity, we additionally assume that can be represented as a permutation matrix , which means that the correspondence originates from two different deformations of the same template. In this case the expression for in (3) can be equivalently rewritten as:(4) 
where . Note that the matrix is now a rank approximation of . The objective of any process converting the functional map back into a pointwise map representation, which is the focus of this work, is to recover the permutation from the sole knowledge of , , and .
Assumptions.
In order to be as general as possible, we do not assume the matching process which generated the given functional map to be known. Additionally, we do not require any particular properties from (e.g., orthogonality), hence allowing to recover maps between shapes undergoing arbitrary deformations. Our only requirement is that the matrix representation is given w.r.t. known bases , . These bases, in turn, may come from diverse optimization processes such as [18, 11, 15].
Mapping indicator functions.
The simplest and most direct way for reconstructing the bijection from the associated functional map consists in mapping indicator functions for each point via , obtaining the image , and then declaring to be the point at which attains the maximum [17]. Such a method, however, suffers from at least two drawbacks. First, it requires constructing and mapping indicator functions for all shape points, which can get easily expensive for large meshes. Second, the lowpass filtering due to the basis truncation has a delocalizing effect on the maximum of (see inset figure), negatively affecting the quality of the final correspondence.
The inverse problem of pointtopoint map recovery.
Considering the problem of recovering from a given according to (4) as a (highly underdetermined) illposed inverse problem, the natural recovery problem to consider is
(5)  
(6) 
for a suitable measure of distance , a regularization function to possibly impose some kind of desired smoothness of the transformation, and a regularization parameter determining a tradeoff between fidelity and regularity.
Unfortunately, the minimization of (5) will be extremely challenging in general. Consider for instance the case of no regularization, , a zero functional map , and the distance measure being the squared Frobenius norm. In this case the minimization problem becomes
for the right hand sides denoting vectorized equations with being formed by the Kronecker product between and . For arbitrary the above problem is a particular reformulation of the quadratic integer programming problem (see equation (8) in [4]). Since the latter was shown to be NPhard we cannot expect to find exact solutions to the above problem with reasonable computational complexity in general. Therefore, we turn our attention to approximations of (5).
Linear assignment problem.
Instead of applying general greedy methods or relaxation techniques to (5), let us recall some observations from [17] regarding the general structure of (4): In particular, note that the indicator function around point has coefficients in the Laplacian basis. This can be seen by observing that the inner product is selecting the column of corresponding to point . Therefore, the image via of all indicator functions centered at points of is simply given by . Recovering the pointtopoint map could then be solved by finding the correspondence between the columns of and the columns of . Measuring the proximity between these columns in the sense gives rise to the linear assignment problem
(7)  
(8) 
Although the problem above admits a polynomial time solution [13], typical values for (in the order of thousands) make computing this solution prohibitively expensive in practice.
Nearest neighbors.
The authors of [17] circumvent the computational costs of the above approach by proposing a nearestneighbor technique for the recovery of the pointtopoint correspondence. In the light of our previous analysis their idea is to consider the matching of every point, i.e., column of , to its nearest neighbor in separately.
Mathematically, the nearestneighbor approach can be seen as a relaxation of problem (7), (8), in which one seeks for the best leftstochastic approximation , i.e.,
(9)  
(10) 
In other words, in comparison to (8) one omits the constraint of . The omission allows to minimize the problem above by independently solving for columns of , one per query. The drawback of such a separable optimization approach, however, is that it may produce onetomany mappings as a result of the recovery process. Moreover, the nearestneighbor approach is an asymmetric method: looking for nearest neighbors from to , or viceversa, will in general produce different results.
Balanced nearest neighbors.
In order to remove the bias, we propose to incorporate additional terms in problem (9), namely minimize
(11) 
where the minimization is performed w.r.t. both and . Note that we incorporated the desired property of being a permutation matrix by a soft constraint, i.e., by penalizing the difference of to . Also note that the limit of leads to a convergence of to a solution of (7) meeting (8).
Instead of solving the minimization problem for increasing values of exactly, we determine an approximate solution by alternating minimization on and . The latter leads to each subproblem being a simple nearest neighbor problem and guarantees the decrease of the objective functions.
Table 1 illustrates the matching accuracy obtained by our symmetrized nearestneighbor method in comparison to the classical nearestneighbor as well as the indicator mapping approaches. As we can see, the symmetrization improves the results of the biased nearestneighbor method by 2–3.
# basis functions  Max  NN  Balanced NN 

25  6.14  30.99  33.24 
50  18.12  43.51  45.65 
75  26.96  52.54  55.07 
A probabilistic model.
The analysis we provided above puts in evidence two major drawbacks, namely: 1) The linear assignment approach, the nearestneighbor search, and our bidirectional variant, all rely on the assumption that the functional map given as input aligns well the columns of with those of in the sense. 2) None of the above approaches incorporates regularity assumptions for the alignment process, i.e., the regularization term in the general inverse problem formulation (5) was omitted.
We propose to cast the pointtopoint map recovery as a probability density estimation problem to obtain both, a better measure of proximity than the
distance and a tool to impose regularity assumptions on the alignment map. Within our model, we interpret the columns of as modes of a continuousprobability distribution defined over (the embedded functional space), while columns of constitute the data, i.e., a discrete sample drawn from the distribution. The task is then to align the modes to the data, such that the pointtopoint mapping can be recovered as the maximum posterior probability.
As a model for the distribution we consider a Gaussian mixture (GMM) with components, having the columns of as centroids in . For simplicity, we assume the components have uniform weight , and equal covariances . With this choice, the GMM density function is:
(12) 
where we write and for to denote the columns of and respectively, and define .
Now let denote the (unknown) transformation aligning the centroid locations to the data points, according to a set of transformation parameters . The alignment problem can then be solved by maximizing the likelihood, or equivalently by minimizing the negative loglikelihood function
(13) 
Note that the argument that minimizes (13) can be also interpreted as the argument that minimizes the KullbackLeibler (KL) divergence between a continuous GMM distribution (represented by ) and a mixture of Dirac distributions (represented by ). Hence, with our probabilistic model we are choosing the distance in Eq. (5) to be the (pseudo)distance .
Given optimal parameters and , the pointtopoint correspondence probability between and can finally be obtained as the posterior probability .
Transformation function.
The above probabilistic model leaves some freedom for the specific choice of a transformation . A simple example is choosing this transformation as a simple rotation, parametrized by . To be more flexible we instead propose to consider the general transformation:
(14) 
where parameters assume the meaning of a displacement field, and . With this definition, the overall behavior of the transformation can be controlled by regularizing .
Assuming that the given functional map represents a reasonably good matching, the refinement procedure should not be allowed to perturb the initial alignment significantly. Here we adopt the Tikhonov regularizer proposed in [28, 14], where is a lowpass operator promoting smoothly changing velocity vectors, hence coherent motion. In the inset we illustrate a smooth velocity field with coherent correspondences (blue) and a mismatch produced by simple nearestneighbors (red). We finally obtain the regularized objective:
(15) 
where is a tradeoff parameter between likelihood and regularity (set to in our experiments).
General alignment problems like (15) have been widely researched in computer vision, and several robust algorithms exist for these tasks [5, 25, 14, 9]. Most of these approaches follow an iterative scheme, optimizing w.r.t. and in an alternating fashion until convergence (EM algorithm [6]). In our experiments, we used publicly available code from [14], which allows to optimize over smooth displacements as in Eq. (15).
In Figure 2 we illustrate a few iterations of the refinement process applied to a pair of nearlyisometric shapes, starting from a functional map obtained by the matching algorithm described in Section 5. As a visual measure of map quality we employ the technique described in [16] (using the top 5 singular vectors), which allows to identify the problematic areas induced by a given functional map. The quality of a map can be judged by the smoothness of the associated plots, with better maps having a more localized behavior.
The output of the EM algorithm is a set of optimal transformation parameters and a leftstochastic correspondence .
4 Pointwise map refinement
As a general representation for shape correspondences, functional maps can be adopted to compactly encode (via Eq. (4)) any dense pointtopoint map obtained by generic matching algorithms. Therefore, one can consider locally refining the functional map with the help of the pointtopoint map recovered with any of the methods described in the previous section. Naturally, one can repeat such a strategy and iterate between updating the pointtopoint correspondence and refining the functional map.
Such an iterative procedure was considered in [17], in which the authors proposed to use the classical Iterative Closest Point (ICP) algorithm [1]. The latter updates according to the nearestneighbor approach (9), followed by a refinement of via
(16)  
(17) 
which is an orthogonal Procrustes problem. Intuitively, this can be seen as a rigid alignment in between point sets and (see Figure 4(b) for an example). The alternating minimization w.r.t. and is repeated until convergence.
Although the ICP approach described above allows to achieve significant improvements in terms of map accuracy, the orthogonality constraints (17) imposed on the functional correspondence limit its applicability to a specific class of transformations, namely the class of volumepreserving isometries (see [17] Theorem 5.1). Therefore, the method cannot be applied to refine maps between shapes undergoing arbitrary deformations.
For a more general refinement procedure, we drop the orthogonality constraint on and consider the problem
(18) 
for a regularization functional which could encourage to correspond to a smooth transformation, or could require to be a small perturbation of the identity. While for the specific example
problem (18) coincides with the rigid alignment problem arising from the constraint (17), a less restrictive choice for the regularization functional makes the method suitable for recovering functional maps for nonisometric shape matching problems.
In our experiments we found that when (18) is combined with our proposed regularized probabilistic model for recovering the pointtopoint correspondence, it is sufficient to simply update in a leastsquares sense: Even without additional regularization, is determined to be a perturbation of the identity as illustrated in Figure 3. The fact that is refined in a nonrigid fashion can improve the refinement results significantly in comparison to orthogonal updates of as illustrated in 3d in Figure 4.
5 Experimental evaluation
We compare our iteratively refined probabilistic pointwise map recovery method with the iterative refinement procedure of Ovsjanikov et al.[17] (denoted as ICP, see Section 4), which is to the best of our knowledge the only existing alternative to date. Both algorithms were implemented in Matlab/MEX and executed (singlecore) on an Intel i73770 3.4GHz CPU with 32GB memory.
As a measure of error, we use the quantitative criterion that was introduced in [10] to evaluate the quality of pointwise maps. The input quantity in our case is a functional map , which is then converted to its pointwise counterpart using the two methods. We plot cumulative curves showing the percent of matches which have geodesic error smaller than a variable threshold.
We evaluate the two methods quantitatively on the FAUST [2] dataset, and qualitatively on the TOSCA [3] and KIDS [21] datasets. The three datasets include isometric as well as nonisometric shapes; in particular, FAUST and KIDS also include pointtopoint ground truth matches between shapes belonging to different classes.
Comparisons.
The functional maps used in the comparisons are constructed by solving a leastsquares system , where matrices and contain the Fourier coefficients of indicator functions for corresponding regions on the two shapes. The region correspondence is established using the ground truth, while the sets of regions are computed using the consensus technique of [22]. This way, we simulate a matching process that provides reasonably good solutions for further refinement.
We show comparisons both in the nearisometric and interclass settings. In the former case, we use 45 pairs of humans in different poses from the intraclass subset of FAUST. All shapes have points, and the functional map is computed with basis functions. Results are reported in Figure 5 (left), where we also included plain nearestneighbors (NN) as a baseline. From the plotted curves we can see that orthogonal ICP is performing slightly better than our method on nearisometric deformations, since the approach is specifically tailored for this particular case. However, initializing our method with the output of ICP allows to achieve further improvement on average.
In Figure 5 (right) we show the same curves for the nonisometric case (interclass matching). In this case, orthogonal refinement cannot be applied due to the different properties of the input functional maps, which relate shapes under non volumepreserving transformations. Additional qualitative comparisons are shown in Figure 6.
Rank.
In a second set of experiments, we evaluate the capability of the different methods to recover pointwise maps from functional maps of increasing rank. In this setting, we assume the input functional map to be as accurate as possible, and for this purpose we construct it explicitly as , where is the groundtruth permutation among the vertices of the two shapes. We do so for a pair of approximately isometric shapes, so that the respective eigenbases and are as repeatable as possible, and further orthogonal refinement is not needed (indeed, applying ICP in this setting actually yielded worse results in our tests).
The results are shown in Figure 7. As the number of basis functions used on the two shapes (i.e., rank of ) increases, so does the amount of exact correspondences recovered by each method. This is also true for the simple indicator mapping approach (Max), since the smoothing effect due to basis truncation is reduced at increasing values of .
Our method allows to recover up to 20% more exact matches than the nearestneighbors approach. In particular, with (a commonly used value in most shape matching pipelines) we are able to perfectly reconstruct 75% of the rows/columns of (a matrix in this example).
Complexity issues.
The time performance of our method depends on two factors: the number of shape points , and the size of the functional map . As we also show in Figure 2, typically a few iterations of the EM algorithm are sufficient to reach accurate solutions, and in practice we used 5 iterations in all our experiments. In the typical case where and , our method takes on average 1 min. to converge, while ICP using efficient search structures (kdtrees) adds up to 3 sec.
It should be noted, however, that while we employed an offtheshelf implementation of the minimization algorithm, this code can be easily parallelized and optimized in several ways. In particular, a GPU implementation of our method remains a practical possibility.
6 Discussion and conclusions
In this paper we considered the problem of pointtopoint map recovery and refinement from its lowrank functional counterpart. Despite the growing success of this representation, the problem has received limited attention to date, and no significant steps have been made in this direction since the framework’s conception. We formulated a general variational recovery approach for the inverse problem of computing pointtopoint correspondences from a given functional map. We demonstrated what simplifications can be used to arrive at the nearest neighbor approach and showed how simply mitigating the asymmetry present in the standard conversion procedure can lead to consistent improvements of 23% in accuracy. We then introduced a probabilistic model for pointwise map recovery and considered a refinement of the functional map that does not rely on the assumption of isometric shapes. The experimental results showed that the proposed approach yields up to 20% accuracy improvements under nonisometric deformations, reaching up to 75% exact pointtopoint matches under good initializations.
The main limitation of our probabilistic method lies in the fact that – similar to the nearestneighbors approach – the optimization procedure is biased towards one of the two shapes, as one can for instance see from the interpretation of minimizing the (nonsymmetric) KullbackLeibler divergence. Extending the ideas of the symmetrized nearestneighbor approach to the probabilistic model for removing this bias represents a possibility. Second, while for simplicity we only considered pairs of shapes related by a bijection, our method can be modified to deal with shapes having different resolutions, as well as partially similar shapes and entire shape collections (joint refinement). We believe these topics to represent exciting directions of future research.
Acknowledgments
We thankfully acknowledge Zorah Lähner and Federico Tombari for useful discussions. E.R. is supported by an Alexander von Humboldt fellowship.
References
 [1] Paul J. Besl and Neil D. McKay. A method for registration of 3d shapes. TPAMI, 14(2):239–256, February 1992.
 [2] Federica Bogo, Javier Romero, Matthew Loper, and Michael J. Black. FAUST: Dataset and evaluation for 3D mesh registration. In Proc. CVPR, June 2014.
 [3] Alexander Bronstein, Michael Bronstein, and Ron Kimmel. Numerical Geometry of NonRigid Shapes. Springer, 2008.

[4]
R.E. Burkard, E. Çela, P.M. Pardalos, and L.S. Pitsoulis.
The quadratic assignment problem.
In
Handbook of Combinatorial Optimization
, pages 241–338. Kluwer Academic Publishers, 1998.  [5] H. Chui and A. Rangarajan. A feature registration framework using mixture models. In Proc. MMBIA, pages 190–197, 2000.
 [6] A. P. Dempster, N. M. Laird, and D. B. Rubin. Maximum likelihood from incomplete data via the EM algorithm. J. R. Stat. Soc. Series B, 39(1):1–38, 1977.
 [7] A. Gasparetto and A. Torsello. A statistical model of Riemannian metric variation for deformable shape analysis. In Proc. CVPR, pages 1219–1228, 2015.
 [8] Qixing Huang, Fan Wang, and Leonidas Guibas. Functional map networks for analyzing and exploring large shape collections. ACM Trans. Graph., 33(4):36:1–36:11, July 2014.

[9]
Bing Jian and B.C. Vemuri.
Robust point set registration using Gaussian mixture models.
TPAMI, 33(8):1633–1645, 2011.  [10] Vladimir G. Kim, Yaron Lipman, and Thomas A. Funkhouser. Blended intrinsic maps. TOG, 30(4):79, 2011.
 [11] A. Kovnatsky, M. Bronstein, A. Bronstein, K. Glashoff, and R. Kimmel. Coupled quasiharmonic bases. Comput. Graph. Forum, 32(2pt4):439–448, 2013.
 [12] A. Kovnatsky, M. M. Bronstein, X. Bresson, and P. Vandergheynst. Functional correspondence by matrix completion. In Proc. CVPR, pages 905–914, 2015.
 [13] Harold W. Kuhn. The Hungarian method for the assignment problem. Naval Research Logistics Quarterly, 2(1–2):83–97, March 1955.
 [14] Andriy Myronenko and Xubo Song. Point set registration: Coherent point drift. TPAMI, 32(12):2262–2275, 2010.
 [15] Thomas Neumann, Kiran Varanasi, Christian Theobalt, Marcus Magnor, and Markus Wacker. Compressed manifold modes for mesh processing. Computer Graphics Forum, 33(5):1–10, 2014.
 [16] M. Ovsjanikov, M. BenChen, F. Chazal, and L. Guibas. Analysis and visualization of maps between shapes. Comput. Graph. Forum, 32(6):135–145, September 2013.
 [17] M. Ovsjanikov, M. BenChen, J. Solomon, A. Butscher, and L. Guibas. Functional maps: a flexible representation of maps between shapes. ACM Trans. Graph., 31(4):30:1–30:11, 2012.
 [18] U. Pinkall and K. Polthier. Computing discrete minimal surfaces and their conjugates. Experimental mathematics, 2(1):15–36, 1993.
 [19] J. Pokrass, A. M. Bronstein, M. M. Bronstein, P. Sprechmann, and G. Sapiro. Sparse modeling of intrinsic correspondences. Computer Graphics Forum, 32(2pt4):459–468, 2013.
 [20] E. Rodolà, L. Cosmo, M. M. Bronstein, A. Torsello, and D. Cremers. Partial functional correspondence. CoRR, abs/1506.05274, 2015.

[21]
E. Rodolà, S. Rota Bulò, T. Windheuser, M. Vestner, and D. Cremers.
Dense nonrigid shape correspondence using random forests.
In Proc. CVPR, pages 4177–4184, 2014.  [22] Emanuele Rodolà, Samuel Rota Bulò, and Daniel Cremers. Robust region detection via consensus segmentation of deformable shapes. Computer Graphics Forum, 33(5), 2014.
 [23] Raif M. Rustamov, Maks Ovsjanikov, Omri Azencot, Mirela BenChen, Frédéric Chazal, and Leonidas Guibas. Mapbased exploration of intrinsic shape differences and variability. ACM Trans. Graph., 32(4):72:1–72:12, July 2013.
 [24] Gabriel Taubin. A signal processing approach to fair surface design. In Proc. SIGGRAPH, pages 351–358. ACM, 1995.
 [25] Yanghai Tsin and Takeo Kanade. A correlationbased approach to robust point set registration. In Proc. ECCV, volume 3023, pages 558–569. Springer Berlin Heidelberg, 2004.
 [26] Oliver Van Kaick, Hao Zhang, Ghassan Hamarneh, and Daniel CohenOr. A survey on shape correspondence. Computer Graphics Forum, 30(6):1681–1707, 2011.
 [27] Fan Wang, Qixing Huang, and Leonidas Guibas. Image cosegmentation via consistent functional maps. In Proc. ICCV, 2013.
 [28] Alan L. Yuille and Norberto M. Grzywacz. A mathematical analysis of the motion coherence theory. International Journal of Computer Vision, 3(2):155–175, 1989.